diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index edd217399d54ab9536e346046abdd6fa146d1785..fbe0e344bb8a0d011f30eabcd84d9406cdc1bee7 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -234,43 +234,6 @@ EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, S
 {
   Internals::slepscComputeForSymmetricMatrix(A, eigenvalues);
 }
-
-void
-EigenvalueSolver::computeExpForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallMatrix<double>& expA)
-{
-  EPS eps;
-
-  EPSCreate(PETSC_COMM_SELF, &eps);
-  EPSSetOperators(eps, A, nullptr);
-  EPSSetProblemType(eps, EPS_HEP);
-
-  Internals::computeAllEigenvaluesOfSymmetricMatrix(eps);
-
-  PetscInt nb_eigenvalues;
-  EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
-
-  SmallArray<double> eigenvalues(nb_eigenvalues);
-  SmallMatrix<double> P(nb_eigenvalues, nb_eigenvalues);
-  SmallMatrix<double> D(nb_eigenvalues, nb_eigenvalues);
-  D    = zero;
-  expA = SmallMatrix<double>(nb_eigenvalues, nb_eigenvalues);
-  Array<double> eigenvector(nb_eigenvalues);
-  for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
-    Vec Vr;
-    VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr);
-    EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr);
-    VecDestroy(&Vr);
-    for (size_t j = 0; j < eigenvector.size(); ++j) {
-      P(j, i) = eigenvector[j];
-    }
-  }
-
-  for (size_t i = 0; i < eigenvalues.size(); ++i) {
-    D(i, i) = std::exp(eigenvalues[i]);
-  }
-  expA = P * D * transpose(P);
-  EPSDestroy(&eps);
-}
 #endif   // PUGS_HAS_SLEPC
 
 template <typename T>
@@ -329,7 +292,7 @@ EigenvalueSolver::findEigenvector(const TinyMatrix<3, 3>& A,
   if (space_size == 3) {
     // std::cout << "space_size = 3 "
     //           << "\n";
-    return eigenvectors0;
+    return eigenvectors;
   }
   TinyMatrix<3, 3> B = A;
 
@@ -485,15 +448,17 @@ EigenvalueSolver::findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eige
 std::tuple<TinyVector<3>, TinyMatrix<3>>
 EigenvalueSolver::findEigen(const TinyMatrix<3> A)
 {
-  const double epsilon = 1.e-6;   // * l2Norm(eigenvalues);
+  constexpr TinyVector<3> eigenvalues0  = zero;
+  constexpr TinyMatrix<3> eigenvectors0 = identity;
+  const double epsilon                  = 1.e-6;   // * l2Norm(eigenvalues);
   if (frobeniusNorm(A) < 1.e-15) {
     // std::cout << " Frobenius norm 0 " << frobeniusNorm(A) << "\n";
     return std::make_tuple(eigenvalues0, eigenvectors0);
   }
   TinyMatrix<3> C = 1. / frobeniusNorm(A) * A;
   if (isDiagonal(C, epsilon)) {
-    // std::cout << "Matrix C " << C << " is diagonal "
-    //           << "\n";
+    std::cout << "Matrix C " << C << " is diagonal "
+              << "\n";
     const TinyVector<3> eigenvalues(A(0, 0), A(1, 1), A(2, 2));
     TinyMatrix<3> Diag = zero;
     for (size_t i = 0; i < 3; ++i) {
@@ -630,7 +595,9 @@ EigenvalueSolver::computeExpForTinyMatrix(const TinyMatrix<3>& A)
   return expA;
 }
 
-EigenvalueSolver::EigenvalueSolver() {}
+#ifdef PUGS_HAS_SLEPC
+
+void
 EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues)
 {
   Internals::slepscComputeForSymmetricMatrix(A, eigenvalues);
diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp
index 1c0024211d89229cc2e19fc34558d0ff744d19d3..98696152215164477c34e4743accb5ad4659e85d 100644
--- a/src/algebra/EigenvalueSolver.hpp
+++ b/src/algebra/EigenvalueSolver.hpp
@@ -16,18 +16,11 @@ class EigenvalueSolver
 {
  private:
   struct Internals;
-  TinyMatrix<3> eigenvectors0{1, 0, 0, 0, 1, 0, 0, 0, 1};
-  TinyVector<3> eigenvalues0{0, 0, 0};
 
   const EigenvalueSolverOptions m_options;
 
   void _slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, SmallArray<double>& eigenvalues);
 
-  void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
-                                 SmallArray<double>& eigenvalues,
-                                 SmallMatrix<double>& P);
-  void computeExpForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallMatrix<double>& expA);
-#endif   // PUGS_HAS_SLEPC
   void _slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues);
 
   void _slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A,
