{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving nonlinear Elasticity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "a rectangle with refinement at corners:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "shape = Rectangle(1,0.1).Face()\n",
    "shape.edges.Max(X).name=\"right\"\n",
    "shape.edges.Min(X).name=\"left\"\n",
    "shape.edges.Max(Y).name=\"top\"\n",
    "shape.edges.Min(Y).name=\"bot\"\n",
    "shape.vertices.Min(X+Y).maxh=0.01\n",
    "shape.vertices.Min(X-Y).maxh=0.01\n",
    "mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.05))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cauchy-Green tensor\n",
    "\n",
    "$$\n",
    "C = F^T F \\qquad \\text{with} \\qquad F = I + \\nabla u\n",
    "$$\n",
    "\n",
    "and hyperelastic energy density\n",
    "\n",
    "$$\n",
    "W : {\\mathbb R}^{d \\times d, sym} \\rightarrow {\\mathbb R}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E, nu = 210, 0.2\n",
    "mu  = E / 2 / (1+nu)\n",
    "lam = E * nu / ((1+nu)*(1-2*nu))\n",
    "\n",
    "def C(u):\n",
    "    F = Id(2) + Grad(u)\n",
    "    return F.trans * F\n",
    "\n",
    "def NeoHooke (C):\n",
    "    return 0.5*mu*(Trace(C-Id(2)) + 2*mu/lam*Det(C)**(-lam/2/mu)-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "stationary point of total energy:\n",
    "\n",
    "$$\n",
    "\\delta \\int W(C(u)) - f u = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factor = Parameter(0)\n",
    "force = CoefficientFunction( (0,factor) )\n",
    "\n",
    "fes = H1(mesh, order=4, dirichlet=\"left\", dim=mesh.dim)\n",
    "u  = fes.TrialFunction()\n",
    "\n",
    "a = BilinearForm(fes, symmetric=True)\n",
    "a += Variation(NeoHooke(C(u)).Compile()*dx)\n",
    "a += Variation((-InnerProduct(force,u)).Compile()*dx)\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "gfu.vec[:] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Variation` function declares that the non-linear form is the derivative of the energy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "a simple Newton solver, using automatic differentiation for residual and tangential stiffness:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveNewton(printrates=False):\n",
    "    for it in range(10):\n",
    "        if (printrates):\n",
    "            print (\"it\", it, \"energy = \", a.Energy(gfu.vec))\n",
    "        res = a.Apply(gfu.vec)\n",
    "        a.AssembleLinearization(gfu.vec)\n",
    "        inv = a.mat.Inverse(fes.FreeDofs() ) \n",
    "        gfu.vec.data -= inv*res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factor.Set(0.4)\n",
    "SolveNewton(printrates=True)\n",
    "scene = Draw (C(gfu)[0,0]-1, mesh, deformation=gfu, min=-0.1, max=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often, we don't have a good starting value for Newton's method. This can be overcome by increasing the load step by step (assuming the solution depends continuously on the loading). The solution of the previous load-step is the initial guess for the next step.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numsteps = 5\n",
    "maxload = 2\n",
    "for ls in range (numsteps):\n",
    "    factor.Set(maxload*(ls+1)/numsteps)\n",
    "    SolveNewton()\n",
    "    Draw (C(gfu)[0,0]-1, mesh, deformation=gfu, min=-0.2, max=0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stress tensor\n",
    "\n",
    "Compute $2^{nd}$ Piola Kirchhoff stress tensor by symbolic differentiation:\n",
    "\n",
    "$$\n",
    "\\Sigma_{i,j} = \\frac{\\partial W}{\\partial C_{i,j}} (C)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C_=C(gfu).MakeVariable()\n",
    "sigma = NeoHooke(C_).Diff(C_)\n",
    "\n",
    "Draw (sigma[0,0], mesh, \"Sxx\", deformation=gfu, min=-10.001, max=10.001); "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The energy functional is represented as an expression tree:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u  = fes.TrialFunction()\n",
    "print (NeoHooke(C(u)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the `Compile` method, the tree is linearized, and common sub-expressions are merged:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (NeoHooke(C(u)).Compile())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (sigma.Compile())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.13.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}