diff --git a/src/scheme/HyperplasticSolverO2.cpp b/src/scheme/HyperplasticSolverO2.cpp
index c99d0100cb902f8d274696acaa2ff179bb996565..1d6174dd139213a263755dd1fb311484a74b4c86 100644
--- a/src/scheme/HyperplasticSolverO2.cpp
+++ b/src/scheme/HyperplasticSolverO2.cpp
@@ -150,7 +150,10 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
     double chi          = 0;
     const double normS2 = frobeniusNorm(Sd) * frobeniusNorm(Sd);
     double limit2       = 2. / 3 * yield * yield;
-    chi                 = std::max(0., std::sqrt(normS2) - std::sqrt(limit2));
+    const double power  = 0.5;
+    if ((normS2 - limit2) > 0) {
+      chi = std::pow((normS2 - limit2), power);
+    }
     return chi;
   }
 
@@ -189,7 +192,7 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
                                                                         symmetry_boundary_descriptor_list);
     DiscreteFunctionDPk<Dimension, Rd> DPk_fh_lim = copy(DPk_fh);
     const auto xj                                 = mesh_data.xj();
-
+    // return DPk_fh_lim;
     const double eps2 = 1E-14;
     double cmin       = 1. - eps;
     double cmax       = 1. + eps;
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 2af32e157475b0fe99b9f2cc9cc6c445e05681b3..4ae9e11f48a4455aef7a017b589b6b1ec719b968 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -25,6 +25,15 @@
 #include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/DiscreteFunctionVariant.hpp>
 
+template <size_t Dimension>
+PUGS_INLINE auto
+symmetrize_matrix(const TinyVector<Dimension>& normal, const TinyMatrix<Dimension>& M)
+{
+  const TinyMatrix<Dimension> Id = identity;
+  const TinyMatrix<Dimension> S  = Id - 2. * tensorProduct(normal, normal);
+  return S * M * S;
+}
+
 template <size_t Dimension>
 PUGS_INLINE auto
 symmetrize_vector(const TinyVector<Dimension>& normal, const TinyVector<Dimension>& u)
@@ -927,12 +936,20 @@ PolynomialReconstruction::_build(
                     } 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");
-                      }
-                      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 {
+                        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);
+                          }
                         }
                       }
                     }