@@ -142,7 +135,7 @@ class EigenvalueSolver
 #ifdef PUGS_HAS_SLEPC
     this->computeExpForSymmetricMatrix(PETScAijMatrixEmbedder{A}, expA);
 #else    // PUGS_HAS_SLEPC
-  throw NotImplementedError("SLEPc is required to solve eigenvalue problems");
+    throw NotImplementedError("SLEPc is required to solve eigenvalue problems");
 #endif   // PUGS_HAS_SLEPC
   }
 
@@ -170,7 +163,6 @@ class EigenvalueSolver
 
   TinyMatrix<3> computeExpForTinyMatrix(const TinyMatrix<3>& A);
 
-  EigenvalueSolver();
   EigenvalueSolver(const EigenvalueSolverOptions& options = EigenvalueSolverOptions::default_options);
   ~EigenvalueSolver() = default;
 };
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 1d5346666727b2e599bf7b44ece73854d2207167..23456fc1da66e86d605cc33694e604af49a3cabf 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -77,7 +77,7 @@ SchemeModule::SchemeModule()
   this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IBoundaryConditionDescriptor>>);
   this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const VariableBCDescriptor>>);
 
-  this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const BasisType>>);
+  //  this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const BasisType>>);
 
   this->_addBuiltinFunction("P0", std::function(
 
@@ -527,30 +527,6 @@ SchemeModule::SchemeModule()
 
                               ));
 
-  this->_addBuiltinFunction("canonicalBasis", std::function(
-
-                                                []() -> std::shared_ptr<const BasisType> {
-                                                  return std::make_shared<BasisType>(BasisType::canonical);
-                                                }
-
-                                                ));
-
-  this->_addBuiltinFunction("taylorBasis", std::function(
-
-                                             []() -> std::shared_ptr<const BasisType> {
-                                               return std::make_shared<BasisType>(BasisType::taylor);
-                                             }
-
-                                             ));
-
-  this->_addBuiltinFunction("lagrangeBasis", std::function(
-
-                                               []() -> std::shared_ptr<const BasisType> {
-                                                 return std::make_shared<BasisType>(BasisType::lagrange);
-                                               }
-
-                                               ));
-
   this->_addBuiltinFunction("apply_acoustic_fluxes",
                             std::function(
 
@@ -1356,6 +1332,8 @@ SchemeModule::SchemeModule()
                                  const std::shared_ptr<const DiscreteFunctionVariant>& sigma,    //
                                  const std::shared_ptr<const ItemValueVariant>& ur,              //
                                  const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,   //
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,   //
                                  const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
@@ -1364,7 +1342,8 @@ SchemeModule::SchemeModule()
                                 return HyperplasticSolverO2Handler{getCommonMesh({rho, u, E, Fe})}   //
                                   .solver()
                                   .apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr,
-                                                HyperplasticSolverO2Handler::RelaxationType::Implicit);
+                                                HyperplasticSolverO2Handler::RelaxationType::Implicit,
+                                                bc_descriptor_list);
                               }
 
                               ));
@@ -1381,6 +1360,8 @@ SchemeModule::SchemeModule()
                                  const std::shared_ptr<const DiscreteFunctionVariant>& sigma,    //
                                  const std::shared_ptr<const ItemValueVariant>& ur,              //
                                  const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,   //
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,   //
                                  const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
@@ -1389,7 +1370,8 @@ SchemeModule::SchemeModule()
                                 return HyperplasticSolverO2Handler{getCommonMesh({rho, u, E, Fe})}   //
                                   .solver()
                                   .apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr,
-                                                HyperplasticSolverO2Handler::RelaxationType::CauchyGreen);
+                                                HyperplasticSolverO2Handler::RelaxationType::CauchyGreen,
+                                                bc_descriptor_list);
                               }
 
                               ));
@@ -1406,6 +1388,8 @@ SchemeModule::SchemeModule()
                                  const std::shared_ptr<const DiscreteFunctionVariant>& sigma,    //
                                  const std::shared_ptr<const ItemValueVariant>& ur,              //
                                  const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,   //
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,   //
                                  const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
@@ -1414,7 +1398,8 @@ SchemeModule::SchemeModule()
                                 return HyperplasticSolverO2Handler{getCommonMesh({rho, u, E, Fe})}   //
                                   .solver()
                                   .apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr,
