Added Hertz contact test

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
+)

Tests/Data/Mechanics/Linear/PythonHertzContact/gen-unit-circle.py

+#!/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"?>
+<?xml-stylesheet type="text/xsl" href="OpenGeoSysGLI.xsl"?>
+
+<OpenGeoSysGLI xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://www.opengeosys.org/images/xsd/OpenGeoSysGLI.xsd" xmlns:ogs="http://www.opengeosys.org">
+    <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>""")

Tests/Data/Mechanics/Linear/PythonHertzContact/hertz_contact.prj

+<?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>

Tests/Data/Mechanics/Linear/PythonHertzContact/hertz_contact_bc.py

+from __future__ import print_function
+
+import OpenGeoSys
+
+
+# enable/disable debug output
+if True:
+    def debug(*args):
+        print(*args)
+else:
+    def debug(*args):
+        pass
+
+
+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
+
+
+    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
+
+
+    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:
+            # there were nodes above the contact area in the previous
+            # iteration => contact area has to grow.
+            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]
+            else:
+                self._a_range[0] = self._a_est
+        else:
+            # contact area has to shrink
+            self._a_range[1] = self._a_est
+
+        self._a_est = 0.5 * sum(self._a_range)
+        print("BC: a_est={:.4f}, ({:.4f}, {:.4f})".format(self._a_est, self._a_range[0], self._a_range[1]))
+
+        self._max_x_with_y_excess = -1.0
+
+
+    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()
+
+        x, y, z = coords
+        ux, uy = primary_vars
+
+        try:
+            # check that we are at the outer boundary
+            assert abs(x**2 + y**2 + z**2 - 1.0) < 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 - 1)
+            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)
+
+        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
+
+
+# instantiate the BC object used by OpenGeoSys
+bc = HertzContactBC()

Tests/Data/Mechanics/Linear/PythonHertzContact/post.py

+#!/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)

web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.md

++++
+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.