{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hertzian contact with a rigid indenter using a penalty approach\n",
    "\n",
    "In this numerical tour, we explore the formulation of frictionless contact between a rigid surface (the indenter) and an elastic domain, representing an infinite half-space for the present case. Contact will be solved using a penalty formulation, that is small interpenetration between the solid and the indenter will be authorized. Prerequisites for this tour are the formulation of an linear elasticity problem and the resolution of a general nonlinear variational problem. See for instance:\n",
    "\n",
    "* https://comet-fenics.readthedocs.io/en/latest/demo/elasticity/2D_elasticity.py.html\n",
    "* https://fenicsproject.org/olddocs/dolfin/latest/python/demos/hyperelasticity/demo_hyperelasticity.py.html"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. figure:: contact_problem.png\n",
    "   :scale: 20%\n",
    "   :align: center\n",
    "   \n",
    "   Hertzian contact problem of a rigid spherical indenter on a semi-infinite domain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The elastic (isotropic $E,\\nu$) solid will be represented by a 3D cubic domain of unit dimension and contact will take place on its top surface $z=0$, centered at $x=y=0$. Symmetry conditions will be applied on the $x=0$ and $y=0$ surfaces whereas the bottom surface $z=-1$ will be fully fixed. If contact appears on a small region of extent $a\\ll 1$, the problem can then be considered to be a good approximation of contact on a semi-infinite domain.\n",
    "\n",
    "The rigid indenter will not be explictly modeled (in particular it will not be meshed) but instead its distance with respect to the solid top surface will be given as an ``Expression``. We will also consider that the indenter radius $R$ is sufficiently large with respect to the contact region characteristic size $a$ so that the spherical surface can be approximated by a parabola. In this case, the distance between such an indenter and the top surface can be written as:\n",
    "\n",
    "$$h(x,y) = h_0 + \\dfrac{1}{2R}(x^2+y^2)$$\n",
    "\n",
    "where $h_0$ is the initial gap between both surfaces at $x=y=0$. Obviously, if $h_0>0$ there is no contact between both surfaces, contact appears only if $h_0=-d<0$ where $d$ will be the indenter depth inside the surface. This classical problem admits the following known analytical solution [[1,2]](#References):\n",
    "\n",
    "* the contact area is of circular shape and radius $a=\\sqrt{Rd}$\n",
    "* the force exerted by the indenter onto the surface is $F=\\dfrac{4}{3}\\dfrac{E}{1-\\nu^2}ad$\n",
    "* the pressure distribution on the contact region is given by $p(r) = p_0\\sqrt{1-(r/a)^2}$ where $p_0=3F/(2\\pi a^2)$ is the maximal pressure\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. figure:: contact_solution.png\n",
    "   :scale: 20%\n",
    "   :align: center\n",
    "   \n",
    "   Indent at depth :math:`d` (left) and contact solution (right)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contact problem formulation and penalty approach\n",
    "\n",
    "The unilateral contact (Signorini) condition on the top surface $\\Gamma$ writes as:\n",
    "\n",
    "$$g \\geq 0, p\\geq 0, g\\cdot p =0 \\text{ on }\\Gamma$$\n",
    "\n",
    "where $g=h-u$ is the gap between the obstacle surface and the solid surface and $p=-\\sigma_{zz}$ is the pressure. One of the most simple way to solve approximately this contact condition consists in replacing the previous complementary conditions by the following penalized condition:\n",
    "\n",
    "$$p = k\\langle -g\\rangle_+ = k\\langle u-h\\rangle_+$$\n",
    "\n",
    "where $\\langle x\\rangle_+ = (|x|+x)/2$ is the positive part (Mackauley bracket) and $k$ is a large penalizing stiffness coefficient. With the previous relation, the pressure will be positive but a small negative gap will be authorized.\n",
    "\n",
    "## Implementation\n",
    "\n",
    "First, a unit cubic mesh is defined and some mapping is applied on the mesh nodes in order to refine the element size near the contact area (around $x=y=z=0$) and change to a $[0;1]\\times [0;1]\\times[-1;0]$ domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "import numpy as np\n",
    "\n",
    "N = 30\n",
    "mesh = UnitCubeMesh.create(N, N, N//2, CellType.Type.hexahedron)\n",
    "# mesh size is smaller near x=y=0\n",
    "mesh.coordinates()[:, :2] = mesh.coordinates()[:, :2]**2\n",
    "# mesh size is smaller near z=0 and mapped to a [-1;0] domain along z\n",
    "mesh.coordinates()[:, 2] = -mesh.coordinates()[:, 2]**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The top surface is defined as a ``SubDomain`` and corresponding exterior facets are marked as 1. Functions for imposition of Dirichlet BC on other boundaries are also defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Top(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        return near(x[2], 0.) and on_boundary\n",
    "def symmetry_x(x, on_boundary):\n",
    "        return near(x[0], 0) and on_boundary\n",
    "def symmetry_y(x, on_boundary):\n",
    "        return near(x[1], 0) and on_boundary\n",
    "def bottom(x, on_boundary):\n",
    "        return near(x[2], -1) and on_boundary\n",
    "    \n",
    "# exterior facets MeshFunction\n",
    "facets = MeshFunction(\"size_t\", mesh, 2)\n",
    "facets.set_all(0)\n",
    "Top().mark(facets, 1)\n",
    "ds = Measure('ds', subdomain_data=facets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The obstacle shape $h(x,y)$ is now defined and a small indentation depth $d$ is considered. Standard function spaces are then defined for the problem formulation and output field exports. In particular, gap and pressure will be saved later. Finally, the ``DirichletBC`` are defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = 0.5\n",
    "d = 0.02\n",
    "obstacle = Expression(\"-d+(pow(x[0],2)+pow(x[1], 2))/2/R\", d=d, R=R, degree=2)\n",
    "\n",
    "V = VectorFunctionSpace(mesh, \"CG\", 1)\n",
    "V2 = FunctionSpace(mesh, \"CG\", 1)\n",
    "V0 = FunctionSpace(mesh, \"DG\", 0)\n",
    "\n",
    "u = Function(V, name=\"Displacement\")\n",
    "du = TrialFunction(V)\n",
    "u_ = TestFunction(V)\n",
    "gap = Function(V2, name=\"Gap\")\n",
    "p = Function(V0, name=\"Contact pressure\")\n",
    "\n",
    "bc =[DirichletBC(V, Constant((0., 0., 0.)), bottom),\n",
    "     DirichletBC(V.sub(0), Constant(0.), symmetry_x),\n",
    "     DirichletBC(V.sub(1), Constant(0.), symmetry_y)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The elastic constants and functions for formulating the weak form are defined. A penalization stiffness ``pen`` is introduced and the penalized weak form\n",
    "\n",
    "$$\\text{Find } u \\text{ such that } \\int_{\\Omega} \\sigma(u):\\varepsilon(v)d\\Omega + k_{pen}\\int_{\\Gamma} \\langle u-h\\rangle_+ v dS = 0$$\n",
    "\n",
    "is defined. The corresponding Jacobian ``J`` is also computed using automatic differentiation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = Constant(10.)\n",
    "nu = Constant(0.3)\n",
    "mu = E/2/(1+nu)\n",
    "lmbda = E*nu/(1+nu)/(1-2*nu)\n",
    "def eps(v):\n",
    "    return sym(grad(v))\n",
    "def sigma(v):\n",
    "    return lmbda*tr(eps(v))*Identity(3) + 2.0*mu*eps(v)\n",
    "def ppos(x):\n",
    "    return (x+abs(x))/2.\n",
    "\n",
    "pen = Constant(1e4)\n",
    "form = inner(sigma(u), eps(u_))*dx + pen*dot(u_[2], ppos(u[2]-obstacle))*ds(1)\n",
    "J = derivative(form, u, du)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A non-linear solver is now defined for solving the previous problem. Due to the 3D nature of the problem with a significant number of elements, an iterative Conjugate-Gradient solver is chosen for the linear solver inside the Newton iterations. Note that choosing a large penalization parameter deteriorates the problem conditioning so that solving time will drastically increase and can eventually fail. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8, True)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "problem = NonlinearVariationalProblem(form, u, bc, J=J)\n",
    "solver = NonlinearVariationalSolver(problem)\n",
    "solver.parameters[\"newton_solver\"][\"linear_solver\"] = \"cg\"\n",
    "solver.parameters[\"newton_solver\"][\"preconditioner\"] = \"ilu\"\n",
    "\n",
    "solver.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As post-processing, the gap is computed (here on the whole mesh) as well as the pressure. The maximal pressure and the total force resultant (note the factor 4 due to the symmetry of the problem) are compared with respect to the analytical solution. Finally, XDMF output is performed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum pressure FE:  1.42763 Hertz:  1.39916\n",
      "Applied force    FE:  0.03206 Hertz:  0.02930\n"
     ]
    }
   ],
   "source": [
    "p.assign(-project(sigma(u)[2, 2], V0))\n",
    "gap.assign(project(obstacle-u[2], V2))\n",
    "\n",
    "a = sqrt(R*d)\n",
    "F = 4/3.*float(E)/(1-float(nu)**2)*a*d\n",
    "p0 = 3*F/(2*pi*a**2)\n",
    "print(\"Maximum pressure FE: {0:8.5f} Hertz: {1:8.5f}\".format(max(np.abs(p.vector().get_local())), p0))\n",
    "print(\"Applied force    FE: {0:8.5f} Hertz: {1:8.5f}\".format(4*assemble(p*ds(1)), F))\n",
    "\n",
    "file_results = XDMFFile(\"contact_penalty_results.xdmf\")\n",
    "file_results.parameters[\"flush_output\"] = True\n",
    "file_results.parameters[\"functions_share_mesh\"] = True\n",
    "file_results.write(u, 0.)\n",
    "file_results.write(gap, 0.)\n",
    "file_results.write(p, 0.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maximum pressure over the contact surface and total applied force by the indenter are computed and compared against the analytical solution. Comparison between both solutions shows good agreement although the quality of the results is quite dependent on the chosen penalty sitffness. Here is what the solution looks like."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. figure:: contact_FE.png\n",
    "   :scale: 50%\n",
    "   :align: center"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "[1] K.L., Johnson, Contact Mechanics, Cambridge University Press, 1987\n",
    "\n",
    "[2] S.P. Timoshenko, J.N. Goodier. Theory of Elasticity, MacGraw Hill International Book, 1982"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}