{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FEniCS for acoustics\n",
    "## Author: Ben Goldsberry\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import FEniCS module\n",
    "To import FEniCS, use:\n",
    "```python\n",
    "from dolfin import *\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import mesh and define function space\n",
    "```python\n",
    "mesh = Mesh('mesh.xml')\n",
    "```\n",
    "For a scalar function (pressure acoustics)\n",
    "```python\n",
    "V = FunctionSpace(mesh,'CG',2)\n",
    "p = TrialFunction(V) #pressure trial function\n",
    "dp = TestFunction(V) #pressure test function\n",
    "```\n",
    "for elasticity:\n",
    "```python\n",
    "V = VectorFunctionSpace(mesh,'CG',2)\n",
    "u = TrialFunction(V) #Displacement trial function\n",
    "du = TestFunction(V) #Displacement test function\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Weak forms\n",
    "## Pressure acoustics\n",
    "The inhomogeneous Helmholtz equation for pressure is:\n",
    "$$ \\nabla^2p+k^2p =f$$\n",
    "The weak form of the above equation is:\n",
    "$$ \\frac{k^2}{\\rho}\\int \\limits_\\Omega  p \\delta p \\, d \\Omega - \\frac{1}{\\rho}\\int \\limits_\\Omega \\nabla p \\cdot \\nabla \\delta p \\,d \\Omega + \\frac{1}{\\rho}\\int \\limits_{\\Gamma_N} \\delta p (\\nabla p \\cdot n) \\,d\\Gamma_N =\\frac{1}{\\rho}\\int \\limits_\\Omega f\\delta p \\, d \\Omega     \\\\\n",
    "p = 0 \\quad \\text{on} \\quad \\Gamma_D$$\n",
    "This is written is FEniCs as\n",
    "```python\n",
    "n = FacetNormal(mesh)\n",
    "StiffnessMatrix = 1/rho*dot(grad(p),grad(dp))*dx\n",
    "MassMatrix = 1/rho*p*dp*dx\n",
    "NeumannBC = 1/rho*dp*dot(p,n)*ds\n",
    "forcingFcn = 1/rho*dp*f*dx\n",
    "WeakForm = -Stiffness + k^2*MassMatrix + NeumannBC\n",
    "```\n",
    "## Linear Elasticity\n",
    "The stiffness tensor is\n",
    "$$ C_{ijkl} = K\\delta_{ij}\\delta_{kl} +\\mu(\\delta_{ik}\\delta_{jl} +\\delta_{il}\\delta_{jk} -\\frac{2}{3}\\delta_{ij}\\delta_{kl})$$\n",
    "In FEniCS, this is implemented as:\n",
    "```python\n",
    "I = identity(2)# Second order identity tensor\n",
    "C = as_tensor(K*I[i,j]*I[k,l] + mu*(I[i,k]*I[j,l] + I[i,l]*I[j,k] - 2./3*I[i,j]*I[k,l],(i,j,k,l))\n",
    "```\n",
    "The weak form of linear elasticity is:\n",
    "$$ \\int \\limits_\\Omega [\\sigma(u)_{ij}\\varepsilon(\\delta u)_{ij} - \\omega^2\\rho \\delta u_i u_i] \\,d\\Omega -\\int \\limits_\\Gamma \\sigma_{ij}n_j\\delta u_i \\,d \\Gamma=0$$\n",
    "This is implemented in FEniCS as:\n",
    "```python\n",
    "StrainU = 0.5*(grad(u)+transpose(grad(u))) = sym(grad(u))\n",
    "StrainDu = 0.5*(grad(du)+transpose(grad(du)))\n",
    "StressU = as_tensor(C[i,j,k,l]*StrainU[k,l],(i,j))\n",
    "WeakForm = inner(StressU,StrainDu)*dx - omega**2*rho*dot(u,du)*dx - p*dot(du,n)*ds\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonlinear Elasticity\n",
    "```python\n",
    "F = grad(u)+I\n",
    "C = F.T*F #Right Cauchy-Green tensor\n",
    "E = 0.5*(C-I) #Finite Strain tensor\n",
    "W = lmbda/2*tr(E)**2+mu*tr(E*E) # St. Venant-Kirchhoff Strain-energy density\n",
    "E = variable(E)\n",
    "S = diff(W,E) #Second Piola-Stress tensor\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subdomains\n",
    "```python\n",
    "class GammaN(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        tol = 1E-14\n",
    "        return on_boundary and near(x[0], 0, tol)\n",
    "        \n",
    "class GammaD(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        tol = 1E-14\n",
    "        return on_boundary and near(x[0], 0, tol)\n",
    "NeumannBoundary = GammaN()\n",
    "DirichletBoundary = GammaD()\n",
    "boundary_markers = FacetFunction('size_t', mesh)\n",
    "boundary_markers.set_all(0)\n",
    "DirichletBoundary.mark(boundary_markers,1)\n",
    "NeumannBoundary.mark(boundary_markers,2)\n",
    "ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)\n",
    "bc = DirichletBC(V,Constant(0.0),DirichletBoundary)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving and plotting\n",
    "```python\n",
    "solve(WeakForm==forcingFcn,u,bc)\n",
    "file = File(\"u.pvd\")\n",
    "file << u\n",
    "plot(u)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}