-                                                HyperplasticSolverO2Handler::RelaxationType::Exponential);
+                                                HyperplasticSolverO2Handler::RelaxationType::Exponential,
+                                                bc_descriptor_list);
                               }
 
                               ));
diff --git a/src/language/modules/SchemeModule.hpp b/src/language/modules/SchemeModule.hpp
index defa9a18fd0146e3da6f71594e981039f3c13d42..6cd0b03427463debf965b9f5662fd21b33f502e3 100644
--- a/src/language/modules/SchemeModule.hpp
+++ b/src/language/modules/SchemeModule.hpp
@@ -3,38 +3,6 @@
 
 #include <analysis/PolynomialBasis.hpp>
 #include <language/modules/BuiltinModule.hpp>
-#include <language/utils/ASTNodeDataTypeTraits.hpp>
-#include <utils/PugsMacros.hpp>
-
-class IBoundaryConditionDescriptor;
-template <>
-inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IBoundaryConditionDescriptor>> =
-  ASTNodeDataType::build<ASTNodeDataType::type_id_t>("boundary_condition");
-
-template <>
-inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const BasisType>> =
-  ASTNodeDataType::build<ASTNodeDataType::type_id_t>("basis_type");
-
-class VariableBCDescriptor;
-
-template <>
-inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const VariableBCDescriptor>> =
-  ASTNodeDataType::build<ASTNodeDataType::type_id_t>("variable_boundary_condition");
-
-class DiscreteFunctionVariant;
-template <>
-inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const DiscreteFunctionVariant>> =
-  ASTNodeDataType::build<ASTNodeDataType::type_id_t>("Vh");
-
-class IDiscreteFunctionDescriptor;
-template <>
-inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IDiscreteFunctionDescriptor>> =
-  ASTNodeDataType::build<ASTNodeDataType::type_id_t>("Vh_type");
-
-class IQuadratureDescriptor;
-template <>
-inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IQuadratureDescriptor>> =
-  ASTNodeDataType::build<ASTNodeDataType::type_id_t>("quadrature");
 
 class SchemeModule : public BuiltinModule
 {
diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp
index 5621a6230ea6898eb361b72b8752f4146c9f6bc4..647298e97ef3ea2333dfeac62a7065234c343d1a 100644
--- a/src/mesh/MeshFlatNodeBoundary.cpp
+++ b/src/mesh/MeshFlatNodeBoundary.cpp
@@ -24,12 +24,12 @@ MeshFlatNodeBoundary<MeshType>::_checkBoundaryIsFlat(const TinyVector<MeshType::
     }
   });
 
-  // if (parallel::allReduceOr(is_bad)) {
-  //   std::ostringstream ost;
-  //   ost << "invalid boundary \"" << rang::fgB::yellow << this->m_ref_node_list.refId() << rang::style::reset
-  //       << "\": boundary is not flat!";
-  //   throw NormalError(ost.str());
-  // }
+  if (parallel::allReduceOr(is_bad)) {
+    std::ostringstream ost;
+    ost << "invalid boundary \"" << rang::fgB::yellow << this->m_ref_node_list.refId() << rang::style::reset
+        << "\": boundary is not flat!";
+    throw NormalError(ost.str());
+  }
 }
 
 template <>
diff --git a/src/scheme/HyperplasticSolverO2.cpp b/src/scheme/HyperplasticSolverO2.cpp
index c4efdd128afab12389b06547370e8c04221fbe7f..9a9eff713b6bebc6b722b87c71b35bb947db25dd 100644
--- a/src/scheme/HyperplasticSolverO2.cpp
+++ b/src/scheme/HyperplasticSolverO2.cpp
@@ -614,6 +614,29 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
     this->_applyVelocityBC(bc_list, Ar, br);
   }
 
+  void
+  _projectOnBoundary(const BoundaryConditionList& bc_list, NodeValue<Rd>& xr) const
+  {
+    for (const auto& boundary_condition : bc_list) {
+      std::visit(
+        [&](auto&& bc) {
+          using T = std::decay_t<decltype(bc)>;
+          if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
+            const Rd& n = bc.outgoingNormal();
+            const Rd& o = bc.origin();
+
+            const Array<const NodeId>& node_list = bc.nodeList();
+            parallel_for(
+              bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
+                const NodeId r = node_list[r_number];
+                Rd& x          = xr[node_list[r]];
+                x -= dot(x - o, n) * n;
+              });
+          }
+        },
+        boundary_condition);
+    }
+  }
   NodeValue<const Rd>
   _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
   {
@@ -740,7 +763,8 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
                const DiscreteScalarFunction& yield,
                const DiscreteTensorFunction3d& sigma,
                const NodeValue<const Rd>& ur,
-               const NodeValuePerCell<const Rd>& Fjr) const
+               const NodeValuePerCell<const Rd>& Fjr,
+               const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
     const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
