interpolation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolation\n",
    "\n",
    "`dolfin.interpolate()` computes the standard (local) nodal interpolant, given by\n",
    "\n",
    "$$ I_K (f) = \\sum_{i = 1}^{n_K} L_i (f) \\phi_i $$\n",
    "\n",
    "for all dofs $L_i$ over cell $K$. This means that `interpolate()` cannot work in Hermite function spaces without extending the interface for `ufc::function` to include differentiation. In particular we need to be able to take partial derivatives along each spatial coordinate (in the notation of the paper $L_{\\alpha_i}, 1 \\leq i \\leq d$) of the function $f$ and use these values as the coefficients for the corresponding shape functions $\\phi_{\\alpha_i}$ in the representation in the global basis of $V_h$.\n",
    "\n",
    "This notebook will contain progress in this matter.\n",
    "\n",
    "## Current status (FEniCS 2017.0.1)\n",
    "\n",
    "The default implementation (in `ffc/evaluatedof.py:_generate_body()`) of code generation for\n",
    "\n",
    "```c++\n",
    "double evaluate_dof(std::size_t i,\n",
    "                      const ufc::function& f,\n",
    "                      const double * coordinate_dofs,\n",
    "                      int cell_orientation,\n",
    "                      const ufc::cell& c) const final override\n",
    "```\n",
    "\n",
    "does the following:\n",
    "\n",
    "1. Get coordinate for dof `i`.\n",
    "2. Call `f.evaluate()` on this coordinate.\n",
    "\n",
    "I.e. it assumes that all dofs are point evaluations! What is the best way of implementing `evaluate_dof()` for derivative evaluations? Do I need to extend the interface of GenericFunctions|whatever to include derivatives? Can I use automatic differentiation in UFL?\n",
    "\n",
    "In this notebook:\n",
    "\n",
    "* Nodal interpolation of `Callable`s into Hermite and DKT `FunctionSpace`s works (actually into any space using only `PointEvaluation` and `PointDerivative`). See [utils](utils.py) for the automatic differentiation.\n",
    "* **FIXME:** $H^1$ norms are [very wrong](#Testing-the-norms). This a problem with `norm()` (Hermite coefficients seem to be ok)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A hack using AD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Boilerplate\n",
    "from dolfin import *\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as pl\n",
    "import autograd.numpy as np\n",
    "from utils import make_derivatives, ExpressionAD\n",
    "from FIAT import functional\n",
    "from ffc import fiatinterface\n",
    "import nbimporter\n",
    "\n",
    "def __nbinit__():\n",
    "    global __all__\n",
    "    __all__ = ['interpolate_hermite', 'make_constant']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to be able to take derivatives. A first attempt might be to [use UFL](#Using-UFL) to compute derivatives of `Expressions` but I can't get this to work. Instead, we use actual python functions and compute their derivatives with `autograd` (see [utils.py](utils.py) for the definition of the decorator `make_derivatives`).\n",
    "\n",
    "Here is how we can replace the evaluation of an already defined `Expression`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = FunctionSpace(UnitSquareMesh(10, 10), 'Hermite', 3)\n",
    "f = Expression(\"x[0]*x[1]*x[1]\", degree=3)\n",
    "# This won't work:\n",
    "#v = interpolate(f, V)\n",
    "\n",
    "# Instead define a function which can (automatically) differentiated\n",
    "@make_derivatives\n",
    "def fun2d(x, y):\n",
    "    return x*(y+2)+x**2\n",
    "\n",
    "# Setup the evaluation on a grid\n",
    "tt = np.linspace(0, 1, 10)\n",
    "xx, yy = np.meshgrid(tt, tt)\n",
    "pp = np.dstack((xx,yy))\n",
    "points_shape = pp.shape[:-1]\n",
    "zz = np.zeros(points_shape)\n",
    "\n",
    "# Evaluate the expression using fun2d.\n",
    "# The original definition of the expression is ignored!\n",
    "for i,j in np.ndindex(points_shape):\n",
    "    zz[i,j] = f.evaluate(pp[i,j], {f: fun2d}, (), None, derivatives=(0,1))\n",
    "\n",
    "_ = pl.imshow(zz, aspect='auto', cmap='bone')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That was nice, but instead we want to subclass `Expression` to be able to reuse the declared expressions both in regular FEniCS code and below in our `interpolate_hermite()`, which requires taking derivatives (again, see [utils.py](utils.py) for the definition). Here's how to use it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The function we are going to wrap\n",
    "def f(x, y):\n",
    "    return x*y\n",
    "\n",
    "# Usage\n",
    "e = ExpressionAD(fun=f, degree=3)\n",
    "\n",
    "# Testing\n",
    "val = np.zeros(1)\n",
    "x = np.array([2., 3.])\n",
    "e.eval(val, x)\n",
    "print(val)\n",
    "e.partial(val, x, (1,1))\n",
    "print(val)\n",
    "e(x, derivatives=(1,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we can take derivatives of our `Expressions`, we can interpolate these onto `FunctionSpaces` made of Lagrange and Hermite elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def interpolate_hermite(f:ExpressionAD, V:FunctionSpace) -> Function:\n",
    "    \"\"\" Nodal interpolation for elements using Hermite or Lagrange DOFs.\n",
    "    \n",
    "    This is just a hack and requires an actual python function\n",
    "          which we can evaluate and differentiate.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "        f: Can be one of:\n",
    "           * ExpressionAD wrapper around a python function of tdim arguments.\n",
    "             Example:\n",
    "                 V = FunctionSpace(some2dmesh, 'Hermite', 3)\n",
    "                 fun = lambda x,y: x*y**2\n",
    "                 f = ExpressionAD(fun=f, degree=3)\n",
    "                 u = interpolate_hermite(f, V)\n",
    "\n",
    "           * Callable with signature (np.ndarray, derivatives=()) where the\n",
    "             array is of dimension equal to the topological dimension of the mesh.\n",
    "             Apply the decorator @make_derivatives to a function of n variables\n",
    "             to use autograd's AD to create the callable.\n",
    "             Example:\n",
    "\n",
    "                @make_derivatives\n",
    "                def fun2d(x, y):\n",
    "                    return x*y\n",
    "                u = interpolate_hermite(fun2d, V)\n",
    "                # Equivalently, with the definition of fun above:\n",
    "                v = interpolate_hermite(make_derivatives(fun), V)\n",
    "\n",
    "            Now fun2d() has signature fun2d(np.ndarray, derivatives=())\n",
    "        V: FunctionSpace using Hermite elements.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "        A Function in V.\n",
    "    \"\"\"\n",
    "    warning(\"WARNING: Nodal interpolation for Hermite dofs. HACK HACK.\")\n",
    "    assert callable(f), \"Need a callable (with signature (x, derivatives=())->f(x))\"\n",
    "    \n",
    "    tdim = V.element().topological_dimension()\n",
    "    e = fiatinterface.create_element(V.ufl_element())\n",
    "    \n",
    "    # This mask filters out the Hermite dofs from the list of dofs of a cell:\n",
    "    # Note that we cannot use list_hermite_dofs() because we need to iterate over\n",
    "    # cells to evaluate the partial derivatives in the right order.\n",
    "    hmask = np.array(list(map(lambda f: isinstance(f, functional.PointDerivative),\n",
    "                              e.dual_basis())))\n",
    "    lmask = np.array(list(map(lambda f: isinstance(f, functional.PointEvaluation),\n",
    "                              e.dual_basis())))\n",
    "    assert ((hmask | lmask) == True).size == len(e.dual_basis()),\\\n",
    "           \"Only Lagrange and Hermite dofs supported.\"\n",
    "\n",
    "    #print(\"hmask = %s\" % hmask)\n",
    "    uv = np.zeros(V.dim())\n",
    "    dm = V.dofmap()\n",
    "    dof_coordinates = V.tabulate_dof_coordinates().reshape((-1, tdim))\n",
    "    if tdim == 1:\n",
    "        derivatives = [(1,)]*2\n",
    "    elif tdim == 2:\n",
    "        derivatives = [(1,0),(0,1)]*3\n",
    "    elif tdim == 3:\n",
    "        derivatives = [(1,0,0),(0,1,0),(0,0,1)]*4\n",
    "    for i in range(V.mesh().num_cells()):\n",
    "        dofs = dm.cell_dofs(i)\n",
    "        hermite_dofs = dofs[hmask]\n",
    "        lagrange_dofs = dofs[lmask]\n",
    "        for dd, hdof in zip(derivatives, hermite_dofs):\n",
    "            uv[hdof] = f(dof_coordinates[hdof], derivatives=dd)\n",
    "            #print(\"hdof=%d, coords=%s, derivs=%s\" % (hdof, dof_coordinates[hdof], dd))\n",
    "        for ldof in lagrange_dofs:\n",
    "            uv[ldof] = f(dof_coordinates[ldof], derivatives=())\n",
    "    u = Function(V)\n",
    "    u.vector().set_local(uv)\n",
    "\n",
    "    return u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is how to use that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = FunctionSpace(UnitSquareMesh(20, 20), 'DKT', 3)\n",
    "def fun2d(x, y):\n",
    "    return x*(y+2)+x**2\n",
    "f = ExpressionAD(fun=fun2d, degree=3)\n",
    "u = interpolate_hermite(f, V)\n",
    "v = project(f, V)\n",
    "w = Function(V, u.vector() - v.vector())\n",
    "\n",
    "pl.figure(figsize=(12,4))\n",
    "pl.subplot(1,3,1)\n",
    "plot(u, title=\"Interpolation into H3\", cmap='bone')\n",
    "pl.subplot(1,3,2)\n",
    "plot(v, title=\"Projection onto H3\", cmap='bone')\n",
    "pl.subplot(1,3,3)\n",
    "_ = plot(w, title=\"Difference (ok?)\", cmap='bone')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The difference is very strange:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wv = u.vector() - v.vector()\n",
    "print((wv.min(), wv.max()))\n",
    "\n",
    "xx = np.linspace(0, 1, 20)\n",
    "for y in xx:\n",
    "    pl.plot(xx, [w(x, y) for x in xx])\n",
    "_ = pl.title(\"$u-v$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Constants\n",
    "\n",
    "Constants are no different and need to be interpolated with `interpolate_hermite()`, since `dolfin.interpolate()` would set all coefficients for the Hermite dofs to the value of the constant. Here is a handy shortcut:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_constant(value:float, V:FunctionSpace):\n",
    "    \"\"\" Returns a constant Function over V.\n",
    "    \n",
    "    This only makes sense when V is built with elements using\n",
    "    Hermite dofs. For other FunctionSpaces this just returns\n",
    "    a Constant.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "        value: a scalar.\n",
    "        V: FunctionSpace of elements using Hermite dofs\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "        A Function in V.\n",
    "    \"\"\"\n",
    "    if V.ufl_element().shortstr().lower()[:3] not in ['her', 'dkt']:\n",
    "        return Constant(value)\n",
    "\n",
    "    def fun(x, derivatives=()):\n",
    "        return value if sum(derivatives) == 0 else 0.\n",
    "\n",
    "    return interpolate_hermite(fun, V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Testing the norms\n",
    "\n",
    "** Warning:**\n",
    "\n",
    "Computing the norm using `dolfin.norm()` is probably an error since *all* coefficients are used, including the Hermite ones, so the value obtained is not an approximation of the real $L^2$ norm of the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import Math\n",
    "print(\"The following norms are wrong!\")\n",
    "Math(\"||u_h - u_p||_{L^2} = %.4f, \\  ||u_h - u_p||_{L^{\\infty}} = %.4f\" % \n",
    "     (norm(w), norm(w.vector(), norm_type='linf')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now compare that to the computation using only the values at the vertices (i.e. the Lagrange dofs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "u_vert = u.compute_vertex_values()\n",
    "v_vert = v.compute_vertex_values()\n",
    "print(\"%.2f%% of vertices where u and v disagree more than 1e-6. However:\" % \n",
    "      (100.*np.nonzero(np.abs(u_vert-v_vert) > 1e-6)[0].size / V.dim()))\n",
    "Math(\"\\\\frac{||u_h - u_p||_{L^2}}{\\\\text{max}(||u_h||_{L^2}, ||u_p||_{L^2})} = %.4f, \\  \"\n",
    "     \"\\\\frac{||u_h - u_p||_{L^{\\infty}}}{\\\\text{max}(||u_h||_{L^{\\infty}},||u_p||_{L^{\\infty}})} = %.4f\" % \n",
    "     (np.linalg.norm(u_vert-v_vert)/np.maximum(np.linalg.norm(u_vert), np.linalg.norm(v_vert)),\n",
    "      np.linalg.norm(u_vert-v_vert, ord=np.inf)/np.max(np.maximum(u_vert, v_vert))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yet another example. Notice how wrong the computation of the $H^1$ norm is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ftest(x, y):\n",
    "    return x**2+y**2\n",
    "f = ExpressionAD(fun=ftest, degree=3)\n",
    "V = FunctionSpace(UnitSquareMesh(20, 20), 'DKT', 3)\n",
    "w1 = interpolate_hermite(f, V)\n",
    "w2 = project(f, V)\n",
    "print(\"           w1\\t\\tw2\\t\\texact \" % ())\n",
    "print(\"L2       : %f, %f, %f\" % (norm(w1, 'L2'), norm(w2, 'L2'), np.sqrt(28/45.)))\n",
    "print(\"H10 norms: %f, %f, %f\" % (norm(w1, 'H10'), norm(w2, 'H10'), np.sqrt(8/3.)))\n",
    "print(\"H1 norms : %f, %f, %f\" % (norm(w1, 'H1'), norm(w2, 'H1'), np.sqrt(8/3. + 28/45.)))\n",
    "#Math(\"\\\\text{Exact: }||f||_{L^2} = %f,\\ ||f||_{H^1_0} = %f,\\ ||f||_{H^1} = %f\" %\n",
    "#      (np.sqrt(28/45.), np.sqrt(8/3.), np.sqrt(8/3. + 28/45.)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Attempts using UFL\n",
    "\n",
    "A couple of tests trying to use AD on expressions in order to implement an `interpolate()` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from ufl.algorithms.apply_derivatives import apply_derivatives\n",
    "\n",
    "V = FunctionSpace(UnitIntervalMesh(2), 'Hermite', 3)\n",
    "x = SpatialCoordinate(V.mesh())\n",
    "f = x[0]**2\n",
    "str(f)\n",
    "#replace(f, {x:0.})  # Doesn't work\n",
    "dfdx = apply_derivatives(diff(f, x))\n",
    "#replace(dfdx, {x: 0.})  # Doesn't work"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A test\n",
    "\n",
    "The first thing to check once `evaluate_dof()` has been implemented for Hermite elements is a constant:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "u = interpolate(Constant(2.0), V)\n",
    "_ = plot(u, cmap='bone')"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "218px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}