{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# 6.2 Contact Problems\n",
    "\n",
    "work in progress ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "from ngsolve import *\n",
    "from ngsolve.solvers import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.webgui import Draw as DrawGeo\n",
    "\n",
    "rect = MoveTo(0,-0.02).Rectangle(0.04, 0.02).Face()\n",
    "ball = MoveTo(0,0).Arc(0.005, 90).Rotate(90).Line(0.005).Close().Face()\n",
    "\n",
    "rect.edges.Min(X).name=\"sym\"\n",
    "rect.edges.Max(X).name=\"free\"\n",
    "rect.edges.Min(Y).name=\"fix\"\n",
    "rect.edges.Max(Y).name=\"contact1\"\n",
    "ball.edges.Min(X).name=\"sym\"\n",
    "ball.edges.Max(Y).name=\"disp\"\n",
    "ball.edges.Max(X-Y).name=\"contact2\"\n",
    "ball.faces.name=\"ball\"\n",
    "rect.faces.name=\"rect\"\n",
    "rect = Glue ([rect, Vertex((0.001,0,0))])\n",
    "rect.vertices.Max(Y-X).maxh=0.0002\n",
    "rect.edges.Max(Y-0.1*X).maxh=0.00002\n",
    "geo = Compound([rect,ball])\n",
    "DrawGeo (geo)\n",
    "mesh = Mesh(OCCGeometry(geo, dim=2).GenerateMesh(maxh=0.002))\n",
    "mesh.Curve(4)\n",
    "Draw(mesh);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.GetBoundaries(), mesh.GetMaterials()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "E = 2.1e11\n",
    "nu = 0.2\n",
    "\n",
    "rho = 7e3\n",
    "mu  = E / 2 / (1+nu)\n",
    "lam = E * nu / ((1+nu)*(1-2*nu))\n",
    "\n",
    "I = Id(mesh.dim)\n",
    "def Pow(a, b):\n",
    "    return a**b  # exp (log(a)*b)\n",
    "\n",
    "def C(u):\n",
    "    # F = Id(2) + Grad(u)\n",
    "    G = Grad(u)\n",
    "    F = Id(3) + CoefficientFunction( (G[0,0], G[1,0], 0,   G[1,0], G[1,1], 0,    0, 0, u[0]/(x+1e-16) ) )\n",
    "    return F.trans*F\n",
    "\n",
    "def NeoHooke (C):\n",
    "    return 0.5 * mu * (Trace(C-Id(C.dims[0])) + 2*mu/lam * Det(C)**(-lam/2/mu) - 1)\n",
    "\n",
    "fes = VectorH1(mesh, order=4, dirichletx=\"sym\", dirichlety=\"fix|disp\")\n",
    "u,v = fes.TnT()\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += Variation((x*NeoHooke(C(u))).Compile()*dx)\n",
    "\n",
    "contact = ContactBoundary(mesh.Boundaries(\"contact1\"), mesh.Boundaries(\"contact2\"))\n",
    "\n",
    "X = CoefficientFunction((x,y))\n",
    "cf = (X + u - (X.Other() + u.Other())) * contact.normal\n",
    "\n",
    "# energy formulation\n",
    "# contact.AddEnergy(IfPos(cf, 1e14*cf*cf, 0))\n",
    "# or integrator\n",
    "contact.AddIntegrator(IfPos(cf, 1e14*cf*(v-v.Other())*contact.normal, 0))\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "\n",
    "disp = -10e-6\n",
    "gfu.Set((0, disp), definedon=mesh.Materials(\"ball\"))\n",
    "\n",
    "sceneu = Draw(gfu, deformation=gfu)\n",
    "\n",
    "C_ = C(gfu).MakeVariable()\n",
    "sigma = NeoHooke(C_).Diff(C_)\n",
    "scenesigma = Draw(sigma[1,1], mesh, \"sigmazz\", min=-2e9, max=2e9)\n",
    "\n",
    "\n",
    "with TaskManager():\n",
    "    contact.Update(gfu, a, 200, 1e-4)\n",
    "\n",
    "    NewtonMinimization(a=a, u=gfu, printing=False, inverse=\"sparsecholesky\")\n",
    "    sceneu.Redraw()\n",
    "    scenesigma.Redraw()\n",
    "\n",
    "Draw(sigma[1,1], mesh, \"sigmazz\", min=-5e9, max=0, deformation=gfu);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "## Dynamic Contact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "from ngsolve import *\n",
    "from ngsolve.solvers import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.webgui import Draw as DrawGeo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "bowl = MoveTo(0,-0.01).Arc(0.2,70).Rotate(90).Line(0.01).Rotate(90) \\\n",
    "    .Arc(0.19,-140).Rotate(90).Line(0.01).Rotate(90).Arc(0.2,70).Face()\n",
    "ball1 = Circle((0,0.1),0.01).Face()\n",
    "ball2 = Circle((0.05,0.1),0.01).Face()\n",
    "ball3 = Circle((-0.05,0.1),0.01).Face()\n",
    "balls = Compound([ball1, ball2, ball3])\n",
    "balls.edges.name = \"balls\"\n",
    "geo = Compound([bowl, balls])\n",
    "bowl.edges.name=\"contact\"\n",
    "bowl.edges.Min(Y+0.01*X).name=\"fix\"\n",
    "bowl.edges.Min(Y-0.01*X).name=\"fix\"\n",
    "DrawGeo(geo);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(OCCGeometry(geo, dim=2).GenerateMesh(maxh=0.005))\n",
    "mesh.Curve(4)\n",
    "Draw(mesh);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.GetBoundaries()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "E, nu = 210e6, 0.2\n",
    "rho = 7e3\n",
    "mu  = E / 2 / (1+nu)\n",
    "lam = E * nu / ((1+nu)*(1-2*nu))\n",
    "\n",
    "I = Id(mesh.dim)\n",
    "\n",
    "def C(u): \n",
    "    F = I+Grad(u)\n",
    "    return F.trans*F\n",
    "def NeoHooke (C):\n",
    "    return 0.5 * mu * (Trace(C-I) + 2*mu/lam * Det(C)**(-lam/2/mu) - 1)\n",
    "\n",
    "fes = VectorH1(mesh, order=3, dirichlet=\"fix\")\n",
    "u,v = fes.TnT()\n",
    "\n",
    "force = CF((0,-9.81*rho))\n",
    "\n",
    "uold = GridFunction(fes)\n",
    "unew = GridFunction(fes)\n",
    "vel = GridFunction(fes)\n",
    "anew = GridFunction(fes)\n",
    "aold = GridFunction(fes)\n",
    "\n",
    "tau = 2e-4\n",
    "\n",
    "bfmstar = BilinearForm(fes)\n",
    "bfmstar += Variation( NeoHooke (C(u)).Compile(False)*dx )\n",
    "bfmstar += Variation( -force*u*dx )\n",
    "bfmstar += Variation( rho/2* 2/tau**2 * (u-uold-tau*vel-tau**2/4*aold)**2 * dx )\n",
    "\n",
    "tend = 1\n",
    "\n",
    "scene = Draw (unew, mesh, \"disp\", deformation=unew)\n",
    "\n",
    "t = 0\n",
    "unew.Set( (0,0) )\n",
    "vel.Set( (0,0) )\n",
    "contact = ContactBoundary(mesh.Boundaries(\"contact|balls\"), mesh.Boundaries(\"contact|balls\"))\n",
    "\n",
    "X = CoefficientFunction((x,y))\n",
    "\n",
    "if True:\n",
    "    cf = (X + u-uold - (X.Other() + u.Other() - uold.Other())) * (-specialcf.normal(2).Other())\n",
    "    # cf = (X + u-uold)*specialcf.normal(2) + \\\n",
    "    #    (X.Other() + u.Other() - uold.Other()) * specialcf.normal(2).Other()\n",
    "    contact.AddEnergy(IfPos(cf, 1e9*cf*cf, 0), deformed=True)\n",
    "else:\n",
    "    cf = (X + u - (X.Other() + u.Other())) * contact.normal\n",
    "    contact.AddEnergy(IfPos(cf, 1e9*cf*cf, 0), deformed=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    while t < tend:\n",
    "        print (\"time\", t)\n",
    "        t += tau\n",
    "        uold.vec.data = unew.vec\n",
    "        aold.vec.data = anew.vec\n",
    "\n",
    "        contact.Update(uold, bfmstar, 5, 0.01)\n",
    "        NewtonMinimization (a=bfmstar, u=unew, printing=False, inverse=\"sparsecholesky\")\n",
    "\n",
    "        anew.vec.data = unew.vec-uold.vec-tau*vel.vec-tau**2/4*aold.vec\n",
    "        anew.vec.data *= 4/tau**2\n",
    "        vel.vec.data += 0.5*tau*aold.vec\n",
    "        vel.vec.data += 0.5*tau*anew.vec\n",
    "\n",
    "        scene.Redraw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "## A 3D example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "from ngsolve import *\n",
    "from ngsolve.solvers import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.webgui import Draw as DrawGeo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "box = Box((0,0,0), (1,1,0.5))\n",
    "rod = Cylinder((0.5, 0.5, -0.5),Z,0.2,1.5)\n",
    "rod.faces.name=\"contact1\"\n",
    "block = box-rod\n",
    "block.faces.Min(X).name=\"fix\"\n",
    "rod.faces.name=\"contact2\"\n",
    "rod.faces.Max(Z).name=\"force\"\n",
    "rod.faces.Min(Z).name=\"free\"\n",
    "geo = Compound([block, rod])\n",
    "DrawGeo (geo);\n",
    "mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.1)).Curve(3)\n",
    "Draw(mesh);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.GetBoundaries()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = VectorH1(mesh, order=3, dirichlet=\"fix\")\n",
    "u = fes.TrialFunction()\n",
    "gfu = GridFunction(fes)\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += Variation(InnerProduct(Sym(Grad(u)), Sym(Grad(u)))*dx)\n",
    "a += Variation(1e-8*u*u*dx)   # regularization\n",
    "a += Variation(-0.01*u[0]*ds(\"force\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "contact = ContactBoundary(mesh.Boundaries(\"contact1\"), mesh.Boundaries(\"contact2\"))\n",
    "X = CoefficientFunction((x,y,z))\n",
    "cf = (X + u - (X.Other() + u.Other())) * specialcf.normal(3)\n",
    "contact.AddEnergy(IfPos(cf, 1e4*cf*cf, 0))\n",
    "contact.Update(gfu, a, 20, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager(pajetrace=10**9):\n",
    "    NewtonMinimization(a=a, u=gfu, printing=True, inverse=\"sparsecholesky\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (gfu, mesh, deformation=True);\n",
    "Draw (Sym(Grad(gfu))[0,0], mesh, draw_surf=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}