@@ -752,7 +776,8 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
     NodeValue<Rd> new_xr = copy(mesh.xr());
     parallel_for(
       mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
-
+    const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
+    this->_projectOnBoundary(bc_list, new_xr);
     std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
 
     CellValue<const double> Vj           = MeshDataManager::instance().getMeshData(mesh).Vj();
@@ -858,7 +883,8 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
                 const DiscreteScalarFunction& yield,
                 const DiscreteTensorFunction3d& sigma,
                 const NodeValue<const Rd>& ur,
-                const NodeValuePerCell<const Rd>& Fjr) const
+                const NodeValuePerCell<const Rd>& Fjr,
+                const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
     const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
@@ -870,6 +896,8 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
     NodeValue<Rd> new_xr = copy(mesh.xr());
     parallel_for(
       mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
+    const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
+    this->_projectOnBoundary(bc_list, new_xr);
 
     std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
 
@@ -948,7 +976,8 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
                  const DiscreteScalarFunction& yield,
                  const DiscreteTensorFunction3d& sigma,
                  const NodeValue<const Rd>& ur,
-                 const NodeValuePerCell<const Rd>& Fjr) const
+                 const NodeValuePerCell<const Rd>& Fjr,
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
     const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
@@ -960,6 +989,8 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
     NodeValue<Rd> new_xr = copy(mesh.xr());
     parallel_for(
       mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
+    const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
+    this->_projectOnBoundary(bc_list, new_xr);
 
     std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
 
@@ -1036,7 +1067,8 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
                const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
                const std::shared_ptr<const ItemValueVariant>& ur,
                const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
-               const RelaxationType& relaxation_type) const
+               const RelaxationType& relaxation_type,
+               const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
     std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
     if (not mesh_v) {
@@ -1058,7 +1090,7 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
                                 yield->get<DiscreteScalarFunction>(),     //
                                 sigma->get<DiscreteTensorFunction3d>(),   //
                                 ur->get<NodeValue<const Rd>>(),           //
-                                Fjr->get<NodeValuePerCell<const Rd>>());
+                                Fjr->get<NodeValuePerCell<const Rd>>(), bc_descriptor_list);
     }
     case RelaxationType::Implicit: {
       return this->apply_fluxes2(dt,                                       //
@@ -1071,7 +1103,7 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
                                  yield->get<DiscreteScalarFunction>(),     //
                                  sigma->get<DiscreteTensorFunction3d>(),   //
                                  ur->get<NodeValue<const Rd>>(),           //
-                                 Fjr->get<NodeValuePerCell<const Rd>>());
+                                 Fjr->get<NodeValuePerCell<const Rd>>(), bc_descriptor_list);
     }
     case RelaxationType::CauchyGreen: {
       return this->apply_fluxes_B(dt,                                       //
@@ -1084,7 +1116,7 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
                                   yield->get<DiscreteScalarFunction>(),     //
                                   sigma->get<DiscreteTensorFunction3d>(),   //
                                   ur->get<NodeValue<const Rd>>(),           //
-                                  Fjr->get<NodeValuePerCell<const Rd>>());
+                                  Fjr->get<NodeValuePerCell<const Rd>>(), bc_descriptor_list);
     }
     default: {
       throw UnexpectedError("invalid relaxation type");
@@ -1298,7 +1330,7 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
     auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, bc_descriptor_list);
-    return apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr, relaxation_type);
+    return apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr, relaxation_type, bc_descriptor_list);
   }
 
   std::tuple<std::shared_ptr<const MeshVariant>,
@@ -1695,6 +1727,12 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2<MeshType>::SymmetryBound
     return m_mesh_flat_node_boundary.outgoingNormal();
   }
 
