{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Tutorial 5: Hyper-elasticity"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "    This tutorial is under construction, but the code below is already functional.\n",
    "\n",
    "\n",
    "## Problem statement"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Gridap\n",
    "using LinearAlgebra: inv, det\n",
    "using LineSearches: BackTracking"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Material parameters"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "const λ = 100.0\n",
    "const μ = 1.0"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Deformation Gradient"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "F(∇u) = one(∇u) + ∇u'\n",
    "\n",
    "J(F) = sqrt(det(C(F)))\n",
    "\n",
    "#Green strain\n",
    "\n",
    "#E(F) = 0.5*( F'*F - one(F) )\n",
    "\n",
    "@law dE(∇du,∇u) = 0.5*( ∇du⋅F(∇u) + (∇du⋅F(∇u))' )"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Right Cauchy-green deformation tensor"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "C(F) = (F')⋅F"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Constitutive law (Neo hookean)"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@law function S(∇u)\n",
    "  Cinv = inv(C(F(∇u)))\n",
    "  μ*(one(∇u)-Cinv) + λ*log(J(F(∇u)))*Cinv\n",
    "end\n",
    "\n",
    "@law function dS(∇du,∇u)\n",
    "  Cinv = inv(C(F(∇u)))\n",
    "  _dE = dE(∇du,∇u)\n",
    "\tλ*(Cinv⊙_dE)*Cinv + 2*(μ-λ*log(J(F(∇u))))*Cinv⋅_dE⋅(Cinv')\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Cauchy stress tensor"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@law σ(∇u) = (1.0/J(F(∇u)))*F(∇u)⋅S(∇u)⋅(F(∇u))'"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Weak form"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "res(u,v) = dE(∇(v),∇(u)) ⊙ S(∇(u))\n",
    "\n",
    "jac_mat(u,du,v) =  dE(∇(v),∇(u)) ⊙ dS(∇(du),∇(u))\n",
    "\n",
    "jac_geo(u,du,v) = ∇(v) ⊙ ( S(∇(u))⋅∇(du) )\n",
    "\n",
    "jac(u,du,v) = jac_mat(u,v,du) + jac_geo(u,v,du)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Model"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "domain = (0,1,0,1)\n",
    "partition = (20,20)\n",
    "model = CartesianDiscreteModel(domain,partition)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Define new boundaries"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "labels = get_face_labeling(model)\n",
    "add_tag_from_tags!(labels,\"diri_0\",[1,3,7])\n",
    "add_tag_from_tags!(labels,\"diri_1\",[2,4,8])"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Construct the FEspace"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "V = TestFESpace(\n",
    "  model=model,valuetype=VectorValue{2,Float64},\n",
    "  reffe=:Lagrangian,conformity=:H1,order=1,\n",
    "  dirichlet_tags = [\"diri_0\", \"diri_1\"])"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Setup integration"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "trian = Triangulation(model)\n",
    "degree = 2\n",
    "quad = CellQuadrature(trian,degree)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Setup weak form terms"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "t_Ω = FETerm(res,jac,trian,quad)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Setup non-linear solver"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "nls = NLSolver(\n",
    "  show_trace=true,\n",
    "  method=:newton,\n",
    "  linesearch=BackTracking())\n",
    "\n",
    "solver = FESolver(nls)\n",
    "\n",
    "function run(x0,disp_x,step,nsteps)\n",
    "\n",
    "  g0 = VectorValue(0.0,0.0)\n",
    "  g1 = VectorValue(disp_x,0.0)\n",
    "  U = TrialFESpace(V,[g0,g1])\n",
    "\n",
    "  #FE problem\n",
    "  op = FEOperator(U,V,t_Ω)\n",
    "\n",
    "  println(\"\\n+++ Solving for disp_x $disp_x in step $step of $nsteps +++\\n\")\n",
    "\n",
    "  uh = FEFunction(U,x0)\n",
    "\n",
    "  uh, = solve!(uh,solver,op)\n",
    "\n",
    "  writevtk(trian,\"results_$(lpad(step,3,'0'))\",cellfields=[\"uh\"=>uh,\"sigma\"=>σ(∇(uh))])\n",
    "\n",
    "  return get_free_values(uh)\n",
    "\n",
    "end\n",
    "\n",
    "function runs()\n",
    "\n",
    " disp_max = 0.75\n",
    " disp_inc = 0.02\n",
    " nsteps = ceil(Int,abs(disp_max)/disp_inc)\n",
    "\n",
    " x0 = zeros(Float64,num_free_dofs(V))\n",
    "\n",
    " for step in 1:nsteps\n",
    "   disp_x = step * disp_max / nsteps\n",
    "   x0 = run(x0,disp_x,step,nsteps)\n",
    " end\n",
    "\n",
    "end\n",
    "\n",
    "#Do the work!\n",
    "runs()"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Picture of the last load step\n",
    "![](../assets/hyperelasticity/neo_hook_2d.png)\n",
    "\n",
    "##  Extension to 3D\n",
    "\n",
    "Extending this tutorial to the 3D case is straightforward. It is leaved as an exercise.\n",
    "\n",
    "![](../assets/hyperelasticity/neo_hook_3d.png)"
   ],
   "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.3.1"
  },
  "kernelspec": {
   "name": "julia-1.3",
   "display_name": "Julia 1.3.1",
   "language": "julia"
  }
 },
 "nbformat": 4
}