diff --git a/Applications/CLI/CMakeLists.txt b/Applications/CLI/CMakeLists.txt
index f15fcb737427d8013770e431caf2a3db3901a992..f4e4a94213a15deeb51a0f63eb7f6b1a620a471d 100644
--- a/Applications/CLI/CMakeLists.txt
+++ b/Applications/CLI/CMakeLists.txt
@@ -26,7 +26,7 @@ if(OGS_USE_PYTHON)
     # appropriate message should be presented. The note is kept for the case
     # that the automatic detection does not work due to whatever reason.
 
-    add_library(ogs_embedded_python STATIC ogs_embedded_python.cpp)
+    add_library(ogs_embedded_python ogs_embedded_python.cpp)
 
     # Performance warning from
     # https://github.com/pybind/pybind11/blob/master/docs/compiling.rst:
@@ -46,7 +46,7 @@ if(OGS_USE_PYTHON)
     if(BUILD_SHARED_LIBS)
         # Add macro definition, because static libs make special handling necessary
         # s.t. the embedded OpenGeoSys Python module won't be removed by the linker.
-        target_compile_definitions(ogs_embedded_python OGS_BUILD_SHARED_LIBS)
+        target_compile_definitions(ogs_embedded_python PRIVATE OGS_BUILD_SHARED_LIBS)
     endif()
 endif()
 
diff --git a/Applications/CLI/ogs_embedded_python.h b/Applications/CLI/ogs_embedded_python.h
index 758718298ba1cb3e00d47338f58901d01ea82074..c215a5f70dfb820fc66fa389ac332b76d63c933d 100644
--- a/Applications/CLI/ogs_embedded_python.h
+++ b/Applications/CLI/ogs_embedded_python.h
@@ -10,11 +10,12 @@
 #pragma once
 
 #include <pybind11/embed.h>
+#include "BaseLib/ExportSymbol.h"
 
 namespace ApplicationsLib
 {
 /// Sets up an embedded Python interpreter and makes sure that the OpenGeoSys
 /// Python module is not removed by the linker.
-pybind11::scoped_interpreter setupEmbeddedPython();
+OGS_EXPORT_SYMBOL pybind11::scoped_interpreter setupEmbeddedPython();
 
 }  // namespace ApplicationsLib
diff --git a/BaseLib/ExportSymbol.h b/BaseLib/ExportSymbol.h
new file mode 100644
index 0000000000000000000000000000000000000000..63e4681ec489307f964c96533bbdb3778a083b5d
--- /dev/null
+++ b/BaseLib/ExportSymbol.h
@@ -0,0 +1,20 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#ifndef OGS_EXPORT_SYMBOL
+#  if defined(WIN32) || defined(_WIN32)
+#    define OGS_EXPORT_SYMBOL __declspec(dllexport)
+#  else
+#    define OGS_EXPORT_SYMBOL __attribute__((visibility("default")))
+#  endif
+#endif  // defined(OGS_EXPORT_SYMBOL)
+
+
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index ba32dd836d7c5b19c54abe5f20dd98f2686609af..832581f3afee4101640df133cfd494d176891d61 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -37,13 +37,11 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
         NumLib::GlobalMatrixProvider::provider.getMatrix(_A_id);
     auto& rhs = NumLib::GlobalVectorProvider::provider.getVector(
         _rhs_id);
-    auto& x_new =
-        NumLib::GlobalVectorProvider::provider.getVector(
-            _x_new_id);
+    auto& x_new = NumLib::GlobalVectorProvider::provider.getVector(_x_new_id);
 
     bool error_norms_met = false;
 
-    LinAlg::copy(x, x_new);  // set initial guess, TODO save the copy
+    LinAlg::copy(x, x_new);  // set initial guess
 
     _convergence_criterion->preFirstIteration();
 
@@ -51,27 +49,34 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
     for (; iteration <= _maxiter;
          ++iteration, _convergence_criterion->reset())
     {
+        BaseLib::RunTime timer_dirichlet;
+        double time_dirichlet = 0.0;
+
         BaseLib::RunTime time_iteration;
         time_iteration.start();
 
-        sys.preIteration(iteration, x);
+        timer_dirichlet.start();
+        sys.computeKnownSolutions(x_new);
+        sys.applyKnownSolutions(x_new);
+        time_dirichlet += timer_dirichlet.elapsed();
+
+        sys.preIteration(iteration, x_new);
 
         BaseLib::RunTime time_assembly;
         time_assembly.start();
-        sys.assemble(x);
+        sys.assemble(x_new);
         sys.getA(A);
         sys.getRhs(rhs);
         INFO("[time] Assembly took %g s.", time_assembly.elapsed());
 
-        BaseLib::RunTime time_dirichlet;
-        time_dirichlet.start();
-        // Here _x_new has to be used and it has to be equal to x!
+        timer_dirichlet.start();
         sys.applyKnownSolutionsPicard(A, rhs, x_new);
-        INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet.elapsed());
+        time_dirichlet += timer_dirichlet.elapsed();
+        INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet);
 
         if (!sys.isLinear() && _convergence_criterion->hasResidualCheck()) {
             GlobalVector res;
-            LinAlg::matMult(A, x_new, res); // res = A * x_new
+            LinAlg::matMult(A, x_new, res);  // res = A * x_new
             LinAlg::axpy(res, -1.0, rhs);   // res -= rhs
             _convergence_criterion->checkResidual(res);
         }
@@ -90,7 +95,7 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
             if (postIterationCallback)
                 postIterationCallback(iteration, x_new);
 
-            switch(sys.postIteration(x_new))
+            switch (sys.postIteration(x_new))
             {
                 case IterationResult::SUCCESS:
                     // Don't copy here. The old x might still be used further
@@ -109,7 +114,8 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
                     INFO(
                         "Picard: The postIteration() hook decided that this "
                         "iteration has to be repeated.");
-                    continue;  // That throws the iteration result away.
+                    LinAlg::copy(x, x_new);  // throw the iteration result away
+                    continue;
             }
         }
 
@@ -184,7 +190,7 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
     bool error_norms_met = false;
 
     // TODO be more efficient
-    // init _minus_delta_x to the right size and 0.0
+    // init minus_delta_x to the right size
     LinAlg::copy(x, minus_delta_x);
 
     _convergence_criterion->preFirstIteration();
@@ -193,9 +199,17 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
     for (; iteration <= _maxiter;
          ++iteration, _convergence_criterion->reset())
     {
+        BaseLib::RunTime timer_dirichlet;
+        double time_dirichlet = 0.0;
+
         BaseLib::RunTime time_iteration;
         time_iteration.start();
 
+        timer_dirichlet.start();
+        sys.computeKnownSolutions(x);
+        sys.applyKnownSolutions(x);
+        time_dirichlet += timer_dirichlet.elapsed();
+
         sys.preIteration(iteration, x);
 
         BaseLib::RunTime time_assembly;
@@ -207,10 +221,10 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
 
         minus_delta_x.setZero();
 
-        BaseLib::RunTime time_dirichlet;
-        time_dirichlet.start();
-        sys.applyKnownSolutionsNewton(J, res, minus_delta_x, x);
-        INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet.elapsed());
+        timer_dirichlet.start();
+        sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
+        time_dirichlet += timer_dirichlet.elapsed();
+        INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet);
 
         if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
             _convergence_criterion->checkResidual(res);
@@ -255,8 +269,6 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
                     continue;  // That throws the iteration result away.
             }
 
-            // TODO could be done via swap. Note: that also requires swapping
-            // the ids. Same for the Picard scheme.
             LinAlg::copy(x_new, x);  // copy new solution to x
             NumLib::GlobalVectorProvider::provider.releaseVector(x_new);
         }
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index d43582d38da3d619e30f1d2a47009a2f5408f708..df992ffe87dd1ef25cd967804c40a51628cc9a1e 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -179,7 +179,7 @@ private:
     std::size_t _A_id = 0u;      //!< ID of the \f$ A \f$ matrix.
     std::size_t _rhs_id = 0u;    //!< ID of the right-hand side vector.
     std::size_t _x_new_id = 0u;  //!< ID of the vector storing the solution of
-                                 //!the linearized equation.
+                                 //! the linearized equation.
 };
 
 /*! Creates a new nonlinear solver from the given configuration.
diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
index b3bf35c59a6a99be5ab17ef1d0f69226b62f9efc..8f78d33151712272f2fe718186e5774627ebbbee 100644
--- a/NumLib/ODESolver/NonlinearSystem.h
+++ b/NumLib/ODESolver/NonlinearSystem.h
@@ -53,14 +53,19 @@ public:
      */
     virtual void getJacobian(GlobalMatrix& Jac) const = 0;
 
+    //! Pre-compute known solutions and possibly store them internally.
+    virtual void computeKnownSolutions(GlobalVector const& x) = 0;
+
     //! Apply known solutions to the solution vector \c x.
+    //! \pre computeKnownSolutions() must have been called before.
     virtual void applyKnownSolutions(GlobalVector& x) const = 0;
 
     //! Apply known solutions to the linearized equation system
     //! \f$ \mathit{Jac} \cdot (-\Delta x) = \mathit{res} \f$.
