Poisson1D.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Poisson 1D\n",
    "\n",
    "This notebook solves the linear problem:\n",
    "\n",
    "$$ -u(x)'' = f(x) \\text{ in } \\Omega = [-1,1]$$ \n",
    "\n",
    "with\n",
    "\n",
    "$$ f(x) = \\frac{\\pi^2}{4} \\cos{\\left(\\frac{\\pi x}{2}\\right)} $$\n",
    "\n",
    "and either homogeneous Dirichlet boundary conditions or mixed Dirichlet and Neumann using Hermite elements. The purpose is to find what the issues are while setting Dirichlet BCs and fix them. See [Manually rebuilding the mass matrix](#Manually-rebuilding-the-mass-matrix)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as pl\n",
    "import FIAT\n",
    "%matplotlib inline\n",
    "import nbimporter\n",
    "from interpolation import make_constant\n",
    "\n",
    "# IMPORTANT! I haven't ported anything to UFLACS yet\n",
    "parameters['form_compiler']['representation'] = 'quadrature'\n",
    "# IMPORTANT! Hermite trafo not implemented in optimisedquadraturetransformer.py\n",
    "parameters['form_compiler']['optimize'] = False\n",
    "\n",
    "# Exact solution with homogeneous Dirichlet BCs:\n",
    "u_exact = lambda x: np.cos(np.pi*x/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def poisson1d(V:FunctionSpace, use_neumann:bool=False):\n",
    "    \"\"\" Solves -u'' = f on [-1,1] with f=... \"\"\"\n",
    "    deg = V.ufl_element().degree()\n",
    "    \n",
    "    LEFT = V.mesh().coordinates().min()\n",
    "    RIGHT = V.mesh().coordinates().max()\n",
    "\n",
    "    u0 = make_constant(0., V) # Work around lack of evaluate_dofs() for Hermite\n",
    "\n",
    "    u = TrialFunction(V)\n",
    "    v = TestFunction(V)\n",
    "    f = Expression(\"pi*pi*cos(pi*x[0]/2)/4\", degree=deg)\n",
    "\n",
    "    a = u.dx(0)*v.dx(0)*dx\n",
    "    L = f*v*dx \n",
    "    \n",
    "    if use_neumann: # Note that this changes the exact solution!\n",
    "        class Neumann(SubDomain):\n",
    "            def inside(self, x, on_boundary):\n",
    "                return on_boundary and near(x[0], RIGHT)\n",
    "        NEUMANN_ID = 1\n",
    "        exterior = FacetFunction('uint', V.mesh())\n",
    "        neu_bdry = Neumann()\n",
    "        neu_bdry.mark(exterior, NEUMANN_ID)\n",
    "        bc = DirichletBC(V, u0, lambda x: near(x[0], LEFT))\n",
    "        L = L - v*ds(subdomain_data=exterior, subdomain_id=NEUMANN_ID)\n",
    "    else:\n",
    "        bc = DirichletBC(V, u0, lambda x: near(x[0], LEFT) or near(x[0], RIGHT))\n",
    "\n",
    "    u = Function(V)\n",
    "    solve(a == L, u, bc)\n",
    "    return u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = IntervalMesh(20, -1, 1)\n",
    "\n",
    "V = FunctionSpace(mesh, \"Lagrange\", 3)\n",
    "u_lag2 = poisson1d(V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = np.linspace(-1, 1, 100)\n",
    "pl.plot(xx, [u_lag2(x) for x in xx], label='$u_{l}$')\n",
    "#pl.plot(xx, [u_exact(x) for x in xx], label='$u_{ex}$')\n",
    "pl.title(\"Solution with CG2 elements\")\n",
    "_ = pl.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Hermite elements\n",
    "\n",
    "Naively applying Dirichlet boundary conditions fails because first derivatives are set to 0 when `DirichletBC` clears all the rows of dofs of nodes at the boundary, since this includes Hermite dofs associated to partial derivatives.\n",
    "\n",
    "## A note about the solutions\n",
    "\n",
    "Note that we always plot the solution by evaluating it at the mesh vertices. If instead we choose to evaluate it at arbitrary points in the domain, the solutions in Hermite spaces show a wiggly behaviour. Is this expected or an error? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = FunctionSpace(mesh, \"Hermite\", 3)\n",
    "u = poisson1d(V)\n",
    "\n",
    "xx = V.mesh().coordinates() #np.linspace(-1, 1, 100)\n",
    "pl.figure(figsize=(8,8))\n",
    "pl.plot(xx, [u(x) for x in xx], label='$u_h$')\n",
    "pl.plot(xx, [u_lag2(x) for x in xx], label='$u_{l}$')\n",
    "_ = pl.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Manually rebuilding the stiffness matrix\n",
    "\n",
    "The solution is to manually restore the relevant rows of the assembled matrix to their values before `DirichletBC.apply()` emptied them. The code has been factored out into [boundary.ipynb](boundary.ipynb).\n",
    "\n",
    "Analogously, the way one applies Neumann boundary conditions has to be adapted in e.g. 4th order problems like the beam equation: in that case Neumann bcs are essential and are included in the function space by setting the relevant dofs to the desired values. See [Euler-Bernoulli.ipynb](Euler-Bernoulli.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nbimporter\n",
    "from boundary import *\n",
    "from utils import fnor\n",
    "from interpolation import interpolate_hermite, make_constant\n",
    "\n",
    "def poisson1d_hermite(V, use_neumann:bool=False):\n",
    "    \"\"\"\n",
    "        V: FunctionSpace.\n",
    "    \"\"\"\n",
    "    assert V.ufl_element().family() == 'Hermite', 'Please pay attention.'\n",
    "\n",
    "    _left = V.mesh().coordinates().min()\n",
    "    _right = V.mesh().coordinates().max()\n",
    "\n",
    "    def left(x):\n",
    "        return near(x[0], _left)\n",
    "    def right(x):\n",
    "        return near(x[0], _right)\n",
    "\n",
    "    u0 = make_constant(0, V)\n",
    "\n",
    "    u = TrialFunction(V)\n",
    "    v = TestFunction(V)\n",
    "    f = Expression(\"pi*pi*cos(pi*x[0]/2)/4\",\n",
    "                   degree=V.ufl_element().degree())\n",
    "\n",
    "    a = u.dx(0)*v.dx(0)*dx\n",
    "    F = f*v*dx \n",
    "\n",
    "    if use_neumann:\n",
    "        dir_bc = DirichletBC(V, u0, left)\n",
    "        class Right(SubDomain):\n",
    "            def inside(self, x, on_boundary):\n",
    "                return on_boundary and near(x[0], _right)\n",
    "\n",
    "        exterior = FacetFunction('uint', V.mesh(), value=0)\n",
    "        neu_bdry = Right()\n",
    "        neu_bdry.mark(exterior, 1)\n",
    "        F = F - 1.0*v*ds(subdomain_data=exterior, subdomain_id=1)\n",
    "    else:\n",
    "        dir_bc = DirichletBC(V, u0, fnor(left, right))\n",
    "\n",
    "    A = assemble(a)\n",
    "    b = assemble(F)\n",
    "    apply_dirichlet_hermite(A, b, dir_bc)    \n",
    "    u = Function(V)\n",
    "    U = u.vector()\n",
    "    solve(A, U, b)\n",
    "\n",
    "    return u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mesh = IntervalMesh(10, -1, 1)\n",
    "V = FunctionSpace(mesh, 'Hermite', 3)\n",
    "u_her = poisson1d_hermite(V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = V.mesh().coordinates() # np.linspace(-1,1,100)\n",
    "pl.figure(figsize=(12,6))\n",
    "pl.plot(xx, [u_her(x) for x in xx], '-', label='$u_h$')\n",
    "pl.plot(xx, [u_lag2(x) for x in xx], '-.', label='$u_{lag2}$')#, alpha=0.7)\n",
    "_ = pl.legend()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}