diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 55a8271a2c6eb547cc50503d0fcf24d4ac59c785..0bc65fc85c5e8f9a650a212c0517a2c0d21b8e6b 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -32,6 +32,14 @@ symmetrize_vector(const TinyVector<Dimension>& normal, const TinyVector<Dimensio
   return u - 2 * dot(u, normal) * normal;
 }
 
+template <size_t Dimension>
+PUGS_INLINE auto
+symmetrize_matrix(const TinyVector<Dimension>& normal, const TinyMatrix<Dimension>& A)
+{
+  const TinyMatrix S = TinyMatrix<Dimension>{identity} - 2 * tensorProduct(normal, normal);
+  return S * A * S;
+}
+
 template <size_t Dimension>
 PUGS_INLINE auto
 symmetrize_coordinates(const TinyVector<Dimension>& origin,
@@ -790,6 +798,44 @@ PolynomialReconstruction::_createMutableDiscreteFunctionDPKVariantList(
   return mutable_discrete_function_dpk_variant_list;
 }
 
+template <MeshConcept MeshType>
+void
+PolynomialReconstruction::_checkDataAndSymmetriesCompatibility(
+  const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+{
+  for (auto&& discrete_function_variant : discrete_function_variant_list) {
+    std::visit(
+      [&](auto&& discrete_function) {
+        using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+        if constexpr (is_discrete_function_P0_v<DiscreteFunctionT> or
+                      is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+          using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
+          if constexpr (is_tiny_vector_v<DataType>) {
+            if constexpr (DataType::Dimension != MeshType::Dimension) {
+              std::stringstream error_msg;
+              error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
+                        << " using a mesh of dimension " << MeshType::Dimension;
+              throw NormalError(error_msg.str());
+            }
+          } else if constexpr (is_tiny_matrix_v<DataType>) {
+            if constexpr ((DataType::NumberOfRows != MeshType::Dimension) or
+                          (DataType::NumberOfColumns != MeshType::Dimension)) {
+              std::stringstream error_msg;
+              error_msg << "cannot symmetrize matrices of dimensions " << DataType::NumberOfRows << 'x'
+                        << DataType::NumberOfColumns << " using a mesh of dimension " << MeshType::Dimension;
+              throw NormalError(error_msg.str());
+            }
+          }
+        } else {
+          // LCOV_EXCL_START
+          throw UnexpectedError("invalid discrete function type");
+          // LCOV_EXCL_STOP
+        }
+      },
+      discrete_function_variant->discreteFunction());
+  }
+}
+
 template <MeshConcept MeshType>
 [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
 PolynomialReconstruction::_build(
@@ -800,6 +846,10 @@ PolynomialReconstruction::_build(
 
   using Rd = TinyVector<MeshType::Dimension>;
 
+  if (m_descriptor.symmetryBoundaryDescriptorList().size() > 0) {
+    this->_checkDataAndSymmetriesCompatibility<MeshType>(discrete_function_variant_list);
+  }
+
   const size_t number_of_columns = this->_getNumberOfColumns(discrete_function_variant_list);
 
   const size_t basis_dimension =
@@ -954,21 +1004,32 @@ PolynomialReconstruction::_build(
                           B(index, kB) = qi_qj[k];
                         }
                       } else {
-                        const DataType& qi_qj = discrete_function[cell_i_id] - qj;
-                        for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
-                          B(index, kB) = qi_qj[k];
-                        }
+                        // 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
                       }
                     } else if constexpr (is_tiny_matrix_v<DataType>) {
                       if constexpr ((DataType::NumberOfColumns == DataType::NumberOfRows) and
                                     (DataType::NumberOfColumns == MeshType::Dimension)) {
-                        throw NotImplementedError("TinyMatrix symmetries for reconstruction of DiscreteFunctionP0");
-                      }
-                      const DataType& qi_qj = discrete_function[cell_i_id] - 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);
+                        const Rd& normal = symmetry_normal_list[i_symmetry];
+
+                        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);
+                          }
                         }
+                      } else {
+                        // LCOV_EXCL_START
+                        std::stringstream error_msg;
+                        error_msg << "cannot symmetrize matrices of dimensions " << DataType::NumberOfRows << 'x'
+                                  << DataType::NumberOfColumns << " using a mesh of dimension " << MeshType::Dimension;
+                        throw UnexpectedError(error_msg.str());
+                        // LCOV_EXCL_STOP
                       }
                     }
                   }
@@ -1031,8 +1092,6 @@ PolynomialReconstruction::_build(
                     }
                   }
                 } else if constexpr (is_tiny_matrix_v<DataType>) {
-                  throw NotImplementedError("NIY TinyMatrix data");
-
                   for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
                        ++i_symmetry) {
                     auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
index b2b855271d1ed4140ebb51d1b21314602caa0dc6..bf6013639f732fd52f9c145ba4982701c0e4bcfc 100644
--- a/src/scheme/PolynomialReconstruction.hpp
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -17,6 +17,10 @@ class PolynomialReconstruction
   size_t _getNumberOfColumns(
     const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
 
+  template <MeshConcept MeshType>
+  void _checkDataAndSymmetriesCompatibility(
+    const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
+
   template <MeshConcept MeshType>
   std::vector<MutableDiscreteFunctionDPkVariant> _createMutableDiscreteFunctionDPKVariantList(
     const std::shared_ptr<const MeshType>& p_mesh,