+  const Rd&
+  origin() const
+  {
+    return m_mesh_flat_node_boundary.origin();
+  }
+
   size_t
   numberOfNodes() const
   {
diff --git a/src/scheme/HyperplasticSolverO2.hpp b/src/scheme/HyperplasticSolverO2.hpp
index ff89780d141f9aaea04221ca6be1ad3f46497085..063dcc7d66226d12847b250a57da530590090902 100644
--- a/src/scheme/HyperplasticSolverO2.hpp
+++ b/src/scheme/HyperplasticSolverO2.hpp
@@ -60,7 +60,8 @@ class HyperplasticSolverO2Handler
                  const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
                  const std::shared_ptr<const ItemValueVariant>& ur,
                  const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
-                 const RelaxationType& relaxation_type) const = 0;
+                 const RelaxationType& relaxation_type,
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
 
     virtual std::tuple<std::shared_ptr<const MeshVariant>,
                        std::shared_ptr<const DiscreteFunctionVariant>,
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 110fa9c077070c3c51f247e5b5209cda662b00ca..1a2b1e072734154516b01001cb9831bfaff8ea73 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -25,15 +25,6 @@
 #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)
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index 26acafc6f02e6c452ff5a94e4d3c6aaa9056eeec..e63f5d2e1b1fa0559074ef996d4888b317e1960e 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -146,33 +146,6 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
           SmallArray<double> eigenvalues;
 
 #ifdef PUGS_HAS_SLEPC
-          EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);
-
-          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
-          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
-          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
-
-          SmallMatrix<double> P{3};
-          SmallMatrix<double> PT{3};
-          SmallMatrix<double> D{3};
-          D = zero;
-
-          for (size_t i = 0; i < 3; ++i) {
-            for (size_t j = 0; j < 3; ++j) {
-              P(i, j)  = eigenvector_list[j][i];
-              PT(i, j) = eigenvector_list[i][j];
-            }
-            D(i, i) = eigenvalue_list[i];
-          }
-          SmallMatrix<double> eigen{3};
-          eigen = zero;
-
-          SmallMatrix PDPT = P * D * PT;
-          for (size_t i = 0; i < 3; ++i) {
-            for (size_t j = 0; j < 3; ++j) {
-              REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
-            }
-          }
           EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues);
           REQUIRE(eigenvalues[0] == Catch::Approx(-1));
           REQUIRE(eigenvalues[1] == Catch::Approx(-1));
@@ -273,32 +246,6 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
           SmallArray<double> eigenvalues;
 
 #ifdef PUGS_HAS_SLEPC
-          EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, P);
-
-          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
-          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
-          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
-
-          SmallMatrix<double> D{3};
-          D = zero;
-          for (size_t i = 0; i < 3; ++i) {
-            D(i, i) = eigenvalue_list[i];
-          }
-          SmallMatrix PT = transpose(P);
-
-          TinyMatrix<3, 3, double> TinyA;
-          for (size_t i = 0; i < 3; ++i) {
-            for (size_t j = 0; j < 3; ++j) {
-              TinyA(i, j) = S(i, j);
-            }
-          }
-
-          SmallMatrix PDPT = P * D * PT;
-          for (size_t i = 0; i < 3; ++i) {
-            for (size_t j = 0; j < 3; ++j) {
-              REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
-            }
-          }
           EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues);
           REQUIRE(eigenvalues[0] == Catch::Approx(-1));
           REQUIRE(eigenvalues[1] == Catch::Approx(-1));
@@ -687,36 +634,8 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
 #endif   // PUGS_HAS_EIGEN3
         }
       }
-      SECTION("exponential of a matrix")
-      {
-        SmallMatrix<double> expA{};
-        SmallMatrix<double> expS{3};
-
-        expS(0, 0) = 1325.074593930307812;
-        expS(0, 1) = 662.353357244568285;
-        expS(0, 2) = 1324.70671448913637;
-        expS(1, 0) = expS(0, 1);
-        expS(1, 1) = 331.544558063455535;
-        expS(1, 2) = 662.353357244568185;
-        expS(2, 0) = expS(0, 2);
-        expS(2, 1) = expS(1, 2);
-        expS(2, 2) = expS(0, 0);
-
-#ifdef PUGS_HAS_SLEPC
-        EigenvalueSolver{}.computeExpForSymmetricMatrix(A, expA);
-
-        for (size_t i = 0; i < 3; ++i) {
-          for (size_t j = 0; j < 3; ++j) {
-            REQUIRE(expA(i, j) - expS(i, j) == Catch::Approx(0).margin(1E-12));
-          }
-        }
-
-#else    //  PUGS_HAS_SLEPC
-        REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeExpForSymmetricMatrix(A, expA),
-                            "not implemented yet: SLEPc is required to solve eigenvalue problems");
-#endif   // PUGS_HAS_SLEPC
-      }
     }
+
     SECTION("symmetric tiny matrix")
     {
       TinyMatrix<3> TestA{3e10, 2e10, 4e10, 2e10, 0, 2e10, 4e10, 2e10, 3e10};