{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Nearly Incompressible Hyperelasticity\n",
    "\n",
    "![](quasi_incompressible_hyperelasticity.gif)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Introduction\n",
    "\n",
    "In this example we study quasi- or nearly-incompressible hyperelasticity using the stable Taylor-Hood approximation. In spirit, this example is the nonlinear analogue of\n",
    "[`incompressible_elasticity`](https://nbviewer.jupyter.org/github/Ferrite-FEM/Ferrite.jl/blob/gh-pages/v0.3.14/examples/incompressible_elasticity.ipynb) and the incompressible analogue of\n",
    "[`hyperelasticity`](https://nbviewer.jupyter.org/github/Ferrite-FEM/Ferrite.jl/blob/gh-pages/v0.3.14/examples/hyperelasticity.ipynb). Much of the code therefore follows from the above two examples.\n",
    "The problem is formulated in the undeformed or reference configuration with the displacement $\\mathbf{u}$ and pressure $p$ being the unknown fields. We now briefly outline\n",
    "the formulation. Consider the standard hyperelasticity problem\n",
    "\n",
    "$$\n",
    "  \\mathbf{u} = \\argmin_{\\mathbf{v}\\in\\mathcal{K}(\\Omega)}\\Pi(\\mathbf{v}),\\quad \\text{where}\\quad \\Pi(\\mathbf{v})  = \\int_\\Omega \\Psi(\\mathbf{v}) \\ \\mathrm{d}\\Omega\\ .\n",
    "$$\n",
    "\n",
    "where $\\mathcal{K}(\\Omega)$ is a suitable function space.\n",
    "\n",
    "For clarity of presentation we ignore any non-zero surface tractions and body forces and instead consider only\n",
    "applied displacements (i.e. non-homogeneous dirichlet boundary conditions). Moreover we stick our attention to the standard Neo-Hookean stored energy density\n",
    "\n",
    "$$\n",
    "    \\Psi(\\mathbf{u}) = \\frac{\\mu}{2}\\left(I_1 - 3 \\right) - \\mu \\log(J) + \\frac{\\lambda}{2}\\left( J - 1\\right){}^2,\n",
    "$$\n",
    "where $I_1 = \\mathrm{tr}(\\mathbf{C}) = \\mathrm{tr}(\\mathbf{F}^\\mathrm{T} \\mathbf{F}) = F_{ij}F_{ij}$ and $J = \\det(\\mathbf{F})$ denote the standard invariants of the deformation gradient tensor $\\mathbf{F} = \\mathbf{I}+\\nabla_{\\mathbf{X}} \\mathbf{u}$.\n",
    "The above problem is ill-posed in the limit of incompressibility (or near-incompressibility), namely when\n",
    "$$\n",
    "    \\lambda/\\mu \\rightarrow +\\infty.\n",
    "$$\n",
    "In order to alleviate the problem, we consider the partial legendre transform of the strain energy density $\\Psi$ with respect to $J = \\det(\\mathbf{F})$, namely\n",
    "$$\n",
    "  \\widehat{\\Psi}(\\mathbf{u}, p) = \\sup_{J} \\left[ p(J - 1) - \\frac{\\mu}{2}\\left(I_1 - 3 \\right) + \\mu \\log(J) - \\frac{\\lambda}{2}\\left( J - 1\\right){}^2 \\right].\n",
    "$$\n",
    "The supremum, say $J^\\star$, can be calculated in closed form by the first order optimailty condition $\\partial\\widehat{\\Psi}/\\partial J = 0$. This gives\n",
    "$$\n",
    "  J^\\star(p) = \\frac{\\lambda + p + \\sqrt{(\\lambda + p){}^2 + 4 \\lambda \\mu }}{(2 \\lambda)}.\n",
    "$$\n",
    "Furthermore, taking the partial legendre transform of $\\widehat{\\Psi}$ once again, gives us back the original problem, i.e.\n",
    "$$\n",
    "    \\Psi(\\mathbf{u}) = \\Psi^\\star(\\mathbf{u}, p) = \\sup_{p} \\left[ p(J - 1) - p(J^\\star - 1) + \\frac{\\mu}{2}\\left(I_1 - 3 \\right) - \\mu \\log(J^\\star) + \\frac{\\lambda}{2}\\left( J^\\star - 1\\right){}^2 \\right].\n",
    "$$\n",
    "Therefore our original hyperelasticity problem can now be reformulated as\n",
    "$$\n",
    "  \\inf_{\\mathbf{u}\\in\\mathcal{K}(\\Omega)}\\sup_{p} \\int_\\Omega\\Psi^{\\star} (\\mathbf{u}, p) \\, \\mathrm{d}\\Omega.\n",
    "$$\n",
    "The total (modified) energy $\\Pi^\\star$ can then be written as\n",
    "$$\n",
    "  \\Pi^\\star(\\mathbf{u}, p) = \\int_\\Omega p (J - J^\\star) \\ \\mathrm{d}\\Omega + \\int_\\Omega \\frac{\\mu}{2} \\left( I_1 - 3\\right) \\ \\mathrm{d}\\Omega - \\int_\\Omega \\mu\\log(J^\\star)\\ \\mathrm{d}\\Omega + \\int_\\Omega \\frac{\\lambda}{2}\\left( J^\\star - 1 \\right){}^2\\ \\mathrm{d}\\Omega\n",
    "$$\n",
    "The Euler-Lagrange equations corresponding to the above energy give us our governing PDEs in the weak form, namely\n",
    "$$\n",
    "  \\int_\\Omega \\frac{\\partial\\Psi^\\star}{\\partial \\mathbf{F}}:\\delta \\mathbf{F} \\ \\mathrm{d}\\Omega = 0\n",
    "$$\n",
    "and\n",
    "$$\n",
    "  \\int_\\Omega \\frac{\\partial \\Psi^\\star}{\\partial p}\\delta p \\ \\mathrm{d}\\Omega = 0,\n",
    "$$\n",
    "where $\\delta \\mathrm{F} = \\delta \\mathrm{grad}_0(\\mathbf{u}) = \\mathrm{grad}_0(\\delta \\mathbf{u})$ and $\\delta \\mathbf{u}$ and $\\delta p$ denote arbitrary variations with respect to displacement and pressure (or the test functions). See the references\n",
    "below for a more detailed explanation of the above mathematical trick. Now, in order to apply Newton's method to the\n",
    "above problem, we further need to linearize the above equations and calculate the respective hessians (or tangents), namely, $\\partial^2\\Psi^\\star/\\partial \\mathbf{F}^2$, $\\partial^2\\Psi^\\star/\\partial p^2$ and $\\partial^2\\Psi^\\star/\\partial \\mathbf{F}\\partial p$\n",
    "which, using `Tensors.jl`, can be determined conveniently using automatic differentiation (see the code below). Hence we only need to define the above potential.\n",
    "The remaineder of the example follows similarly.\n",
    "## References\n",
    "1. [A paradigm for higher-order polygonal elements in finite elasticity using a gradient correction scheme, CMAME 2016, 306, 216–251](http://pamies.cee.illinois.edu/Publications_files/CMAME_2016.pdf)\n",
    "2. [Approximation of incompressible large deformation elastic problems: some unresolved issues, Computational Mechanics, 2013](https://link.springer.com/content/pdf/10.1007/s00466-013-0869-0.pdf)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Implementation\n",
    "We now get to the actual code. First, we import the respective packages"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Ferrite, Tensors, ProgressMeter\n",
    "using BlockArrays, SparseArrays, LinearAlgebra"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "and the corresponding `struct` to store our material properties."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "struct NeoHooke\n",
    "    μ::Float64\n",
    "    λ::Float64\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "We then create a function to generate a simple test mesh on which to compute FE solution. We also mark the boundaries\n",
    "to later assign Dirichlet boundary conditions"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function importTestGrid()\n",
    "    grid = generate_grid(Tetrahedron, (5, 5, 5), zero(Vec{3}), ones(Vec{3}));\n",
    "    addfaceset!(grid, \"myBottom\", x -> norm(x[2]) ≈ 0.0);\n",
    "    addfaceset!(grid, \"myBack\", x -> norm(x[3]) ≈ 0.0);\n",
    "    addfaceset!(grid, \"myRight\", x -> norm(x[1]) ≈ 1.0);\n",
    "    addfaceset!(grid, \"myLeft\", x -> norm(x[1]) ≈ 0.0);\n",
    "    return grid\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "The function to create corresponding cellvalues for the displacement field `u` and pressure `p`\n",
    "follows in a similar fashion from the `incompressible_elasticity` example"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function create_values(interpolation_u, interpolation_p)\n",
    "    # quadrature rules\n",
    "    qr      = QuadratureRule{3,RefTetrahedron}(4)\n",
    "    face_qr = QuadratureRule{2,RefTetrahedron}(4)\n",
    "\n",
    "    # geometric interpolation\n",
    "    interpolation_geom = Lagrange{3,RefTetrahedron,1}()\n",
    "\n",
    "    # cell and facevalues for u\n",
    "    cellvalues_u = CellVectorValues(qr, interpolation_u, interpolation_geom)\n",
    "    facevalues_u = FaceVectorValues(face_qr, interpolation_u, interpolation_geom)\n",
    "\n",
    "    # cellvalues for p\n",
    "    cellvalues_p = CellScalarValues(qr, interpolation_p, interpolation_geom)\n",
    "\n",
    "    return cellvalues_u, cellvalues_p, facevalues_u\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "We now create the function for Ψ*"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function Ψ(F, p, mp::NeoHooke)\n",
    "    μ = mp.μ\n",
    "    λ = mp.λ\n",
    "    Ic = tr(tdot(F))\n",
    "    J = det(F)\n",
    "    Js = (λ + p + sqrt((λ + p)^2. + 4. * λ * μ ))/(2. * λ)\n",
    "    return p * (Js - J) + μ / 2 * (Ic - 3) - μ * log(Js) + λ / 2 * (Js - 1)^2\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "and it's derivatives (required in computing the jacobian and hessian respectively)"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function constitutive_driver(F, p, mp::NeoHooke)\n",
    "    # Compute all derivatives in one function call\n",
    "    ∂²Ψ∂F², ∂Ψ∂F = Tensors.hessian(y -> Ψ(y, p, mp), F, :all)\n",
    "    ∂²Ψ∂p², ∂Ψ∂p = Tensors.hessian(y -> Ψ(F, y, mp), p, :all)\n",
    "    ∂²Ψ∂F∂p = Tensors.gradient(q -> Tensors.gradient(y -> Ψ(y, q, mp), F), p)\n",
    "    return ∂Ψ∂F, ∂²Ψ∂F², ∂Ψ∂p, ∂²Ψ∂p², ∂²Ψ∂F∂p\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "source": [
    "The functions to create the `DofHandler` and `ConstraintHandler` (to assign corresponding boundary conditions) follow\n",
    "likewise from the incompressible elasticity example, namely"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function create_dofhandler(grid, ipu, ipp)\n",
    "    dh = DofHandler(grid)\n",
    "    add!(dh, :u, 3, ipu) # displacement dim = 3\n",
    "    add!(dh, :p, 1, ipp) # pressure dim = 1\n",
    "    close!(dh)\n",
    "    return dh\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 7
  },
  {
   "cell_type": "markdown",
   "source": [
    "We are simulating a uniaxial tensile loading of a unit cube. Hence we apply a displacement field (`:u`) in `x` direction on the right face.\n",
    "The left, bottom and back faces are fixed in the `x`, `y` and `z` components of the displacement so as to emulate the uniaxial nature\n",
    "of the loading."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function create_bc(dh)\n",
    "    dbc = ConstraintHandler(dh)\n",
    "    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, \"myLeft\"), (x,t) -> zero(Vec{1}), [1]))\n",
    "    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, \"myBottom\"), (x,t) -> zero(Vec{1}), [2]))\n",
    "    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, \"myBack\"), (x,t) -> zero(Vec{1}), [3]))\n",
    "    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, \"myRight\"), (x,t) -> t*ones(Vec{1}), [1]))\n",
    "    close!(dbc)\n",
    "    Ferrite.update!(dbc, 0.0)\n",
    "    return dbc\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 8
  },
  {
   "cell_type": "markdown",
   "source": [
    "Also, since we are considering incompressible hyperelasticity, an interesting quantity that we can compute is the deformed volume of the solid.\n",
    "It is easy to show that this is equal to ∫J*dΩ where J=det(F). This can be done at the level of each element (cell)"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function calculate_element_volume(cell, cellvalues_u, ue)\n",
    "    reinit!(cellvalues_u, cell)\n",
    "    evol::Float64=0.0;\n",
    "    for qp in 1:getnquadpoints(cellvalues_u)\n",
    "        dΩ = getdetJdV(cellvalues_u, qp)\n",
    "        ∇u = function_gradient(cellvalues_u, qp, ue)\n",
    "        F = one(∇u) + ∇u\n",
    "        J = det(F)\n",
    "        evol += J * dΩ\n",
    "    end\n",
    "    return evol\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 9
  },
  {
   "cell_type": "markdown",
   "source": [
    "and then assembled over all the cells (elements)"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function calculate_volume_deformed_mesh(w, dh::DofHandler, cellvalues_u)\n",
    "    evol::Float64 = 0.0;\n",
    "    for cell in CellIterator(dh)\n",
    "        global_dofs = celldofs(cell)\n",
    "        nu = getnbasefunctions(cellvalues_u)\n",
    "        global_dofs_u = global_dofs[1:nu]\n",
    "        ue = w[global_dofs_u]\n",
    "        δevol = calculate_element_volume(cell, cellvalues_u, ue)\n",
    "        evol += δevol;\n",
    "    end\n",
    "    return evol\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 10
  },
  {
   "cell_type": "markdown",
   "source": [
    "The function to assemble the element stiffness matrix for each element in the mesh now has a block structure like in\n",
    "`incompressible_elasticity`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function assemble_element!(Ke, fe, cell, cellvalues_u, cellvalues_p, mp, ue, pe)\n",
    "    # Reinitialize cell values, and reset output arrays\n",
    "    ublock, pblock = 1, 2\n",
    "    reinit!(cellvalues_u, cell)\n",
    "    reinit!(cellvalues_p, cell)\n",
    "    fill!(Ke, 0.0)\n",
    "    fill!(fe, 0.0)\n",
    "\n",
    "    n_basefuncs_u = getnbasefunctions(cellvalues_u)\n",
    "    n_basefuncs_p = getnbasefunctions(cellvalues_p)\n",
    "\n",
    "    for qp in 1:getnquadpoints(cellvalues_u)\n",
    "        dΩ = getdetJdV(cellvalues_u, qp)\n",
    "        # Compute deformation gradient F\n",
    "        ∇u = function_gradient(cellvalues_u, qp, ue)\n",
    "        p = function_value(cellvalues_p, qp, pe)\n",
    "        F = one(∇u) + ∇u\n",
    "\n",
    "        # Compute first Piola-Kirchhoff stress and tangent modulus\n",
    "        ∂Ψ∂F, ∂²Ψ∂F², ∂Ψ∂p, ∂²Ψ∂p², ∂²Ψ∂F∂p = constitutive_driver(F, p, mp)\n",
    "\n",
    "        # Loop over the `u`-test functions to calculate the `u`-`u` and `u`-`p` blocks\n",
    "        for i in 1:n_basefuncs_u\n",
    "            # gradient of the test function\n",
    "            ∇δui = shape_gradient(cellvalues_u, qp, i)\n",
    "            # Add contribution to the residual from this test function\n",
    "            fe[BlockIndex((ublock), (i))] += ( ∇δui ⊡ ∂Ψ∂F) * dΩ\n",
    "\n",
    "            ∇δui∂S∂F = ∇δui ⊡ ∂²Ψ∂F²\n",
    "            for j in 1:n_basefuncs_u\n",
    "                ∇δuj = shape_gradient(cellvalues_u, qp, j)\n",
    "\n",
    "                # Add contribution to the tangent\n",
    "                Ke[BlockIndex((ublock, ublock), (i, j))] += ( ∇δui∂S∂F ⊡ ∇δuj ) * dΩ\n",
    "            end\n",
    "            # Loop over the `p`-test functions\n",
    "            for j in 1:n_basefuncs_p\n",
    "                δp = shape_value(cellvalues_p, qp, j)\n",
    "                # Add contribution to the tangent\n",
    "                Ke[BlockIndex((ublock, pblock), (i, j))] += ( ∂²Ψ∂F∂p ⊡ ∇δui ) * δp * dΩ\n",
    "            end\n",
    "        end\n",
    "        # Loop over the `p`-test functions to calculate the `p-`u` and `p`-`p` blocks\n",
    "        for i in 1:n_basefuncs_p\n",
    "            δp = shape_value(cellvalues_p, qp, i)\n",
    "            fe[BlockIndex((pblock), (i))] += ( δp * ∂Ψ∂p) * dΩ\n",
    "\n",
    "            for j in 1:n_basefuncs_u\n",
    "                ∇δuj = shape_gradient(cellvalues_u, qp, j)\n",
    "                Ke[BlockIndex((pblock, ublock), (i, j))] += ∇δuj ⊡ ∂²Ψ∂F∂p * δp * dΩ\n",
    "            end\n",
    "            for j in 1:n_basefuncs_p\n",
    "                δp = shape_value(cellvalues_p, qp, j)\n",
    "                Ke[BlockIndex((pblock, pblock), (i, j))] += δp * ∂²Ψ∂p² * δp * dΩ\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 11
  },
  {
   "cell_type": "markdown",
   "source": [
    "The only thing that changes in the assembly of the global stiffness matrix is slicing the corresponding element\n",
    "dofs for the displacement (see `global_dofsu`) and pressure (`global_dofsp`)."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function assemble_global!(K::SparseMatrixCSC, f, cellvalues_u::CellVectorValues{dim},\n",
    "                         cellvalues_p::CellScalarValues{dim}, dh::DofHandler, mp::NeoHooke, w) where {dim}\n",
    "    nu = getnbasefunctions(cellvalues_u)\n",
    "    np = getnbasefunctions(cellvalues_p)\n",
    "\n",
    "    # start_assemble resets K and f\n",
    "    fe = PseudoBlockArray(zeros(nu + np), [nu, np]) # local force vector\n",
    "    ke = PseudoBlockArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # local stiffness matrix\n",
    "\n",
    "    assembler = start_assemble(K, f)\n",
    "    # Loop over all cells in the grid\n",
    "    for cell in CellIterator(dh)\n",
    "        global_dofs = celldofs(cell)\n",
    "        global_dofsu = global_dofs[1:nu]; # first nu dofs are displacement\n",
    "        global_dofsp = global_dofs[nu + 1:end]; # last np dofs are pressure\n",
    "        @assert size(global_dofs, 1) == nu + np # sanity check\n",
    "        ue = w[global_dofsu] # displacement dofs for the current cell\n",
    "        pe = w[global_dofsp] # pressure dofs for the current cell\n",
    "        assemble_element!(ke, fe, cell, cellvalues_u, cellvalues_p, mp, ue, pe)\n",
    "        assemble!(assembler, global_dofs, fe, ke)\n",
    "    end\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 12
  },
  {
   "cell_type": "markdown",
   "source": [
    "We now define a main function `solve`. For nonlinear quasistatic problems we often like to parameterize the\n",
    "solution in terms of a pseudo time like parameter, which in this case is used to gradually apply the boundary\n",
    "displacement on the right face. Also for definitenessm we consider λ/μ = 10⁴"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function solve(interpolation_u, interpolation_p)\n",
    "\n",
    "    # import the mesh\n",
    "    grid = importTestGrid()\n",
    "\n",
    "    # Material parameters\n",
    "    μ = 1.\n",
    "    λ = 1.E4 * μ\n",
    "    mp = NeoHooke(μ, λ)\n",
    "\n",
    "    # Create the DofHandler and CellValues\n",
    "    dh = create_dofhandler(grid, interpolation_u, interpolation_p)\n",
    "    cellvalues_u, cellvalues_p, facevalues_u = create_values(interpolation_u, interpolation_p)\n",
    "\n",
    "    # Create the DirichletBCs\n",
    "    dbc = create_bc(dh)\n",
    "\n",
    "    # Pre-allocation of vectors for the solution and Newton increments\n",
    "    _ndofs = ndofs(dh)\n",
    "    w  = zeros(_ndofs)\n",
    "    ΔΔw = zeros(_ndofs)\n",
    "    apply!(w, dbc)\n",
    "\n",
    "    # Create the sparse matrix and residual vector\n",
    "    K = create_sparsity_pattern(dh)\n",
    "    f = zeros(_ndofs)\n",
    "\n",
    "    # We run the simulation parameterized by a time like parameter. `Tf` denotes the final value\n",
    "    # of this parameter, and Δt denotes its increment in each step\n",
    "    Tf = 2.0;\n",
    "    Δt = 0.1;\n",
    "    NEWTON_TOL = 1e-8\n",
    "\n",
    "    pvd = paraview_collection(\"hyperelasticity_incomp_mixed.pvd\");\n",
    "    for t ∈ 0.0:Δt:Tf\n",
    "        # Perform Newton iterations\n",
    "        Ferrite.update!(dbc, t)\n",
    "        apply!(w, dbc)\n",
    "        newton_itr = -1\n",
    "        prog = ProgressMeter.ProgressThresh(NEWTON_TOL, \"Solving @ time $t of $Tf;\")\n",
    "        fill!(ΔΔw, 0.0);\n",
    "        while true; newton_itr += 1\n",
    "            assemble_global!(K, f, cellvalues_u, cellvalues_p, dh, mp, w)\n",
    "            norm_res = norm(f[Ferrite.free_dofs(dbc)])\n",
    "            apply_zero!(K, f, dbc)\n",
    "            # Only display output at specific load steps\n",
    "            if t%(5*Δt) == 0\n",
    "                ProgressMeter.update!(prog, norm_res; showvalues = [(:iter, newton_itr)])\n",
    "            end\n",
    "            if norm_res < NEWTON_TOL\n",
    "                break\n",
    "            elseif newton_itr > 30\n",
    "                error(\"Reached maximum Newton iterations, aborting\")\n",
    "            end\n",
    "            # Compute the incremental `dof`-vector (both displacement and pressure)\n",
    "            ΔΔw .= K\\f;\n",
    "\n",
    "            apply_zero!(ΔΔw, dbc)\n",
    "            w .-= ΔΔw\n",
    "        end;\n",
    "\n",
    "        # Save the solution fields\n",
    "        vtk_grid(\"hyperelasticity_incomp_mixed_$t.vtu\", dh) do vtkfile\n",
    "            vtk_point_data(vtkfile, dh, w)\n",
    "            vtk_save(vtkfile)\n",
    "            pvd[t] = vtkfile\n",
    "        end\n",
    "    end;\n",
    "    vtk_save(pvd);\n",
    "    vol_def = calculate_volume_deformed_mesh(w, dh, cellvalues_u);\n",
    "    print(\"Deformed volume is $vol_def\")\n",
    "    return vol_def;\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 13
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can now test the solution using the Taylor-Hood approximation"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\rSolving @ time 0.0 of 2.0; (thresh = 1e-08, value = 4.78011e-05)\u001b[K\r\n",
      "  iter:  1\u001b[K\r\u001b[A\n",
      "\r\u001b[K\u001b[A\rSolving @ time 0.0 of 2.0; Time: 0:00:00 (3 iterations)\u001b[K\r\n",
      "  iter:  2\u001b[K\n",
      "\rSolving @ time 0.5 of 2.0; (thresh = 1e-08, value = 0.00164904)\u001b[K\r\n",
      "  iter:  1\u001b[K\r\u001b[A\n",
      "\r\u001b[K\u001b[A\rSolving @ time 0.5 of 2.0; (thresh = 1e-08, value = 5.74317e-05)\u001b[K\r\n",
      "  iter:  2\u001b[K\r\u001b[A\n",
      "\r\u001b[K\u001b[A\rSolving @ time 0.5 of 2.0; (thresh = 1e-08, value = 1.07293e-08)\u001b[K\r\n",
      "  iter:  3\u001b[K\r\u001b[A\n",
      "\r\u001b[K\u001b[A\rSolving @ time 0.5 of 2.0; Time: 0:00:01 (5 iterations)\u001b[K\r\n",
      "  iter:  4\u001b[K\n",
      "\rSolving @ time 1.0 of 2.0; (thresh = 1e-08, value = 0.000819569)\u001b[K\r\n",
      "  iter:  1\u001b[K\r\u001b[A\n",
      "\r\u001b[K\u001b[A\rSolving @ time 1.0 of 2.0; (thresh = 1e-08, value = 2.60479e-05)\u001b[K\r\n",
      "  iter:  2\u001b[K\r\u001b[A\n",
      "\r\u001b[K\u001b[A\rSolving @ time 1.0 of 2.0; Time: 0:00:00 (4 iterations)\u001b[K\r\n",
      "  iter:  3\u001b[K\n",
      "\rSolving @ time 1.5 of 2.0; (thresh = 1e-08, value = 0.000500654)\u001b[K\r\n",
      "  iter:  1\u001b[K\r\u001b[A\n",
      "\r\u001b[K\u001b[A\rSolving @ time 1.5 of 2.0; (thresh = 1e-08, value = 1.30232e-05)\u001b[K\r\n",
      "  iter:  2\u001b[K\r\u001b[A\n",
      "\r\u001b[K\u001b[A\rSolving @ time 1.5 of 2.0; Time: 0:00:00 (4 iterations)\u001b[K\r\n",
      "  iter:  3\u001b[K\n",
      "\rSolving @ time 2.0 of 2.0; (thresh = 1e-08, value = 0.000339016)\u001b[K\r\n",
      "  iter:  1\u001b[K\r\u001b[A\n",
      "\r\u001b[K\u001b[A\rSolving @ time 2.0 of 2.0; (thresh = 1e-08, value = 6.94365e-06)\u001b[K\r\n",
      "  iter:  2\u001b[K\r\u001b[A\n",
      "\r\u001b[K\u001b[A\rSolving @ time 2.0 of 2.0; Time: 0:00:00 (4 iterations)\u001b[K\r\n",
      "  iter:  3\u001b[K\n",
      "Deformed volume is 1.0001999771699892"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "1.0001999771699892"
     },
     "metadata": {},
     "execution_count": 14
    }
   ],
   "cell_type": "code",
   "source": [
    "quadratic = Lagrange{3, RefTetrahedron, 2}()\n",
    "linear = Lagrange{3, RefTetrahedron, 1}()\n",
    "vol_def = solve(quadratic, linear)"
   ],
   "metadata": {},
   "execution_count": 14
  },
  {
   "cell_type": "markdown",
   "source": [
    "The deformed volume is indeed close to 1 (as should be for a nearly incompressible material)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "---\n",
    "\n",
    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.5"
  },
  "kernelspec": {
   "name": "julia-1.8",
   "display_name": "Julia 1.8.5",
   "language": "julia"
  }
 },
 "nbformat": 4
}