{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# 6.3 Elasto-Plasticity\n",
    "\n",
    "\n",
    "This tutorial considers Hencky-type [1] plasticity with isotropic hardening in 2d (plane strain).\n",
    "It showcases the classes `IntegrationRuleSpace`, `NewtonCF` and `MinimizationCF`.\n",
    "\n",
    "We begin with a classical formulation with an explicit yield surface that is probably more familiar to engineers. \n",
    "In a second stage, we consider a (constrained) minimization formulation which leads to a more streamlined implementation. Both formulations correspond to the isotropic hardening model discussed in [2].\n",
    "\n",
    "The essential model ingredients are given in terms of the dissipation potential $\\Phi$ and the stored energy density contribution $\\Psi$, see eg, [3]. For convenience, we first introduce the strain energy density $\\Psi^\\text{e}$\n",
    "\n",
    "\\begin{align*}\n",
    "\\Psi^\\text{e}(\\epsilon) = \\frac{1}{2}\\| \\epsilon\\|^2_{M} &= \\frac{E}{2(1+\\nu)(1-2\\nu)}\\left((1-2\\nu)\\epsilon:\\epsilon + \\nu \\,\\text{tr}(\\epsilon)^2\\right) \\\\\n",
    "&= \\mu\\left(\\epsilon:\\epsilon\\right) + \\frac{\\lambda}{2}\\,\\text{tr}(\\epsilon)^2 \\\\\n",
    "&= \\frac{1}{2} \\epsilon : \\mathbf{C} : \\epsilon.\n",
    "\\end{align*}\n",
    "\n",
    "Moreover, we consider the full strain $\\epsilon$ to be the sum of the elastic $\\epsilon^e$ and the plastic strain contribution $p$.\n",
    "Conversely, the elastic strain is the given as \n",
    "\\begin{equation} \\epsilon^\\text{e} = \\epsilon - p ,\\end{equation}\n",
    "which is the strain quantity that will be actually employed as argument to $\\Psi^\\text{e}$.\n",
    "\n",
    "The isotropic hardening variable will be denoted by $\\alpha$ throughout this tutorial. The hardening potential\n",
    "density is simply chosen as\n",
    "$\\Psi^\\text{h} = \\frac{1}{2}\\alpha^2$ such that the total energy density reads\n",
    "\\begin{equation}\n",
    "\\Psi(\\epsilon, p, \\alpha) = \\Psi^\\text{e}(\\epsilon-p) + \\Psi^\\text{h}(\\alpha) = (\\epsilon - p) : \\mathbf{C} : (\\epsilon - p) + \\frac{1}{2} \\alpha^2.\n",
    "\\end{equation}\n",
    "\n",
    "The yield function $f$ that represents the elastic limit is given as\n",
    "\n",
    "\\begin{equation}\n",
    "f(\\sigma, \\beta) = \\left|\\mathrm{dev} \\sigma\\right| - \\sqrt{2/3}\\sigma_Y(1 - H \\beta),\n",
    "\\end{equation}\n",
    "\n",
    "where we note that the hardening parameter $H$ enters in $f$ rather than in $\\Psi$ for consistency with [2]. Deviating from [2], we scale the yield stress $\\sigma_Y$ by $\\sqrt{2/3}$ to obtain the von-Mises yield surface radius. We remark that by virtue of $\\mathrm{dev}\\sigma$ appearing in $f$, the plastic strain (rate) $p$ ($\\dot p$) will be deviatoric.\n",
    "\n",
    "The boundary value problem solved in this tutorial is in plain strains. However, this does not mean that the plastic\n",
    "strains are plain and neither are the elastic ones. Only their sum has vanishing out-of-plane components.\n",
    "Therefore, for simplicity, the material model is implemented in 3d whereas displacements are in 2d.\n",
    "\n",
    "\n",
    "## Formulation with explicit yield surface\n",
    "\n",
    "Based on the principle of maximum dissipation we can abstractly formulate the dissipation potential as\n",
    "\n",
    "\\begin{align*}\n",
    "\\Phi(\\dot{p}, \\dot{\\alpha}) = \n",
    "  \\sup_{\\sigma} \\sup_{\\beta} \\inf_{\\Lambda \\geq 0} \n",
    "    \\: \\sigma : \\dot{p} + \\beta \\dot{\\alpha} - \\Lambda f(\\sigma, \\beta)\n",
    "\\end{align*}\n",
    "\n",
    "The stationary conditions of the above functional read as\n",
    "\n",
    "\\begin{align*} \n",
    "  &\\dot{p} = \\Lambda \\frac{\\partial \\left|\\mathrm{dev} \\sigma\\right|}{\\partial \\sigma} \\\\\n",
    "  &\\dot{\\alpha} = \\Lambda\\,\\sqrt{2/3}\\sigma_YH \\\\\n",
    "  &\\{\\Lambda \\geq 0 \\: \\land \\: (\\left|\\mathrm{dev} \\sigma\\right| - \\sqrt{2/3}\\sigma_Y(1 - H\\beta)) \\leq 0 \\: \\land \\:\n",
    "  \\Lambda(\\left|\\mathrm{dev} \\sigma\\right| - \\sqrt{2/3}\\sigma_Y(1 - H\\beta)) = 0\\}\n",
    "\\end{align*}\n",
    "\n",
    "which are regarded as evolution equations for $p$ and $\\alpha$. The actual values of $\\sigma$ and $\\beta$ are\n",
    "obtained through the constitutive relations provided by the framework of generalized standard materials, we have\n",
    "\n",
    "\\begin{align*}\n",
    "&\\frac{\\partial\\Psi((\\epsilon - p), \\alpha)}{\\partial p} + \\frac{\\partial\\Phi}{\\partial \\dot{p}} = 0 \\quad \\Rightarrow \\quad \\sigma = \\mathbf{C} : (\\epsilon - p) \\\\\n",
    "&\\frac{\\partial\\Psi((\\epsilon - p), \\alpha)}{\\partial \\alpha} + \\frac{\\partial\\Phi}{\\partial \\dot{\\alpha}} = 0\n",
    "\\quad \\Rightarrow \\quad \\beta = - \\alpha .\n",
    "\\end{align*}\n",
    "\n",
    "The time-discrete form is simply obtained by replacing the time rates with the approximations \n",
    "$\\dot p \\approx \\frac{p_{k+1} - p_k}{\\Delta t}$ and $\\dot \\alpha \\approx \\frac{\\alpha_{k+1} - \\alpha_k}{\\Delta t}$,\n",
    "whereby the value of $\\Delta t$ has no effect for the present rate-independent model.\n",
    "\n",
    "\n",
    "\n",
    "## Minimization formulation\n",
    "\n",
    "\n",
    "A semi-implicit time-discrete formulation of the model including hardening reads [2]:  \n",
    "For the given state $(u^k,p^k, \\alpha^k)$ find a new solution $(u^{k+1},p^{k+1})$ such that\n",
    "\\begin{align*}\n",
    "\\text{1.} & \\: && \\int\\Psi(\\epsilon(u^{k+1}), p^{k+1},p^{k}, \\alpha^{k})+\\Phi(p^{k+1},p^k)- f^{k+1}\\cdot \n",
    "u^{k+1}\\,dx \\rightarrow \\min \\\\\n",
    "\\text{2.} & \\: && \\alpha^{k+1}=\\alpha^k+\\sqrt{2/3}\\sigma_y H \\|p^{k+1}-p^k\\| ,\n",
    "\\end{align*}\n",
    "with\n",
    "\\begin{align*}\n",
    "&\\Psi(\\epsilon(u), p, p^k, \\alpha^k) = \\frac{1}{2} \\|\\epsilon(u)-p\\|_{M}^2 \n",
    "    + \\frac{1}{2} \\left(\\alpha^k + \\sqrt{2/3}\\sigma_y H \\|p-p^k\\|_{\\varepsilon}\\right)^2 \\quad\\text{and}\\quad\n",
    "\\\\\n",
    "&\\Phi(p,p^k) = \\sqrt{2/3}\\sigma_Y \\|p-p^k\\|_{\\varepsilon},\n",
    "\\end{align*}\n",
    "\n",
    "where $\\Psi$ denotes the augmented stored energy density and $\\|\\bullet \\|_{\\varepsilon}$ indicates a perturbed norm in order to avoid divisions by zero in the evaluation of the time-discrete evolution equation derived from an incremental variational principle.\n",
    "\n",
    "\n",
    "## Notes in the spatial discretization\n",
    "\n",
    "The state $u$ is spatially discretized by continuous Lagrange elements, the internal states $p$ and $\\alpha$ as well as the multiplier $\\Lambda$ reside at quadrature points, for which NGSolve has the notion of an integration rule space. Thus, the problem above can be decomposed into a \"global\" problem for the displacements and local problems for the internal variables. The latter can be solved individually for each quadrature point. However, there is a coupling between the local problems and displacements, which has to be accounted for by what is called \"algorithmically consistent linearization\". The latter is automatically obtained by using an appropriate\n",
    "`BilinearForm` object representing the boundary value problem and the internal evolution problems as demonstrated for both implementations.\n",
    "\n",
    "\n",
    "**References**\n",
    " 1. [W. Han and B.D. Reddy: Mathematical Theory and Numerical Analysis, Springer 1999](https://www.springer.com/gp/book/9781461459392)\n",
    " 2. [C. Carstensen, Domain Decomposition for a Non-smooth Convex Minimization Problem and its Application in Plasticity, 1997](https://doi.org/10.1002/(SICI)1099-1506(199705/06)4:3<177::AID-NLA106>3.0.CO;2-B)\n",
    " 3. [K. Hackl, F.D. Fischer, On the relation between the principle of maximum dissipation and inelastic evolution given by dissipation potentials, PRSA, 2008](https://doi.org/10.1098/rspa.2007.0086)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "\n",
    "\n",
    "### General setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from ngsolve import *\n",
    "from ngsolve.comp import IntegrationRuleSpace\n",
    "from ngsolve.fem import MinimizationCF, NewtonCF\n",
    "from netgen.geom2d import CSG2d, Circle, Rectangle\n",
    "from ngsolve.webgui import Draw\n",
    "SetNumThreads(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "The problem domain is the classical plate with hole. Thanks to symmetry of the problem to be solved, we actually use only one quarter of the domain together with additional symmetry boundary conditions.\n",
    "Parameters for the geometry and the material employed are taken from [4].\n",
    "\n",
    "\n",
    "4. [A. Düster and E. Rank: The p-version of the fnite element method compared to an\n",
    "adaptive h-version for the deformation theory of plasticity, CMAME 190 (2001) 1925-1935](https://doi.org/10.1016/S0045-7825(00)00215-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# expression (jit) compilation flag\n",
    "realcompile = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# polynomial order for geometry and displacements\n",
    "order = 3\n",
    "\n",
    "geo = CSG2d()\n",
    "circle = Circle(center=(100,100), radius=10.0, bc=\"curve\").Maxh(1)\n",
    "rect = Rectangle(pmin=(0,100), pmax=(100,200), bottom=\"bottom\", left=\"left\", top=\"top\", right=\"right\")\n",
    "geo.Add(rect-circle)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=5))\n",
    "mesh.Curve(order)\n",
    "Draw(mesh)\n",
    "\n",
    "# points of interest for post-processing\n",
    "node_A = mesh(100,200)\n",
    "node_B = mesh(0,200)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "### Strain energy density and generic helper functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "#### Model parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Elasticity parameters correspond to those in Ref. 4, ie. mu = 80193.8 and kappa = 164206\n",
    "\n",
    "# Young's modulus\n",
    "Emod = 206900\n",
    "\n",
    "# Poisson's ratio\n",
    "nu = 0.29\n",
    "\n",
    "print(\"E =\", Emod, \"\\n\\nnu =\", nu)\n",
    "\n",
    "mu = Parameter(Emod / (2 * (1 + nu)))\n",
    "lmbda = Parameter(Emod * nu / ((1 + nu) * (1 - 2 * nu)))\n",
    "\n",
    "print()\n",
    "print(\"mu =\", mu.Get())\n",
    "print(\"lambda =\", lmbda.Get())\n",
    "\n",
    "# just for reference:\n",
    "kappa = Emod / (3 * (1 - 2 * nu))\n",
    "\n",
    "print(\"kappa =\", kappa)\n",
    "\n",
    "\n",
    "# Hardening parameter\n",
    "H = Parameter(1.0)\n",
    "\n",
    "# Yield stress\n",
    "sigma_Y = Parameter(450)\n",
    "\n",
    "# Perturbation parameter for tensor norms: this value is decisive for the accuracy of the result!\n",
    "#  for higher values of H it must be sufficently small\n",
    "norm_pert = 1e-16"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "#### Function definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "def strain(u):\n",
    "    # compute symmetric gradient (here a 2x2 block) and \"embed\" it in a 3x3 matrix via \"ExtendDimension\"\n",
    "    return Sym(Grad(u)).ExtendDimension((3, 3), pos=(0, 0))\n",
    "\n",
    "\n",
    "def elastic_strain_energy_mu_lambda(eps, mu, lmbda):\n",
    "    return mu * InnerProduct(eps, eps) + lmbda/2 * Trace(eps)**2\n",
    "\n",
    "\n",
    "def elastic_strain_energy(eps):\n",
    "    return elastic_strain_energy_mu_lambda(eps, mu, lmbda)\n",
    "\n",
    "\n",
    "def tensor_norm(a, pert=0):\n",
    "    return sqrt(InnerProduct(a, a) + pert)\n",
    "\n",
    "\n",
    "def dev(a):\n",
    "    return a - Trace(a) * Id(3) / 3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "### Common FE spaces\n",
    "\n",
    "We employ a $H^1$ space for the displacements and a \"integration rule space\" (IR space) as basis space for all internal variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes_u = VectorH1(mesh, order=order, dirichletx=\"right\", dirichlety=\"bottom\")\n",
    "fes_ir = IntegrationRuleSpace(mesh, order=order-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "Scalar internal states can be defined based on `fes_ir` as needed.\n",
    "\n",
    "The \"measure\" corresponding to `fes_ir`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract integration rules from IntegrationRuleSpace\n",
    "irs_dx = dx(intrules=fes_ir.GetIntegrationRules())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "### Load definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "force = 450\n",
    "# For loadsteps\n",
    "loadfactor = Parameter(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17",
   "metadata": {},
   "source": [
    "### Some data structures for postprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for postprocessing\n",
    "drawfes = MatrixValued(L2(mesh, order=order-1), dim=3, symmetric=True, deviatoric=False)\n",
    "pd, qd = drawfes.TnT()\n",
    "pdraw = GridFunction(drawfes)\n",
    "ad = BilinearForm(InnerProduct(pd, qd)*irs_dx, symmetric=True).Assemble()\n",
    "invad = ad.mat.Inverse()\n",
    "\n",
    "drawafes = L2(mesh, order=order-1)\n",
    "pda, qda = drawafes.TnT()\n",
    "adraw = GridFunction(drawafes)\n",
    "ada = BilinearForm(InnerProduct(pda, qda)*irs_dx, symmetric=True).Assemble()\n",
    "invada = ada.mat.Inverse()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19",
   "metadata": {},
   "source": [
    "### Implementation with explicit yield surface\n",
    "\n",
    "In this approach we employ a generic symmetric-matrix-valued space for $p$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes_p = MatrixValued(fes_ir, symmetric=True, deviatoric=False, dim=3)\n",
    "\n",
    "fes_int = fes_p * fes_ir * fes_ir # p x alpha x Lambda\n",
    "\n",
    "# Trial and test functions\n",
    "u, u_test = fes_u.TnT()\n",
    "int_trial, int_test = fes_int.TnT()\n",
    "p, alpha, Lambda = int_trial\n",
    "\n",
    "# GridFunction for \"external\" state\n",
    "gfu = GridFunction(fes_u)\n",
    "\n",
    "# GridFunction for internal states\n",
    "gfint = GridFunction(fes_int)\n",
    "gfp, gfalpha, gfLambda = gfint.components\n",
    "\n",
    "# For history states\n",
    "gfhist = GridFunction(fes_int)\n",
    "gfp_k, gfalpha_k, gfLambda_k = gfhist.components\n",
    "\n",
    "# \"Trial\" states for plastic evolution (these are not trial functions in FE parlance but elastic predictors)\n",
    "gftrial = GridFunction(fes_int)\n",
    "gfsigma_trial, gfbeta_trial, _ = gftrial.components\n",
    "\n",
    "# Time increment\n",
    "Delta_t = 1.0\n",
    "# Not really needed for a rate-independent model. Value could be used for scaling the equations though.\n",
    "\n",
    "# Short-hand for storing history variables. Note that the history value of Lambda is not relevant and thus omitted.\n",
    "def store_internal():\n",
    "    gfp_k.Interpolate(gfp)\n",
    "    gfalpha_k.Interpolate(gfalpha)\n",
    "    gfLambda_k.Interpolate(gfLambda)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21",
   "metadata": {},
   "source": [
    "#### Functions for the energy density $\\Psi$, the dissipation potential $\\Phi$ and the evolution equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Energy density\n",
    "def Psi(strain, p, alpha):\n",
    "    strain_energy = elastic_strain_energy(strain-p)\n",
    "    hardening     = 1/2 * alpha**2\n",
    "    return strain_energy + hardening\n",
    "\n",
    "\n",
    "# The yield condition: result > 0 ==> plastic evolution\n",
    "def yield_condition(sigma, beta, pert=0):\n",
    "    return tensor_norm(dev(sigma), pert=pert) - sqrt(2/3) * sigma_Y * (1 - H*beta)\n",
    "\n",
    "\n",
    "# The objective function defining the dissipation potential\n",
    "def Phi_obj(sigma, beta, p_dot, alpha_dot, Lambda):\n",
    "    return InnerProduct(sigma, p_dot) + InnerProduct(beta, alpha_dot) - \\\n",
    "            Lambda * yield_condition(sigma, beta, pert=norm_pert)\n",
    "\n",
    "\n",
    "# The evolution equations\n",
    "def evolution_eqs(strain, p, alpha, Lambda, evol_tol):\n",
    "    p_dot = (p - gfp_k) / Delta_t\n",
    "    alpha_dot = (alpha - gfalpha_k) / Delta_t\n",
    "    sigma = -Psi(strain, p, alpha).Diff(p)\n",
    "    sigma = sigma.MakeVariable()\n",
    "    beta = -Psi(strain, p, alpha).Diff(alpha)\n",
    "    beta = beta.MakeVariable()\n",
    "    _Phi_obj = Phi_obj(sigma, beta, p_dot, alpha_dot, Lambda)\n",
    "    \n",
    "    dsigma = _Phi_obj.Diff(sigma)\n",
    "    \n",
    "    # only active when there is hardening\n",
    "    dbeta = IfPos(H - 1e-16, _Phi_obj.Diff(beta), alpha_dot)\n",
    "    \n",
    "    dLambda = _Phi_obj.Diff(Lambda)\n",
    "\n",
    "    # Only solve the evolution equations if the elastic trial state is beyond the yield surface.\n",
    "    # This essentially leads to the fulfillment of Lambda >= 0\n",
    "    return IfPos(yield_condition(gfsigma_trial, gfbeta_trial) - evol_tol, \n",
    "                 CF((dsigma, dbeta, dLambda)),\n",
    "                 CF((p_dot, alpha_dot, Lambda - gfLambda_k)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23",
   "metadata": {},
   "source": [
    "#### Handling of the evolution of internal variables. \n",
    "\n",
    "The evolution equations are a system of nonlinear equations per quadrature point.\n",
    "`NewtonCF` solves such equations with given starting points, tolerance and maximum iterations.\n",
    "On failure, it returns `NaN`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters for evolution\n",
    "tol = Parameter(1e-10)\n",
    "maxiter = 20\n",
    "evol_tol = tol\n",
    "\n",
    "eqs = evolution_eqs(strain(gfu), p, alpha, Lambda, evol_tol).Compile(realcompile=realcompile)\n",
    "evolution = NewtonCF(eqs, gfint, tol=tol.Get(), maxiter=maxiter)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25",
   "metadata": {},
   "source": [
    "The \"trial state\" that indicates whether there will be a plastic evolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26",
   "metadata": {},
   "outputs": [],
   "source": [
    "Psi_trial = Psi(strain(gfu), *gfhist.components[:2])\n",
    "def compute_trial_state():\n",
    "    gftrial.components[0].Interpolate(-Psi_trial.Diff(gfhist.components[0]))\n",
    "    gftrial.components[1].Interpolate(-Psi_trial.Diff(gfhist.components[1]))\n",
    "    # no need to interpolate Lambda"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27",
   "metadata": {},
   "source": [
    "The actual solution of the equations happens during an again exact (interpolating) projection that sets `gfint`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A Linearform for which the evaluation triggers the solution of the evolution equations \n",
    "testfunc = CoefficientFunction(tuple(fes_int.TestFunction()))\n",
    "L_evol = LinearForm(fes_int)\n",
    "L_evol += InnerProduct(CacheCF(evolution), testfunc) * irs_dx\n",
    "# NOTE: CacheCF(evolution) ensures that the evolution equations are solved only once per evaluation\n",
    "\n",
    "# Mass-Matrix for the internal variable space as generic interpolation tool\n",
    "M_evol = BilinearForm(fes_int, symmetric=True, diagonal=True)\n",
    "for _trial, _test in zip(*fes_int.TnT()):\n",
    "    M_evol += InnerProduct(_trial,  _test).Compile() * irs_dx\n",
    "    \n",
    "M_evol.Assemble()\n",
    "M_evol_inv = M_evol.mat.Inverse()\n",
    "\n",
    "\n",
    "def evolve_internal():\n",
    "    compute_trial_state()\n",
    "    L_evol.Assemble()\n",
    "    gfint.vec.data = M_evol_inv * L_evol.vec"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29",
   "metadata": {},
   "source": [
    "#### The `LinearForm` for the global boundary value problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a strain variable for partial differentiation\n",
    "strain_gf = strain(gfu).MakeVariable()\n",
    "\n",
    "# stress\n",
    "sigma = Psi(strain_gf, gfp, gfalpha).Diff(strain_gf)\n",
    "\n",
    "L = LinearForm(fes_u)\n",
    "L += InnerProduct(sigma, strain(u_test)).Compile(realcompile=realcompile) * irs_dx\n",
    "L += -u_test[1] * force * loadfactor * ds(\"top\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31",
   "metadata": {},
   "source": [
    "#### The \"elastic\" part of the  `BilinearForm `"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32",
   "metadata": {},
   "outputs": [],
   "source": [
    "# elasticity tensor\n",
    "cc = sigma.Diff(strain_gf).Reshape((strain_gf.dim, strain_gf.dim))\n",
    "\n",
    "# reshape test and trial strians\n",
    "delta_eps = strain(u_test).Reshape((strain_gf.dim, ))\n",
    "Delta_eps = strain(u).Reshape((strain_gf.dim, ))\n",
    "\n",
    "a = BilinearForm(fes_u, symmetric=False)\n",
    "\n",
    "# \"elastic contribution\"\n",
    "a += InnerProduct(delta_eps, cc * Delta_eps).Compile(realcompile=realcompile) * irs_dx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33",
   "metadata": {},
   "source": [
    "#### The \"algorithmic\" part of the  `BilinearForm `"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34",
   "metadata": {},
   "source": [
    "An algorithmically consistent linearization also needs to account for the linearization of internal states.\n",
    "\n",
    "The first ingredients are the derivatives of the stress with respect to internal states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Partial derivates w.r.t. internal states; need to be reshaped in order to put them together \n",
    "#  in a single expression below.\n",
    "\n",
    "dim_sigma = sigma.dim\n",
    "dp_sigma = sigma.Diff(gfp).Reshape((dim_sigma, gfp.dim))\n",
    "dalpha_sigma = sigma.Diff(gfalpha).Reshape((dim_sigma, 1))\n",
    "dLambda_sigma = sigma.Diff(gfLambda).Reshape((dim_sigma, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36",
   "metadata": {},
   "source": [
    "Next, we collect these derivative in a single object with dimensions\n",
    "`(dim_sigma, dim_int)`. Using the method `ExtendDimension(...)` we can define\n",
    "individual contributions/blocks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37",
   "metadata": {},
   "outputs": [],
   "source": [
    "dim_int = gfp.dim + gfalpha.dim + gfLambda.dim\n",
    "\n",
    "dint_sigma = dp_sigma.ExtendDimension((dim_sigma, dim_int), pos=(0,0)) \\\n",
    " + dalpha_sigma.ExtendDimension((dim_sigma, dim_int), pos=(0, gfp.dim)) \\\n",
    " + dLambda_sigma.ExtendDimension((dim_sigma, dim_int), pos=(0, gfp.dim + gfalpha.dim))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38",
   "metadata": {},
   "source": [
    "In order to obtain the sensitivities of the internal states with respect to the strain, \n",
    "we need to define the partial deriatives of the evolution equations with respect to\n",
    "\"external\" (strain) and \"internal\" (p, alpha, Lamda) states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39",
   "metadata": {},
   "outputs": [],
   "source": [
    "# instance of evolution equations based only on grid functions for direct evaluation\n",
    "eqs_gf = evolution_eqs(strain_gf, gfp, gfalpha, gfLambda, evol_tol)\n",
    "\n",
    "dim_eqs = eqs_gf.dim\n",
    "\n",
    "# Partial derivate w.r.t. external (strain) state\n",
    "dext_eqs_gf = eqs_gf.Diff(strain_gf).Reshape((dim_eqs, strain_gf.dim))\n",
    "\n",
    "# Partial derivates w.r.t. internal states; need to be reshaped in order to put them together \n",
    "#  in a single expression below.\n",
    "dp_eqs_gf = eqs_gf.Diff(gfp).Reshape((dim_eqs, gfp.dim))\n",
    "dalpha_eqs_gf = eqs_gf.Diff(gfalpha).Reshape((dim_eqs, 1))\n",
    "dLambda_eqs_gf = eqs_gf.Diff(gfLambda).Reshape((dim_eqs, 1))\n",
    "\n",
    "dint_eqs_gf = \\\n",
    "   dp_eqs_gf.ExtendDimension((dim_eqs, dim_int), pos=(0,0)) \\\n",
    " + dalpha_eqs_gf.ExtendDimension((dim_eqs, dim_int), pos=(0, gfp.dim)) \\\n",
    " + dLambda_eqs_gf.ExtendDimension((dim_eqs, dim_int), pos=(0, gfp.dim + gfalpha.dim))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40",
   "metadata": {},
   "source": [
    "The matrix represented by `dint_eqs` is singular, because `p`\n",
    "is a symmetric deviatoric (trace-less) tensor, it has a \n",
    "non-trivial embedding into a vector space for its *independent* components.\n",
    "Therefore, we need to create such an embedding for the space of internal\n",
    "variables as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41",
   "metadata": {},
   "outputs": [],
   "source": [
    "# VS embedding\n",
    "\n",
    "# Retreive non-trivial embedding for p\n",
    "VS_p = gfp.space.VSEmbedding().NumPy()\n",
    "\n",
    "# Create a common embedding for all internal states\n",
    "VS_all = np.zeros((dim_eqs, VS_p.shape[1] + 2))\n",
    "VS_all[:VS_p.shape[0], :VS_p.shape[1]] = VS_p\n",
    "VS_all[VS_p.shape[0], VS_p.shape[1]] = 1\n",
    "VS_all[VS_p.shape[0] + 1, VS_p.shape[1] + 1] = 1\n",
    "\n",
    "# Create a CoefficienFunction from the array\n",
    "VS_cf = CF(tuple(VS_all.flatten().tolist()), dims=VS_all.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42",
   "metadata": {},
   "source": [
    "With this embedding we project `dint_eqs_gf` to a reduced space and obtain the invertible\n",
    "matrix `A_red`. After inversion we can compute the sensitivity of the internal states with respect to the external (strain) state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43",
   "metadata": {},
   "outputs": [],
   "source": [
    "A_red = VS_cf.trans * dint_eqs_gf * VS_cf\n",
    "A_red_inv = Inv(A_red)\n",
    "A_inv = VS_cf * A_red_inv * VS_cf.trans\n",
    "\n",
    "# sensitivity of internal states wrt to external states, i.e. strain\n",
    "dext_int = -A_inv * dext_eqs_gf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44",
   "metadata": {},
   "source": [
    "Finally, the complete linearization term is given as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45",
   "metadata": {},
   "outputs": [],
   "source": [
    "linearization_term = delta_eps * (dint_sigma * (dext_int * Delta_eps))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46",
   "metadata": {},
   "outputs": [],
   "source": [
    "a += linearization_term.Compile(realcompile=realcompile) * irs_dx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47",
   "metadata": {},
   "source": [
    "#### Solution procedure"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48",
   "metadata": {},
   "source": [
    "The nonlinear PDE is solved with a Newton-Raphson scheme. Thereby, in each Newton iteration\n",
    "we solve a linearized PDE for the displacements and a nonlinear evolution problem at each quadrature point via *`NewtonCF`*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49",
   "metadata": {},
   "source": [
    "In order to run the problem, we need to define vectors that store respective solution data, define loadsteps\n",
    "and setup a loop that iterates through the latter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A set of functions/vectors for keeping data\n",
    "w = gfu.vec.CreateVector()\n",
    "rhs = gfu.vec.CreateVector()\n",
    "\n",
    "# Postprocess internal states\n",
    "fd = LinearForm(InnerProduct(gfp, qd)*irs_dx)\n",
    "fda = LinearForm(InnerProduct(gfalpha_k, qda)*irs_dx)\n",
    "\n",
    "# Load steps (chosen based on experience)\n",
    "loadsteps = [0.1,0.3,0.5,0.7,0.8,0.9,0.95,1.0]\n",
    "\n",
    "# Set solution to zero initially\n",
    "gfu.vec[:] = 0\n",
    "gfint.vec[:] = 0\n",
    "gfhist.vec[:] = 0\n",
    "gftrial.vec[:] = 0\n",
    "\n",
    "\n",
    "# Iterate through load steps\n",
    "for ls in loadsteps:\n",
    "    loadfactor.Set(ls)\n",
    "        \n",
    "    # Update old solution at time t = t_k\n",
    "    store_internal()\n",
    "    \n",
    "    with TaskManager():\n",
    "        for i in range(20):\n",
    "            \n",
    "            a.Assemble()\n",
    "            L.Assemble()\n",
    "            \n",
    "            w.data = a.mat.Inverse(freedofs=fes_u.FreeDofs(False), inverse=\"umfpack\") * L.vec\n",
    "            gfu.vec.data -= w\n",
    "            \n",
    "            # Solve the evolution equation (quadrature-point-wise)\n",
    "            evolve_internal()\n",
    "                        \n",
    "            if np.isnan(np.max(np.abs(gfp.vec.FV().NumPy()))):\n",
    "                raise Exception(\"Evolution solver failed\")\n",
    "            \n",
    "            # NOTE: This problem is not a minimization problem. \n",
    "            # Therefore, we compute || |W| . |R| || instead of || W . R ||.\n",
    "            err = np.linalg.norm(np.abs(w.FV().NumPy()) * np.abs(L.vec.FV().NumPy()))\n",
    "            print(\"step \", i, \"err = \", err)\n",
    "            \n",
    "            # Check convergence\n",
    "            if err < 1e-6: break\n",
    "        \n",
    "    print(\"force = \", ls * force, \", uy_A =\", gfu(node_A)[1], \", ux_B =\", gfu(node_B)[0],\\\n",
    "          \", int u2 =\", Integrate(gfu[1] * ds(\"top\"),mesh))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51",
   "metadata": {},
   "outputs": [],
   "source": [
    "fd.Assemble()\n",
    "pdraw.vec.data = invad * fd.vec\n",
    "Draw(Norm(pdraw), mesh, \"p\")\n",
    "Draw(gfu,mesh, \"u\")\n",
    "\n",
    "if H.Get() > 1e-16:\n",
    "    fda.Assemble()\n",
    "    adraw.vec.data = invada * fda.vec\n",
    "    Draw(Norm(adraw), mesh, \"alpha\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52",
   "metadata": {},
   "outputs": [],
   "source": [
    "# confirm that the trace vanishes\n",
    "p1 = np.array(pdraw(mesh(90, 100, 0))).reshape((3,3))\n",
    "np.einsum('ij,ij', p1, np.eye(3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "54",
   "metadata": {},
   "source": [
    "### Implementation of the minimization formulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55",
   "metadata": {},
   "source": [
    "We again start with the (coupled) function spaces. This time, we won't need a `TrialFunction` \"variable\" for $\\alpha$ but only a `GridFunction` to store its value. Also, there is no need for $\\Lambda$ in this formulation.\n",
    "\n",
    "Moreover, we use a *deviatoric* space for $p$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes_p = MatrixValued(fes_ir, symmetric=True, deviatoric=True, dim=3)\n",
    "\n",
    "fes = fes_u * fes_p\n",
    "u, p = fes.TrialFunction()\n",
    "\n",
    "# GridFunction for solution\n",
    "gfsol = GridFunction(fes)\n",
    "gfu, gfp = gfsol.components\n",
    "\n",
    "# Save previous solution\n",
    "gfsol_k = GridFunction(fes)\n",
    "gfu_k, gfp_k = gfsol_k.components\n",
    "\n",
    "# alpha_0 = 0\n",
    "alpha_k = GridFunction(fes_ir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57",
   "metadata": {},
   "source": [
    "The energy density $\\Psi$ and the dissipation potential $\\Phi$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The full energy density (incl. hardening)\n",
    "def Psi(p, strain, p_k, alpha_k):\n",
    "    strain_energy = elastic_strain_energy(strain-p)\n",
    "    alpha = alpha_k + IfPos(H - 1e-16, sqrt(2/3) * sigma_Y * H * tensor_norm(p - p_k, pert=norm_pert), 0)\n",
    "    hardening    = 1/2 * alpha**2\n",
    "    return strain_energy + hardening\n",
    "\n",
    "\n",
    "# Dissipation potential: For rate-independent plasticity (homogeneous of degree 1 in the rate of p) \n",
    "# this is the dissipation! Also, because of rate independence, the parameter Delta_t is not needed for\n",
    "# the evolution.\n",
    "def Phi(p, p_k):\n",
    "    return sqrt(2/3) * sigma_Y * tensor_norm(p - p_k, pert=norm_pert)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59",
   "metadata": {},
   "source": [
    "Finally, we have all ingredients for the variational formulation:\n",
    "\n",
    "Find for given $(u^k,p^k)$ a solution $(u^{k+1},p^{k+1})$ such that\n",
    "\\begin{align*}\n",
    "\\Psi(\\epsilon(u^{k+1}), p^{k+1},p^{k}) + \\Phi(p^{k+1},p^k)-\\int f^{k+1}\\cdot u^{k+1}\\,dx \\rightarrow \\min!\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60",
   "metadata": {},
   "outputs": [],
   "source": [
    "realcompile = False\n",
    "\n",
    "# Possible optimization: use static condensation for p\n",
    "a = BilinearForm(fes, symmetric=True)\n",
    "a += Variation( Psi(p, strain(u), gfp_k, alpha_k) * irs_dx ).Compile(realcompile=realcompile)\n",
    "a += Variation( Phi(p, gfp_k) * irs_dx ).Compile(realcompile=realcompile)\n",
    "a += Variation( -u[1] * force * loadfactor * ds(\"top\") ).Compile(realcompile=realcompile)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61",
   "metadata": {},
   "source": [
    "Since Phi(p, gfp_k) is included above, we automatically obtain an \"algorithmically\n",
    "consistent linearization\", corresponding to the local evolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62",
   "metadata": {},
   "outputs": [],
   "source": [
    "evolution_objective = (Psi(p, strain(gfu), gfp_k, alpha_k) + Phi(p,gfp_k)).Compile(realcompile=realcompile)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63",
   "metadata": {},
   "source": [
    "which will be used to solve for `p` in a nested way, ie. for each \"guess\" for `u^{k+1}`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64",
   "metadata": {},
   "source": [
    "In order to run the problem, we need to define vectors that store respective solution data, define loadsteps\n",
    "and setup a loop that iterates through the latter.\n",
    "\n",
    "The nonlinear PDE is solved with a Newton-Raphson scheme. Thereby, in each Newton iteration\n",
    "we solve a linearized PDE for the displacements and a nonlinear evolution problem at each quadrature point via *`MinimizationCF`*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65",
   "metadata": {},
   "outputs": [],
   "source": [
    "# a set of vectors for keeping data\n",
    "vec_k = gfsol.vec.CreateVector()\n",
    "w      = gfsol.vec.CreateVector()\n",
    "rhs    = gfsol.vec.CreateVector()\n",
    "\n",
    "# postprocess internal states\n",
    "fd = LinearForm(InnerProduct(gfp, qd)*irs_dx)\n",
    "fda = LinearForm(InnerProduct(gfalpha_k, qda)*irs_dx)\n",
    "\n",
    "# load steps (chosen based on experience)\n",
    "loadsteps = [0.1,0.3,0.5,0.7,0.8,0.9,0.95,1.0]\n",
    "\n",
    "# set solution to zero initially\n",
    "gfsol.vec[:] = 0\n",
    "alpha_k.vec[:] = 0\n",
    "\n",
    "# evolution params (default values of MinimizationCF at the time of writing)\n",
    "tol = 1e-6\n",
    "maxiter = 20\n",
    "\n",
    "# iterate through load steps\n",
    "for ls in loadsteps:\n",
    "    loadfactor.Set(ls)\n",
    "    \n",
    "    # update old solution at time t = t_k\n",
    "    gfsol_k.vec.data = gfsol.vec\n",
    "    with TaskManager():\n",
    "        for i in range(20):\n",
    "            energy_k = a.Energy(gfsol.vec)\n",
    "            #print (\"energy(t_k) = \", energy_k)\n",
    "            \n",
    "            a.AssembleLinearization(gfsol.vec)\n",
    "            a.Apply(gfsol.vec, rhs)\n",
    "            \n",
    "            # static condensation could be employed to speedup solving\n",
    "            w.data = a.mat.Inverse(freedofs=fes.FreeDofs(False), inverse=\"sparsecholesky\") * rhs\n",
    "\n",
    "            # linesearch ( with damping)\n",
    "            vec_k.data = gfsol.vec\n",
    "            scale = 1\n",
    "            while scale > 1e-7:\n",
    "                gfsol.vec.data -= scale * w\n",
    "                energy1 = a.Energy(gfsol.vec)\n",
    "                \n",
    "                # Here we solve the evolution equation at each quadrature point through evaluation of \n",
    "                # the MinimizationCF. This function takes the \"objective\", the initial guess (most \n",
    "                # likely the previous solution), tolerances and the maximum number of iterations.\n",
    "                gfp.Interpolate(MinimizationCF(evolution_objective, gfp, tol=tol, maxiter=maxiter))\n",
    "                \n",
    "                energy2 = a.Energy(gfsol.vec)\n",
    "                \n",
    "                if energy2 < energy_k + 1e-12: \n",
    "                    break\n",
    "                scale *= 0.5\n",
    "                gfsol.vec.data = vec_k\n",
    "                #print(scale)\n",
    "                \n",
    "            err = sqrt(InnerProduct(w, rhs))\n",
    "            print(\"step \", i, \"err = \", err)\n",
    "            \n",
    "            # check convergence\n",
    "            if err < 1e-6: break\n",
    "    \n",
    "    alpha_k.Interpolate(alpha_k + sqrt(2/3) * sigma_Y * H * sqrt(InnerProduct(gfp - gfp_k, gfp - gfp_k)))\n",
    "    \n",
    "    print(\"force = \", ls * force, \", uy_A =\", gfu(node_A)[1], \", ux_B =\", gfu(node_B)[0],\\\n",
    "          \", int u2 =\", Integrate(gfu[1] * ds(\"top\"),mesh))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66",
   "metadata": {},
   "outputs": [],
   "source": [
    "fd.Assemble()\n",
    "pdraw.vec.data = invad * fd.vec\n",
    "Draw(Norm(pdraw), mesh, \"p\")\n",
    "Draw(gfu,mesh, \"u\")\n",
    "\n",
    "if H.Get() > 1e-16:\n",
    "    fda.Assemble()\n",
    "    adraw.vec.data = invada * fda.vec\n",
    "    Draw(Norm(adraw), mesh, \"alpha\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67",
   "metadata": {},
   "outputs": [],
   "source": [
    "# confirm that the trace vanishes\n",
    "p1 = np.array(pdraw(mesh(90, 100, 0))).reshape((3,3))\n",
    "np.einsum('ij,ij', p1, np.eye(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68",
   "metadata": {},
   "source": [
    "#### A note on performance\n",
    "\n",
    "For the (rather small) problem at hand, the minimization implementation is considerable faster than the\n",
    "\"classical\" implementation with explicit yield surface, dispite the fact that the former method solves\n",
    "a larger linear system than the latter. \n",
    "The main bottleneck is the assembly of the consistent linearization term. However, for larger problems, the effort for solving\n",
    "the linear system grows faster than that for the assembly. Hence, the picture will change for larger problems, in particular in three dimensions.\n",
    "A possibility to reduce the system size in case of the minimization implementation is static condensation,\n",
    "which corresponds to \"consistent linearization\" computed at the level of element matrices\n",
    "instead of quadrature points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69",
   "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.2"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "384px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}