-    virtual void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
-                                           GlobalVector& minus_delta_x,
-                                           GlobalVector& x) = 0;
+    //! \pre computeKnownSolutions() must have been called before.
+    virtual void applyKnownSolutionsNewton(
+        GlobalMatrix& Jac, GlobalVector& res,
+        GlobalVector& minus_delta_x) const = 0;
 };
 
 /*! A System of nonlinear equations to be solved with the Picard fixpoint
@@ -86,13 +91,18 @@ public:
     //! \pre assemble() must have been called before.
     virtual void getRhs(GlobalVector& rhs) const = 0;
 
+    //! Pre-compute known solutions and possibly store them internally.
+    virtual void computeKnownSolutions(GlobalVector const& x) = 0;
+
     //! Apply known solutions to the solution vector \c x.
+    //! \pre computeKnownSolutions() must have been called before.
     virtual void applyKnownSolutions(GlobalVector& x) const = 0;
 
     //! Apply known solutions to the linearized equation system
     //! \f$ A \cdot x = \mathit{rhs} \f$.
+    //! \pre computeKnownSolutions() must have been called before.
     virtual void applyKnownSolutionsPicard(GlobalMatrix& A, GlobalVector& rhs,
-                                           GlobalVector& x) = 0;
+                                           GlobalVector& x) const = 0;
 };
 
 //! @}
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index 617c722a65fe4f41124d7b7712fcf457b35bb1e7..0fb5bd62d0a57942f1377a7c62150fee9d595f1a 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -120,28 +120,31 @@ void TimeDiscretizedODESystem<
     _mat_trans->computeJacobian(*_Jac, Jac);
 }
 
+void TimeDiscretizedODESystem<
+    ODESystemTag::FirstOrderImplicitQuasilinear,
+    NonlinearSolverTag::Newton>::computeKnownSolutions(GlobalVector const& x)
+{
+    _known_solutions = _ode.getKnownSolutions(_time_disc.getCurrentTime(), x);
+}
+
 void TimeDiscretizedODESystem<
     ODESystemTag::FirstOrderImplicitQuasilinear,
     NonlinearSolverTag::Newton>::applyKnownSolutions(GlobalVector& x) const
 {
-    ::detail::applyKnownSolutions(
-        _ode.getKnownSolutions(_time_disc.getCurrentTime(), x), x);
+    ::detail::applyKnownSolutions(_known_solutions, x);
 }
 
 void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
                               NonlinearSolverTag::Newton>::
     applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
-                              GlobalVector& minus_delta_x, GlobalVector& x)
+                              GlobalVector& minus_delta_x) const
 {
-    auto const* known_solutions =
-        _ode.getKnownSolutions(_time_disc.getCurrentTime(), x);
-
-    if (!known_solutions || known_solutions->empty())
+    if (!_known_solutions)
         return;
 
     using IndexType = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
     std::vector<IndexType> ids;
-    for (auto const& bc : *known_solutions)
+    for (auto const& bc : *_known_solutions)
     {
         std::copy(bc.ids.cbegin(), bc.ids.cend(), std::back_inserter(ids));
     }
@@ -149,8 +152,6 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
     // For the Newton method the values must be zero
     std::vector<double> values(ids.size(), 0);
     MathLib::applyKnownSolution(Jac, res, minus_delta_x, ids, values);
-
-    ::detail::applyKnownSolutions(known_solutions, x);
 }
 
 TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
@@ -201,34 +202,37 @@ void TimeDiscretizedODESystem<
 
 void TimeDiscretizedODESystem<
     ODESystemTag::FirstOrderImplicitQuasilinear,
-    NonlinearSolverTag::Picard>::applyKnownSolutions(GlobalVector& x) const
+    NonlinearSolverTag::Picard>::computeKnownSolutions(GlobalVector const& x)
 {
-    ::detail::applyKnownSolutions(
-        _ode.getKnownSolutions(_time_disc.getCurrentTime(), x), x);
+    _known_solutions = _ode.getKnownSolutions(_time_disc.getCurrentTime(), x);
 }
 
 void TimeDiscretizedODESystem<
     ODESystemTag::FirstOrderImplicitQuasilinear,
-    NonlinearSolverTag::Picard>::applyKnownSolutionsPicard(GlobalMatrix& A,
-                                                           GlobalVector& rhs,
-                                                           GlobalVector& x)
+    NonlinearSolverTag::Picard>::applyKnownSolutions(GlobalVector& x) const
 {
-    auto const* known_solutions =
-        _ode.getKnownSolutions(_time_disc.getCurrentTime(), x);
+    ::detail::applyKnownSolutions(_known_solutions, x);
+}
+
+void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
+                              NonlinearSolverTag::Picard>::
+    applyKnownSolutionsPicard(GlobalMatrix& A,
+                              GlobalVector& rhs,
+                              GlobalVector& x) const
+{
+    if (!_known_solutions)
+        return;
 
-    if (known_solutions)
+    using IndexType = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
+    std::vector<IndexType> ids;
+    std::vector<double> values;
+    for (auto const& bc : *_known_solutions)
     {
-        using IndexType = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
-        std::vector<IndexType> ids;
-        std::vector<double> values;
-        for (auto const& bc : *known_solutions)
-        {
-            std::copy(bc.ids.cbegin(), bc.ids.cend(), std::back_inserter(ids));
-            std::copy(bc.values.cbegin(), bc.values.cend(),
-                      std::back_inserter(values));
-        }
-        MathLib::applyKnownSolution(A, rhs, x, ids, values);
+        std::copy(bc.ids.cbegin(), bc.ids.cend(), std::back_inserter(ids));
+        std::copy(bc.values.cbegin(), bc.values.cend(),
+                  std::back_inserter(values));
     }
+    MathLib::applyKnownSolution(A, rhs, x, ids, values);
 }
 
 }  // namespace NumLib
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index 59ead507782483e9247ee37c7fd2d91f13e20287..498bb8ef808e90072ffbaee6455ef1053d00d60c 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -89,11 +89,12 @@ public:
 
     void getJacobian(GlobalMatrix& Jac) const override;
 
+    void computeKnownSolutions(GlobalVector const& x) override;
+
     void applyKnownSolutions(GlobalVector& x) const override;
 
     void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
-                                   GlobalVector& minus_delta_x,
-                                   GlobalVector& x) override;
+                                   GlobalVector& minus_delta_x) const override;
 
     bool isLinear() const override
     {
@@ -129,6 +130,10 @@ private:
     //! the object used to compute the matrix/vector for the nonlinear solver
     std::unique_ptr<MatTrans> _mat_trans;
 
+    using Index = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
+    std::vector<NumLib::IndexValueVector<Index>> const* _known_solutions =
+        nullptr;  //!< stores precomputed values for known solutions
+
     GlobalMatrix* _Jac;  //!< the Jacobian of the residual
     GlobalMatrix* _M;    //!< Matrix \f$ M \f$.
     GlobalMatrix* _K;    //!< Matrix \f$ K \f$.
@@ -185,10 +190,12 @@ public:
         _mat_trans->computeRhs(*_M, *_K, *_b, rhs);
     }
 
+    void computeKnownSolutions(GlobalVector const& x) override;
+
     void applyKnownSolutions(GlobalVector& x) const override;
 
     void applyKnownSolutionsPicard(GlobalMatrix& A, GlobalVector& rhs,
-                                   GlobalVector& x) override;
+                                   GlobalVector& x) const override;
 
     bool isLinear() const override
     {
@@ -224,6 +231,10 @@ private:
     //! the object used to compute the matrix/vector for the nonlinear solver
     std::unique_ptr<MatTrans> _mat_trans;
 
+    using Index = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
+    std::vector<NumLib::IndexValueVector<Index>> const* _known_solutions =
+        nullptr;  //!< stores precomputed values for known solutions
+
     GlobalMatrix* _M;  //!< Matrix \f$ M \f$.
     GlobalMatrix* _K;  //!< Matrix \f$ K \f$.
     GlobalVector* _b;  //!< Matrix \f$ b \f$.
diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp
index b7e640956248ee3a06ffc426d4aeec1e031c99fe..3df5a5473aebaf702bc91f4a713cdc4aa9d1effa 100644
--- a/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp
@@ -235,8 +235,9 @@ std::unique_ptr<PythonBoundaryCondition> createPythonBoundaryCondition(
     // parameters are not read and will cause an error.
     // TODO (naumov): Add a function to ConfigTree for skipping the tags of the
     // subtree and move the code up in createBoundaryCondition().
-    if (bc_mesh.getDimension() == 0 && bc_mesh.getNumberOfNodes() == 0 &&
-        bc_mesh.getNumberOfElements() == 0)
+    if (boundary_mesh.getDimension() == 0 &&
+        boundary_mesh.getNumberOfNodes() == 0 &&
+        boundary_mesh.getNumberOfElements() == 0)
     {
         return nullptr;
     }
diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionModule.h b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionModule.h
index 82c141a4bc3aa73734b85f315aeb0576874b6cea..c52f4f9e7ca1f97a59892a7483b476258b7db1ec 100644
--- a/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionModule.h
+++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionModule.h
@@ -10,9 +10,10 @@
 #pragma once
 
 #include <pybind11/pybind11.h>
+#include "BaseLib/ExportSymbol.h"
 
 namespace ProcessLib
 {
 //! Creates Python bindings for the Python BC class.
-void pythonBindBoundaryCondition(pybind11::module& m);
+OGS_EXPORT_SYMBOL void pythonBindBoundaryCondition(pybind11::module& m);
 }  // namespace ProcessLib
diff --git a/ProcessLib/GroundwaterFlow/Tests.cmake b/ProcessLib/GroundwaterFlow/Tests.cmake
index 5e113dd90cf600ef6e57aeedec6add6d0de76fda..77a293d5b08e4b6bee026178ce00641f903d07bc 100644
--- a/ProcessLib/GroundwaterFlow/Tests.cmake
+++ b/ProcessLib/GroundwaterFlow/Tests.cmake
@@ -682,3 +682,15 @@ AddTest(
     VIS square_1e6__nodal_sources_expected_pcs_0_ts_1_t_1.000000.vtu
 )
 
+AddTest(
+    NAME PythonBCGroundWaterFlowProcessLaplaceEqDirichletNeumann
+    PATH Elliptic/square_1x1_GroundWaterFlow_Python
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS square_1e3_laplace_eq.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS OGS_USE_PYTHON AND NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    DIFF_DATA
+    python_laplace_eq_ref.vtu square_1e3_neumann_pcs_0_ts_1_t_1.000000.vtu pressure_expected pressure 4e-4 1e-16
+)
+
diff --git a/ProcessLib/SmallDeformation/Tests.cmake b/ProcessLib/SmallDeformation/Tests.cmake
index 704d03ba4abcd8d6a29ac9e87d4bf97d304ec2de..1787f24f1e42cd063841776e3d7a8858c2ab471a 100644
--- a/ProcessLib/SmallDeformation/Tests.cmake
+++ b/ProcessLib/SmallDeformation/Tests.cmake
@@ -316,3 +316,20 @@ AddTest(
     ref_piston_pcs_0_ts_10_t_10.000000.vtu piston_pcs_0_ts_10_t_10.000000.vtu epsilon epsilon 1e-13 0
     ref_piston_pcs_0_ts_10_t_10.000000.vtu piston_pcs_0_ts_10_t_10.000000.vtu sigma sigma 1e-7 0
 )
+
+AddTest(
+    NAME LARGE_PythonBCSmallDeformationHertzContact
+    PATH Mechanics/Linear/PythonHertzContact
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS hertz_contact.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS OGS_USE_PYTHON AND NOT OGS_USE_MPI
+    DIFF_DATA
+    ref_hertz_contact_ts_5.vtu  hertz_pcs_0_ts_5_t_5.000000.vtu  displacement displacement 1e-16 0
+    ref_hertz_contact_ts_5.vtu  hertz_pcs_0_ts_5_t_5.000000.vtu  epsilon epsilon 1e-16 0
+    ref_hertz_contact_ts_5.vtu  hertz_pcs_0_ts_5_t_5.000000.vtu  sigma sigma 1e-16 0
+    ref_hertz_contact_ts_10.vtu hertz_pcs_0_ts_10_t_10.000000.vtu displacement displacement 1e-16 0
+    ref_hertz_contact_ts_10.vtu hertz_pcs_0_ts_10_t_10.000000.vtu epsilon epsilon 1e-16 0
+    ref_hertz_contact_ts_10.vtu hertz_pcs_0_ts_10_t_10.000000.vtu sigma sigma 1e-16 0
+)
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index 9ece9ab2b24f5e6249147d191702c1e5febdc889..19450c9ea720cbb8a64d332d6de21ec29fd73ae8 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -58,37 +58,6 @@ static void setEquationSystem(NumLib::NonlinearSolverBase& nonlinear_solver,
     }
 }
 
-//! Applies known solutions to the solution vector \c x, transparently
-//! for equation systems linearized with either the Picard or Newton method.
-template <NumLib::NonlinearSolverTag NLTag>
-static void applyKnownSolutions(NumLib::EquationSystem const& eq_sys,
-                                GlobalVector& x)
-{
-    using EqSys = NumLib::NonlinearSystem<NLTag>;
-    assert(dynamic_cast<EqSys const*>(&eq_sys) != nullptr);
-    auto& eq_sys_ = static_cast<EqSys const&>(eq_sys);
-
-    eq_sys_.applyKnownSolutions(x);
-}
-
-//! Applies known solutions to the solution vector \c x, transparently
-//! for equation systems linearized with either the Picard or Newton method.
-static void applyKnownSolutions(NumLib::EquationSystem const& eq_sys,
-                                NumLib::NonlinearSolverTag const nl_tag,
-                                GlobalVector& x)
-{
-    using Tag = NumLib::NonlinearSolverTag;
-    switch (nl_tag)
-    {
-        case Tag::Picard:
-            applyKnownSolutions<Tag::Picard>(eq_sys, x);
-            break;
-        case Tag::Newton:
-            applyKnownSolutions<Tag::Newton>(eq_sys, x);
-            break;
-    }
-}
-
 namespace ProcessLib
 {
 template <NumLib::ODESystemTag ODETag>
@@ -286,8 +255,6 @@ bool solveOneTimeStepOneProcess(int const process_id, GlobalVector& x,
 
     time_disc.nextTimestep(t, delta_t);
 
-    applyKnownSolutions(ode_sys, nl_tag, x);
-
     auto const post_iteration_callback = [&](unsigned iteration,
                                              GlobalVector const& x) {
         output_control.doOutputNonlinearIteration(process, process_id,
diff --git a/Tests/Data/Mechanics/Linear/PythonHertzContact/gen-unit-circle.py b/Tests/Data/Mechanics/Linear/PythonHertzContact/gen-unit-circle.py
new file mode 100755
index 0000000000000000000000000000000000000000..ef4bf9072aea53d301399e2d51cb61ba338548da
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonHertzContact/gen-unit-circle.py
@@ -0,0 +1,124 @@
+#!/usr/bin/vtkpython
+
+import sys
+import numpy as np
+from vtk.util.numpy_support import vtk_to_numpy, numpy_to_vtk
+import vtk
+
+in_grid, out_grid, out_geom = sys.argv[1:]
+
+
+def distribute_points_evenly(c2):
+    assert np.sqrt(c2.shape[0]).is_integer()
+
+    CELLS_PER_DIRECTION = int(np.sqrt(c2.shape[0])) - 1
+    r2 = np.sqrt(c2[:,0]**2 + c2[:,1]**2)
+    alpha2 = np.arctan2(c2[:,1], c2[:,0])
+
+    bins = [ [] for _ in range(CELLS_PER_DIRECTION+1) ]
+    nbin = np.round(r2 * CELLS_PER_DIRECTION).astype(int)
+    for node, b in enumerate(nbin):
+        bins[b].append(node)
+
+    for b in bins:
+        b.sort(key=lambda n: alpha2[n])
+        # print(len(b))
+
+    c3 = np.zeros_like(c2)
+
+    for node, r in enumerate(r2):
+        b = nbin[node]
+        i = bins[b].index(node)
+        if len(bins[b]) == 1:
+            phi = 0.0
+        else:
+            phi = np.pi * 0.5 / (len(bins[b]) - 1) * i
+
+        c3[node, 0] = r * np.cos(phi)
+        c3[node, 1] = r * np.sin(phi)
+
+    return c3
+
+
+reader = vtk.vtkXMLUnstructuredGridReader()
+reader.SetFileName(in_grid)
+reader.Update()
+
+grid = reader.GetOutput()
+assert grid.GetBounds() == (0.0, 1.0, 0.0, 1.0, 0.0, 0.0)
+
+coords = vtk_to_numpy(grid.GetPoints().GetData())
+
+a = np.empty(coords.shape[0])
+
+for i, (x, y, z) in enumerate(coords):
+    if x > y:
+        R = np.sqrt(1 + (y/x)**2)
+        a[i] = 1.0 / R
+    elif x < y:
+        R = np.sqrt(1 + (x/y)**2)
+        a[i] = 1.0 / R
+    else:
+        a[i] = np.sqrt(0.5)
+
+# scale coordinates
+new_coords = np.multiply(coords.T, a).T
+
+if False:
+    # If the input is a regular mesh, there is the possibility to make bins of
+    # points with equal radius. For each radius those points can than be
+    # distributed with even spacing at the quarter circles.
+    new_coords = distribute_points_evenly(new_coords)
+
+new_coords_vtk = numpy_to_vtk(new_coords, 1)
+pts = vtk.vtkPoints()
+pts.SetData(new_coords_vtk)
+grid.SetPoints(pts)
+
+writer = vtk.vtkXMLUnstructuredGridWriter()
+writer.SetFileName(out_grid)
+writer.SetInputData(grid)
+writer.Write()
+
+# extract boundary of unit circle
+R_squared = new_coords[:,0]**2 + new_coords[:,1]**2
+phi = np.arctan2(new_coords[:,1], new_coords[:,0])
+
+idcs = np.where(abs(R_squared - 1) < 1e-8)[0]
+idcs = idcs[np.argsort(phi[idcs])]  # sorted with ascending polar angle
+
+
+with open(out_geom, "w") as fh:
+    fh.write("""<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysGLI>
+    <name>geom</name>
+    <points>
+        <point id="0" x="0" y="0" z="0" name="center"/>
+        <point id="1" x="0" y="1" z="0" name="top"/>
+        <point id="2" x="1" y="0" z="0"/>
+""")
+
+    for i, (x, y, z) in enumerate(new_coords[idcs]):
+        fh.write('        <point id="{}" x="{}" y="{}" z="{}" />\n'.format(
+            i+3, x, y, z))
+
+    fh.write("""
+    </points>
+
+    <polylines>
+        <polyline id="0" name="left">
+            <pnt>0</pnt>
+            <pnt>1</pnt>
+        </polyline>
+        <polyline id="1" name="bottom">
+            <pnt>0</pnt>
+            <pnt>2</pnt>
+        </polyline>
+        <polyline id="2" name="outer">\n""")
+
+    for i in range(len(idcs)):
+        fh.write("            <pnt>{}</pnt>\n".format(i+3))
+
+    fh.write("""</polyline>
+    </polylines>
+</OpenGeoSysGLI>""")
diff --git a/Tests/Data/Mechanics/Linear/PythonHertzContact/hertz_contact.prj b/Tests/Data/Mechanics/Linear/PythonHertzContact/hertz_contact.prj
new file mode 100644
index 0000000000000000000000000000000000000000..8cbeeb41991bf27ac1c70922025bff33b701ec23
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonHertzContact/hertz_contact.prj
@@ -0,0 +1,154 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <mesh axially_symmetric="true">unit-circle.vtu</mesh>
+    <geometry>unit-circle.gml</geometry>
+    <python_script>hertz_contact_bc.py</python_script>
+    <processes>
+        <process>
+            <name>SD</name>
+            <type>SMALL_DEFORMATION</type>
+            <integration_order>2</integration_order>
+            <constitutive_relation>
+                <type>LinearElasticIsotropic</type>
+                <youngs_modulus>E</youngs_modulus>
+                <poissons_ratio>nu</poissons_ratio>
+            </constitutive_relation>
+            <solid_density>rho_sr</solid_density>
+            <specific_body_force>0 0</specific_body_force>
+            <process_variables>
+                <process_variable>displacement</process_variable>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+                <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+            </secondary_variables>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="SD">
+                <nonlinear_solver>basic_newton</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>INFINITY_N</norm_type>
+                    <abstol>1e-10</abstol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <output>
+                    <variables>
+                        <variable>displacement</variable>
+                        <variable>sigma</variable>
+                        <variable>epsilon</variable>
+                    </variables>
+                </output>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0</t_initial>
+                    <t_end>10</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>1</repeat>
+                            <delta_t>1</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>hertz</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>1</repeat>
+                    <each_steps>1</each_steps>
+                </pair>
+            </timesteps>
+            <output_iteration_results>false</output_iteration_results>
+        </output>
+    </time_loop>
+    <parameters>
+        <parameter>
+            <name>E</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>nu</name>
+            <type>Constant</type>
+            <value>.3</value>
+        </parameter>
+        <parameter>
+            <name>rho_sr</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>displacement0</name>
+            <type>Constant</type>
+            <values>0 0</values>
+        </parameter>
+        <parameter>
+            <name>zero</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>2</components>
+            <order>1</order>
+            <initial_condition>displacement0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>geom</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>zero</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>geom</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>zero</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>geom</geometrical_set>
+                    <geometry>outer</geometry>
+                    <type>Python</type>
+                    <component>1</component>
+                    <bc_object>bc</bc_object>
+                    <flush_stdout>false</flush_stdout> <!-- for debugging only -->
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_newton</name>
+            <type>Newton</type>
+            <max_iter>50</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+            <eigen>
+                <solver_type>CG</solver_type>
+                <precon_type>DIAGONAL</precon_type>
+                <max_iteration_step>10000</max_iteration_step>
+                <error_tolerance>1e-16</error_tolerance>
+            </eigen>
+            <petsc>
+                <prefix>sd</prefix>
+                <parameters>-sd_ksp_type cg -sd_pc_type bjacobi -sd_ksp_rtol 1e-16 -sd_ksp_max_it 10000</parameters>
+            </petsc>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Mechanics/Linear/PythonHertzContact/hertz_contact_bc.py b/Tests/Data/Mechanics/Linear/PythonHertzContact/hertz_contact_bc.py
new file mode 100644
index 0000000000000000000000000000000000000000..4abfa50b91b36e48b425cadf5606e150420e5e29
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonHertzContact/hertz_contact_bc.py
@@ -0,0 +1,150 @@
+from __future__ import print_function
+
+import OpenGeoSys
+
+
+SPHERE_RADIUS = 1.0
+START_TIME = 0.0
+
+
+class HertzContactBC(OpenGeoSys.BoundaryCondition):
+    def __init__(self):
+        super(HertzContactBC, self).__init__()
+
+        self._first_node = None         # ID of the first node of this BC's geometry
+        self._t_old = START_TIME - 1.0  # time of previous invocation of this BC
+
+        self._boundary_x_coords = []    # the x coordinates of all boundary nodes
+
+
+    def _init_timestep(self):
+        """Initializes the internal state at the beginning of a new timestep."""
+
+        self._a_range = [ 0.0, SPHERE_RADIUS ]  # range of possible contact radii
+        self._a_est = SPHERE_RADIUS             # estimated contact radius
+
+        self._max_x_with_y_excess = -1.0        # maximum x value where a node is above the contact line
+
+        self._tendency_is_up = True             # whether the contact area is supposed to grow in the current iteration
+        self._a_curr = 0.0                      # the radius of the contact area in this iteration (continuously updated)
+        self._a_prev = 0.0                      # the radius of the contact area that was prescripted in the previous iteration
+
+        self._iteration = 0                     # the nonlinear solver iteration number
+
+
+    def _init_iteration(self):
+        """Initializes the internal state at the beginning of a new nonlinear solver iteration."""
+
+        # variant of a bisection algorithm
+        if self._max_x_with_y_excess >= 0.0:
+            if self._max_x_with_y_excess < self._a_est:
+                # This is an ad-hoc variation of the bisection algorithm for
+                # finding the radius of the contact area. The variation is
+                # necessary, because we obtain the value of
+                # self._max_x_with_y_excess only after we already applied the
+                # BCs for the current nonlinear solver iteration, i.e., the
+                # applied BCs and the information based on which the BCs are
+                # applied are off by one iteration.
+                self._a_range[1] = self._max_x_with_y_excess
+                assert self._a_range[0] <= self._a_range[1]
+                self._tendency_is_up = False
+            else:
+                # there were nodes above the contact area in the previous
+                # iteration => contact area has to grow.
+                self._a_range[0] = self._a_est
+                self._tendency_is_up = True
+        else:
+            # contact area has to shrink
+            self._a_range[1] = self._a_est
+            self._tendency_is_up = False
+
+        self._a_est = 0.5 * sum(self._a_range)
+
+        self._max_x_with_y_excess = -1.0
+        self._a_prev = self._a_curr
+        self._a_curr = -1.0
+
+        self._iteration += 1
+
+        print("BC: a_est={:.4f}, a_prev={:.4f} ({:.4f}, {:.4f})".format(
+            self._a_est, self._a_prev, self._a_range[0], self._a_range[1]))
+
+
+    def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+        if self._t_old < t:
+            self._t_old = t
+            self._init_timestep()
+
+        # detect nonlinear iteration, assumes that nodes are always processed in
+        # the same order
+        if self._first_node is None:
+            self._first_node = node_id
+
+        if self._first_node == node_id:
+            self._init_iteration()
+            if self._iteration == 2:
+                self._boundary_x_coords.sort()
+
+        x, y, z = coords
+        ux, uy = primary_vars
+
+        if not self._boundary_x_coords_are_initialized(t):
+            self._boundary_x_coords.append(x)
+
+        try:
+            # check that we are at the outer boundary
+            assert abs(x**2 + y**2 + z**2 - SPHERE_RADIUS**2) < 1e-15
+        except:
+            print("assert abs(x**2 + y**2 + z**2 - 1.0) < 1e-15",
+                    x, y, z, x**2 + y**2 + z**2 - SPHERE_RADIUS**2)
+            raise
+
+        y_deformed = y + uy
+        y_top = HertzContactBC._get_y_top(t)
+
+        res = (False, 0.0)
+
+        if y_deformed >= y_top:
+            self._max_x_with_y_excess = max(self._max_x_with_y_excess, x)
+
+            if x <= self._a_est:
+                res = (True, y_top - y)
+                self._a_curr = max(x, self._a_curr)
+
+            elif self._tendency_is_up:
+                # This branch catches some corner-cases where Dirichlet BCs are
+                # set on too few nodes.
+                #
+                # The existence of these corner cases and the fact that the
+                # bisection algorithm for finding the correct contact radius is
+                # rather ad-hoc show that in general for similar problems it
+                # will be necessary to have more flexible convergence criteria
+                # or nonlinear solvers.  In particular said similar problems are
+                # those where the surfaces where Dirichlet or Neumann BCs are
+                # set change for changing primary variables.
+
+                if x <= self._a_prev:
+                    assert False  # this case shouldn't happen
+                    res = (True, y_top - y)
+                elif self._boundary_x_coords_are_initialized(t):
+                    idx = self._boundary_x_coords.index(x)
+                    if idx != 0 and self._boundary_x_coords[idx-1] == self._a_prev:
+                        res = (True, y_top - y)
+
+
+        return res
+
+
+    @staticmethod
+    def _get_y_top(t):
+        """Returns the y-position of the contact area depending on the load step t."""
+
+        return 1.0 - 0.005 * t
+
+
+    def _boundary_x_coords_are_initialized(self, t):
+        return self._iteration >= 2
+
+
+# instantiate the BC object used by OpenGeoSys
+bc = HertzContactBC()
diff --git a/Tests/Data/Mechanics/Linear/PythonHertzContact/post.py b/Tests/Data/Mechanics/Linear/PythonHertzContact/post.py
new file mode 100755
index 0000000000000000000000000000000000000000..8dc81bf1acb607fce9a338edc6f4bb267aff5f07
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonHertzContact/post.py
@@ -0,0 +1,393 @@
+#!/usr/bin/vtkpython
+
+from __future__ import print_function
+import vtk
+from vtk.util.numpy_support import vtk_to_numpy, numpy_to_vtk
+import numpy as np
+from scipy.interpolate import interp1d
+
+import matplotlib.pyplot as plt
+
+pvd_file = "out/hertz_pcs_0.pvd"
+
+
+R12 = 1.0
+R = R12 / 2.0
+nu_12 = 0.3
+E_12 = 1.0
+
+lambda_ = E_12 * nu_12 / (1 + nu_12) / (1 - 2 * nu_12)
+mu = E_12 / 2.0 / (1 + nu_12)
+G_12 = mu
+
+kappa = 0.5 * G_12 / R / (1-nu_12)
+print("kappa:", kappa)
+
+C = lambda_ * np.matrix([
+    1, 1, 1, 0,
+    1, 1, 1, 0,
+    1, 1, 1, 0,
+    0, 0, 0, 0 ]).reshape(4, 4) \
+            + 2 * mu * np.identity(4)
+
+def p_contact(r, r_contact):
+    return kappa * np.sqrt(r_contact**2 - r**2)
+
+
+### helpers ##############################################
+
+import os
+try:
+    import xml.etree.cElementTree as ET
+except:
+    import xml.etree.ElementTree as ET
+
+
+def relpathfrom(origin, relpath):
+    if os.path.isabs(relpath):
+        return relpath
+    return os.path.join(origin, relpath)
+
+def read_pvd_file(fn):
+    try:
+        path = fn.name
+    except AttributeError:
+        path = fn
+    pathroot = os.path.dirname(path)
+    pvdtree = ET.parse(fn)
+    node = pvdtree.getroot()
+    if node.tag != "VTKFile": return None, None
+    children = list(node)
+    if len(children) != 1: return None, None
+    node = children[0]
+    if node.tag != "Collection": return None, None
+
+    ts = []
+    fs = []
+
+    for child in node:
+        if child.tag != "DataSet": return None, None
+        ts.append(float(child.get("timestep")))
+        fs.append(relpathfrom(pathroot, child.get("file")))
+
+    return ts, fs
+
+### helpers end ##########################################
+
+
+def get_y_top(t):
+    return 1.0 - 0.005 * t
+
+ts, fns = read_pvd_file(pvd_file)
+
+reader = vtk.vtkXMLUnstructuredGridReader()
+
+destroyTopology = vtk.vtkShrinkFilter()
+destroyTopology.SetShrinkFactor(1.0)
+# destroyTopology.SetInputConnection(strainFilter.GetOutputPort())
+# destroyTopology.SetInputConnection(reader.GetOutputPort())
+
+# strainFilter = vtk.vtkCellDerivatives()
+# strainFilter.SetVectorModeToPassVectors()
+# strainFilter.SetTensorModeToComputeStrain()
+# strainFilter.SetOutputPointsPrecision(vtk.vtkAlgorithm.DOUBLE_PRECISION)
+
+warpVector = vtk.vtkWarpVector()
+# warpVector.SetInputConnection(strainFilter.GetOutputPort())
+
+# cell2point = vtk.vtkCellDataToPointData()
+# cell2point.SetInputConnection(destroyTopology.GetOutputPort())
+# cell2point.PassCellDataOff()
+# cell2point.SetOutputPointsPrecision(vtk.vtkAlgorithm.DOUBLE_PRECISION)
+
+plane = vtk.vtkPlane()
+# plane.SetOrigin(0, 0, 0)
+plane.SetNormal(0, 1, 0)
+
+cutter = vtk.vtkCutter()
+cutter.SetCutFunction(plane)
+cutter.SetInputConnection(warpVector.GetOutputPort())
+cutter.SetOutputPointsPrecision(vtk.vtkAlgorithm.DOUBLE_PRECISION)
+
+writer = vtk.vtkXMLUnstructuredGridWriter()
+
+ws = []
+rs_contact = []
+Fs = []
+
+fig, ax = plt.subplots()
+
+fig.subplots_adjust(right=0.75)
+ax.set_xlabel(r"$\xi$ / m")
+ax.set_ylabel(r"$\bar p$ / Pa")
+add_leg = True
+
+
+for t, fn in zip(ts, fns):
+    print("###### time", t)
+    reader.SetFileName(fn)
+    reader.Update()
+    grid = reader.GetOutput()
+
+    disp_2d = vtk_to_numpy(grid.GetPointData().GetArray("displacement"))
+    disp_3d = np.zeros((disp_2d.shape[0], 3))
+    disp_3d[:,(0,1)] = disp_2d
+    disp_3d_vtk = numpy_to_vtk(disp_3d, 1)
+    disp_3d_vtk.SetName("u")
+
+    grid.GetPointData().AddArray(disp_3d_vtk)
+    # grid.GetPointData().SetActiveVectors("u")
+
+    if False:
+        # compute strain
+        def strain_triangle_axi(cell, point_data, strain_data):
+            cell_pts = np.matrix(vtk_to_numpy(cell.GetPoints().GetData())[:,:-1])
+            assert cell_pts.shape[0] == 3  # triangles
+            assert point_data.shape[1] == 2  # 2D vector field
+            node_ids = [ cell.GetPointId(i) for i in range(cell.GetNumberOfPoints()) ]
+
+            # interpolation using barycentric coordinates on linear triangles
+            T = np.matrix(np.empty((2,2)))
+            T[:,0] = (cell_pts[0,:] - cell_pts[2,:]).T
+            T[:,1] = (cell_pts[1,:] - cell_pts[2,:]).T
+            T_inv = np.linalg.inv(T)
+
+            dl1 = T_inv[0,:]  # row 0
+            dl2 = T_inv[1,:]  # row 1
+
+            for node in range(3):
+                l1, l2 = T_inv * (cell_pts[node, :].T - cell_pts[2, :].T)
+                assert -1e-15 < l1 and 1 + 1e-15 > l1 \
+                        and -1e-15 < l2 and 1+ 1e-15 > l2
+
+            grad = np.empty((2,2))
+            for comp in range(2):
+                nodal_values = point_data[node_ids, comp]
+                # nodal_values = cell_pts[:, comp].flat
+                # if t > 0 and cell_pts[0,1] > 0.95 and comp == 1:
+                #     print(nodal_values[0])
+                grad[comp,:] = dl1 * nodal_values[0] \
+                        + dl2 * nodal_values[1] \
+                        - (dl1 + dl2) * nodal_values[2]
+
+            # if t > 0 and cell_pts[0,1] > 0.95:
+            #     print(grad)
+
+            strain = 0.5 * (grad + grad.T)  # rr, rz, zr,zz components
+
+            for node in range(3):
+                r = cell_pts[node, 0]
+                node_id = node_ids[node]
+
+                if r == 0:
+                    dvdr = grad[0,0]
+                    v_over_r = dvdr
+                else:
+                    v_over_r = point_data[node_id, 0] / r
+
+                strain_kelvin = np.array([
+                    strain[0,0], strain[1,1], v_over_r,
+                    strain[0,1] * np.sqrt(2.0)
+                    ])
+                strain_data[node_id, :] = strain_kelvin
+
+        def computeStrain(grid):
+            destroyTopology.SetInputData(grid)
+            destroyTopology.Update()
+            grid = destroyTopology.GetOutput()
+
+            disp_2d = vtk_to_numpy(grid.GetPointData().GetArray("displacement"))
+            strain_kelvin = np.empty((disp_2d.shape[0], 4))
+
+            n_cells = grid.GetNumberOfCells()
+            for c in xrange(n_cells):
+                cell = grid.GetCell(c)
+                strain_triangle_axi(cell, disp_2d, strain_kelvin)
+
+            strain_kelvin_vtk = numpy_to_vtk(strain_kelvin, 1)
+            strain_kelvin_vtk.SetName("strain_post_kelvin")
+            grid.GetPointData().AddArray(strain_kelvin_vtk)
+
+            strain = strain_kelvin.copy()
+            strain[:,3] /= np.sqrt(2.0)
+            strain_vtk = numpy_to_vtk(strain, 1)
+            strain_vtk.SetName("strain_post")
+            grid.GetPointData().AddArray(strain_vtk)
+
+            # ( (4 x 4) * (nodes x 4).T ).T
+            stress_kelvin = (C * strain_kelvin.T).T
+            # stress_kelv = np.empty_like(strain_kelv)
+            # for c, eps in enumerate(strain_kelv):
+            #     stress_kelv[c, :] = (C * np.atleast_2d(eps).T).flat
+
+            stress_kelvin_vtk = numpy_to_vtk(stress_kelvin, 1)
+            stress_kelvin_vtk.SetName("stress_post_kelvin")
+
+            stress_symm_tensor = stress_kelvin.copy()
+            stress_symm_tensor[:,3] /= np.sqrt(2.0)
+
+            stress_symm_tensor_vtk = numpy_to_vtk(stress_symm_tensor, 1)
+            stress_symm_tensor_vtk.SetName("stress_post")
+            grid.GetPointData().AddArray(stress_symm_tensor_vtk)
+
+            writer.SetInputData(grid)
+            writer.SetFileName(os.path.join(
+                os.path.dirname(fn), "post_{:.0f}.vtu".format(t)))
+            writer.Write()
+
+            return grid
+
+        grid = computeStrain(grid)
+
+    grid.GetPointData().SetActiveVectors("u")
+    warpVector.SetInputData(grid)
+    warpVector.Update()
+    grid = warpVector.GetOutput()
+
+    xmin, xmax, ymin, ymax, zmin, zmax = grid.GetBounds()
+    y_top = get_y_top(t)
+    try:
+        assert abs(ymax - y_top) < 1e-7
+    except:
+        print(ymax, y_top, ymax - y_top)
+        raise
+    ws.append(2 * (1.0 - y_top))
+
+    # determine top boundary
+    coords = vtk_to_numpy(grid.GetPoints().GetData())
+    assert abs(min(coords[:,0])) < 1e-8
+    idcs_top_boundary = np.where(coords[:,1] > y_top - 1e-7)[0]
+    # print(idcs_top_boundary)
+    assert len(idcs_top_boundary) != 0
+    xs_top_boundary = coords[idcs_top_boundary, 0]
+    idx_max = np.argmax(xs_top_boundary)
+    r_contact = max(xs_top_boundary[idx_max], 0.0)
+    y_at_r_contact = coords[idcs_top_boundary[idx_max], 1]
+    print("radius of contact area:", r_contact, "at y =", y_at_r_contact)
+
+    rs_contact.append(r_contact)
+
+    def average_stress(rs, stress):
+        r_contact = max(rs)
+        rs_int = np.linspace(min(rs), max(rs)-1e-8, max(len(rs), 200))
+        stress_int = interp1d(rs, stress, bounds_error=False, fill_value=0.0)
+        avg_stress = np.empty_like(rs_int)
+
+        for i, r in enumerate(rs_int):
+            rho_max = np.sqrt(r_contact**2 - r**2)
+            rhos = np.linspace(0, rho_max, 100)
+            xis = np.sqrt(rhos**2 + r**2)
+            try:
+                assert max(xis) <= r_contact + 1e-8
+            except:
+                print(max(xis), r_contact)
+                raise
+            avg_stress[i] = 1.0 / rho_max * np.trapz(x=rhos, y=stress_int(xis))
+        # avg_stress[-1] = 0.0
+
+        return rs_int, avg_stress
+
+    def total_force(rs, stress):
+        assert all(rs > -1e-6)
+        rs_int = np.linspace(min(rs), max(rs), max(len(rs), 200))
+        stress_int = interp1d(rs, stress, bounds_error=False, fill_value=0.0)
+
+        F = 2.0 * np.pi * np.trapz(x=rs_int, y=rs_int * stress_int(rs_int))
+
+        return F
+
+
+    def stress_at_contact_area():
+        global add_leg
+
+        plane.SetOrigin(0, y_at_r_contact, 0)
+        cutter.Update()
+        grid = cutter.GetOutput()
+
+        # for a in range(grid.GetPointData().GetNumberOfArrays()):
+        #     print(grid.GetPointData().GetArrayName(a))
+
+        xs = vtk_to_numpy(grid.GetPoints().GetData())[:, 0]
+        try:
+            assert abs(max(xs) - r_contact) < 1e-7
+            assert abs(min(xs)) < 1e-8
+        except:
+            print(min(xs), max(xs), r_contact, max(xs) - r_contact)
+            raise
+
+        w_0 = 2.0 * (1.0 - y_top)
+
+        rs = np.linspace(0, r_contact, 200)
+        if add_leg: ax.plot([], [], color="k", ls=":", label="ref")
+        h, = ax.plot(rs, -p_contact(rs, r_contact), ls=":")
+
+        if False:
+            r_contact_ana = np.sqrt(w_0 * R)
+
+            rs2 = np.linspace(0, r_contact_ana, 200)
+            if add_leg: ax.plot([], [], color="k", ls="--", label="ref2")
+            ax.plot(rs2, -p_contact(rs, r_contact_ana), ls="--", color=h.get_color())
+
+        stress_yy = vtk_to_numpy(grid.GetPointData().GetArray("sigma"))[:,1]
+        rs, avg_stress_yy = average_stress(xs, stress_yy)
+        if add_leg: ax.plot([], [], color="k", ls="-", label="ogs")
+        ax.plot(rs, avg_stress_yy, color=h.get_color(), ls="-",
+                label=r"$w_0 = {}$".format(w_0))
+
+        if False:
+            stress = vtk_to_numpy(grid.GetPointData().GetArray("stress_post"))
+            rs, avg_stress_yy = average_stress(xs, stress[:,1])
+            if add_leg: ax.plot([], [], color="k", label="post")
+            ax.plot(rs, avg_stress_yy, color=h.get_color(),
+                    label=r"$w_0 = {}$".format(2*(1.0 - y_top)))
+
+        ax.scatter([r_contact], [0], color=h.get_color())
+        ax.legend(loc="center left", bbox_to_anchor=(1.0, 0.5))
+
+        fig.savefig("stress_at_contact.png")
+        add_leg = False
+
+        Fs.append(-total_force(xs, stress_yy))
+
+    if t > 0:
+        stress_at_contact_area()
+
+fig.savefig("stress_at_contact.png")
+plt.close(fig)
+
+
+fig, ax = plt.subplots()
+ax.scatter(ws, rs_contact, label="ogs")
+
+ws_ref = np.linspace(0, max(ws), 200)
+rs_ref = np.sqrt(ws_ref * R)
+
+ax.plot(ws_ref, rs_ref, label="ref")
+
+ax.legend()
+
+ax.set_xlabel("displacement of sphere centers $w_0$ / m")
+ax.set_ylabel("radius of contact area $a$ / m")
+
+fig.subplots_adjust(left=0.15, right=0.98)
+fig.savefig("contact_radii.png")
+plt.close(fig)
+
+
+fig, ax = plt.subplots()
+ax.scatter(rs_contact[1:], Fs, label="ogs")
+
+rs_ref = np.linspace(0, max(rs_contact), 200)
+Fs_ref = 8. * rs_ref**3 * kappa / 3.0
+
+ax.plot(rs_ref, Fs_ref, label="ref")
+
+l = ax.legend()
+l.get_frame().set_facecolor('white')
+
+ax.set_xlabel("radius of contact area $a$ / m")
+ax.set_ylabel("applied force $F$ / N")
+
+fig.subplots_adjust(left=0.15, right=0.98)
+fig.savefig("total_force.png")
+plt.close(fig)
diff --git a/Tests/Data/Mechanics/Linear/PythonHertzContact/ref_hertz_contact_ts_10.vtu b/Tests/Data/Mechanics/Linear/PythonHertzContact/ref_hertz_contact_ts_10.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..7c76d78d096efd3e098be1fc73a906fc313da710
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonHertzContact/ref_hertz_contact_ts_10.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:71bc3a358b09c749174bb760d050cc76ab5b6ca748de2b67c54421980b041e00
+size 3679599
diff --git a/Tests/Data/Mechanics/Linear/PythonHertzContact/ref_hertz_contact_ts_5.vtu b/Tests/Data/Mechanics/Linear/PythonHertzContact/ref_hertz_contact_ts_5.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..fc948ec69f524baebd69398aa396a2875caae011
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonHertzContact/ref_hertz_contact_ts_5.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d0837cb953245e61369946031c78d6e8c1abfe183ed7fc1bb62655be91a173ed
+size 3688217
diff --git a/Tests/Data/Mechanics/Linear/PythonHertzContact/unit-circle.gml b/Tests/Data/Mechanics/Linear/PythonHertzContact/unit-circle.gml
new file mode 100644
index 0000000000000000000000000000000000000000..7f452341f1e87acf5880730aa0d48bda13302227
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonHertzContact/unit-circle.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fc1fa83a8c78e8850ac8536908d409c67cdddc7cc6eeb13e8d7ded8a7e52c447
+size 20412
diff --git a/Tests/Data/Mechanics/Linear/PythonHertzContact/unit-circle.vtu b/Tests/Data/Mechanics/Linear/PythonHertzContact/unit-circle.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..7d98625cf830548a032f6f54a98aba6ed117243e
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonHertzContact/unit-circle.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:aa46da4364e860a71edfbd753adce21ca60a5c75a8890df25d847e28e1d7036a
+size 326511
diff --git a/web/content/docs/benchmarks/python-bc/hertz-contact/contact_radii.png b/web/content/docs/benchmarks/python-bc/hertz-contact/contact_radii.png
new file mode 100644
index 0000000000000000000000000000000000000000..92b7664b73b9b299284ba979ef6f0e69384dc7eb
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/hertz-contact/contact_radii.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7f45567cca7e784dba58159935b5370cd3b4c28d882fb486f6110f926f0e2ef2
+size 21649
diff --git a/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.pandoc b/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.pandoc
new file mode 100644
index 0000000000000000000000000000000000000000..57af17b0ac0cbf26b3488ca802410361817dc772
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.pandoc
@@ -0,0 +1,100 @@
++++
+project = "https://github.com/ufz/ogs-data/blob/master/Mechanics/Linear/Python/hertz-contact.prj"
+author = "Christoph Lehmann"
+date = "2018-08-06T11:41:00+02:00"
+title = "Hertz Contact using Python Boundary Conditions"
+weight = 4
+
+[menu]
+  [menu.benchmarks]
+    parent = "python-bc"
+
++++
+
+{{< data-link >}}
+
+## Motivation of this test case
+
+The aim of this test is:
+
+* to show that it is possible to apply Dirichlet BCs at a boundary that changes over the course of the simulation
+* to give an advanced use-case of the Python BC.
+  Here essentially an iterative procedure for a contact problem is implemented
+  within the Python BC.
+
+## Problem description
+
+Two elastic spheres of same radius $R$ are brought into contact.
+The sphere centers are displaced towards each other by $w_0$, with increasing
+values in every load step.
+Due to symmetry reasons a flat circular contact area of radius $a$ forms.
+
+{{< img src="../hertz-contact.png" >}}
+
+The contact between the two spheres is modelled as a Dirichlet BC
+on a varying boundary. The exact boundary and Dirichlet values for the
+$y$ displacements are determined in a Python script.
+Compared to the sketch above the prescribed $y$ displacements correspond
+to $w_0/2$, because due to symmetry only half of the system (a section of the
+lower sphere) is simulated.
+
+
+## Analytical solution
+
+The derivation of the formulae below can be found, e.g.,
+in [this book (in German)](http://www.uni-magdeburg.de/ifme/l-festigkeit/pdf/Bertram-Gluege_Festkoerpermechanik2012.pdf).
+
+The radius of the contact area is given by
+$$
+\begin{equation}
+a = \sqrt{\frac{w_0 R}{2}}
+\end{equation}
+$$
+
+The average pressure $\bar p$ over a the secant with distance $\xi$ to the
+center of the contact area (cf. vertical dashed line in the sketch above) is assumed to be
+$$
+\begin{equation}
+\bar p(\xi) = \kappa \sqrt{a^2 - \xi^2}
+\label{eq:bar-p}
+\end{equation}
+$$
+with the prefactor $\kappa$ given by
+$$
+\begin{equation}
+\kappa = \frac{G}{R \cdot (1-\nu)}
+\,.
+\end{equation}
+$$
+
+The total force $F$ exerted on the contact area is given by
+$$
+\begin{equation}
+F = \frac{8 a^3}{3\kappa}
+\,.
+\end{equation}
+$$
+
+
+## Results
+
+Contact radii:
+
+{{<img src="../contact_radii.png">}}
+
+Average pressure $\bar{p}$:
+
+{{<img src="../stress_at_contact.png">}}
+
+Total force $F$:
+
+{{<img src="../total_force.png">}}
+
+The simulation results for contact radii and total force reproduce the
+analytical square root and cubic laws, respectively.
+For the average pressure $\bar p$ the analytical form of
+Eq.&nbsp;$(\ref{eq:bar-p})$ is reproduced, too.
+But for $\bar p$ there are rather strong deviations between the numerical
+and analytical results, which might be due to deviations in the
+contact radii&nbsp;$a$, due to unsufficient mesh resolution or due to
+the chosen linear order of FEM ansatz functions.
diff --git a/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.png b/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.png
new file mode 100644
index 0000000000000000000000000000000000000000..336d096192a3f3a931729057fb2e2e15eebe3df4
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d129535c6bbe64e46936b8760f71720d293ee206ce617bff6a06160c1b25b5ef
+size 21644
diff --git a/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.svg b/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.svg
new file mode 100644
index 0000000000000000000000000000000000000000..07dc91e4480d299fa11260be2401b008b71d9a92
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.svg
@@ -0,0 +1,382 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!-- Created with Inkscape (http://www.inkscape.org/) -->
+
+<svg
+   xmlns:dc="http://purl.org/dc/elements/1.1/"
+   xmlns:cc="http://creativecommons.org/ns#"
+   xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
+   xmlns:svg="http://www.w3.org/2000/svg"
+   xmlns="http://www.w3.org/2000/svg"
+   xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
+   xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
+   width="241.4953"
+   height="199.97981"
+   viewBox="0 0 63.895632 52.911329"
+   version="1.1"
+   id="svg8"
+   inkscape:version="0.92.2 2405546, 2018-03-11"
+   sodipodi:docname="hertz-contact.svg"
+   inkscape:export-filename="/home/lehmannc/prog/ogs6/github-chleh-PRs/src/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.png"
+   inkscape:export-xdpi="200"
+   inkscape:export-ydpi="200">
+  <defs
+     id="defs2">
+    <marker
+       inkscape:stockid="Arrow1Mend"
+       orient="auto"
+       refY="0"
+       refX="0"
+       id="marker5094"
+       style="overflow:visible"
+       inkscape:isstock="true">
+      <path
+         inkscape:connector-curvature="0"
+         id="path5092"
+         d="M 0,0 5,-5 -12.5,0 5,5 Z"
+         style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.00000003pt;stroke-opacity:1"
+         transform="matrix(-0.4,0,0,-0.4,-4,0)" />
+    </marker>
+    <marker
+       inkscape:isstock="true"
+       style="overflow:visible"
+       id="marker4996"
+       refX="0"
+       refY="0"
+       orient="auto"
+       inkscape:stockid="Arrow1Mend"
+       inkscape:collect="always">
+      <path
+         transform="matrix(-0.4,0,0,-0.4,-4,0)"
+         style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.00000003pt;stroke-opacity:1"
+         d="M 0,0 5,-5 -12.5,0 5,5 Z"
+         id="path4994"
+         inkscape:connector-curvature="0" />
+    </marker>
+    <marker
+       inkscape:stockid="Arrow1Mend"
+       orient="auto"
+       refY="0"
+       refX="0"
+       id="marker4544"
+       style="overflow:visible"
+       inkscape:isstock="true"
+       inkscape:collect="always">
+      <path
+         id="path4542"
+         d="M 0,0 5,-5 -12.5,0 5,5 Z"
+         style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.00000003pt;stroke-opacity:1"
+         transform="matrix(-0.4,0,0,-0.4,-4,0)"
+         inkscape:connector-curvature="0" />
+    </marker>
+    <marker
+       inkscape:isstock="true"
+       style="overflow:visible"
+       id="marker4136"
+       refX="0"
+       refY="0"
+       orient="auto"
+       inkscape:stockid="Arrow1Mend">
+      <path
+         transform="matrix(-0.4,0,0,-0.4,-4,0)"
+         style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.00000003pt;stroke-opacity:1"
+         d="M 0,0 5,-5 -12.5,0 5,5 Z"
+         id="path4134"
+         inkscape:connector-curvature="0" />
+    </marker>
+    <marker
+       inkscape:stockid="Arrow1Mend"
+       orient="auto"
+       refY="0"
+       refX="0"
+       id="marker3552"
+       style="overflow:visible"
+       inkscape:isstock="true"
+       inkscape:collect="always">
+      <path
+         id="path3550"
+         d="M 0,0 5,-5 -12.5,0 5,5 Z"
+         style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.00000003pt;stroke-opacity:1"
+         transform="matrix(-0.4,0,0,-0.4,-4,0)"
+         inkscape:connector-curvature="0" />
+    </marker>
+    <marker
+       inkscape:isstock="true"
+       style="overflow:visible"
+       id="marker3052"
+       refX="0"
+       refY="0"
+       orient="auto"
+       inkscape:stockid="Arrow1Mend"
+       inkscape:collect="always">
+      <path
+         transform="matrix(-0.4,0,0,-0.4,-4,0)"
+         style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.00000003pt;stroke-opacity:1"
+         d="M 0,0 5,-5 -12.5,0 5,5 Z"
+         id="path3050"
+         inkscape:connector-curvature="0" />
+    </marker>
+    <marker
+       inkscape:stockid="Arrow1Mend"
+       orient="auto"
+       refY="0"
+       refX="0"
+       id="Arrow1Mend"
+       style="overflow:visible"
+       inkscape:isstock="true"
+       inkscape:collect="always">
+      <path
+         id="path833"
+         d="M 0,0 5,-5 -12.5,0 5,5 Z"
+         style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.00000003pt;stroke-opacity:1"
+         transform="matrix(-0.4,0,0,-0.4,-4,0)"
+         inkscape:connector-curvature="0" />
+    </marker>
+    <marker
+       inkscape:stockid="Arrow1Mstart"
+       orient="auto"
+       refY="0"
+       refX="0"
+       id="Arrow1Mstart"
+       style="overflow:visible"
+       inkscape:isstock="true">
+      <path
+         id="path830"
+         d="M 0,0 5,-5 -12.5,0 5,5 Z"
+         style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.00000003pt;stroke-opacity:1"
+         transform="matrix(0.4,0,0,0.4,4,0)"
+         inkscape:connector-curvature="0" />
+    </marker>
+    <marker
+       inkscape:stockid="Arrow1Send"
+       orient="auto"
+       refY="0"
+       refX="0"
+       id="Arrow1Send"
+       style="overflow:visible"
+       inkscape:isstock="true">
+      <path
+         id="path839"
+         d="M 0,0 5,-5 -12.5,0 5,5 Z"
+         style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.00000003pt;stroke-opacity:1"
+         transform="matrix(-0.2,0,0,-0.2,-1.2,0)"
+         inkscape:connector-curvature="0" />
+    </marker>
+    <marker
+       inkscape:stockid="Arrow1Sstart"
+       orient="auto"
+       refY="0"
+       refX="0"
+       id="Arrow1Sstart"
+       style="overflow:visible"
+       inkscape:isstock="true">
+      <path
+         id="path836"
+         d="M 0,0 5,-5 -12.5,0 5,5 Z"
+         style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.00000003pt;stroke-opacity:1"
+         transform="matrix(0.2,0,0,0.2,1.2,0)"
+         inkscape:connector-curvature="0" />
+    </marker>
+  </defs>
+  <sodipodi:namedview
+     id="base"
+     pagecolor="#ffffff"
+     bordercolor="#666666"
+     borderopacity="1.0"
+     inkscape:pageopacity="0.0"
+     inkscape:pageshadow="2"
+     inkscape:zoom="1.979899"
+     inkscape:cx="159.92526"
+     inkscape:cy="104.36074"
+     inkscape:document-units="mm"
+     inkscape:current-layer="layer1"
+     showgrid="false"
+     inkscape:window-width="1600"
+     inkscape:window-height="875"
+     inkscape:window-x="0"
+     inkscape:window-y="25"
+     inkscape:window-maximized="1"
+     units="px"
+     inkscape:snap-global="false"
+     fit-margin-top="0"
+     fit-margin-left="0"
+     fit-margin-right="0"
+     fit-margin-bottom="0" />
+  <metadata
+     id="metadata5">
+    <rdf:RDF>
+      <cc:Work
+         rdf:about="">
+        <dc:format>image/svg+xml</dc:format>
+        <dc:type
+           rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
+        <dc:title></dc:title>
+      </cc:Work>
+    </rdf:RDF>
+  </metadata>
+  <g
+     inkscape:label="Ebene 1"
+     inkscape:groupmode="layer"
+     id="layer1"
+     transform="translate(-38.044515,-206.59195)">
+    <path
+       d="m 53.67722,228.37355 a 15.497022,15.497022 0 0 1 15.497021,15.49702 H 53.677219 Z"
+       sodipodi:end="0"
+       sodipodi:start="4.712389"
+       sodipodi:ry="15.497022"
+       sodipodi:rx="15.497022"
+       sodipodi:cy="243.87057"
+       sodipodi:cx="53.677219"
+       sodipodi:type="arc"
+       id="path4075"
+       style="opacity:1;fill:#cccccc;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1" />
+    <circle
+       style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+       id="path815"
+       cx="-85.251534"
+       cy="212.12242"
+       r="15.497022"
+       transform="rotate(-35.474654)" />
+    <path
+       style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:0.52916665, 0.52916665;stroke-dashoffset:0;stroke-opacity:1"
+       d="M 42.472232,233.04379 H 64.774705"
+       id="path818"
+       inkscape:connector-curvature="0" />
+    <circle
+       r="15.497022"
+       cy="50.393448"
+       cx="244.57024"
+       id="circle820"
+       style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+       transform="rotate(65.944031)" />
+    <path
+       style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-start:url(#Arrow1Mend);marker-end:url(#Arrow1Mstart)"
+       d="m 53.677221,228.20874 v 9.61452"
+       id="path1360"
+       inkscape:connector-curvature="0" />
+    <text
+       xml:space="preserve"
+       style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;line-height:6.61458302px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans';font-variant-ligatures:normal;font-variant-position:normal;font-variant-caps:normal;font-variant-numeric:normal;font-variant-alternates:normal;font-feature-settings:normal;text-indent:0;text-align:start;text-decoration:none;text-decoration-line:none;text-decoration-style:solid;text-decoration-color:#000000;letter-spacing:0px;word-spacing:0px;text-transform:none;writing-mode:lr-tb;direction:ltr;text-orientation:mixed;dominant-baseline:auto;baseline-shift:baseline;text-anchor:start;white-space:normal;shape-padding:0;opacity:1;vector-effect:none;fill:#000000;fill-opacity:1;stroke:none;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+       x="54.004654"
+       y="230.85077"
+       id="text3030"><tspan
+         sodipodi:role="line"
+         id="tspan3028"
+         x="54.004654"
+         y="230.85077"
+         style="stroke-width:0.26458332px"><tspan
+           style="font-style:oblique;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans, Oblique';font-variant-ligatures:normal;font-variant-caps:normal;font-variant-numeric:normal;font-feature-settings:normal;text-align:start;writing-mode:lr-tb;text-anchor:start"
+           id="tspan3032">w</tspan><tspan
+           style="font-size:64.99999762%;baseline-shift:sub"
+           id="tspan3034">0</tspan></tspan></text>
+    <circle
+       r="11.018708"
+       cy="238.72971"
+       cx="-37.353424"
+       id="circle3040"
+       style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+       transform="rotate(-30.961967)" />
+    <path
+       style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#marker3052)"
+       d="m 90.887678,223.87032 9.015421,-5.40887"
+       id="path3042"
+       inkscape:connector-curvature="0" />
+    <text
+       xml:space="preserve"
+       style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;line-height:6.61458302px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans';font-variant-ligatures:normal;font-variant-position:normal;font-variant-caps:normal;font-variant-numeric:normal;font-variant-alternates:normal;font-feature-settings:normal;text-indent:0;text-align:start;text-decoration:none;text-decoration-line:none;text-decoration-style:solid;text-decoration-color:#000000;letter-spacing:0px;word-spacing:0px;text-transform:none;writing-mode:lr-tb;direction:ltr;text-orientation:mixed;dominant-baseline:auto;baseline-shift:baseline;text-anchor:start;white-space:normal;shape-padding:0;opacity:1;vector-effect:none;fill:#000000;fill-opacity:1;stroke:none;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+       x="94.475052"
+       y="219.86319"
+       id="text3098"><tspan
+         sodipodi:role="line"
+         id="tspan3096"
+         x="94.475052"
+         y="219.86319"
+         style="font-style:oblique;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans, Oblique';font-variant-ligatures:normal;font-variant-caps:normal;font-variant-numeric:normal;font-feature-settings:normal;text-align:start;writing-mode:lr-tb;text-anchor:start;stroke-width:0.26458332px">a</tspan></text>
+    <text
+       xml:space="preserve"
+       style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;line-height:6.61458302px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans';font-variant-ligatures:normal;font-variant-position:normal;font-variant-caps:normal;font-variant-numeric:normal;font-variant-alternates:normal;font-feature-settings:normal;text-indent:0;text-align:start;text-decoration:none;text-decoration-line:none;text-decoration-style:solid;text-decoration-color:#000000;letter-spacing:0px;word-spacing:0px;text-transform:none;writing-mode:lr-tb;direction:ltr;text-orientation:mixed;dominant-baseline:auto;baseline-shift:baseline;text-anchor:start;white-space:normal;shape-padding:0;opacity:1;vector-effect:none;fill:#000000;fill-opacity:1;stroke:none;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+       x="85.289948"
+       y="227.97792"
+       id="text3540"><tspan
+         sodipodi:role="line"
+         id="tspan3538"
+         x="85.289948"
+         y="227.97792"
+         style="stroke-width:0.26458332px">contact area</tspan></text>
+    <path
+       inkscape:connector-curvature="0"
+       id="path3548"
+       d="m 60.050003,232.8548 19.1589,-5.6054"
+       style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#marker3552)"
+       sodipodi:nodetypes="cc" />
+    <path
+       sodipodi:nodetypes="cc"
+       style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#marker4136)"
+       d="m 60.317273,241.40742 16.35257,2.67996"
+       id="path4132"
+       inkscape:connector-curvature="0" />
+    <text
+       xml:space="preserve"
+       style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;line-height:6.61458302px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans';font-variant-ligatures:normal;font-variant-position:normal;font-variant-caps:normal;font-variant-numeric:normal;font-variant-alternates:normal;font-feature-settings:normal;text-indent:0;text-align:start;text-decoration:none;text-decoration-line:none;text-decoration-style:solid;text-decoration-color:#000000;letter-spacing:0px;word-spacing:0px;text-transform:none;writing-mode:lr-tb;direction:ltr;text-orientation:mixed;dominant-baseline:auto;baseline-shift:baseline;text-anchor:start;white-space:normal;shape-padding:0;opacity:1;vector-effect:none;fill:#000000;fill-opacity:1;stroke:none;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+       x="78.176338"
+       y="244.88242"
+       id="text4930"><tspan
+         sodipodi:role="line"
+         id="tspan4928"
+         x="78.176338"
+         y="244.88242"
+         style="stroke-width:0.26458332px">simulation domain</tspan></text>
+    <path
+       style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#marker4544)"
+       d="m 53.461273,222.37563 12.519405,-8.92165"
+       id="path4932"
+       inkscape:connector-curvature="0" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path4992"
+       d="m 53.569132,243.62843 6.266508,14.03789"
+       style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#marker4996)" />
+    <text
+       xml:space="preserve"
+       style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;line-height:6.61458302px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans';font-variant-ligatures:normal;font-variant-position:normal;font-variant-caps:normal;font-variant-numeric:normal;font-variant-alternates:normal;font-feature-settings:normal;text-indent:0;text-align:start;text-decoration:none;text-decoration-line:none;text-decoration-style:solid;text-decoration-color:#000000;letter-spacing:0px;word-spacing:0px;text-transform:none;writing-mode:lr-tb;direction:ltr;text-orientation:mixed;dominant-baseline:auto;baseline-shift:baseline;text-anchor:start;white-space:normal;shape-padding:0;opacity:1;vector-effect:none;fill:#000000;fill-opacity:1;stroke:none;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+       x="55.7257"
+       y="218.15546"
+       id="text5078"><tspan
+         sodipodi:role="line"
+         id="tspan5076"
+         x="55.7257"
+         y="218.15546"
+         style="font-style:oblique;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans, Oblique';font-variant-ligatures:normal;font-variant-caps:normal;font-variant-numeric:normal;font-feature-settings:normal;text-align:start;writing-mode:lr-tb;text-anchor:start;stroke-width:0.26458332px">R</tspan></text>
+    <text
+       id="text5082"
+       y="252.09869"
+       x="53.854813"
+       style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;line-height:6.61458302px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans';font-variant-ligatures:normal;font-variant-position:normal;font-variant-caps:normal;font-variant-numeric:normal;font-variant-alternates:normal;font-feature-settings:normal;text-indent:0;text-align:start;text-decoration:none;text-decoration-line:none;text-decoration-style:solid;text-decoration-color:#000000;letter-spacing:0px;word-spacing:0px;text-transform:none;writing-mode:lr-tb;direction:ltr;text-orientation:mixed;dominant-baseline:auto;baseline-shift:baseline;text-anchor:start;white-space:normal;shape-padding:0;opacity:1;vector-effect:none;fill:#000000;fill-opacity:1;stroke:none;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+       xml:space="preserve"><tspan
+         style="font-style:oblique;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans, Oblique';font-variant-ligatures:normal;font-variant-caps:normal;font-variant-numeric:normal;font-feature-settings:normal;text-align:start;writing-mode:lr-tb;text-anchor:start;stroke-width:0.26458332px"
+         y="252.09869"
+         x="53.854813"
+         id="tspan5080"
+         sodipodi:role="line">R</tspan></text>
+    <path
+       style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#marker5094)"
+       d="M 90.865278,223.93005 H 83.581564"
+       id="path5084"
+       inkscape:connector-curvature="0" />
+    <text
+       xml:space="preserve"
+       style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;line-height:6.61458302px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans';font-variant-ligatures:normal;font-variant-position:normal;font-variant-caps:normal;font-variant-numeric:normal;font-variant-alternates:normal;font-feature-settings:normal;text-indent:0;text-align:start;text-decoration:none;text-decoration-line:none;text-decoration-style:solid;text-decoration-color:#000000;letter-spacing:0px;word-spacing:0px;text-transform:none;writing-mode:lr-tb;direction:ltr;text-orientation:mixed;dominant-baseline:auto;baseline-shift:baseline;text-anchor:start;white-space:normal;shape-padding:0;opacity:1;vector-effect:none;fill:#000000;fill-opacity:1;stroke:none;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+       x="87.559395"
+       y="222.74509"
+       id="text5200"><tspan
+         sodipodi:role="line"
+         id="tspan5198"
+         x="87.559395"
+         y="222.74509"
+         style="font-style:oblique;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans, Oblique';font-variant-ligatures:normal;font-variant-caps:normal;font-variant-numeric:normal;font-feature-settings:normal;text-align:start;writing-mode:lr-tb;text-anchor:start;stroke-width:0.26458332px">ξ</tspan></text>
+    <path
+       style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:0.52916661, 0.52916661;stroke-dashoffset:0;stroke-opacity:1"
+       d="m 83.023678,216.03601 v 15.63875"
+       id="path5316"
+       inkscape:connector-curvature="0" />
+  </g>
+</svg>
diff --git a/web/content/docs/benchmarks/python-bc/hertz-contact/stress_at_contact.png b/web/content/docs/benchmarks/python-bc/hertz-contact/stress_at_contact.png
new file mode 100644
index 0000000000000000000000000000000000000000..3aa8f0146621819024074f347e147de509822b28
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/hertz-contact/stress_at_contact.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8ec5b3f1339f960f105b99b7bd092ce00bc985bb52624df51df9182b7d4b6e65
+size 86699
diff --git a/web/content/docs/benchmarks/python-bc/hertz-contact/total_force.png b/web/content/docs/benchmarks/python-bc/hertz-contact/total_force.png
new file mode 100644
index 0000000000000000000000000000000000000000..2cb423ffa402b68fb81373ba505361c9ccf0c662
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/hertz-contact/total_force.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4284b5904a527b7e5228fb45b7d211064f9a373764f4096ed0e2e7a13fb225a3
+size 19421