{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Quick overview of the finite element method\n",
"<div id=\"ch:overview\"></div>\n",
"\n",
"<!-- dom:FIGURE: [fig/dolfin_mesh.png, width=500 frac=0.8] Example on a complicated domain for solving PDEs. <div id=\"overview:meshex\"></div> -->\n",
"<!-- begin figure -->\n",
"<div id=\"overview:meshex\"></div>\n",
"\n",
"<p>Example on a complicated domain for solving PDEs.</p>\n",
"<img src=\"fig/dolfin_mesh.png\" width=500>\n",
"\n",
"<!-- end figure -->\n",
"\n",
"\n",
"The finite element method is a rich and versatile approach to construct\n",
"computational schemes to solve any partial differential equation on\n",
"any domain in any dimension. The method may at first glance appear\n",
"cumbersome and even unnatural as it relies on variational formulations\n",
"and polynomial spaces.\n",
"\n",
"Let us start by outlining the concepts briefly.\n",
"Consider the following PDE in 2D:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"-\\nabla^2 u = -u_{xx} - u_{yy} = f,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"equipped with suitable boundary conditions.\n",
"A finite difference scheme to solve the current PDE\n",
"would in the simplest case be described by the stencil"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"overview:2d:fdm0\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\label{overview:2d:fdm0} \\tag{1}\n",
"-\\frac{u_{i-1,j} - 2 u_{i,j} + u_{i+1,j}}{h^2}\n",
"-\\frac{u_{i,j-1} - 2 u_{i,j} + u_{i,j+1}}{h^2}\n",
" = f_{i}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"or reordered to the more recognizable"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"overview:2d:fdm\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\label{overview:2d:fdm} \\tag{2}\n",
"\\frac{-u_{i-1,j} -u_{i,j-1} + 4 u_{i,j} - u_{i+1,j} -u_{i,j+1}}{h^2} = f_{i}\n",
"{\\thinspace .}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On a structured mesh, the stencil appears natural and\n",
"is convenient to implement.\n",
"However, for a unstructured, \"complicated\" domain\n",
"as shown in [Figure](#overview:meshex),\n",
"we would need to be careful when placing\n",
"points and evaluating stencils and functions. \n",
"In particular, \n",
"it will be difficult to evaluate the stencil near the dolphin in \n",
"[Figure](#overview:meshex) because some points will be on the inside and some outside on the outside of the dolphin. \n",
"Both accuracy and efficiency\n",
"may easily be sacrificed by a reckless implementation.\n",
"\n",
"In general, a domain like the one represented in [Figure](#overview:meshex) will be represented by a triangulation. The\n",
"finite element method (and the finite volume method which often is a\n",
"special case of the finite element method) is a methodology for\n",
"creating stencils in a structured manner\n",
"that adapt to the underlying triangulation.\n",
"\n",
"The triangulation in [Figure](#overview:meshex) is a mesh that\n",
"consists of cells that are connected and defined in terms of\n",
"vertices. The fundamental idea of the finite element method is\n",
"to construct a procedure to compute a stencil on a general element and\n",
"then apply this procedure to each element of the mesh. Let\n",
"us therefore denote the mesh as $\\Omega$ while $\\Omega_e$ is the domain\n",
"of a generic element such that $\\Omega=\\cup_e \\Omega_e$.\n",
"\n",
"This is exactly the point where the challenges of the finite element\n",
"method start and where we need some new concepts. The basic question\n",
"is: How should we create a stencil for a\n",
"general element and a general PDE that has the maximal accuracy and\n",
"minimal computational complexity at the current triangulation? The\n",
"two basic building blocks of the finite element method are\n",
"\n",
"1. the solution is represented in terms of a polynomial expression on the\n",
" given general element, and\n",
"\n",
"2. a variational formulation of the PDE\n",
" where element-wise integration enables the PDE to be transformed to a\n",
" stencil.\n",
"\n",
"Step 1 is, as will be explained later, conveniently represented\n",
"both implementation-wise and mathematically as a solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"overview:u:fem\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\label{overview:u:fem} \\tag{3}\n",
"u = \\sum_{i=0}^N c_i {\\psi}_i(x,y),\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $\\{c_i\\}$ are the coefficients to be determined\n",
"(often called the degrees of freedom)\n",
"and ${\\psi}_i(x,y)$ are prescribed polynomials.\n",
"The basis functions ${\\psi}_i(x,y)$ used to express the solution\n",
"is often called the trial functions. \n",
"\n",
"\n",
"The next step is the variational formulation. This step\n",
"may seem like a magic trick or a cumbersome\n",
"mathematical exercise at first glance.\n",
"We take the PDE and multiply by a function $v$ (usually called \n",
"the test function)\n",
"and integrate over an element $\\Omega_e$ and obtain the expression"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"overview:poisson\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\label{overview:poisson} \\tag{4}\n",
"\\int_{\\Omega_e} -\\nabla^2 u \\, v {\\, \\mathrm{d}x} = \\int_{\\Omega_e} f \\, v {\\, \\mathrm{d}x}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A perfectly natural question at this point is: Why multiply\n",
"with a test function $v$? The simple answer is that\n",
"there are $N+1$ unknowns that need to be determined in $u$\n",
"in ([3](#overview:u:fem))\n",
"and for this we need $N+1$ equations. The equations are\n",
"obtained by using $N+1$ different test functions which when used\n",
"in ([5](#overview:fem:a))\n",
"give rise to $N+1$ linearly independent equations.\n",
"\n",
"\n",
"While ([4](#overview:poisson)) is a variational formulation of\n",
"our PDE problem, it is not the most common form.\n",
"It is common to re-write"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"overview:fem:a\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\label{overview:fem:a} \\tag{5}\n",
"\\int_{\\Omega_e} -\\nabla^2 u \\, v {\\, \\mathrm{d}x}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"to weaken the requirement of the polynomial space used for the\n",
"trial functions (that here needs to be twice differentiable)\n",
"and write this term in its corresponding weak form.\n",
"That\n",
"is, the term is rewritten in terms of first-derivatives only (of\n",
"both the trial and the test function) with the aid of Gauss-Green's lemma:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"overview:fem:a:weak\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\label{overview:fem:a:weak} \\tag{6}\n",
"\\int_{\\Omega_e} -\\nabla^2 u \\, v {\\, \\mathrm{d}x} =\n",
"\\int_{\\Omega_e} \\nabla u \\cdot \\nabla v {\\, \\mathrm{d}x} - \\int_{\\partial \\Omega_e} \\frac{\\partial u}{\\partial n} \\, v \\, dS\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The reasons behind this alternative formulation are rather mathematical and will\n",
"not be a major subject of this book as they are well described elsewhere.\n",
"In fact, a precise explanation would need tools from functional analysis.\n",
"\n",
"With the above rewrite and assuming now that the boundary term vanishes due to\n",
"boundary conditions (why this is possible will be dealt with in detail\n",
"later) the\n",
"stencil, corresponding to ([2](#overview:2d:fdm)), is represented by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_{\\Omega_e} \\nabla u \\cdot \\nabla v {\\, \\mathrm{d}x}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $u$ is called the *trial function*, $v$ is called a *test function*,\n",
"and $\\Omega$ is an element of\n",
"a triangulated mesh. The idea of software like FEniCS is that this\n",
"piece of mathematics can be directly expressed in terms of Python code as"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# DO NOT RUN THIS CELL, THIS IS A DEMONSTRATION\n",
"mesh = Mesh(\"some_file\")\n",
"V = FunctionSpace(mesh, \"some polynomial\")\n",
"u = TrialFunction(V)\n",
"v = TestFunction(V)\n",
"a = dot(grad(u), grad(v))*dx"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The methodology and code in this example is not tied to a particular\n",
"equation, except the formula for `a`, holding the derivatives of our\n",
"sample PDE, but any other PDE terms could be expressed via `u`, `v`,\n",
"`grad`, and other symbolic operators in this line of code. In fact,\n",
"finite element packages like FEniCS are typically structured as\n",
"general toolboxes that can be adapted to any PDE as soon as the\n",
"derivation of variational formulations is mastered. The main obstacle\n",
"here for a novice FEM user is then to understand the concept of trial\n",
"functions and test functions realized in terms of polynomial spaces.\n",
"\n",
"\n",
"Hence, a finite element formulation (or a weak formulation) of\n",
"the Poisson problem that works on any mesh $\\Omega$ can be written\n",
"in terms of solving the problem:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_\\Omega\\nabla u\\cdot\\nabla vd{\\, \\mathrm{d}x} = \\int_\\Omega fv{\\, \\mathrm{d}x}{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"By varying the trial and test spaces we obtain different stencils,\n",
"some of which will be identical to finite difference schemes on\n",
"particular meshes. We will now show a complete FEniCS program to\n",
"illustrate how a typical finite element code may be structured"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# DO NOT RUN THIS CELL, THIS IS A DEMONSTRATION\n",
"mesh = Mesh(\"some_file\")\n",
"V = FunctionSpace(mesh, \"some polynomial\")\n",
"u = TrialFunction(V)\n",
"v = TestFunction(V)\n",
"a = dot(grad(u), grad(v))*dx\n",
"L = f*v*dx\n",
"\n",
"bc = DirichletBC(V, \"some_function\", \"some_domain\")\n",
"solution = Function(V) # unknown FEM function\n",
"solve(a == L, solution, bc)\n",
"plot(solution)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- # -->\n",
"<!-- # -->\n",
"<!-- # -->\n",
"<!-- # -->\n",
"\n",
"\n",
"While the finite element method is versatile and may be adapted to any\n",
"PDE on any domain in any dimension, the different methods that are\n",
"derived by using different trial and test functions may vary\n",
"significantly in terms of accuracy and efficiency. In fact, a bad choice of polynomial space may in some cases lead to a\n",
"completely wrong result. This is particularly the case for complicated\n",
"PDEs. For this reason, it is dangerous to regard the method as a black\n",
"box and not do proper verification of the method for a particular\n",
"application.\n",
"\n",
"In our view, there\n",
"are three important tests that should be frequently employed\n",
"during verification:\n",
"\n",
"1. reducing the model problem to 1D and carefully check the calculations involved in the variational formulation on a small 1D mesh\n",
"\n",
"2. perform the calculation involved on one general or random element\n",
"\n",
"3. test whether convergences is obtained and to what order the method converge by refining the mesh\n",
"\n",
"The two first tasks here should ideally be performed by independent calculations\n",
"outside the framework used for the simulations. In our view `sympy` is a\n",
"convenient tool that can be used to assist hand calculations.\n",
"\n",
"So far, we have outlined how the finite element method handles derivatives\n",
"in a PDE, but we also had a right-hand side function $f$. This term is multiplied\n",
"by the test function $v$ as well, such that the entire Poisson equation\n",
"is transformed to"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_\\Omega\\nabla u\\cdot\\nabla vd{\\, \\mathrm{d}x} = \\int_\\Omega fv{\\, \\mathrm{d}x}{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This statement is assumed valid for all test functions $v$ in some\n",
"function space $V$ of polynomials. The right-hand side expression is\n",
"coded in FEniCS as"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# DO NOT RUN THIS CELL. RUN THE FULL PROGRAM (IN THE END) INSTEAD.\n",
"L = f*v*dx"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and the problem is then solved by the statements"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# DO NOT RUN THIS CELL. RUN THE FULL PROGRAM (IN THE END) INSTEAD.\n",
"u = Function(V) # unknown FEM function\n",
"solve(a == L, u, bc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where `bc` holds information about boundary conditions. This information\n",
"is connected to information about the triangulation, the *mesh*.\n",
"Assuming $u=0$ on the boundary, we can in FEniCS generate a triangular\n",
"mesh over a rectangular domain $[-1,-1]\\times [-1,1]$ as follows:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# DO NOT RUN THIS CELL. RUN THE FULL PROGRAM (IN THE END) INSTEAD.\n",
"mesh = RectangleMesh(Point(-1, -1), Point(1, 1), 10, 10)\n",
"bc = DirichletBC(V, 0, 'on_boundary')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Mathematically, the finite element method transforms our PDE to\n",
"a sparse linear system. The `solve` step performs two tasks:\n",
"construction of the linear system based on the given information about\n",
"the domain and its elements, and then solution of the linear system by\n",
"either an iterative or direct method.\n",
"\n",
"We are now in a position to summarize all the parts of a FEniCS program\n",
"that solves the Poisson equation by the finite element method:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"from fenics import *\n",
"mesh = RectangleMesh(Point(-1, -1), Point(1, 1), 10, 10)\n",
"V = FunctionSpace(mesh, 'P', 2) # quadratic polynomials\n",
"bc = DirichletBC(V, 0, 'on_boundary')\n",
"u = TrialFunction(V)\n",
"v = TestFunction(V)\n",
"a = dot(grad(u), grad(v))*dx\n",
"L = f*v*dx\n",
"u = Function(V) # unknown FEM function to be computed\n",
"solve(a == L, u, bc)\n",
"vtkfile = File('poisson.pvd'); vtkfile << u # store solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solving a different PDE is a matter of changing `a` and `L`. \n",
"\n",
"Although we assert here that the finite element method is a tool that\n",
"can solve any PDE problem on any domain of any complexity, the\n",
"fundamental ideas of the method are in fact even more general. \n",
"We will therefore start the book by variational \n",
"methods for approximation in general, then consider the finite\n",
"element in a wide range of applications.\n",
"\n"
]
}
],
"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.7.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}