JuliaFEM.jl/docs/tutorials/2015-06-25-elasticity-solver-example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving elasticity problems using JuliaFEM\n",
    "\n",
    "Author(s): Jukka Aho\n",
    "\n",
    "**Abstract**: Solving elasticity equations using JuliaFEM.\n",
    "\n",
    "### Weak form\n",
    "\n",
    "Given function spaces\n",
    "\\begin{align}\n",
    "\\boldsymbol{\\mathcal{U}} & =\\left\\{ \\boldsymbol{u}\\in H^{1}\\left(\\Omega\\right)|\\boldsymbol{u}\\left(\\boldsymbol{X},t\\right)=\\hat{\\boldsymbol{u}}\\left(\\boldsymbol{X},t\\right)\\text{ on }\\Gamma_{\\mathrm{u}}\\right\\} ,\\\\\n",
    "\\boldsymbol{\\mathcal{V}} & =\\left\\{ \\delta\\boldsymbol{u}\\in H^{1}\\left(\\Omega\\right)|\\delta\\boldsymbol{u}\\left(\\boldsymbol{X}\\right)=0\\text{ on }\\Gamma_{\\mathrm{u}}\\right\\} ,\n",
    "\\end{align}\n",
    "find $\\boldsymbol{u}\\in\\boldsymbol{\\mathcal{U}}$ such that\n",
    "\\begin{equation}\n",
    "\\delta\\mathcal{W}:=\\int_{\\Omega_{0}}\\rho_{0}\\ddot{\\boldsymbol{u}}\\cdot\\delta\\boldsymbol{u}\\,\\mathrm{d}V_{0}+\\int_{\\Omega_{0}}\\boldsymbol{S}:\\delta\\boldsymbol{E}\\,\\mathrm{d}V_{0}-\\int_{\\Omega_{0}}\\hat{\\boldsymbol{b}}_{0}\\cdot\\delta\\boldsymbol{u}\\,\\mathrm{d}V_{0}-\\int_{\\Gamma_{\\sigma}}\\hat{\\boldsymbol{t}}_{0}\\cdot\\delta\\boldsymbol{u}\\,\\mathrm{d}A_{0} =0 \\qquad\\forall\\delta\\boldsymbol{u}\\in\\boldsymbol{\\mathcal{V}}\n",
    "\\end{equation}\n",
    "\n",
    "### Some formulas\n",
    "\\begin{align}\n",
    "J & =\\det\\left(F\\right)\\\\\n",
    "I_{c} & =\\mbox{tr}\\left(C\\right)\\\\\n",
    "\\mathbf{C} & =\\mathbf{F}^{\\mathrm{T}}\\mathbf{F}\\\\\n",
    "\\mathbf{F} & =\\mathbf{I}+\\nabla\\mathbf{u}\\\\\n",
    "\\mathbf{E} & =\\frac{1}{2}\\left(\\mathbf{F}^{\\mathrm{T}}\\mathbf{F}-\\mathbf{I}\\right)\n",
    "\\end{align}\n",
    "\n",
    "### Potential energy\n",
    "\n",
    "\\begin{equation}\n",
    "\\underset{u\\in\\boldsymbol{\\mathcal{U}}}{\\min}\\Pi\\left(\\mathbf{u}\\right)\n",
    "\\end{equation}\n",
    "\\begin{equation}\n",
    "\\Pi\\left(\\mathbf{u}\\right)=\\int_{\\Omega}\\psi\\left(\\mathbf{u}\\right)-\\int_{\\Omega}\\hat{\\mathbf{b}}_{0}\\cdot\\mathbf{u}-\\int_{\\Gamma_{\\sigma}}\\hat{\\mathbf{t}}_{0}\\cdot\\mathbf{u}\\,\\mathrm{d}A_{0}\n",
    "\\end{equation}\n",
    "\n",
    "### Material models\n",
    "\n",
    "https://en.wikipedia.org/wiki/Strain_energy_density_function\n",
    "\n",
    "Saint-Venant-Kirchhoff model https://en.wikipedia.org/wiki/Hyperelastic_material\n",
    "\\begin{equation}\n",
    "\\psi\\left(\\mathbf{E}\\right)=\\frac{\\lambda}{2}\\left[\\mbox{tr}\\left(\\mathbf{E}\\right)\\right]^{2}+\\mu\\mbox{tr}\\left(\\mathbf{E}^2\\right)\n",
    "\\end{equation}\n",
    "\n",
    "neo-Hookean material https://en.wikipedia.org/wiki/Neo-Hookean_solid\n",
    "\\begin{equation}\n",
    "\\psi=\\frac{\\mu}{2}\\left(I_{c}-3\\right)-\\mu\\ln\\left(J\\right)+\\frac{\\lambda}{2}\\ln\\left(J\\right)^{2}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using JuliaFEM.API: Model, Material, add_material!, LoadCase, add_element_set!, ForceBC, DisplacementBC\n",
    "using JuliaFEM.API: add_boundary_condition!, add_load_case!, add_solver!\n",
    "using JuliaFEM.Interfaces: solve!\n",
    "using JuliaFEM.Preprocess: parse_abaqus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Parsing nodes\n",
      "INFO: Parsing elements. Type: Tet10\n",
      "INFO: Parsing elements. Type: Tri6\n",
      "INFO: Creating elset BC1\n",
      "INFO: Creating elset BC2\n",
      "INFO: model loaded.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  4.209444 seconds (19.12 M allocations: 565.872 MB, 11.55% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time begin\n",
    "    # Linear models\n",
    "    #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_8789_P1.inp\")\n",
    "    #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_16436_P1.inp\")\n",
    "    #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_27343_P1.inp\")\n",
    "    #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_45510_P1.inp\")\n",
    "    #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_75470_P1.inp\")\n",
    "    #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/wrench/wrench_128903_P1.inp\")\n",
    "\n",
    "    # Quadratic models\n",
    "    #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_19611_P2.inp\")\n",
    "    #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_55950_P2.inp\")\n",
    "    #model = open(JuliaFEM.Core.parse_abaqus, \"/Temp/piston_107168_P2.inp\")\n",
    "    #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_345757_P2.inp\")\n",
    "\n",
    "    mesh = open(parse_abaqus, \"/Temp/piston_107168_P2.inp\")\n",
    "    model = Model(\"piston model\", mesh)\n",
    "    info(\"model loaded.\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "JuliaFEM.API.Material(\"steel\",Dict(\"poissons ratio\"=>0.3,\"youngs modulus\"=>210000.0))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mat = Material(\"steel\")\n",
    "mat[\"youngs modulus\"] = 210.0e3\n",
    "mat[\"poissons ratio\"] = 0.3\n",
    "add_material!(model, \"PISTON\", mat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1-element Array{ASCIIString,1}:\n",
       " \"PISTON\""
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simulation = LoadCase(:ElasticityProblem)\n",
    "\n",
    "# simulation = Simulation(:ElasticityProblem, \"elasticity equations\")\n",
    "# add_element_set!(simulation, \"PISTON\")\n",
    "# add_boundary_condition!(simulation, bc1, bc2)\n",
    "# add_simulation!(model, simulation)\n",
    "# add_solver!(simulation, :DirectSolver)\n",
    "# solve!(model, \"elasticity equations\", 0.0)\n",
    "# solve!(model, \"heat equations\", 0.0)\n",
    "# solve!(model, \"solve equations 3\", 0.0)\n",
    "\n",
    "add_element_set!(simulation, \"PISTON\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1-element Array{JuliaFEM.API.DirichletBC,1}:\n",
       " JuliaFEM.API.DirichletBC(\"BC2\",\"displacement\"=>0.0)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "traction = Vector{Float64}[[0.0, 0.0, -10.0] for i in 1:6]\n",
    "bc1 = ForceBC(\"BC1\", \"displacement traction force\" => traction)\n",
    "bc2 = DisplacementBC(\"BC2\", \"displacement\" => 0.0)\n",
    "\n",
    "# JuliaFEM.API.add_boundary_condition!(prob, bc1, bc2)\n",
    "add_boundary_condition!(simulation, bc1)\n",
    "add_boundary_condition!(simulation, bc2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "JuliaFEM.API.LoadCase(:ElasticityProblem,[JuliaFEM.API.NeumannBC(\"BC1\",\"displacement traction force\"=>[[0.0,0.0,-10.0],[0.0,0.0,-10.0],[0.0,0.0,-10.0],[0.0,0.0,-10.0],[0.0,0.0,-10.0],[0.0,0.0,-10.0]])],[JuliaFEM.API.DirichletBC(\"BC2\",\"displacement\"=>0.0)],:DirectSolver,ASCIIString[\"PISTON\"])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "add_solver!(simulation, :DirectSolver)\n",
    "add_load_case!(model, \"elasticity equations\", simulation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "solver = JuliaFEM.Interfaces.get_solver(model, \"elasticity equations\", 0.0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "JuliaFEM.Core.Element{JuliaFEM.Core.Tet10}([1,2,3,4,5,6,7,8,9,10],Dict{ASCIIString,JuliaFEM.Core.Field{A<:Union{JuliaFEM.Core.Continuous,JuliaFEM.Core.Discrete},B<:Union{JuliaFEM.Core.Constant,JuliaFEM.Core.Variable},C<:Union{JuliaFEM.Core.TimeInvariant,JuliaFEM.Core.TimeVariant}}}(\"poissons ratio\"=>JuliaFEM.Core.Field{JuliaFEM.Core.Discrete,JuliaFEM.Core.Constant,JuliaFEM.Core.TimeInvariant}(0.3),\"geometry\"=>JuliaFEM.Core.Field{JuliaFEM.Core.Discrete,JuliaFEM.Core.Variable,JuliaFEM.Core.TimeInvariant}([[13.96894,1.54855,-2.99382],[15.71724,2.88794,-2.79243],[14.76115,1.1419,-1.23125],[15.93754,1.0623,-2.82067],[14.84309,2.21825,-2.89313],[15.2392,2.01492,-2.01184],[14.36505,1.34523,-2.11254],[14.95324,1.30543,-2.90725],[15.82739,1.97512,-2.80655],[15.34935,1.1021,-2.02596]]),\"youngs modulus\"=>JuliaFEM.Core.Field{JuliaFEM.Core.Discrete,JuliaFEM.Core.Constant,JuliaFEM.Core.TimeInvariant}(210000.0)))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "celements = solver.field_problems[1].elements;\n",
    "celements[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: 0 surface elements\n"
     ]
    }
   ],
   "source": [
    "surface_elements = filter((e) -> isa(e, JuliaFEM.Core.Element{JuliaFEM.Core.Tri6}), celements)\n",
    "info(\"$(length(surface_elements)) surface elements\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "62454"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "length(celements)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: # of field problems: 1\n",
      "INFO: # of boundary problems: 1\n",
      "INFO: Starting iteration 1\n",
      "INFO: Assembling field problems...\n",
      "INFO: Assembling body 1...\n",
      "INFO: Assembly: 10.0 % done. \n",
      "INFO: Assembly: 20.0 % done. \n",
      "INFO: Assembly: 30.0 % done. \n",
      "INFO: Assembly: 40.0 % done. \n",
      "INFO: Assembly: 50.0 % done. \n",
      "INFO: Assembly: 60.0 % done. \n",
      "INFO: Assembly: 70.0 % done. \n",
      "INFO: Assembly: 80.0 % done. \n",
      "INFO: Assembly: 90.0 % done. \n",
      "INFO: Assembly: 100.0 % done. \n",
      "INFO: Assembling boundary problems...\n",
      "INFO: Assembling boundary 1...\n",
      "INFO: Assembly: 10.0 % done. \n",
      "INFO: Assembly: 20.0 % done. \n",
      "INFO: Assembly: 30.0 % done. \n",
      "INFO: Assembly: 40.0 % done. \n",
      "INFO: Assembly: 50.0 % done. \n",
      "INFO: Assembly: 60.0 % done. \n",
      "INFO: Assembly: 70.0 % done. \n",
      "INFO: Assembly: 80.0 % done. \n",
      "INFO: Assembly: 90.0 % done. \n",
      "INFO: Assembly: 100.0 % done. \n",
      "INFO: Solving system\n",
      "INFO: all dofs = 321504\n",
      "INFO: interior dofs = 312936\n",
      "INFO: boundary dofs = 8568\n",
      "INFO: preparation in 3.385999917984009 seconds\n",
      "INFO: displacement on boundary solved.\n",
      "INFO: norm[u_boundary_dofs] = 0.0\n",
      "INFO: homogeneous dirichlet boundary\n",
      "INFO: solve boundary = 0.42100000381469727\n",
      "INFO: factorizations done in 14.133000135421753 seconds\n",
      "INFO: solved interior in 0.2969999313354492 seconds. norm = 0.0\n",
      "INFO: timing info for non-linear iteration:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(1,true)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: boundary assembly       : 2.1679999828338623\n",
      "INFO: field assembly          : 279.9069998264313\n",
      "INFO: dump matrices to disk   : 0.0\n",
      "INFO: solve problem           : 19.37600016593933\n",
      "INFO: update element data     : 1.6999998092651367\n",
      "INFO: non-linear iteration    : 303.1509997844696\n",
      "INFO: solver finished in 304.5089998245239 seconds.\n"
     ]
    }
   ],
   "source": [
    "solve!(model, \"elasticity equations\", 0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saving results to file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Grid Name=\"Grid\">\n",
       "  <Time Value=\"0\"/>\n",
       "</Grid>\n"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xdoc, xmodel = JuliaFEM.Postprocess.xdmf_new_model()\n",
    "temporal_collection = JuliaFEM.Postprocess.xdmf_new_temporal_collection(xmodel)\n",
    "grid = JuliaFEM.Postprocess.xdmf_new_grid(temporal_collection; time=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save geometry to xdmf file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: number of nodes in model: 8789\n",
      "INFO: Number of elements in model: 31437\n"
     ]
    }
   ],
   "source": [
    "nnodes = length(model[\"nodes\"])\n",
    "info(\"number of nodes in model: $nnodes\")\n",
    "\n",
    "X = zeros(3, nnodes)\n",
    "for nid in keys(model[\"nodes\"])\n",
    "    X[:, perm[nid]] = model[\"nodes\"][nid]\n",
    "end\n",
    "\n",
    "nelements = length(model[\"elsets\"][\"PISTON\"])\n",
    "info(\"Number of elements in model: $nelements\")\n",
    "elmap = zeros(Int64, 5, nelements)\n",
    "#elmap[1,:] = 0x0026\n",
    "elmap[1,:] = 0x6\n",
    "for (i, elid) in enumerate(model[\"elsets\"][\"PISTON\"])\n",
    "    elmap[2:end,i] = Int64[perm[nid] for nid in model[\"elements\"][elid]]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "JuliaFEM.Postprocess.xdmf_new_mesh(grid, X, elmap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save nodal data to model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3,8789)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using JuliaFEM.Core: get_connectivity\n",
    "\n",
    "u = zeros(3, nnodes)\n",
    "\n",
    "for element in values(elements)\n",
    "    isa(element, Element{Tet4}) || continue\n",
    "    connectivity = get_connectivity(element)\n",
    "    field = element[\"displacement\"](0.0)\n",
    "    for (i, nid) in enumerate(connectivity)\n",
    "        u[:, nid] = field[i]\n",
    "    end\n",
    "end\n",
    "size(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "JuliaFEM.Postprocess.xdmf_new_field(grid, \"Displacement\", \"nodes\", u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1525505"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "JuliaFEM.Postprocess.xdmf_save_model(xdoc, \"/tmp/piston_8789_P1.xmf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3x10 Array{Float64,2}:\n",
       " 0.0582911   0.339694  0.431768    …   0.124351  0.346492    0.000617283\n",
       " 0.0900323   0.22174   0.433814        0.140503  0.275036    0.0044649  \n",
       " 0.0838928  -0.148224  0.00262483     -0.11171   0.0843743  -0.0175412  "
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u[:, 1:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1x3 Array{Float64,2}:\n",
       " -0.0471355  -0.468786  -0.222571"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "minimum(u, 2)'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1x3 Array{Float64,2}:\n",
       " 0.511168  0.563004  0.66328"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "maximum(u, 2)'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "38.52232721814439"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "piston_8789_P1_solution_norm = 38.52232721814439"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@assert isapprox(norm(vec(u)), piston_8789_P1_solution_norm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.1",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}