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/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..e7a4d0110886d4d7b6b3e0338146acfcdb0cadd0
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonHertzContact/gen-unit-circle.py
@@ -0,0 +1,126 @@
+#!/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>""")
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..5351190b431c9cac449abd5c8b1600783368e83b
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonHertzContact/hertz_contact_bc.py
@@ -0,0 +1,111 @@
+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()
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..390741dc42747d99437b5199aae8549464030286
--- /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:4c557c6b2f2dea1cc8d7be0a258a344ea5f14172996789fdfe3f43bfee738938
+size 20651
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.md b/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.md
new file mode 100644
index 0000000000000000000000000000000000000000..c249033c4300d8dc7cbc86cf44b540de77f519c9
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.md
@@ -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. $(\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 $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