Google Summer of Code 2017 overview. https://summerofcode.withgoogle.com/projects/#5847002177863680

gsoc_overview.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Develop assembly of finite element forms on quadrilateral and hexahedral meshes in FEniCS",
    "\n",
    "### by Ivan Yashchuk" 
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook contains the overview of the work done for the Google Summer of Code 2017.\n",
    "\n",
    "This project involved implementing [Lagrange finite elements](http://femwiki.wikidot.com/elements:lagrange-elements)\n",
    "for quadrilaterals and hexahedrons. The project idea aimed at being able to assemble and solve the simplest PDE, a Poisson's equation, in 2D (quadrilateral mesh) and 3D (hexahedral mesh) in FEniCS. FEniCS Project consists of [several repositories](https://bitbucket.org/fenics-project/), main changes were added to [FIAT](https://bitbucket.org/fenics-project/fiat) and [FFC](https://bitbucket.org/fenics-project/ffc).\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### List of the contributions for the GSoC project:\n",
    "\n",
    "*DOLFIN*  \n",
    " • Merged - https://bitbucket.org/fenics-project/dolfin/pull-requests/363  \n",
    " • Rebased and Merged - https://bitbucket.org/fenics-project/dolfin/pull-requests/371  \n",
    " • Merged - https://bitbucket.org/fenics-project/dolfin/pull-requests/378  \n",
    " • Merged - https://bitbucket.org/fenics-project/dolfin/pull-requests/381  \n",
    " • Merged - https://bitbucket.org/fenics-project/dolfin/pull-requests/383  \n",
    " • Merged - https://bitbucket.org/fenics-project/dolfin/pull-requests/385  \n",
    "*FIAT*  \n",
    " • Closed - https://bitbucket.org/fenics-project/fiat/pull-requests/38  \n",
    " • Merged - https://bitbucket.org/fenics-project/fiat/pull-requests/39  \n",
    "*FFC*  \n",
    " • Closed - https://bitbucket.org/fenics-project/ffc/pull-requests/73  \n",
    " • Merged - https://bitbucket.org/fenics-project/ffc/pull-requests/77  \n",
    " • Closed -  https://bitbucket.org/fenics-project/ffc/pull-requests/82  \n",
    " • Merged - https://bitbucket.org/fenics-project/ffc/pull-requests/83  \n",
    " • Merged - https://bitbucket.org/fenics-project/ffc/pull-requests/84  \n",
    " • Merged - https://bitbucket.org/fenics-project/ffc/pull-requests/85  \n",
    " • Open - https://bitbucket.org/fenics-project/ffc/pull-requests/87  \n",
    "*UFL*  \n",
    " • Merged - https://bitbucket.org/fenics-project/ufl/pull-requests/77\n",
    " ___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Description"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many parts to assemble and solve on quad/hex meshes were already in FEniCS,\n",
    "but there were missing links to get the pieces working as a whole."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It was decided to construct Lagrange finite elements on quadrilateral and hexahedron as a tensor product of 1D Lagrange elements on interval, which is [mathematically correct](http://hplgit.github.io/num-methods-for-PDEs/doc/pub/approx/html/._approx009.html).\n",
    "\n",
    "FIAT has already had a class for constructing [`TensorProductElements`](http://fenics.readthedocs.io/projects/fiat/en/latest/api-doc/FIAT.html#FIAT.tensor_product.TensorProductElement)\n",
    "which creates a finite element that is the tensor product of two existing finite elements.\n",
    "\n",
    "It was possible to create quad and hex element in FIAT:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import FIAT\n",
    "\n",
    "# Create interval reference cell\n",
    "interval_ref_el = FIAT.reference_element.UFCInterval()\n",
    "# Define polynomial degree for the element\n",
    "degree = 1\n",
    "# Create 1D Lagrange element\n",
    "interval_element = FIAT.lagrange.Lagrange(interval_ref_el, degree)\n",
    "# Create quadrilateral element\n",
    "quad_element = FIAT.tensor_product.TensorProductElement(interval_element, interval_element)\n",
    "# Create hexahedral element\n",
    "hex_element = FIAT.tensor_product.TensorProductElement(quad_element, interval_element)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reference cell of the `quad_element` and `hex_element` clearly has the shape of quadrilateral and hexahedron."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((0.0, 0.0), (0.0, 1.0), (1.0, 0.0), (1.0, 1.0))\n"
     ]
    }
   ],
   "source": [
    "# Print coordinates of reference cell vertices\n",
    "print(quad_element.ref_el.get_vertices())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((0.0, 0.0, 0.0), (0.0, 0.0, 1.0), (0.0, 1.0, 0.0), (0.0, 1.0, 1.0), (1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (1.0, 1.0, 0.0), (1.0, 1.0, 1.0))\n"
     ]
    }
   ],
   "source": [
    "# Print coordinates of reference cell vertices\n",
    "print(hex_element.ref_el.get_vertices())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, FFC automatic code generation did not work for these elements. \n",
    "\n",
    "The first step was to build an appropriate interface between FFC and FIAT so that C++ code could be generated by FFC\n",
    "to evaluate the tensor product finite element basis functions on the quad/hex cell geometry.\n",
    "\n",
    "First rough implementation that took care of getting the data from FIAT and computing intermediate representation for the code generation for quadhex case in FFC was submitted in [FFC PR #73](https://bitbucket.org/fenics-project/ffc/pull-requests/73) and [FIAT PR #38](https://bitbucket.org/fenics-project/fiat/pull-requests/38). Both declined as the changes were moved to other places in the later pull requests."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These changes made possible assembling the simplest form on quadhex mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n"
     ]
    }
   ],
   "source": [
    "# Importing FEniCS components\n",
    "from fenics import *\n",
    "# Creating unit square domain devided into 4 x 6 cells\n",
    "mesh = UnitQuadMesh.create(4, 6)\n",
    "# Assembling the expression. It should return 1 for unit square domain\n",
    "print(assemble(1.0 * dx(mesh)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assembling functions defined as `Expression` and interpolation was fixed with [FFC PR #77](https://bitbucket.org/fenics-project/ffc/pull-requests/77/).\n",
    "\n",
    "Topological dimension of TensorProductElement's reference cell is stored as the tuple of dimensions in different “directions”, it comes from the tensor product operation, FFC can work only if topological dimension is an integer. This issue is described in my blogpost [here](https://ivanyashchuk.github.io/post/numbering-of-mesh-entities/). [A new wrapper class](http://fenics.readthedocs.io/projects/fiat/en/latest/api-doc/FIAT.html#FIAT.tensor_product.FlattenedDimensions) that reconstructs a FIAT element defined on a [`TensorProductCell`](http://fenics.readthedocs.io/projects/fiat/en/latest/api-doc/FIAT.html#FIAT.reference_element.TensorProductCell) to be defined on [`Quadrilateral`](http://fenics.readthedocs.io/projects/fiat/en/latest/api-doc/FIAT.html#FIAT.reference_element.UFCQuadrilateral) or [`Hexahedron`](http://fenics.readthedocs.io/projects/fiat/en/latest/api-doc/FIAT.html#FIAT.reference_element.UFCHexahedron) cells was added with [FIAT PR #39](https://bitbucket.org/fenics-project/fiat/pull-requests/39). Short description of the class is in [my blogpost](https://ivanyashchuk.github.io/post/meet-new-class-for-FIAT/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Changes in FFC to support new FIAT class were added in [FFC PR #83](https://bitbucket.org/fenics-project/ffc/pull-requests/83).\n",
    "\n",
    "The issue with [`DirichletBC`](https://ivanyashchuk.github.io/post/dirichlet-bc/) was fixed with [DOLFIN PR #371](https://bitbucket.org/fenics-project/dolfin/pull-requests/371) and [FFC PR #84](https://bitbucket.org/fenics-project/ffc/pull-requests/84). The problem was that definition of vertices and facets for quadrilateral and hexahedron cells was different in different parts of FEniCS.\n",
    "\n",
    "With the above introduced changes it is possible to solve the Poisson's equation from the [demo](http://fenics.readthedocs.io/projects/dolfin/en/latest/demos/poisson/python/demo_poisson.py.html) just by changing the mesh to one consisting of quadrilaterals or hexahedrons.\n",
    "\n",
    "Projection of `Constant` function onto the Lagrange FunctionSpace for quads and hexes was fixed in [FFC PR #85](https://bitbucket.org/fenics-project/ffc/pull-requests/85) allowing to solve Stokes or Hyperelasticity type of problems correctly, where it is common to have constant loading."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the result of this GSoC project, `Q` and `dQ` elements (see [Periodic Table of the Finite Elements](https://femtable.org/)) were implemented in FEniCS. All needed changes are in master branch already. One can try it with any of [demos](https://bitbucket.org/fenics-project/dolfin/src/0c347b1f22e0/demo/?at=master) which uses `CG` or `DG` FunctionSpace, by changing the mesh to one consisting of quadrilaterals instead of triangles or hexahedrons instead of tetrahedrons.\n",
    "```python\n",
    "quadrilateral_mesh = UnitQuadMesh.create(n, m)\n",
    "\n",
    "hexahedral_mesh = UnitHexMesh.create(n, m , k)\n",
    "```\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assembling and solving the PDE on quad and hex mesh also work in parallel!\n",
    "\n",
    "Weak scaling perfomance with ~640000 DOFs per process was tested for hex mesh and the results are [here](https://ivanyashchuk.github.io/post/weak-scaling/).\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Further improvements:\n",
    "\n",
    " • Function evaluation at a point does not work\n",
    "\n",
    "```python\n",
    "from fenics import *\n",
    "mesh = UnitQuadMesh.create(3, 3)\n",
    "V = FunctionSpace(mesh, 'P', 1)\n",
    "u = Function(V)\n",
    "u(0.1, 0.1) # Returns \"not implemented\" error\n",
    "```\n",
    "It requires `evaluate_basis` method of the finite element class, which is not implemented. [`CiarletElement`](http://fenics.readthedocs.io/projects/fiat/en/latest/api-doc/FIAT.html#FIAT.finite_element.CiarletElement) interface shall be implemented for elements on quadrilateral and hexahedron cells. Changes are needed in `FIAT/expansions.py`, `FIAT/polynomial_set.py`, `ffc/uflacs/backends/ufc/jacobian.py`.\n",
    "\n",
    " • Discontinuous Galerkin Poisson [demo](https://bitbucket.org/fenics-project/dolfin/src/0c347b1f22e0f3ca818135d7e0d2cff4f456bb75/demo/undocumented/dg-poisson/python/demo_dg-poisson.py?at=master&fileviewer=file-view-default) does not work as computation of `Circumradius` is not implemented for quad and hex cells. Therefore `h = CellSize(mesh)` does not work.\n",
    "\n",
    " • Many parts of DOLFIN that interact with quadrilateral and hexahedral mesh is not implemented yet. Collision detection, mesh refinement.\n",
    "\n",
    " • Changing [`TensorProductElement`](http://fenics.readthedocs.io/projects/fiat/en/latest/api-doc/FIAT.html#FIAT.tensor_product.TensorProductElement) class to accept list of FiniteElements will simplify the code to generate tensor product element of more than two elements. \n",
    "\n",
    "\n",
    "To conclude, the original goals of the project are met and exceeded, so I consider GSoC as successfully completed. I hope the community will appreciate my work and contributions :)\n",
    "___\n",
    "Some additional information about my journey through the GSoC 2017 can be found in [my blog](https://ivanyashchuk.github.io/).\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Demonstration\n",
    "\n",
    "\n",
    "[Cahn–Hilliard equation](https://en.wikipedia.org/wiki/Cahn%E2%80%93Hilliard_equation) solved on a quadrilateral mesh:\n",
    "\n",
    "![alt text](https://drive.google.com/uc?export=view&id=0B7_NVEEiR5wUV00tSXhVTWVXeGs \"Cahn–Hilliard equation gif\")\n",
    "\n",
    "[Hyperelasticity problem](http://fenics.readthedocs.io/projects/dolfin/en/latest/demos/hyperelasticity/python/demo_hyperelasticity.py.html) solved on a hexahedral mesh:\n",
    "![alt text](https://drive.google.com/uc?export=view&id=0B7_NVEEiR5wUakhCTGhkSENiOEk \"Twisted cube\")\n",
    "\n",
    "Lid driven cavity Stokes flow solved on a quadrilateral mesh:\n",
    "![alt text](https://drive.google.com/uc?export=view&id=0B7_NVEEiR5wUQ3psTHdEeC1xVWs \"Cavity flow streamlines\")\n",
    "\n",
    "And [Poisson's equation](http://fenics.readthedocs.io/projects/dolfin/en/latest/demos/poisson/python/demo_poisson.py.html) solved on a quad mesh:\n",
    "![alt text](https://drive.google.com/uc?export=view&id=0B7_NVEEiR5wUSUlqSXY0NkJjV0k \"Poisson solution\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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": 2
}