5- varform-fe.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variational formulations with finite elements\n",
    "<div id=\"ch:varform:fe\"></div>\n",
    "\n",
    "\n",
    "We shall now take the ideas from previous chapters and put together such that\n",
    "we can solve PDEs using the flexible finite element basis functions.\n",
    "This is quite a machinery with many details, but the chapter is mostly an\n",
    "assembly of concepts and details we have already met.\n",
    "\n",
    "\n",
    "# Computing with finite elements\n",
    "<div id=\"fem:deq:1D:fem1\"></div>\n",
    "\n",
    "The purpose of this section is to demonstrate in detail how\n",
    "the finite element method can then be applied to the model problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "-u''(x) = 2,\\quad x\\in (0,L),\\ u(0)=u(L)=0,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with variational formulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(u',v') = (2,v)\\quad\\forall v\\in V{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any $v\\in V$ must obey $v(0)=v(L)=0$ because of the Dirichlet conditions\n",
    "on $u$.\n",
    "\n",
    "## Finite element mesh and basis functions\n",
    "\n",
    "We introduce a finite element mesh with $N_e$ cells, all\n",
    "with length $h$, and number\n",
    "the cells from left to right.\n",
    "Choosing P1 elements, there are two\n",
    "nodes per cell, and the coordinates of the nodes become"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "x_{i} = i h,\\quad h=L/N_e,\\quad i=0,\\ldots,N_n-1=N_e{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any node $i$ is associated with a finite element basis function\n",
    "${\\varphi}_i(x)$.  When approximating a given function $f$ by a finite\n",
    "element function $u$, we expand $u$ using finite element basis\n",
    "functions associated with *all* nodes in the mesh. The parameter\n",
    "$N$, which counts the unknowns from 0 to $N$, is then equal to\n",
    "$N_n-1$ such that the total number of unknowns, $N+1$, is the\n",
    "total number of nodes.\n",
    "However, when solving differential equations we will often have\n",
    "$N<N_n-1$ because of Dirichlet boundary conditions. The reason is\n",
    "simple: we know what $u$ are at some (here two) nodes, and the number of\n",
    "unknown parameters is naturally reduced.\n",
    "\n",
    "In our case with homogeneous Dirichlet boundary conditions, we do not\n",
    "need any boundary function $B(x)$, so we can work with the expansion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:fem1:ex:u\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u(x) = \\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_j(x){\\thinspace .}\n",
    "\\label{fem:deq:1D:fem1:ex:u} \\tag{64}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because of the boundary conditions, we must demand\n",
    "${\\psi}_i(0)={\\psi}_i(L)=0$, $i\\in{\\mathcal{I}_s}$. When ${\\psi}_i$ for all\n",
    "$i=0,\\ldots,N$ is to be selected among the finite element basis\n",
    "functions ${\\varphi}_j$, $j=0,\\ldots,N_n-1$, we have to avoid using\n",
    "${\\varphi}_j$ functions that do not vanish at $x_{0}=0$ and\n",
    "$x_{N_n-1}=L$. However, all ${\\varphi}_j$ vanish at these two nodes for\n",
    "$j=1,\\ldots,N_n-2$.  Only basis functions associated with the end nodes,\n",
    "${\\varphi}_0$ and ${\\varphi}_{N_n-1}$, violate the boundary conditions of\n",
    "our differential equation. Therefore, we select the basis functions\n",
    "${\\varphi}_i$ to be the set of finite element basis functions associated\n",
    "with all the interior nodes in the mesh:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "{\\psi}_i={\\varphi}_{i+1},\\quad i=0,\\ldots,N{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $i$ index runs over all the unknowns $c_i$ in the expansion\n",
    "for $u$, and in this case $N=N_n-3$.\n",
    "\n",
    "In the general case, and in particular on domains in higher dimensions,\n",
    "the nodes are not necessarily numbered from left\n",
    "to right, so we introduce a mapping from the node numbering, or more\n",
    "precisely the degree of freedom numbering, to the numbering of\n",
    "the unknowns in the final equation system. These unknowns take on\n",
    "the numbers $0,\\ldots,N$. Unknown number $j$ in the linear system\n",
    "corresponds to degree of freedom number $\\nu (j)$, $j\\in{\\mathcal{I}_s}$.\n",
    "We can then write"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "{\\psi}_i={\\varphi}_{\\nu(i)},\\quad i=0,\\ldots,N{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a regular numbering as in the present example,\n",
    "$\\nu(j) = j+1$, $j=0,\\ldots,N=N_n-3$.\n",
    "\n",
    "\n",
    "## Computation in the global physical domain\n",
    "<div id=\"fem:deq:1D:comp:global\"></div>\n",
    "\n",
    "\n",
    "We shall first perform a computation in the $x$\n",
    "coordinate system because the integrals can be easily computed\n",
    "here by simple, visual,\n",
    "geometric considerations. This is called a global approach\n",
    "since we work in the $x$ coordinate system and compute integrals on\n",
    "the global domain $[0,L]$.\n",
    "\n",
    "The entries in the coefficient matrix and right-hand side are"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A_{i,j}=\\int_0^L{\\psi}_i'(x){\\psi}_j'(x) {\\, \\mathrm{d}x},\\quad\n",
    "b_i=\\int_0^L2{\\psi}_i(x) {\\, \\mathrm{d}x}, \\quad i,j\\in{\\mathcal{I}_s}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expressed in terms of finite element basis functions ${\\varphi}_i$ we\n",
    "get the alternative expressions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A_{i,j}=\\int_0^L{\\varphi}_{i+1}'(x){\\varphi}_{j+1}'(x) {\\, \\mathrm{d}x},\\quad\n",
    "b_i=\\int_0^L2{\\varphi}_{i+1}(x) {\\, \\mathrm{d}x},\\quad i,j\\in{\\mathcal{I}_s}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the following calculations the subscripts on the finite\n",
    "element basis functions are more conveniently written as\n",
    "$i$ and $j$ instead of $i+1$ and $j+1$, so our notation becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A_{i-1,j-1}=\\int_0^L{\\varphi}_{i}'(x){\\varphi}_{j}'(x) {\\, \\mathrm{d}x},\\quad\n",
    "b_{i-1}=\\int_0^L2{\\varphi}_{i}(x) {\\, \\mathrm{d}x},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where the $i$ and $j$ indices run as $i,j=1,\\ldots,N+1=N_n-2$.\n",
    "\n",
    "The ${\\varphi}_i(x)$ function is a hat function with peak at $x=x_{i}$\n",
    "and a linear variation in $[x_{i-1},x_{i}]$ and\n",
    "$[x_{i},x_{i+1}]$.\n",
    "The derivative is $1/h$ to the left of $x_{i}$ and $-1/h$ to\n",
    "the right, or more formally,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:approx:fe:Dphi:1:formula2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "{\\varphi}_i'(x) = \\left\\lbrace\\begin{array}{ll}\n",
    "0, & x < x_{i-1},\\\\ \n",
    "h^{-1},\n",
    "& x_{i-1} \\leq x < x_{i},\\\\ \n",
    "-h^{-1},\n",
    "& x_{i} \\leq x < x_{i+1},\\\\ \n",
    "0, & x\\geq x_{i+1}\n",
    "\\end{array}\n",
    "\\right.\n",
    "\\label{fem:approx:fe:Dphi:1:formula2} \\tag{65}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Figure](#fem:approx:fe:fig:dP1:2) shows ${\\varphi}_2'(x)$ and ${\\varphi}_3'(x)$.\n",
    "\n",
    "\n",
    "<!-- dom:FIGURE: [fig/fe_mesh1D_dphi_2_3.png, width=400]  Illustration of the derivative of piecewise linear basis functions associated with nodes in cell 2.  <div id=\"fem:approx:fe:fig:dP1:2\"></div> -->\n",
    "<!-- begin figure -->\n",
    "<div id=\"fem:approx:fe:fig:dP1:2\"></div>\n",
    "\n",
    "<p>Illustration of the derivative of piecewise linear basis functions associated with nodes in cell 2.</p>\n",
    "<img src=\"fig/fe_mesh1D_dphi_2_3.png\" width=400>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "We realize that ${\\varphi}_i'$ and ${\\varphi}_j'$ has no overlap, and hence their\n",
    "product vanishes, unless $i$ and $j$ are nodes belonging to the same\n",
    "cell. The only nonzero contributions to the coefficient matrix are\n",
    "therefore"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "A_{i-1,i-2} &=\\int_0^L{\\varphi}_i'(x) {\\varphi}_{i-1}'(x) {\\, \\mathrm{d}x},\\\\ \n",
    "A_{i-1,i-1}&=\\int_0^L{\\varphi}_{i}'(x)^2 {\\, \\mathrm{d}x}, \\\\ \n",
    "A_{i-1,i}&=\\int_0^L{\\varphi}_{i}'(x){\\varphi}_{i+1}'(x) {\\, \\mathrm{d}x},\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for $i=1,\\ldots,N+1$, but for $i=1$, $A_{i-1,i-2}$ is not defined,\n",
    "and for $i=N+1$, $A_{i-1,i}$ is not defined.\n",
    "\n",
    "From [Figure](#fem:approx:fe:fig:dP1:2),\n",
    "we see that ${\\varphi}_{i-1}'(x)$ and ${\\varphi}_i'(x)$ have overlap of one\n",
    "cell $\\Omega^{(i-1)}=[x_{i-1},x_{i}]$ and that their product\n",
    "then is $-1/h^{2}$. The integrand is constant and therefore\n",
    "$A_{i-1,i-2}=-h^{-2}h=-h^{-1}$.\n",
    "A similar reasoning can be applied to\n",
    "$A_{i-1,i}$, which also becomes $-h^{-1}$. The integral of\n",
    "${\\varphi}_i'(x)^2$ gets contributions from two cells,\n",
    "$\\Omega^{(i-1)}=[x_{i-1},x_{i}]$ and\n",
    "$\\Omega^{(i)}=[x_{i},x_{i+1}]$, but ${\\varphi}_i'(x)^2=h^{-2}$ in\n",
    "both cells, and the length of the integration interval is $2h$ so\n",
    "we get\n",
    "$A_{i-1,i-1}=2h^{-1}$.\n",
    "\n",
    "The right-hand side involves an integral of $2{\\varphi}_i(x)$,\n",
    "$i=1,\\ldots,N_n-2$,\n",
    "which is just the area under a hat function of height 1 and width\n",
    "$2h$, i.e., equal to $h$. Hence, $b_{i-1}=2h$.\n",
    "\n",
    "To summarize the linear system, we switch from $i$ to $i+1$ such that\n",
    "we can write"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A_{i,i-1}=A_{i,i+1}=-h^{-1},\\quad A_{i,i}=2h^{-1},\\quad\n",
    "b_i = 2h{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equation system to be solved only involves the unknowns\n",
    "$c_i$ for $i\\in{\\mathcal{I}_s}$. With our numbering of unknowns and\n",
    "nodes, we have that $c_i$ equals $u(x_{i+1})$.\n",
    "The complete matrix system then takes the following form:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:ex1:Ab:glob\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{1}{h}\\left(\n",
    "\\begin{array}{ccccccccc}\n",
    "1 & -1 & 0 &\\cdots & \\cdots & \\cdots & \\cdots & \\cdots & 0 \\\\ \n",
    "-1 & 2 & -1 & \\ddots &   & &  & &  \\vdots \\\\ \n",
    "0 & -1 & 2 & -1 &\n",
    "\\ddots & &  &  & \\vdots \\\\ \n",
    "\\vdots & \\ddots &  & \\ddots & \\ddots & 0 &  & & \\vdots \\\\ \n",
    "\\vdots &  & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & & \\vdots \\\\ \n",
    "\\vdots & &  & 0 & -1 & 2 & -1 & \\ddots & \\vdots \\\\ \n",
    "\\vdots & & &  & \\ddots & \\ddots & \\ddots &\\ddots  & 0 \\\\ \n",
    "\\vdots & & & &  &\\ddots  & \\ddots &\\ddots  & -1 \\\\ \n",
    "0 &\\cdots & \\cdots &\\cdots & \\cdots & \\cdots  & 0 & -1 & 1\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "c_0 \\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots\\\\ \n",
    "c_{N}\n",
    "\\end{array}\n",
    "\\right)\n",
    "=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "2h \\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots\\\\ \n",
    "2h\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\label{fem:deq:1D:ex1:Ab:glob} \\tag{66}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison with a finite difference discretization\n",
    "<div id=\"fem:deq:1D:fdm_vs_fem\"></div>\n",
    "\n",
    "A typical row in the matrix system ([66](#fem:deq:1D:ex1:Ab:glob))\n",
    "can be written as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:fem:ex1:c\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "-\\frac{1}{h}c_{i-1} + \\frac{2}{h}c_{i} - \\frac{1}{h}c_{i+1} = 2h{\\thinspace .}\n",
    "\\label{fem:deq:1D:fem:ex1:c} \\tag{67}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us introduce the notation $u_j$ for the value of $u$ at node $j$:\n",
    "$u_j=u(x_{j})$, since we have the interpretation\n",
    "$u(x_{j})=\\sum_jc_j{\\varphi}(x_{j})=\\sum_j c_j\\delta_{ij}=c_j$.\n",
    "The unknowns $c_0,\\ldots,c_N$ are $u_1,\\ldots,u_{N_n-2}$.\n",
    "Shifting $i$ with $i+1$ in ([67](#fem:deq:1D:fem:ex1:c)) and inserting\n",
    "$u_i = c_{i-1}$, we get"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:fem:ex1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "-\\frac{1}{h}u_{i-1} + \\frac{2}{h}u_{i} - \\frac{1}{h}u_{i+1} = 2h,\n",
    "\\label{fem:deq:1D:fem:ex1} \\tag{68}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A finite difference discretization of $-u''(x)=2$ by a centered,\n",
    "second-order finite difference approximation $u''(x_i)\\approx [D_x D_x u]_i$\n",
    "with $\\Delta x = h$\n",
    "yields"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto35\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "-\\frac{u_{i-1} - 2u_{i} + u_{i+1}}{h^2} = 2,\n",
    "\\label{_auto35} \\tag{69}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which is, in fact, equal to ([68](#fem:deq:1D:fem:ex1)) if\n",
    "([68](#fem:deq:1D:fem:ex1)) is divided by $h$.\n",
    "Therefore, the finite difference and the finite element method are\n",
    "equivalent in this simple test problem.\n",
    "\n",
    "Sometimes a finite element method generates the finite difference\n",
    "equations on a uniform mesh, and sometimes the finite element method\n",
    "generates equations that are different.  The differences are modest,\n",
    "but may influence the numerical quality of the solution significantly,\n",
    "especially in time-dependent problems.\n",
    "It depends on the problem at hand\n",
    "whether a finite element discretization is more or less accurate than\n",
    "a corresponding finite difference discretization.\n",
    "<!-- There will be many examples illustrating this point. -->\n",
    "\n",
    "## Cellwise computations\n",
    "<div id=\"fem:deq:1D:comp:elmwise\"></div>\n",
    "\n",
    "Software for finite element computations normally employs\n",
    "the cell by cell computational procedure where\n",
    "an element matrix and vector are calculated for each cell and\n",
    "assembled in the global linear system.\n",
    "Let us go through the details of this type of algorithm.\n",
    "\n",
    "All integrals are mapped to the local reference coordinate system\n",
    "$X\\in [-1,1]$.\n",
    "In the present case, the matrix entries contain derivatives\n",
    "with respect to $x$,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A_{i-1,j-1}^{(e)}=\\int_{\\Omega^{(e)}} {\\varphi}_i'(x){\\varphi}_j'(x) {\\, \\mathrm{d}x}\n",
    "= \\int_{-1}^1 \\frac{d}{dx}{\\tilde{\\varphi}}_r(X)\\frac{d}{dx}{\\tilde{\\varphi}}_s(X)\n",
    "\\frac{h}{2} {\\, \\mathrm{d}X},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where the global degree of freedom $i$ is related to the local\n",
    "degree of freedom $r$ through $i=q(e,r)$. Similarly,\n",
    "$j=q(e,s)$. The local degrees of freedom run as $r,s=0,1$ for a P1\n",
    "element.\n",
    "\n",
    "### The integral for the element matrix\n",
    "\n",
    "There are simple formulas for the basis functions ${\\tilde{\\varphi}}_r(X)$ as\n",
    "functions of $X$.\n",
    "However, we now\n",
    "need to find the derivative of ${\\tilde{\\varphi}}_r(X)$ with respect to $x$.\n",
    "Given"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "{\\tilde{\\varphi}}_0(X)=\\frac{1}{2}(1-X),\\quad{\\tilde{\\varphi}}_1(X)=\\frac{1}{2}(1+X),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we can easily compute $d{\\tilde{\\varphi}}_r/ dX$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{d{\\tilde{\\varphi}}_0}{dX} = -\\frac{1}{2},\\quad  \\frac{d{\\tilde{\\varphi}}_1}{dX} = \\frac{1}{2}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the chain rule,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto36\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{d{\\tilde{\\varphi}}_r}{dx} = \\frac{d{\\tilde{\\varphi}}_r}{dX}\\frac{dX}{dx}\n",
    "= \\frac{2}{h}\\frac{d{\\tilde{\\varphi}}_r}{dX}{\\thinspace .}   \\label{_auto36} \\tag{70}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The transformed integral is then"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A_{i-1,j-1}^{(e)}=\\int_{\\Omega^{(e)}} {\\varphi}_i'(x){\\varphi}_j'(x) {\\, \\mathrm{d}x}\n",
    "= \\int_{-1}^1 \\frac{2}{h}\\frac{d{\\tilde{\\varphi}}_r}{dX}\\frac{2}{h}\\frac{d{\\tilde{\\varphi}}_s}{dX}\n",
    "\\frac{h}{2} {\\, \\mathrm{d}X}\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The integral for the element vector\n",
    "\n",
    "The right-hand side is transformed according to"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "b_{i-1}^{(e)} = \\int_{\\Omega^{(e)}} 2{\\varphi}_i(x) {\\, \\mathrm{d}x} =\n",
    "\\int_{-1}^12{\\tilde{\\varphi}}_r(X)\\frac{h}{2} {\\, \\mathrm{d}X},\\quad i=q(e,r),\\ r=0,1\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Detailed calculations of the element matrix and vector\n",
    "\n",
    "Specifically for P1 elements we arrive at the following calculations for\n",
    "the element matrix entries:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\tilde A_{0,0}^{(e)} &= \\int_{-1}^1\\frac{2}{h}\\left(-\\frac{1}{2}\\right)\n",
    "\\frac{2}{h}\\left(-\\frac{1}{2}\\right)\\frac{h}{2} {\\, \\mathrm{d}X} = \\frac{1}{h}\\\\ \n",
    "\\tilde A_{0,1}^{(e)} &= \\int_{-1}^1\\frac{2}{h}\\left(-\\frac{1}{2}\\right)\n",
    "\\frac{2}{h}\\left(\\frac{1}{2}\\right)\\frac{h}{2} {\\, \\mathrm{d}X} = -\\frac{1}{h}\\\\ \n",
    "\\tilde A_{1,0}^{(e)} &= \\int_{-1}^1\\frac{2}{h}\\left(\\frac{1}{2}\\right)\n",
    "\\frac{2}{h}\\left(-\\frac{1}{2}\\right)\\frac{h}{2} {\\, \\mathrm{d}X} = -\\frac{1}{h}\\\\ \n",
    "\\tilde A_{1,1}^{(e)} &= \\int_{-1}^1\\frac{2}{h}\\left(\\frac{1}{2}\\right)\n",
    "\\frac{2}{h}\\left(\\frac{1}{2}\\right)\\frac{h}{2} {\\, \\mathrm{d}X} = \\frac{1}{h}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The element vector entries become"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\tilde b_0^{(e)} &= \\int_{-1}^12\\frac{1}{2}(1-X)\\frac{h}{2} {\\, \\mathrm{d}X} = h\\\\ \n",
    "\\tilde b_1^{(e)} &= \\int_{-1}^12\\frac{1}{2}(1+X)\\frac{h}{2} {\\, \\mathrm{d}X} = h{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expressing these entries in matrix and vector notation, we have"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:ex1:Ab:elm\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\tilde A^{(e)} =\\frac{1}{h}\\left(\\begin{array}{rr}\n",
    "1 & -1\\\\ \n",
    "-1 & 1\n",
    "\\end{array}\\right),\\quad\n",
    "\\tilde b^{(e)} = h\\left(\\begin{array}{c}\n",
    "1\\\\ \n",
    "1\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "\\label{fem:deq:1D:ex1:Ab:elm} \\tag{71}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Contributions from the first and last cell\n",
    "\n",
    "The first and last cell involve only one unknown and one basis function\n",
    "because of the Dirichlet boundary conditions at the first and last\n",
    "node.\n",
    "The element matrix therefore becomes a $1\\times 1$ matrix and there\n",
    "is only one entry in the element vector. On cell 0, only ${\\psi}_0={\\varphi}_1$\n",
    "is involved, corresponding to integration with ${\\tilde{\\varphi}}_1$. On cell $N_e$,\n",
    "only ${\\psi}_N={\\varphi}_{N_n-2}$ is involved, corresponding to\n",
    "integration with ${\\tilde{\\varphi}}_0$.\n",
    "We then get the special end-cell contributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:ex1:Ab:elm:ends\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\tilde A^{(e)} =\\frac{1}{h}\\left(\\begin{array}{r}\n",
    "1\n",
    "\\end{array}\\right),\\quad\n",
    "\\tilde b^{(e)} = h\\left(\\begin{array}{c}\n",
    "1\n",
    "\\end{array}\\right),\n",
    "\\label{fem:deq:1D:ex1:Ab:elm:ends} \\tag{72}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for $e=0$ and $e=N_e$. In these cells, we have only one degree of\n",
    "freedom, not two as in the interior cells.\n",
    "\n",
    "### Assembly\n",
    "\n",
    "The next step is to assemble the contributions from the various cells.\n",
    "The assembly of an element matrix and vector into the global matrix\n",
    "and right-hand side can be expressed as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A_{q(e,r),q(e,s)} = A_{q(e,r),q(e,s)} + \\tilde A^{(e)}_{r,s},\\quad\n",
    "b_{q(e,r)} = b_{q(e,r)} + \\tilde b^{(e)}_{r},\\quad\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for $r$ and $s$ running over all local degrees of freedom in cell $e$.\n",
    "\n",
    "To make the assembly algorithm more precise, it is convenient to set up\n",
    "Python data structures and a code snippet for carrying out all details\n",
    "of the algorithm.\n",
    "For a mesh of four equal-sized P1 elements and $L=2$ we have"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vertices = [0, 0.5, 1, 1.5, 2]\n",
    "cells = [[0, 1], [1, 2], [2, 3], [3, 4]]\n",
    "dof_map = [[0], [0, 1], [1, 2], [2]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The total number of degrees of freedom is 3, being the function\n",
    "values at the internal 3 nodes where $u$ is unknown.\n",
    "In cell 0 we have global degree of freedom 0, the next\n",
    "cell has $u$ unknown at its two nodes, which become\n",
    "global degrees of freedom 0 and 1, and so forth according to\n",
    "the `dof_map` list. The mathematical $q(e,r)$ quantity is nothing\n",
    "but the `dof_map` list.\n",
    "\n",
    "Assume all element matrices are stored in a list `Ae` such that\n",
    "`Ae[e][i,j]` is $\\tilde A_{i,j}^{(e)}$. A corresponding list\n",
    "for the element vectors is named `be`, where `be[e][r]` is\n",
    "$\\tilde b_r^{(e)}$.\n",
    "A Python code snippet\n",
    "illustrates all details of the assembly algorithm:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "# A[i,j]: coefficient matrix, b[i]: right-hand side\n",
    "for e in range(len(Ae)):\n",
    "    for r in range(Ae[e].shape[0]):\n",
    "        for s in range(Ae[e].shape[1]):\n",
    "            A[dof_map[e,r],dof_map[e,s]] += Ae[e][i,j]\n",
    "        b[dof_map[e,r]] += be[e][i,j]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The general case with `N_e` P1 elements of length `h` has"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "N_n = N_e + 1\n",
    "vertices = [i*h for i in range(N_n)]\n",
    "cells = [[e, e+1] for e in range(N_e)]\n",
    "dof_map = [[0]] + [[e-1, e] for i in range(1, N_e)] + [[N_n-2]]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Carrying out the assembly results in a linear system that is identical\n",
    "to ([66](#fem:deq:1D:ex1:Ab:glob)), which is not surprising, since\n",
    "the procedures is mathematically equivalent to the calculations\n",
    "in the physical domain.\n",
    "\n",
    "So far, our technique for computing the matrix system have assumed\n",
    "that $u(0)=u(L)=0$. The next section deals with the extension to\n",
    "nonzero Dirichlet conditions.\n",
    "\n",
    "# Boundary conditions: specified nonzero value\n",
    "<div id=\"fem:deq:1D:essBC\"></div>\n",
    "\n",
    "We have to take special actions to incorporate nonzero\n",
    "Dirichlet conditions,\n",
    "such as $u(L)=D$, into the computational procedures. The present\n",
    "section outlines alternative, yet mathematically equivalent, methods.\n",
    "\n",
    "\n",
    "## General construction of a boundary function\n",
    "<div id=\"fem:deq:1D:fem:essBC:Bfunc\"></div>\n",
    "\n",
    "In previous sections, we introduced a boundary function $B(x)$\n",
    "to deal with nonzero Dirichlet boundary conditions for $u$. The\n",
    "construction of such a function is not always trivial, especially not\n",
    "in multiple dimensions. However, a simple and general constructive idea\n",
    "exists when the\n",
    "basis functions have the property"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "{\\varphi}_i(x_{j}) = \\delta_{ij},\\quad\n",
    "\\delta_{ij} = \\left\\lbrace\\begin{array}{ll}\n",
    "1, & i=j,\\\\ \n",
    "0, & i\\neq j,\n",
    "\\end{array}\\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $x_{j}$ is a boundary point. Examples on such\n",
    "functions are the Lagrange interpolating polynomials and finite\n",
    "element functions.\n",
    "\n",
    "Suppose now that $u$ has Dirichlet boundary conditions at nodes\n",
    "with numbers $i\\in{I_b}$. For example, ${I_b} = \\{0,N_n-1\\}$ in a 1D\n",
    "mesh with node numbering from left to right and Dirichlet conditions\n",
    "at the end nodes $i=0$ and $i=N_n-1$.\n",
    "Let $U_i$ be the corresponding prescribed values of $u(x_{i})$.\n",
    "We can then, in general, use"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto37\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "B(x) = \\sum_{j\\in{I_b}} U_j{\\varphi}_j(x){\\thinspace .}\n",
    "\\label{_auto37} \\tag{73}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is easy to verify that\n",
    "$B(x_{i})= \\sum_{j\\in{I_b}} U_j{\\varphi}_j(x_{i})=U_i$.\n",
    "\n",
    "\n",
    "The unknown function can then be written as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto38\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u(x) = \\sum_{j\\in{I_b}} U_j{\\varphi}_j(x) + \\sum_{j\\in{\\mathcal{I}_s}}c_j{\\varphi}_{\\nu(j)},\n",
    "\\label{_auto38} \\tag{74}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $\\nu(j)$ maps unknown number $j$ in the equation system to\n",
    "node $\\nu(j)$, ${I_b}$ is the set of indices corresponding to basis functions\n",
    "associated with nodes where Dirichlet conditions apply, and ${\\mathcal{I}_s}$ is the\n",
    "set of indices used to number the unknowns from zero to $N$.\n",
    "We can easily show that with this $u$, a Dirichlet\n",
    "condition $u(x_{k})=U_k$ is fulfilled:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u(x_{k}) = \\sum_{j\\in{I_b}} U_j\\underbrace{{\\varphi}_j(x_k)}_{\\neq 0\\,\n",
    "\\Leftrightarrow\\,j=k} +\n",
    "\\sum_{j\\in{\\mathcal{I}_s}} c_j\\underbrace{{\\varphi}_{\\nu(j)}(x_{k})}_{=0,\\ k\\not\\in{\\mathcal{I}_s}}\n",
    "= U_k\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some examples will further clarify the notation. With a regular\n",
    "left-to-right numbering of nodes in a mesh with P1 elements,\n",
    "and Dirichlet conditions at $x=0$, we use finite element basis\n",
    "functions associated with the nodes $1, 2, \\ldots, N_n-1$, implying\n",
    "that  $\\nu(j)=j+1$, $j=0,\\ldots,N$, where $N=N_n-2$. Consider\n",
    "a particular mesh:\n",
    "\n",
    "<!-- dom:FIGURE: [fig/fe_mesh1D_P1.png, width=500 frac=0.8] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/fe_mesh1D_P1.png\" width=500>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "The expansion associated with this mesh becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u(x) = U_0{\\varphi}_0(x) + c_0{\\varphi}_1(x) +\n",
    "c_1{\\varphi}_2(x) + \\cdots + c_4{\\varphi}_5(x){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Switching to the more standard case of left-to-right numbering and\n",
    "boundary conditions $u(0)=C$, $u(L)=D$, we have $N=N_n-3$ and"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "u(x) &= C{\\varphi}_0 + D{\\varphi}_{N_n-1} + \\sum_{j\\in{\\mathcal{I}_s}} c_j{\\varphi}_{j+1}\\\\ \n",
    "&= C{\\varphi}_0 + D{\\varphi}_{N_n} + c_0{\\varphi}_1 + c_1{\\varphi}_2 +\\cdots\n",
    "+ c_N{\\varphi}_{N_n-2}{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finite element meshes in non-trivial 2D and 3D geometries usually leads\n",
    "to an irregular cell and node numbering. Let us therefore take a look\n",
    "at an irregular numbering in 1D:\n",
    "\n",
    "<!-- dom:FIGURE: [fig/fe_mesh1D_random_numbering.png, width=500 frac=0.8] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/fe_mesh1D_random_numbering.png\" width=500>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "Say we in this mesh have Dirichlet conditions on the left-most and\n",
    "right-most node, with numbers 3 and 1, respectively.  We can number\n",
    "the unknowns at the interior nodes as we want, e.g., from left to\n",
    "right, resulting in $\\nu(0)=0$, $\\nu(1)=4$, $\\nu(2)=5$, $\\nu(3)=2$.\n",
    "This gives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "B(x) = U_3{\\varphi}_3(x) + U_1{\\varphi}_1(x),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u(x) = B(x) + \\sum_{j=0}^3 c_j{\\varphi}_{\\nu(j)}\n",
    "= U_3{\\varphi}_3 + U_1{\\varphi}_1 + c_0{\\varphi}_0 + c_1{\\varphi}_4\n",
    "+ c_2{\\varphi}_5 + c_3{\\varphi}_2{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idea of constructing $B$, described here, generalizes almost\n",
    "trivially to 2D and 3D problems: $B=\\sum_{j\\in{I_b}}U_j{\\varphi}_j$,\n",
    "where ${I_b}$ is the index set containing the numbers of all the\n",
    "nodes on the boundaries where Dirichlet values are prescribed.\n",
    "\n",
    "## Example on computing with a finite element-based boundary function\n",
    "\n",
    "Let us see how the model problem $-u''=2$, $u(0)=C$, $u(L)=D$,\n",
    "is affected by a $B(x)$ to incorporate boundary values.\n",
    "Inserting the expression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u(x) = B(x) + \\sum_{j\\in{\\mathcal{I}_s}}c_j{\\psi}_j(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "in $-(u'',{\\psi}_i)=(f,{\\psi}_i)$ and\n",
    "integrating by parts results in a linear system with"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A_{i,j} = \\int_0^L {\\psi}_i'(x){\\psi}_j'(x) {\\, \\mathrm{d}x},\\quad\n",
    "b_i = \\int_0^L (f(x){\\psi}_i(x) - B'(x){\\psi}_i'(x)) {\\, \\mathrm{d}x}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We choose ${\\psi}_i={\\varphi}_{i+1}$, $i=0,\\ldots,N=N_n-3$\n",
    "if the node numbering is from left\n",
    "to right. (Later we also need the assumption that cells too\n",
    "are numbered from left to right.)\n",
    "The boundary function becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "B(x) = C{\\varphi}_0(x) + D{\\varphi}_{N_n-1}(x){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The expansion for $u(x)$ is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u(x)  = B(x) + \\sum_{j\\in{\\mathcal{I}_s}} c_j{\\varphi}_{j+1}(x){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write the matrix and right-hand side entries as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "A_{i-1,j-1} &= \\int_0^L {\\varphi}_i'(x){\\varphi}_j'(x) {\\, \\mathrm{d}x},\\\\ \n",
    "b_{i-1} &=\n",
    "\\int_0^L (f(x){\\varphi}_i'(x) - (C{\\varphi}_{0}'(x) + D{\\varphi}_{N_n-1}'(x)){\\varphi}_i'(x) ){\\, \\mathrm{d}x},\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for $i,j = 1,\\ldots,N+1=N_n-2$. Note that we have here used\n",
    "$B'=C{\\varphi}_0' + D{\\varphi}_{N_n-1}'$.\n",
    "\n",
    "### Computations in physical coordinates\n",
    "\n",
    "Most of the terms in the linear system have already been computed\n",
    "so we concentrate on the new contribution from the boundary function.\n",
    "The integral $C\\int_0^L {\\varphi}_{0}'(x)){\\varphi}_i'(x) {\\, \\mathrm{d}x}$, associated\n",
    "with the Dirichlet condition in $x=0$,  can only get\n",
    "a nonzero contribution from the first cell,\n",
    "$\\Omega^{(0)}=[x_{0},x_{1}]$\n",
    "since ${\\varphi}_{0}'(x)=0$ on all other cells. Moreover,\n",
    "${\\varphi}_{0}'(x){\\varphi}_i'(x) {\\, \\mathrm{d}x} \\neq 0$ only for $i=0$ and $i=1$\n",
    "(but node $i=0$ is excluded from the formulation),\n",
    "since ${\\varphi}_{i}=0$ on the first cell if $i>1$.\n",
    "With a similar reasoning we realize that\n",
    "$D\\int_0^L {\\varphi}_{N_n-1}'(x)){\\varphi}_i'(x) {\\, \\mathrm{d}x}$ can only get\n",
    "a nonzero contribution from the last cell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\int_0^L {\\varphi}_{0}'(x){\\varphi}_{1}'(x) {\\, \\mathrm{d}x} &=\n",
    "(-\\frac{1}{h})\\cdot\\frac{1}{h}\\cdot h = -\\frac{1}{h},\\\\ \n",
    "\\int_0^L {\\varphi}_{N_n-1}'(x){\\varphi}_{N_n-2}'(x) {\\, \\mathrm{d}x} &=\n",
    "\\frac{1}{h}\\cdot(-\\frac{1}{h})\\cdot h = -\\frac{1}{h}{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With these expressions we get"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "b_0 = \\int_0^Lf(x){\\varphi}_1{\\, \\mathrm{d}x} - C(-\\frac{1}{h}),\\quad\n",
    "b_N = \\int_0^L f(x){\\varphi}_{N_n-2}{\\, \\mathrm{d}x} - D(-\\frac{1}{h}){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cellwise computations on the reference element\n",
    "\n",
    "As an equivalent alternative, we now turn to cellwise computations.\n",
    "The element matrices and vectors are calculated as in the section [Cellwise computations](#fem:deq:1D:comp:elmwise), so we concentrate on the impact of\n",
    "the new term involving $B(x)$. This new term,\n",
    "$B'=C{\\varphi}_0' + D{\\varphi}_{N_n-1}'$, vanishes on all cells except\n",
    "for $e=0$ and $e=N_e$. Over the first cell ($e=0$) the $B'(x)$ function\n",
    "in local coordinates reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{dB}{dx} = C\\frac{2}{h}\\frac{d{\\tilde{\\varphi}}_0}{dX},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "while over the last cell ($e=N_e$) it looks like"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{dB}{dx} = D\\frac{2}{h}\\frac{d{\\tilde{\\varphi}}_1}{dX}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For an arbitrary interior cell, we have the formula"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\tilde b_r^{(e)} = \\int_{-1}^1 f(x(X)){\\tilde{\\varphi}}_r(X)\\frac{h}{2}{\\, \\mathrm{d}X},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for an entry in the local element vector.\n",
    "In the first cell, the value at local node 0 is known so only the value\n",
    "at local node 1 is unknown. The associated element vector entry becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\tilde b_0^{(1)} &= \\int_{-1}^1 \\left(f{\\tilde{\\varphi}}_1 -\n",
    "C\\frac{2}{h}\\frac{d{\\tilde{\\varphi}}_0}{dX}\\frac{2}{h}\\frac{d{\\tilde{\\varphi}}_1}{dX}\\right)\n",
    "\\frac{h}{2} {\\, \\mathrm{d}X} \\\\ \n",
    "& = \\frac{h}{2} 2\\int_{-1}^1 {\\tilde{\\varphi}}_1  {\\, \\mathrm{d}X}\n",
    "- C\\frac{2}{h}(-\\frac{1}{2})\\frac{2}{h}\\frac{1}{2}\\frac{h}{2}\\cdot 2\n",
    "= h + C\\frac{1}{h}{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The value at local node 1 in the last cell is known so the\n",
    "element vector here is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\tilde b_0^{N_e} &= \\int_{-1}^1 \\left(f{\\tilde{\\varphi}}_0 -\n",
    "D\\frac{2}{h}\\frac{d{\\tilde{\\varphi}}_1}{dX}\\frac{2}{h}\\frac{d{\\tilde{\\varphi}}_0}{dX}\\right)\n",
    "\\frac{h}{2} {\\, \\mathrm{d}X}\\\\ \n",
    "& = \\frac{h}{2} 2\\int_{-1}^1 {\\tilde{\\varphi}}_0  {\\, \\mathrm{d}X}\n",
    "- D\\frac{2}{h}\\frac{1}{2}\\frac{2}{h}(-\\frac{1}{2})\\frac{h}{2}\\cdot 2\n",
    "= h + D\\frac{1}{h}{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The contributions from the $B(x)$ function to the global right-hand side\n",
    "vector becomes $C/h$ for $b_0$ and $D/h$ for $b_N$, exactly as we\n",
    "computed in the physical domain.\n",
    "\n",
    "## Modification of the linear system\n",
    "<div id=\"fem:deq:1D:fem:essBC:Bfunc:modsys\"></div>\n",
    "\n",
    "From an implementational point of view, there is a convenient alternative\n",
    "to adding the $B(x)$ function and using only the basis functions associated\n",
    "with nodes where $u$ is truly unknown.\n",
    "Instead of seeking"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:fem:essBC:Bfunc:modsys:utrad\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u(x) = \\sum_{j\\in{I_b}} U_j{\\varphi}_j(x)\n",
    "+ \\sum_{j\\in{\\mathcal{I}_s}}c_j{\\varphi}_{\\nu(j)}(x),\n",
    "\\label{fem:deq:1D:fem:essBC:Bfunc:modsys:utrad} \\tag{75}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we use the sum over all degrees of freedom, including the known boundary\n",
    "values:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:fem:essBC:Bfunc:modsys:uall\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u(x) = \\sum_{j\\in{\\mathcal{I}_s}}c_j{\\varphi}_j(x){\\thinspace .}\n",
    "\\label{fem:deq:1D:fem:essBC:Bfunc:modsys:uall} \\tag{76}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the collections of unknowns\n",
    "$\\left\\{ {c}_i \\right\\}_{i\\in{\\mathcal{I}_s}}$ in ([75](#fem:deq:1D:fem:essBC:Bfunc:modsys:utrad))\n",
    "and ([76](#fem:deq:1D:fem:essBC:Bfunc:modsys:uall)) are different.\n",
    "The index set ${\\mathcal{I}_s}=\\{0,\\ldots,N\\}$ always goes to $N$, and the\n",
    "number of unknowns is $N+1$, but\n",
    "in ([75](#fem:deq:1D:fem:essBC:Bfunc:modsys:utrad)) the unknowns\n",
    "correspond to nodes where $u$ is not known, while\n",
    "in ([76](#fem:deq:1D:fem:essBC:Bfunc:modsys:uall)) the unknowns\n",
    "cover $u$ values at all the nodes. So, if the index set\n",
    "${I_b}$ contains $N_b$ node numbers where $u$ is prescribed,\n",
    "we have that $N=N_n-N_b$ in\n",
    "([75](#fem:deq:1D:fem:essBC:Bfunc:modsys:utrad)) and\n",
    "$N=N_n$ in ([76](#fem:deq:1D:fem:essBC:Bfunc:modsys:uall)).\n",
    "\n",
    "The idea is to compute the entries in the linear system as if no\n",
    "Dirichlet values are prescribed. Afterwards, we modify the linear system\n",
    "to ensure that the known $c_j$ values are incorporated.\n",
    "\n",
    "A potential problem arises for the boundary term $[u'v]_0^L$ from the\n",
    "integration by parts: imagining no Dirichlet conditions means that we\n",
    "no longer require $v=0$ at Dirichlet points, and the boundary term is\n",
    "then nonzero at these points. However, when we modify the linear\n",
    "system, we will erase whatever the contribution from $[u'v]_0^L$\n",
    "should be at the Dirichlet points in the right-hand side of the linear\n",
    "system. We can therefore safely forget $[u'v]_0^L$ at any point where\n",
    "a Dirichlet condition applies.\n",
    "\n",
    "\n",
    "\n",
    "### Computations in the physical system\n",
    "\n",
    "Let us redo the computations in the example in\n",
    "the section [General construction of a boundary function](#fem:deq:1D:fem:essBC:Bfunc). We solve $-u''=2$ with\n",
    "$u(0)=0$ and $u(L)=D$. The expressions for $A_{i,j}$ and $b_i$\n",
    "are the same, but the numbering is different as the numbering of\n",
    "unknowns and nodes now coincide:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A_{i,j} = \\int_0^L {\\varphi}_i'(x){\\varphi}_j'(x) {\\, \\mathrm{d}x},\\quad\n",
    "b_{i} = \\int_0^L f(x){\\varphi}_i(x) {\\, \\mathrm{d}x},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for $i,j = 0,\\ldots,N=N_n-1$.\n",
    "The integrals involving basis functions\n",
    "corresponding to interior mesh nodes, $i,j=1,\\ldots,N_n-2$, are\n",
    "obviously the same as before. We concentrate on the contributions\n",
    "from ${\\varphi}_0$ and ${\\varphi}_{N_n-1}$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "A_{0,0} &= \\int_0^L ({\\varphi}_0')^2{\\, \\mathrm{d}x} = \\int_{0}^{x_{1}}\n",
    "= ({\\varphi}_0')^2{\\, \\mathrm{d}x} \\frac{1}{h},\\\\ \n",
    "A_{0,1} &= \\int_0^L {\\varphi}_0'{\\varphi}_1'{\\, \\mathrm{d}x}\n",
    "= \\int_{0}^{x_{1}} {\\varphi}_0'{\\varphi}_1'{\\, \\mathrm{d}x} = -\\frac{1}{h},\\\\ \n",
    "A_{N,N} &= \\int_0^L ({\\varphi}_N')^2{\\, \\mathrm{d}x}\n",
    "= \\int_{x_{N_n-2}}^{x_{N_n-1}} ({\\varphi}_N')^2{\\, \\mathrm{d}x} = \\frac{1}{h},\\\\ \n",
    "A_{N,N-1} &= \\int_0^L {\\varphi}_N'{\\varphi}_{N-1}'{\\, \\mathrm{d}x}\n",
    "=\\int_{x_{N_n-2}}^{x_{N_n-1}} {\\varphi}_N'{\\varphi}_{N-1}'{\\, \\mathrm{d}x} = -\\frac{1}{h}{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The new terms on the right-hand side are also those involving\n",
    "${\\varphi}_0$ and ${\\varphi}_{N_n-1}$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "b_0 &= \\int_0^L 2{\\varphi}_0(x) {\\, \\mathrm{d}x} = \\int_0^{x_{1}} 2{\\varphi}_0(x){\\, \\mathrm{d}x} = h,\\\\ \n",
    "b_N &=  \\int_0^L 2{\\varphi}_{N_n-1}{\\, \\mathrm{d}x} =\n",
    "\\int_{x_{N_n-2}}^{x_{N_n-1}} 2{\\varphi}_{N_n-1}{\\, \\mathrm{d}x} = h{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The complete matrix system, involving all degrees of freedom, takes the form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:ex1:Ab:glob2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{1}{h}\\left(\n",
    "\\begin{array}{ccccccccc}\n",
    "1 & -1 & 0\n",
    "&\\cdots &\n",
    "\\cdots & \\cdots & \\cdots &\n",
    "\\cdots & 0 \\\\ \n",
    "-1 & 2 & -1 & \\ddots &   & &  & &  \\vdots \\\\ \n",
    "0 & -1 & 2 & -1 &\n",
    "\\ddots & &  &  & \\vdots \\\\ \n",
    "\\vdots & \\ddots &  & \\ddots & \\ddots & 0 &  & & \\vdots \\\\ \n",
    "\\vdots &  & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & & \\vdots \\\\ \n",
    "\\vdots & &  & 0 & -1 & 2 & -1 & \\ddots & \\vdots \\\\ \n",
    "\\vdots & & &  & \\ddots & \\ddots & \\ddots &\\ddots  & 0 \\\\ \n",
    "\\vdots & & & &  &\\ddots  & \\ddots &\\ddots  & -1 \\\\ \n",
    "0 &\\cdots & \\cdots &\\cdots & \\cdots & \\cdots  & 0 & -1 & 1\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "c_0 \\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots\\\\ \n",
    "c_{N}\n",
    "\\end{array}\n",
    "\\right)\n",
    "=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "h \\\\ \n",
    "2h\\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "2h\\\\ \n",
    "h\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\label{fem:deq:1D:ex1:Ab:glob2} \\tag{77}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Incorporation of Dirichlet values can now be done by replacing\n",
    "the first and last equation by\n",
    "the very simple equations $c_0=0$ and $c_N=D$, respectively.\n",
    "Note that the factor $1/h$ in front of the matrix then requires\n",
    "a factor $h$ to be introduce appropriately on the diagonal in the first and last\n",
    "row of the matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:ex1:Ab:glob3\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{1}{h}\\left(\n",
    "\\begin{array}{ccccccccc}\n",
    "h & 0 & 0\n",
    "&\\cdots &\n",
    "\\cdots & \\cdots & \\cdots &\n",
    "\\cdots & 0 \\\\ \n",
    "-1 & 2 & -1 & \\ddots &   & &  & &  \\vdots \\\\ \n",
    "0 & -1 & 2 & -1 &\n",
    "\\ddots & &  &  & \\vdots \\\\ \n",
    "\\vdots & \\ddots &  & \\ddots & \\ddots & 0 &  & & \\vdots \\\\ \n",
    "\\vdots &  & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & & \\vdots \\\\ \n",
    "\\vdots & &  & 0 & -1 & 2 & -1 & \\ddots & \\vdots \\\\ \n",
    "\\vdots & & &  & \\ddots & \\ddots & \\ddots &\\ddots  & 0 \\\\ \n",
    "\\vdots & & & &  &\\ddots  & \\ddots &\\ddots  & -1 \\\\ \n",
    "0 &\\cdots & \\cdots &\\cdots & \\cdots & \\cdots  & 0 & 0 & h\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "c_0 \\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots\\\\ \n",
    "c_{N}\n",
    "\\end{array}\n",
    "\\right)\n",
    "=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "0 \\\\ \n",
    "2h\\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "2h\\\\ \n",
    "D\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\label{fem:deq:1D:ex1:Ab:glob3} \\tag{78}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that because we do not require ${\\varphi}_i(0)=0$ and\n",
    "${\\varphi}_i(L)=0$, $i\\in{\\mathcal{I}_s}$, the boundary term $[u'v]_0^L$,\n",
    "in principle, gives\n",
    "contributions $u'(0){\\varphi}_0(0)$ to $b_0$ and\n",
    "$u'(L){\\varphi}_N(L)$ to $b_N$ ($u'{\\varphi}_i$ vanishes for $x=0$ or\n",
    "$x=L$ for $i=1,\\ldots,N-1$).  Nevertheless, we erase these\n",
    "contributions in $b_0$ and $b_N$ and insert boundary values instead. This\n",
    "argument shows why we can drop computing $[u'v]_0^L$ at Dirichlet\n",
    "nodes when we implement the Dirichlet values by modifying the linear\n",
    "system.\n",
    "\n",
    "\n",
    "## Symmetric modification of the linear system\n",
    "<div id=\"fem:deq:1D:fem:essBC:Bfunc:modsys:symm\"></div>\n",
    "\n",
    "The original matrix system ([66](#fem:deq:1D:ex1:Ab:glob)) is symmetric,\n",
    "but the modifications in ([78](#fem:deq:1D:ex1:Ab:glob3)) destroy this\n",
    "symmetry. Our described modification will in general destroy an\n",
    "initial symmetry in the matrix system. This is not a particular\n",
    "computational disadvantage for tridiagonal systems arising in 1D\n",
    "problems, but may be more serious in 2D and 3D problems when the\n",
    "systems are large and exploiting symmetry can be important for halving\n",
    "the storage demands and speeding up computations. Methods for solving\n",
    "symmetric matrix are also usually more stable and efficient than those\n",
    "for non-symmetric systems.  Therefore, an alternative modification\n",
    "which preserves symmetry is attractive.\n",
    "\n",
    "One can formulate a general algorithm for incorporating a Dirichlet\n",
    "condition in a symmetric way.\n",
    "Let $c_k$ be a coefficient corresponding to a known value\n",
    "$u(x_{k}) = U_k$.\n",
    "We want to replace equation $k$ in the system by $c_k=U_k$, i.e.,\n",
    "insert zeroes in row number $k$ in the coefficient matrix,\n",
    "set 1 on the diagonal, and replace $b_k$ by $U_k$.\n",
    "A symmetry-preserving modification consists in first\n",
    "subtracting column number $k$ in the coefficient matrix, i.e., $A_{i,k}$\n",
    "for $i\\in{\\mathcal{I}_s}$, times the boundary value $U_k$, from the\n",
    "right-hand side: $b_i \\leftarrow b_i - A_{i,k}U_k$,\n",
    "$i=0,\\ldots,N$. Then we put\n",
    "zeroes in both row number $k$ *and* column number $k$ in the coefficient matrix,\n",
    "and finally set $b_k=U_k$. The steps in algorithmic form becomes\n",
    "\n",
    "1. $b_i \\leftarrow b_i - A_{i,k}U_k$ for $i\\in{\\mathcal{I}_s}$\n",
    "\n",
    "2. $A_{i,k} = A_{k,i} = 0$ for $i\\in{\\mathcal{I}_s}$\n",
    "\n",
    "3. $A_{k,k}=1$\n",
    "\n",
    "4. $b_i = U_k$\n",
    "\n",
    "This modification goes as follows for the specific linear system\n",
    "written out in ([77](#fem:deq:1D:ex1:Ab:glob2)) in\n",
    "the section [Modification of the linear system](#fem:deq:1D:fem:essBC:Bfunc:modsys). First we\n",
    "subtract the first column in the coefficient matrix, times the boundary\n",
    "value, from the right-hand side. Because $c_0=0$, this subtraction\n",
    "has no effect. Then we subtract the last column, times the boundary value $D$,\n",
    "from the right-hand side. This action results in $b_{N-1}=2h+D/h$ and\n",
    "$b_N=h-2D/h$. Thereafter, we place zeros in the first and last row and\n",
    "column in the coefficient matrix and 1 on the two corresponding diagonal\n",
    "entries. Finally, we set $b_0=0$ and $b_N=D$. The result becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:ex1:Ab:glob3:symm\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{1}{h}\\left(\n",
    "\\begin{array}{ccccccccc}\n",
    "h & 0 & 0\n",
    "&\\cdots &\n",
    "\\cdots & \\cdots & \\cdots &\n",
    "\\cdots & 0 \\\\ \n",
    "0 & 2 & -1 & \\ddots &   & &  & &  \\vdots \\\\ \n",
    "0 & -1 & 2 & -1 &\n",
    "\\ddots & &  &  & \\vdots \\\\ \n",
    "\\vdots & \\ddots &  & \\ddots & \\ddots & 0 &  & & \\vdots \\\\ \n",
    "\\vdots &  & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & & \\vdots \\\\ \n",
    "\\vdots & &  & 0 & -1 & 2 & -1 & \\ddots & \\vdots \\\\ \n",
    "\\vdots & & &  & \\ddots & \\ddots & \\ddots &\\ddots  & 0 \\\\ \n",
    "\\vdots & & & &  &\\ddots  & \\ddots &\\ddots  & 0 \\\\ \n",
    "0 &\\cdots & \\cdots &\\cdots & \\cdots & \\cdots  & 0 & 0 & h\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "c_0 \\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots\\\\ \n",
    "c_{N}\n",
    "\\end{array}\n",
    "\\right)\n",
    "=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "0 \\\\ \n",
    "2h\\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "2h +D/h\\\\ \n",
    "D\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\label{fem:deq:1D:ex1:Ab:glob3:symm} \\tag{79}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modification of the element matrix and vector\n",
    "<div id=\"fem:bc:elmat:mod\"></div>\n",
    "\n",
    "The modifications of the global linear system can alternatively\n",
    "be done for the element matrix and vector. Let us perform the\n",
    "associated calculations in the computational example where\n",
    "the element matrix and vector is given by\n",
    "([71](#fem:deq:1D:ex1:Ab:elm)). The modifications are needed in\n",
    "cells where one of the degrees of freedom is known. In the\n",
    "present example, this means\n",
    "the first and last cell. We compute the element matrix\n",
    "and vector as if there were no Dirichlet conditions. The boundary term\n",
    "$[u'v]_0^L$ is simply forgotten at nodes that have Dirichlet conditions\n",
    "because the modification of the element vector will anyway erase the\n",
    "contribution from the boundary term. In the first cell,\n",
    "local degree of freedom number 0\n",
    "is known and the modification becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:ex1:Ab:elm:bc:0\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\tilde A^{(0)} =\n",
    "A = \\frac{1}{h}\\left(\\begin{array}{rr}\n",
    "h & 0\\\\ \n",
    "-1 & 1\n",
    "\\end{array}\\right),\\quad\n",
    "\\tilde b^{(0)} = \\left(\\begin{array}{c}\n",
    "0\\\\ \n",
    "h\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "\\label{fem:deq:1D:ex1:Ab:elm:bc:0} \\tag{80}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the last cell we set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:ex1:Ab:elm:bc:N\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\tilde A^{(N_e)} =\n",
    "A = \\frac{1}{h}\\left(\\begin{array}{rr}\n",
    "1 & -1\\\\ \n",
    "0 & h\n",
    "\\end{array}\\right),\\quad\n",
    "\\tilde b^{(N_e)} = \\left(\\begin{array}{c}\n",
    "h\\\\ \n",
    "D\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "\\label{fem:deq:1D:ex1:Ab:elm:bc:N} \\tag{81}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also perform the symmetric modification. This operation affects\n",
    "only the last cell with a nonzero Dirichlet condition. The algorithm\n",
    "is the same as for the global linear system, resulting in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:ex1:Ab:elm:bc:N:symm\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\tilde A^{(N_e)} =\n",
    "A = \\frac{1}{h}\\left(\\begin{array}{rr}\n",
    "1 & 0\\\\ \n",
    "0 & h\n",
    "\\end{array}\\right),\\quad\n",
    "\\tilde b^{(N_e)} = \\left(\\begin{array}{c}\n",
    "h + D/h\\\\ \n",
    "D\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "\\label{fem:deq:1D:ex1:Ab:elm:bc:N:symm} \\tag{82}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reader is encouraged to assemble the element matrices and vectors and\n",
    "check that the result coincides with the system\n",
    "([79](#fem:deq:1D:ex1:Ab:glob3:symm)).\n",
    "\n",
    "<!-- As a final remark, we repeat that Dirichlet conditions are referred -->\n",
    "<!-- to as essential boundary conditions because they require -->\n",
    "<!-- quite some work with modifying -->\n",
    "<!-- either the linear system or the expansion formula -->\n",
    "<!-- for $u$. Boundary conditions for the derivative are much easier to -->\n",
    "<!-- implement, as shown next, and therefore deserve the name *natural -->\n",
    "<!-- boundary condition*. -->\n",
    "\n",
    "# Boundary conditions: specified derivative\n",
    "<div id=\"fem:deq:1D:BC:nat\"></div>\n",
    "\n",
    "Suppose our model problem $-u''(x)=f(x)$ features\n",
    "the boundary conditions $u'(0)=C$ and $u(L)=D$.\n",
    "As already indicated,\n",
    "the former condition can be incorporated through the boundary term\n",
    "that arises from integration by parts. The details of this method will now be\n",
    "illustrated in the context of finite element basis functions.\n",
    "\n",
    "## The variational formulation\n",
    "\n",
    "Starting with the Galerkin method,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\int_0^L(u''(x)+f(x)){\\psi}_i(x) {\\, \\mathrm{d}x} = 0,\\quad i\\in{\\mathcal{I}_s},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "integrating $u''{\\psi}_i$ by parts results in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\int_0^Lu'(x)'{\\psi}_i'(x) {\\, \\mathrm{d}x} -(u'(L){\\psi}_i(L) - u'(0){\\psi}_i(0)) =\n",
    "\\int_0^L f(x){\\psi}_i(x) {\\, \\mathrm{d}x}, \\quad i\\in{\\mathcal{I}_s}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first boundary term, $u'(L){\\psi}_i(L)$,\n",
    "vanishes because $u(L)=D$.\n",
    "The second boundary\n",
    "term, $u'(0){\\psi}_i(0)$, can be used to implement the condition $u'(0)=C$,\n",
    "provided ${\\psi}_i(0)\\neq 0$ for some $i$ (but with finite elements\n",
    "we fortunately have ${\\psi}_0(0)=1$).\n",
    "The variational form of the differential equation then becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:BC:nat:varform\"></div>\n",
    "\n",
    "$$\n",
    "\\int_0^Lu'(x){\\varphi}_i'(x) {\\, \\mathrm{d}x} + C{\\varphi}_i(0) =\n",
    "\\int_0^L f(x){\\varphi}_i(x) {\\, \\mathrm{d}x},\\quad i\\in{\\mathcal{I}_s}{\\thinspace .}\n",
    "\\label{fem:deq:1D:BC:nat:varform} \\tag{83}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boundary term vanishes because of the test functions\n",
    "<div id=\"fem:deq:1D:BC:nat:uLtest\"></div>\n",
    "\n",
    "At points where $u$ is known we may require ${\\psi}_i$ to vanish.\n",
    "Here, $u(L)=D$ and then ${\\psi}_i(L)=0$, $i\\in{\\mathcal{I}_s}$. Obviously,\n",
    "the boundary term $u'(L){\\psi}_i(L)$ then vanishes.\n",
    "\n",
    "The set of basis functions $\\left\\{ {{\\psi}}_i \\right\\}_{i\\in{\\mathcal{I}_s}}$ contains, in this\n",
    "case, all the finite element basis functions on the mesh, except\n",
    "the one that is 1 at $x=L$. The basis function that is left out is\n",
    "used in a boundary function $B(x)$ instead.\n",
    "With a left-to-right numbering,\n",
    "${\\psi}_i = {\\varphi}_i$, $i=0,\\ldots,N_n-2$, and $B(x)=D{\\varphi}_{N_n-1}$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u(x) = D{\\varphi}_{N_n-1}(x) + \\sum_{j=0}^{N=N_n-2} c_j{\\varphi}_j(x){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inserting this expansion for $u$ in the variational form\n",
    "([83](#fem:deq:1D:BC:nat:varform)) leads to the linear system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:natBC\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\sum_{j=0}^{N}\\left(\n",
    "\\int_0^L {\\varphi}_i'(x){\\varphi}_j'(x) {\\, \\mathrm{d}x} \\right)c_j =\n",
    "\\int_0^L\\left(f(x){\\varphi}_i(x) -D{\\varphi}_{N_n-1}'(x){\\varphi}_i'(x)\\right) {\\, \\mathrm{d}x}\n",
    " - C{\\varphi}_i(0),\n",
    "\\label{fem:deq:1D:natBC} \\tag{84}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for $i=0,\\ldots,N=N_n-2$.\n",
    "\n",
    "\n",
    "## Boundary term vanishes because of linear system modifications\n",
    "<div id=\"fem:deq:1D:BC:nat:uLmod\"></div>\n",
    "\n",
    "We may, as an alternative to the approach in the previous section, use\n",
    "a basis $\\left\\{ {{\\psi}}_i \\right\\}_{i\\in{\\mathcal{I}_s}}$ which contains all the finite element\n",
    "functions on the mesh: ${\\psi}_i={\\varphi}_i$, $i=0,\\ldots,N_n-1=N$.  In\n",
    "this case, $u'(L){\\psi}_i(L)=u'(L){\\varphi}_i(L)\\neq 0$ for the $i$\n",
    "corresponding to the boundary node at $x=L$ (where ${\\varphi}_i=1$).\n",
    "The number of this node is $i=N_n-1=N$ if a left-to-right numbering of\n",
    "nodes is utilized.\n",
    "\n",
    "However, even though $u'(L){\\varphi}_{N_n-1}(L)\\neq 0$, we do not need to\n",
    "compute this term.  For $i<N_n-1$ we realize that ${\\varphi}_i(L)=0$.  The\n",
    "only nonzero contribution to the right-hand side comes from $i=N$\n",
    "($b_N$). Without a boundary function we must implement the\n",
    "condition $u(L)=D$ by the equivalent statement $c_N=D$ and modify the\n",
    "linear system accordingly. This modification will erase the last\n",
    "row and replace $b_N$ by another value. Any attempt to compute\n",
    "the boundary term $u'(L){\\varphi}_{N_n-1}(L)$ and store it in $b_N$ will be\n",
    "lost. Therefore, we can safely forget about boundary terms\n",
    "corresponding to Dirichlet boundary conditions also when we use\n",
    "the methods from the section [Modification of the linear system](#fem:deq:1D:fem:essBC:Bfunc:modsys)\n",
    "or the section [Symmetric modification of the linear system](#fem:deq:1D:fem:essBC:Bfunc:modsys:symm).\n",
    "\n",
    "The expansion for $u$ reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u(x) = \\sum_{j\\in{\\mathcal{I}_s}} c_j{\\varphi}_j(x){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Insertion in the variational form\n",
    "([83](#fem:deq:1D:BC:nat:varform)) leads to\n",
    "the linear system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:natBC2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\sum_{j\\in{\\mathcal{I}_s}}\\left(\n",
    "\\int_0^L {\\varphi}_i'(x){\\varphi}_j'(x) {\\, \\mathrm{d}x} \\right)c_j =\n",
    "\\int_0^L\\left(f(x){\\varphi}_i(x)\\right) {\\, \\mathrm{d}x}\n",
    " - C{\\varphi}_i(0),\\quad i\\in{\\mathcal{I}_s}\n",
    "{\\thinspace .}\n",
    "\\label{fem:deq:1D:natBC2} \\tag{85}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After having computed the system, we replace the last row by\n",
    "$c_N=D$, either straightforwardly as in\n",
    "the section [Modification of the linear system](#fem:deq:1D:fem:essBC:Bfunc:modsys) or in a symmetric\n",
    "fashion as in the section [Symmetric modification of the linear system](#fem:deq:1D:fem:essBC:Bfunc:modsys:symm).\n",
    "These modifications can also be performed in the element matrix and\n",
    "vector for the right-most cell.\n",
    "\n",
    "\n",
    "## Direct computation of the global linear system\n",
    "<div id=\"fem:deq:1D:BC:nat:Aub\"></div>\n",
    "\n",
    "We now turn to actual computations with P1 finite elements.\n",
    "The focus is on how the linear system and\n",
    "the element matrices and vectors are modified by the\n",
    "condition $u'(0)=C$.\n",
    "\n",
    "Consider first the approach where Dirichlet conditions are incorporated\n",
    "by a $B(x)$ function and the known degree of freedom\n",
    "$C_{N_n-1}$ is left out of the linear system\n",
    "(see the section [Boundary term vanishes because of the test functions](#fem:deq:1D:BC:nat:uLtest)).\n",
    "The relevant formula for the linear\n",
    "system is given by ([84](#fem:deq:1D:natBC)).\n",
    "There are three differences compared to the extensively\n",
    "computed case where $u(0)=0$ in the sections [Computation in the global physical domain](#fem:deq:1D:comp:global) and [Cellwise computations](#fem:deq:1D:comp:elmwise).\n",
    "First, because we do not have a Dirichlet\n",
    "condition at the left boundary, we need to extend the linear system\n",
    "([66](#fem:deq:1D:ex1:Ab:glob)) with an equation associated with the node\n",
    "$x_{0}=0$.\n",
    "According to the section [Modification of the linear system](#fem:deq:1D:fem:essBC:Bfunc:modsys), this\n",
    "extension consists of including $A_{0,0}=1/h$, $A_{0,1}=-1/h$, and $b_0=h$.\n",
    "For $i>0$ we have $A_{i,i}=2/h$, $A_{i-1,i}=A_{i,i+1}=-1/h$.\n",
    "Second, we need to include\n",
    "the extra term\n",
    "$-C{\\varphi}_i(0)$ on the right-hand side. Since all ${\\varphi}_i(0)=0$\n",
    "for $i=1,\\ldots,N$, this term reduces to $-C{\\varphi}_0(0)=-C$ and\n",
    "affects only the first equation ($i=0$). We simply add $-C$ to $b_0$\n",
    "such that $b_0=h - C$.\n",
    "Third, the boundary term $-\\int_0^L D{\\varphi}_{N_n-1}(x){\\varphi}_i{\\, \\mathrm{d}x}$\n",
    "must be computed. Since $i=0,\\ldots,N=N_n-2$, this integral can only\n",
    "get a nonzero contribution with $i=N_n-2$ over the last cell.\n",
    "The result becomes $-Dh/6$.\n",
    "The resulting linear system can be summarized in the form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:BC:nat:Aub:system\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{1}{h}\\left(\n",
    "\\begin{array}{ccccccccc}\n",
    "1 & -1 & 0 &\\cdots & \\cdots & \\cdots & \\cdots & \\cdots & 0 \\\\ \n",
    "-1 & 2 & -1 & \\ddots &   & &  & &  \\vdots \\\\ \n",
    "0 & -1 & 2 & -1 &\n",
    "\\ddots & &  &  & \\vdots \\\\ \n",
    "\\vdots & \\ddots &  & \\ddots & \\ddots & 0 &  & & \\vdots \\\\ \n",
    "\\vdots &  & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & & \\vdots \\\\ \n",
    "\\vdots & &  & 0 & -1 & 2 & -1 & \\ddots & \\vdots \\\\ \n",
    "\\vdots & & &  & \\ddots & \\ddots & \\ddots &\\ddots  & 0 \\\\ \n",
    "\\vdots & & & &  &\\ddots  & \\ddots &\\ddots  & -1 \\\\ \n",
    "0 &\\cdots & \\cdots &\\cdots & \\cdots & \\cdots  & 0 & -1 & 2\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "c_0 \\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots\\\\ \n",
    "c_{N}\n",
    "\\end{array}\n",
    "\\right)\n",
    "=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "h - C \\\\ \n",
    "2h\\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots\\\\ \n",
    "2h - Dh/6\n",
    "\\end{array}\n",
    "\\right){\\thinspace .}\n",
    "\\label{fem:deq:1D:BC:nat:Aub:system} \\tag{86}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we consider the technique where we modify the linear system to\n",
    "incorporate Dirichlet conditions\n",
    "(see the section [Boundary term vanishes because of linear system modifications](#fem:deq:1D:BC:nat:uLmod)). Now $N=N_n-1$.\n",
    "The two differences from the\n",
    "case above is that the $-\\int_0^LD{\\varphi}_{N_n-1}{\\varphi}_i{\\, \\mathrm{d}x}$ term is\n",
    "left out of the right-hand side and an extra last row associated\n",
    "with the node $x_{N_n-1}=L$ where the Dirichlet condition applies\n",
    "is appended to the system.\n",
    "This last row is anyway replaced by the condition $c_N=D$ or this\n",
    "condition can be incorporated in a symmetric fashion. Using the simplest,\n",
    "former approach gives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:BC:nat:Aub:system:mod\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{1}{h}\\left(\n",
    "\\begin{array}{ccccccccc}\n",
    "1 & -1 & 0 &\\cdots & \\cdots & \\cdots & \\cdots & \\cdots & 0 \\\\ \n",
    "-1 & 2 & -1 & \\ddots &   & &  & &  \\vdots \\\\ \n",
    "0 & -1 & 2 & -1 &\n",
    "\\ddots & &  &  & \\vdots \\\\ \n",
    "\\vdots & \\ddots &  & \\ddots & \\ddots & 0 &  & & \\vdots \\\\ \n",
    "\\vdots &  & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & & \\vdots \\\\ \n",
    "\\vdots & &  & 0 & -1 & 2 & -1 & \\ddots & \\vdots \\\\ \n",
    "\\vdots & & &  & \\ddots & \\ddots & \\ddots &\\ddots  & 0 \\\\ \n",
    "\\vdots & & & &  &\\ddots  & -1 & 2  & -1 \\\\ \n",
    "0 &\\cdots & \\cdots &\\cdots & \\cdots & \\cdots  & 0 & 0 & h\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "c_0 \\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots\\\\ \n",
    "c_{N}\n",
    "\\end{array}\n",
    "\\right)\n",
    "=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "h - C \\\\ \n",
    "2h\\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "2h \\\\ \n",
    "D\n",
    "\\end{array}\n",
    "\\right){\\thinspace .}\n",
    "\\label{fem:deq:1D:BC:nat:Aub:system:mod} \\tag{87}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Note: show equivalence between the two systems!!! -->\n",
    "\n",
    "## Cellwise computations\n",
    "\n",
    "Now we compute with one element at a time, working in the reference\n",
    "coordinate system $X\\in [-1,1]$.\n",
    "We need to see how the\n",
    "$u'(0)=C$ condition affects the element matrix and vector.\n",
    "The extra term $-C{\\varphi}_i(0)$ in the variational formulation\n",
    "only affects the element vector in the first cell.\n",
    "On the reference cell, $-C{\\varphi}_i(0)$ is transformed to\n",
    "$-C{\\tilde{\\varphi}}_r(-1)$, where $r$ counts local degrees of freedom.\n",
    "We have ${\\tilde{\\varphi}}_0(-1)=1$ and ${\\tilde{\\varphi}}_1(-1)=0$ so\n",
    "we are left with the contribution\n",
    "$-C{\\tilde{\\varphi}}_0(-1)=-C$ to $\\tilde b^{(0)}_0$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:ex1:Ab:elm:bc:nat\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\tilde A^{(0)} =\n",
    "A = \\frac{1}{h}\\left(\\begin{array}{rr}\n",
    "1 & 1\\\\ \n",
    "-1 & 1\n",
    "\\end{array}\\right),\\quad\n",
    "\\tilde b^{(0)} = \\left(\\begin{array}{c}\n",
    "h - C\\\\ \n",
    "h\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "\\label{fem:deq:1D:ex1:Ab:elm:bc:nat} \\tag{88}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No other element matrices or vectors are affected by the $-C{\\varphi}_i(0)$\n",
    "boundary term.\n",
    "\n",
    "There are two alternative ways of incorporating the Dirichlet condition.\n",
    "Following the section [Boundary term vanishes because of the test functions](#fem:deq:1D:BC:nat:uLtest), we get\n",
    "a $1\\times 1$ element matrix in the last cell and\n",
    "an element vector with an extra term containing $D$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:ex1:Ab:elm:ends2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\tilde A^{(e)} =\\frac{1}{h}\\left(\\begin{array}{r}\n",
    "1\n",
    "\\end{array}\\right),\\quad\n",
    "\\tilde b^{(e)} = h\\left(\\begin{array}{c}\n",
    "1 - D/6\n",
    "\\end{array}\\right),\n",
    "\\label{fem:deq:1D:ex1:Ab:elm:ends2} \\tag{89}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we include the degree of freedom at the node with\n",
    "$u$ specified. The element matrix and vector must then be modified\n",
    "to constrain the $\\tilde c_1 = c_N$ value at local node $r=1$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:ex1:Ab:elm:bc:nat:mod\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\tilde A^{(N_e)} =\n",
    "A = \\frac{1}{h}\\left(\\begin{array}{rr}\n",
    "1 & 1\\\\ \n",
    "0 & h\n",
    "\\end{array}\\right),\\quad\n",
    "\\tilde b^{(N_e)} = \\left(\\begin{array}{c}\n",
    "h\\\\ \n",
    "D\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "\\label{fem:deq:1D:ex1:Ab:elm:bc:nat:mod} \\tag{90}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Assemble and show that it is correct -->\n",
    "\n",
    "\n",
    "# Implementation of finite element algorithms\n",
    "<div id=\"fem:deq:1D:code\"></div>\n",
    "\n",
    "At this point, it is sensible to create a program with symbolic\n",
    "calculations to perform all the steps in the computational machinery,\n",
    "both for automating the work and for documenting the complete\n",
    "algorithms.  As we have seen, there are quite many details involved\n",
    "with finite element computations and incorporation of boundary\n",
    "conditions.  An implementation will also act as a structured summary\n",
    "of all these details.\n",
    "\n",
    "\n",
    "## Extensions of the code for approximation\n",
    "<div id=\"fem:deq:1D:code:fe\"></div>\n",
    "\n",
    "Implementation of the finite element algorithms for differential\n",
    "equations follows closely the algorithm for approximation of functions.\n",
    "The new additional ingredients are\n",
    "\n",
    "1. other types of integrands (as implied by the variational formulation)\n",
    "\n",
    "2. additional boundary terms in the variational formulation for\n",
    "   Neumann boundary conditions\n",
    "\n",
    "3. modification of element matrices and vectors due to Dirichlet\n",
    "   boundary conditions\n",
    "\n",
    "Point 1 and 2 can be taken care of by letting the user supply\n",
    "functions defining the integrands and boundary terms on the\n",
    "left- and right-hand side of the equation system:\n",
    "\n",
    " * Integrand on the left-hand side: `ilhs(e, phi, r, s, X, x, h)`\n",
    "\n",
    " * Integrand on the right-hand side: `irhs(e, phi, r, X, x, h)`\n",
    "\n",
    " * Boundary term on the left-hand side: `blhs (e, phi, r, s, X, x, h)`\n",
    "\n",
    " * Boundary term on the right-hand side: `brhs (e, phi, r, s, X, x, h)`\n",
    "\n",
    "Here, `phi` is a dictionary where `phi[q]` holds a list of the\n",
    "derivatives of order `q` of the basis functions with respect to the\n",
    "physical coordinate $x$.  The derivatives are available as Python\n",
    "functions of `X`.  For example, `phi[0][r](X)` means ${\\tilde{\\varphi}}_r(X)$,\n",
    "and `phi[1][s](X, h)` means $d{\\tilde{\\varphi}}_s (X)/dx$ (we refer to\n",
    "the file [fe1D.py](src/fe1D.py) for details\n",
    "regarding the function `basis` that computes the `phi`\n",
    "dictionary).  The `r` and `s`\n",
    "arguments in the above functions correspond to the index in the\n",
    "integrand contribution from an integration point to $\\tilde\n",
    "A^{(e)}_{r,s}$ and $\\tilde b^{(e)}_r$. The variables `e` and `h` are\n",
    "the current element number and the length of the cell, respectively.\n",
    "Specific examples below will make it clear how to construct these\n",
    "Python functions.\n",
    "\n",
    "Given a mesh represented by `vertices`, `cells`, and `dof_map` as\n",
    "explained before, we can write a pseudo Python code to list all\n",
    "the steps in the computational algorithm for finite element solution\n",
    "of a differential equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "<Declare global matrix and rhs: A, b>\n",
    "\n",
    "for e in range(len(cells)):\n",
    "\n",
    "    # Compute element matrix and vector\n",
    "    n = len(dof_map[e])  # no of dofs in this element\n",
    "    h = vertices[cells[e][1]] - vertices[cells[e][1]]\n",
    "    <Initialize element matrix and vector: A_e, b_e>\n",
    "\n",
    "    # Integrate over the reference cell\n",
    "    points, weights = <numerical integration rule>\n",
    "    for X, w in zip(points, weights):\n",
    "        phi = <basis functions and derivatives at X>\n",
    "        detJ = h/2\n",
    "        dX = detJ*w\n",
    "\n",
    "        x = <affine mapping from X>\n",
    "        for r in range(n):\n",
    "            for s in range(n):\n",
    "                A_e[r,s] += ilhs(e, phi, r, s, X, x, h)*dX\n",
    "            b_e[r] += irhs(e, phi, r, X, x, h)*dX\n",
    "\n",
    "    # Add boundary terms\n",
    "    for r in range(n):\n",
    "        for s in range(n):\n",
    "            A_e[r,s] += blhs(e, phi, r, s, X, x)*dX\n",
    "        b_e[r] += brhs(e, phi, r, X, x, h)*dX\n",
    "\n",
    "    # Incorporate essential boundary conditions\n",
    "    for r in range(n):\n",
    "        global_dof = dof_map[e][r]\n",
    "        if global_dof in essbc:\n",
    "            # local dof r is subject to an essential condition\n",
    "            value = essbc[global_dof]\n",
    "            # Symmetric modification\n",
    "            b_e -= value*A_e[:,r]\n",
    "            A_e[r,:] = 0\n",
    "            A_e[:,r] = 0\n",
    "            A_e[r,r] = 1\n",
    "            b_e[r] = value\n",
    "\n",
    "    # Assemble\n",
    "    for r in range(n):\n",
    "        for s in range(n):\n",
    "            A[dof_map[e][r], dof_map[e][s]] += A_e[r,s]\n",
    "        b[dof_map[e][r] += b_e[r]\n",
    "\n",
    "<solve linear system>\n",
    "```          "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having implemented this function, the user only has supply the `ilhs`,\n",
    "`irhs`, `blhs`, and `brhs` functions that specify the PDE and boundary\n",
    "conditions. The rest of the implementation forms a generic\n",
    "computational engine. The big finite element packages Diffpack and FEniCS\n",
    "are structured exactly the same way. A sample implementation of `ilhs`\n",
    "for the 1D Poisson problem is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "   def integrand_lhs(phi, i, j):\n",
    "        return phi[1][i]*phi[1][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which returns  $d{\\tilde{\\varphi}}_i (X)/dx d{\\tilde{\\varphi}}_j (X)/dx$.  Reducing the\n",
    "amount of code the user has to supply and making the code as close as\n",
    "possible to the mathematical formulation makes the software user-friendly and\n",
    "easy to debug.\n",
    "\n",
    "A complete function `finite_element1D_naive` for the 1D finite\n",
    "algorithm above, is found in the file [fe1D.py](src/fe1D.py). The term \"naive\" refers to a version of the\n",
    "algorithm where we use a standard dense square matrix as global matrix\n",
    "`A`. The implementation also has a verbose mode for printing out the\n",
    "element matrices and vectors as they are computed.  Below is the\n",
    "complete function without the print statements.  You should study in\n",
    "detail since it contains all the steps in the finite element\n",
    "algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def finite_element1D_naive(\n",
    "    vertices, cells, dof_map,     # mesh\n",
    "    essbc,                        # essbc[globdof]=value\n",
    "    ilhs,                         # integrand left-hand side\n",
    "    irhs,                         # integrand right-hand side\n",
    "    blhs=lambda e, phi, r, s, X, x, h: 0,\n",
    "    brhs=lambda e, phi, r, X, x, h: 0,\n",
    "    intrule='GaussLegendre',      # integration rule class\n",
    "    verbose=False,                # print intermediate results?\n",
    "    ):\n",
    "    N_e = len(cells)\n",
    "    N_n = np.array(dof_map).max() + 1\n",
    "\n",
    "    A = np.zeros((N_n, N_n))\n",
    "    b = np.zeros(N_n)\n",
    "\n",
    "    for e in range(N_e):\n",
    "        Omega_e = [vertices[cells[e][0]], vertices[cells[e][1]]]\n",
    "        h = Omega_e[1] - Omega_e[0]\n",
    "\n",
    "        d = len(dof_map[e]) - 1  # Polynomial degree\n",
    "        # Compute all element basis functions and their derivatives\n",
    "        phi = basis(d)\n",
    "\n",
    "        # Element matrix and vector\n",
    "        n = d+1  # No of dofs per element\n",
    "        A_e = np.zeros((n, n))\n",
    "        b_e = np.zeros(n)\n",
    "\n",
    "        # Integrate over the reference cell\n",
    "        if intrule == 'GaussLegendre':\n",
    "            points, weights = GaussLegendre(d+1)\n",
    "        elif intrule == 'NewtonCotes':\n",
    "            points, weights = NewtonCotes(d+1)\n",
    "\n",
    "        for X, w in zip(points, weights):\n",
    "            detJ = h/2\n",
    "            x = affine_mapping(X, Omega_e)\n",
    "            dX = detJ*w\n",
    "\n",
    "            # Compute contribution to element matrix and vector\n",
    "            for r in range(n):\n",
    "                for s in range(n):\n",
    "                    A_e[r,s] += ilhs(phi, r, s, X, x, h)*dX\n",
    "                b_e[r] += irhs(phi, r, X, x, h)*dX\n",
    "\n",
    "        # Add boundary terms\n",
    "        for r in range(n):\n",
    "            for s in range(n):\n",
    "                A_e[r,s] += blhs(phi, r, s, X, x, h)\n",
    "            b_e[r] += brhs(phi, r, X, x, h)\n",
    "\n",
    "        # Incorporate essential boundary conditions\n",
    "        modified = False\n",
    "        for r in range(n):\n",
    "            global_dof = dof_map[e][r]\n",
    "            if global_dof in essbc:\n",
    "                # dof r is subject to an essential condition\n",
    "                value = essbc[global_dof]\n",
    "                # Symmetric modification\n",
    "                b_e -= value*A_e[:,r]\n",
    "                A_e[r,:] = 0\n",
    "                A_e[:,r] = 0\n",
    "                A_e[r,r] = 1\n",
    "                b_e[r] = value\n",
    "                modified = True\n",
    "\n",
    "        # Assemble\n",
    "        for r in range(n):\n",
    "            for s in range(n):\n",
    "                A[dof_map[e][r], dof_map[e][s]] += A_e[r,s]\n",
    "            b[dof_map[e][r]] += b_e[r]\n",
    "\n",
    "    c = np.linalg.solve(A, b)\n",
    "    return c, A, b, timing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `timing` object is a dictionary holding the CPU spent on computing\n",
    "`A` and the CPU time spent on solving the linear system. (We have\n",
    "left out the timing statements.)\n",
    "\n",
    "## Utilizing a sparse matrix\n",
    "<div id=\"fem:deq:1D:code:fe_sparse\"></div>\n",
    "\n",
    "A potential efficiency problem with the `finite_element1D_naive` function\n",
    "is that it uses dense $(N+1)\\times (N+1)$ matrices, while we know that\n",
    "only $2d+1$ diagonals around the main diagonal are different from zero.\n",
    "Switching to a sparse matrix is very easy. Using the DOK (dictionary of\n",
    "keys) format, we declare `A` as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "import scipy.sparse\n",
    "A = scipy.sparse.dok_matrix((N_n, N_n))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assignments or in-place arithmetics are done as for a dense matrix,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        A[i,j] += term\n",
    "        A[i,j]  = term\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "but only the index pairs `(i,j)` we have used in assignments or\n",
    "in-place arithmetics are actually stored.\n",
    "A tailored solution algorithm is needed. The most reliable is\n",
    "sparse Gaussian elimination. SciPy gives access to the\n",
    "[UMFPACK](https://en.wikipedia.org/wiki/UMFPACK) algorithm\n",
    "for this purpose:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "import scipy.sparse.linalg\n",
    "c = scipy.sparse.linalg.spsolve(A.tocsr(), b, use_umfpack=True)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The declaration of `A` and the solve statement are the only\n",
    "changes needed in the `finite_element1D_naive` to utilize\n",
    "sparse matrices. The resulting modification is found in the\n",
    "function `finite_element1D`.\n",
    "\n",
    "## Application to our model problem\n",
    "\n",
    "Let us demonstrate the finite element software on"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "-u''(x)=f(x),\\quad x\\in (0,L),\\quad u'(0)=C,\\ u(L)=D{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This problem can be analytically solved by the\n",
    "`model2` function from the section \"Simple model problems and their solutions\" in previous chapter.\n",
    "Let $f(x)=x^2$. Calling `model2(x**2, L, C, D)` gives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u(x) = D + C(x-L) + \\frac{1}{12}(L^4 - x^4)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variational formulation reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(u', v) = (x^2, v) - Cv(0){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The entries in the element matrix and vector,\n",
    "which we need to set up the `ilhs`, `irhs`,\n",
    "`blhs`, and `brhs` functions, becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "A^{(e)}_{r,s} &= \\int_{-1}^1 \\frac{d{\\tilde{\\varphi}}_r}{dx}\\frac{{\\tilde{\\varphi}}_s}{dx}(\\det J{\\, \\mathrm{d}X}),\\\\ \n",
    "b^{(e)} &= \\int_{-1}^1 x^2{\\tilde{\\varphi}}_r\\det J{\\, \\mathrm{d}X} - C{\\tilde{\\varphi}}_r(-1)I(e,0),\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $I(e)$ is an indicator function: $I(e,q)=1$ if $e=q$, otherwise $I(e)=0$.\n",
    "We use this indicator function to formulate that the boundary term\n",
    "$Cv(0)$, which in the local element coordinate system becomes $C{\\tilde{\\varphi}}_r(-1)$,\n",
    "is only included for the element $e=0$.\n",
    "\n",
    "The functions for specifying the element matrix and vector entries\n",
    "must contain the integrand, but without the $\\det J{\\, \\mathrm{d}X}$ term\n",
    "as this term is taken care of by the quadrature loop, and\n",
    "the derivatives $d{\\tilde{\\varphi}}_r(X)/dx$\n",
    "with respect to the physical $x$ coordinates are\n",
    "contained in `phi[1][r](X)`, computed by the function `basis`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ilhs(e, phi, r, s, X, x, h):\n",
    "    return phi[1][r](X, h)*phi[1][s](X, h)\n",
    "\n",
    "def irhs(e, phi, r, X, x, h):\n",
    "    return x**2*phi[0][r](X)\n",
    "\n",
    "def blhs(e, phi, r, s, X, x, h):\n",
    "    return 0\n",
    "\n",
    "def brhs(e, phi, r, X, x, h):\n",
    "    return -C*phi[0][r](-1) if e == 0 else 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then make the call to `finite_element1D_naive` or `finite_element1D`\n",
    "to solve the problem with two P1 elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.fe1D import finite_element1D_naive, mesh_uniform\n",
    "C = 5;  D = 2;  L = 4\n",
    "d = 1\n",
    "\n",
    "vertices, cells, dof_map = mesh_uniform(\n",
    "    N_e=2, d=d, Omega=[0,L], symbolic=False)\n",
    "essbc = {}\n",
    "essbc[dof_map[-1][-1]] = D\n",
    "\n",
    "c, A, b, timing = finite_element1D_naive(\n",
    "    vertices, cells, dof_map, essbc,\n",
    "    ilhs=ilhs, irhs=irhs, blhs=blhs, brhs=brhs,\n",
    "    intrule='GaussLegendre')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It remains to plot the solution (with high resolution in each element).\n",
    "To this end, we use the `u_glob` function imported from\n",
    "`fe1D`, which imports it from `fe_approx1D_numit` (the\n",
    "`u_glob` function in `fe_approx1D.py`\n",
    "works with `elements` and `nodes`, while `u_glob` in\n",
    "`fe_approx1D_numint` works with `cells`, `vertices`,\n",
    "and `dof_map`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_exact = lambda x: D + C*(x-L) + (1./6)*(L**3 - x**3)\n",
    "from src.fe1D import u_glob\n",
    "x, u, nodes = u_glob(c, cells, vertices, dof_map)\n",
    "u_e = u_exact(x)\n",
    "# print((u_exact(nodes) - c))  # difference at the nodes\n",
    "\n",
    "# import matplotlib.pyplot as plt\n",
    "# plt.plot(x, u, 'b-', x, u_e, 'r--')\n",
    "# plt.legend(['finite elements, d=%d' %d, 'exact'], loc='upper left')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is shown in [Figure](#fem:deq:1D:code:fe:fig1). We see\n",
    "that the solution using P1 elements is exact at the nodes, but\n",
    "feature considerable discrepancy between the nodes.\n",
    "[Exercise 10: Investigate exact finite element solutions](#fem:deq:exer:1D:exact_numerics) asks you to explore\n",
    "this problem further using other $m$ and $d$ values.\n",
    "\n",
    "<!-- dom:FIGURE: [fig/uxx_x2.png, width=500 frac=0.8] Finite element and exact solution using two cells. <div id=\"fem:deq:1D:code:fe:fig1\"></div> -->\n",
    "<!-- begin figure -->\n",
    "<div id=\"fem:deq:1D:code:fe:fig1\"></div>\n",
    "\n",
    "<p>Finite element and exact solution using two cells.</p>\n",
    "<img src=\"fig/uxx_x2.png\" width=500>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "\n",
    "# Variational formulations in 2D and 3D\n",
    "<div id=\"fem:deq:2D:varform\"></div>\n",
    "\n",
    "The major difference between deriving variational formulations in 2D\n",
    "and 3D compared to 1D is the rule for integrating by parts.  The cells\n",
    "have shapes different from an interval, so basis functions look a bit\n",
    "different, and there is a technical difference in actually calculating\n",
    "the integrals over cells. Otherwise, going to 2D and 3D is not a big\n",
    "step from 1D. All the fundamental ideas still apply.\n",
    "\n",
    "## Integration by parts\n",
    "\n",
    "A typical second-order term in a PDE may be written in dimension-independent\n",
    "notation as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\nabla^2 u \\quad\\hbox{or}\\quad \\nabla\\cdot\\left( {\\alpha}(\\boldsymbol{x})\\nabla u\\right)\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The explicit forms in a 2D problem become"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\nabla^2 u = \\nabla\\cdot\\nabla u =\n",
    "\\frac{\\partial^2 u}{\\partial x^2} +\n",
    "\\frac{\\partial^2 u}{\\partial y^2},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\nabla\\cdot\\left( a(\\boldsymbol{x})\\nabla u\\right) =\n",
    "\\frac{\\partial}{\\partial x}\\left( {\\alpha}(x,y)\\frac{\\partial u}{\\partial x}\\right) +\n",
    "\\frac{\\partial}{\\partial y}\\left( {\\alpha}(x,y)\\frac{\\partial u}{\\partial y}\\right)\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We shall continue with the latter operator as the former arises from\n",
    "just setting ${\\alpha} =1$.\n",
    "\n",
    "**The integration by parts formula for $\\int\\nabla\\cdot({\\alpha}\\nabla)$.**\n",
    "\n",
    "The general rule for integrating by parts is often referred to as\n",
    "[Green's first identity](http://en.wikipedia.org/wiki/Green's_identities):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:2D:int:by:parts\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "-\\int_{\\Omega} \\nabla\\cdot ({\\alpha}(\\boldsymbol{x})\\nabla u) v{\\, \\mathrm{d}x} =\n",
    "\\int_{\\Omega} {\\alpha}(\\boldsymbol{x})\\nabla u\\cdot\\nabla v {\\, \\mathrm{d}x} -\n",
    "\\int_{\\partial\\Omega} a\\frac{\\partial u}{\\partial n} v {\\, \\mathrm{d}s},\n",
    "\\label{fem:deq:2D:int:by:parts} \\tag{91}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $\\partial\\Omega$ is the boundary of $\\Omega$ and\n",
    "$\\partial u/\\partial n = \\boldsymbol{n}\\cdot\\nabla u$ is the derivative\n",
    "of $u$ in the outward normal direction, $\\boldsymbol{n}$ being an outward\n",
    "unit normal to $\\partial\\Omega$. The integrals $\\int_\\Omega (){\\, \\mathrm{d}x}$ are\n",
    "area integrals in 2D and volume integrals in 3D, while\n",
    "$\\int_{\\partial\\Omega} (){\\, \\mathrm{d}s}$ is a line integral in 2D and a surface\n",
    "integral in 3D.\n",
    "\n",
    "\n",
    "\n",
    "It will be convenient to divide the boundary into two parts:\n",
    "\n",
    " * $\\partial\\Omega_N$, where we have Neumann conditions\n",
    "   $-a\\frac{\\partial u}{\\partial n} = g$, and\n",
    "\n",
    " * $\\partial\\Omega_D$, where we have Dirichlet conditions\n",
    "   $u = u_0$.\n",
    "\n",
    "The test functions $v$ are (as usual) required to vanish on\n",
    "$\\partial\\Omega_D$.\n",
    "\n",
    "## Example on a multi-dimensional variational problem\n",
    "<div id=\"sec:varform:general:convdiff\"></div>\n",
    "\n",
    "Here is a quite general, stationary, linear PDE arising in many problems:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"varform:conv:diff:pde:pre\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\label{varform:conv:diff:pde:pre} \\tag{92}\n",
    "\\boldsymbol{v}\\cdot\\nabla u + \\beta u = \\nabla\\cdot\\left( {\\alpha}\\nabla u\\right) + f,\n",
    "\\quad\\boldsymbol{x}\\in\\Omega,\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"varform:conv:diff:bc1:pre\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "\\label{varform:conv:diff:bc1:pre} \\tag{93}\n",
    "u = u_0,\\quad\\boldsymbol{x}\\in\\partial\\Omega_D,\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"varform:conv:diff:bc2:pre\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "\\label{varform:conv:diff:bc2:pre} \\tag{94}\n",
    "-{\\alpha}\\frac{\\partial u}{\\partial n} = g,\\quad\\boldsymbol{x}\\in\\partial\\Omega_N\n",
    "{\\thinspace .}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vector field $\\boldsymbol{v}$ and the scalar functions $a$, $\\alpha$, $f$, $u_0$, and\n",
    "$g$ may vary with the spatial coordinate $\\boldsymbol{x}$ and must be known.\n",
    "\n",
    "Such a second-order PDE needs exactly one boundary condition at each\n",
    "point of the boundary, so $\\partial\\Omega_N\\cup\\partial\\Omega_D$\n",
    "must be the complete boundary $\\partial\\Omega$.\n",
    "\n",
    "Assume that the boundary function $u_0(\\boldsymbol{x})$ is defined for all $\\boldsymbol{x}\\in\\Omega$.\n",
    "The unknown function can then be expanded as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u = B + \\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_j,\\quad B = u_0 {\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As long as any ${\\psi}_j=0$ on $\\partial\\Omega_D$, we realize that $u=u_0$\n",
    "on $\\partial\\Omega_D$.\n",
    "\n",
    "The variational formula is obtained from Galerkin's method, which\n",
    "technically means multiplying the PDE by a test\n",
    "function $v$ and integrating over $\\Omega$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\int_{\\Omega} (\\boldsymbol{v}\\cdot\\nabla u + \\beta u)v{\\, \\mathrm{d}x} =\n",
    "\\int_{\\Omega} \\nabla\\cdot\\left( {\\alpha}\\nabla u\\right){\\, \\mathrm{d}x} + \\int_{\\Omega}fv {\\, \\mathrm{d}x}\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second-order term is integrated by parts, according to the formula\n",
    "([91](#fem:deq:2D:int:by:parts)):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\int_{\\Omega} \\nabla\\cdot\\left( {\\alpha}\\nabla u\\right)v {\\, \\mathrm{d}x} =\n",
    "-\\int_{\\Omega} {\\alpha}\\nabla u\\cdot\\nabla v{\\, \\mathrm{d}x}\n",
    "+ \\int_{\\partial\\Omega} {\\alpha}\\frac{\\partial u}{\\partial n} v{\\, \\mathrm{d}s}\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Galerkin's method therefore leads to"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\int_{\\Omega} (\\boldsymbol{v}\\cdot\\nabla u + \\beta u)v{\\, \\mathrm{d}x} =\n",
    "-\\int_{\\Omega} {\\alpha}\\nabla u\\cdot\\nabla v{\\, \\mathrm{d}x}\n",
    "+ \\int_{\\partial\\Omega} {\\alpha}\\frac{\\partial u}{\\partial n} v{\\, \\mathrm{d}s}\n",
    "+ \\int_{\\Omega} fv {\\, \\mathrm{d}x}\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The boundary term can be developed further by noticing that $v\\neq 0$\n",
    "only on $\\partial\\Omega_N$,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\int_{\\partial\\Omega} {\\alpha}\\frac{\\partial u}{\\partial n} v{\\, \\mathrm{d}s}\n",
    "= \\int_{\\partial\\Omega_N} {\\alpha}\\frac{\\partial u}{\\partial n} v{\\, \\mathrm{d}s},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and that on $\\partial\\Omega_N$, we have the condition\n",
    "$a\\frac{\\partial u}{\\partial n}=-g$, so the term becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "-\\int_{\\partial\\Omega_N} gv{\\, \\mathrm{d}s}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final variational form is then"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\int_{\\Omega} (\\boldsymbol{v}\\cdot\\nabla u + \\beta u)v{\\, \\mathrm{d}x} =\n",
    "-\\int_{\\Omega} {\\alpha}\\nabla u\\cdot\\nabla v {\\, \\mathrm{d}x}\n",
    "- \\int_{\\partial\\Omega_N} g v{\\, \\mathrm{d}s}\n",
    "+ \\int_{\\Omega} fv {\\, \\mathrm{d}x}\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of using the integral signs, we may use the inner product\n",
    "notation:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(\\boldsymbol{v}\\cdot\\nabla u, v) + (\\beta u,v) =\n",
    "- ({\\alpha}\\nabla u,\\nabla v) - (g,v)_{N} + (f,v)\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The subscript $\\,{}_N$ in $(g,v)_{N}$ is a notation for a line or surface\n",
    "integral over $\\partial\\Omega_N$, while $(\\cdot,\\cdot)$ is the area/volume\n",
    "integral over $\\Omega$.\n",
    "\n",
    "We can derive explicit expressions for the linear system for $\\left\\{ {c}_j \\right\\}_{j\\in{\\mathcal{I}_s}}$\n",
    "that arises from the variational formulation.\n",
    "Inserting the $u$ expansion results in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\sum_{j\\in{\\mathcal{I}_s}} ((\\boldsymbol{v}\\cdot\\nabla {\\psi}_j, {\\psi}_i) &+ (\\beta {\\psi}_j ,{\\psi}_i) + ({\\alpha}\\nabla {\\psi}_j,\\nabla {\\psi}_i))c_j = \\\\ \n",
    "& (g,{\\psi}_i)_{N} + (f,{\\psi}_i) -\n",
    "(\\boldsymbol{v}\\cdot\\nabla u_0, {\\psi}_i) + (\\beta u_0 ,{\\psi}_i) +\n",
    "({\\alpha}\\nabla u_0,\\nabla {\\psi}_i)\n",
    "{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a linear system with matrix entries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A_{i,j} = (\\boldsymbol{v}\\cdot\\nabla {\\psi}_j, {\\psi}_i) + (\\beta {\\psi}_j ,{\\psi}_i) + ({\\alpha}\\nabla {\\psi}_j,\\nabla {\\psi}_i)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and right-hand side entries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "b_i = (g,{\\psi}_i)_{N} + (f,{\\psi}_i) -\n",
    "(\\boldsymbol{v}\\cdot\\nabla u_0, {\\psi}_i) + (\\beta u_0 ,{\\psi}_i) +\n",
    "({\\alpha}\\nabla u_0,\\nabla {\\psi}_i),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for $i,j\\in{\\mathcal{I}_s}$.\n",
    "\n",
    "In the finite element method, we usually express $u_0$ in terms of\n",
    "basis functions and restrict $i$ and $j$ to run over the degrees of\n",
    "freedom that are not prescribed as Dirichlet conditions.\n",
    "However, we can also keep all the $\\left\\{ {c}_j \\right\\}_{j\\in{\\mathcal{I}_s}}$ as unknowns,\n",
    "drop the $u_0$ in the expansion for $u$, and incorporate all the\n",
    "known $c_j$ values in the linear system. This has been explained\n",
    "in detail in the 1D case, and the technique is the same for 2D and\n",
    "3D problems.\n",
    "\n",
    "## Transformation to a reference cell in 2D and 3D\n",
    "\n",
    "The real power of the finite element method first becomes evident when\n",
    "we want to solve partial differential equations posed on two- and\n",
    "three-dimensional domains of non-trivial geometric shape.  As in 1D,\n",
    "the domain $\\Omega$ is divided into $N_e$ non-overlapping cells. The\n",
    "elements have simple shapes: triangles and quadrilaterals are popular\n",
    "in 2D, while tetrahedra and box-shapes elements dominate in 3D.  The\n",
    "finite element basis functions ${\\varphi}_i$ are, as in 1D, polynomials\n",
    "over each cell.  The integrals in the variational formulation are, as\n",
    "in 1D, split into contributions from each cell, and these\n",
    "contributions are calculated by mapping a physical cell, expressed in\n",
    "physical coordinates $\\boldsymbol{x}$, to a reference cell in a local coordinate\n",
    "system $\\boldsymbol{X}$. This mapping will now be explained in detail.\n",
    "\n",
    "\n",
    "We consider an integral of the type"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto39\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\int_{{\\Omega}^{(e)}} {\\alpha}(\\boldsymbol{x})\\nabla{\\varphi}_i\\cdot\\nabla{\\varphi}_j{\\, \\mathrm{d}x},\n",
    "\\label{_auto39} \\tag{95}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where the ${\\varphi}_i$ functions are finite element basis functions in\n",
    "2D or 3D, defined in the physical domain.\n",
    "Suppose we want to calculate this integral over a reference cell,\n",
    "denoted by $\\tilde\\Omega^r$, in a coordinate system with coordinates\n",
    "$\\boldsymbol{X} = (X_0, X_1)$ (2D) or $\\boldsymbol{X} = (X_0, X_1, X_2)$ (3D).\n",
    "The mapping between a point $\\boldsymbol{X}$ in the reference coordinate system  and\n",
    "the corresponding point $\\boldsymbol{x}$ in the physical coordinate system is\n",
    "given by a vector relation $\\boldsymbol{x}(\\boldsymbol{X})$.\n",
    "The corresponding Jacobian, $J$, of this mapping has entries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "J_{i,j}=\\frac{\\partial x_j}{\\partial X_i}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The change of variables requires ${\\, \\mathrm{d}x}$ to be replaced by $\\det J{\\, \\mathrm{d}X}$.\n",
    "The derivatives in the $\\nabla$ operator in the variational form are\n",
    "with respect to $\\boldsymbol{x}$, which we may denote by $\\nabla_{\\boldsymbol{x}}$.\n",
    "The ${\\varphi}_i(\\boldsymbol{x})$ functions in the integral\n",
    "are replaced by local basis functions ${\\tilde{\\varphi}}_r(\\boldsymbol{X})$ so\n",
    "the integral features $\\nabla_{\\boldsymbol{x}}{\\tilde{\\varphi}}_r(\\boldsymbol{X})$. We readily have\n",
    "$\\nabla_{\\boldsymbol{X}}{\\tilde{\\varphi}}_r(\\boldsymbol{X})$ from formulas for the basis functions in\n",
    "the reference cell, but\n",
    "the desired quantity $\\nabla_{\\boldsymbol{x}}{\\tilde{\\varphi}}_r(\\boldsymbol{X})$ requires some efforts\n",
    "to compute. All the details are provided below.\n",
    "\n",
    "Let $i=q(e,r)$ and consider two space dimensions. By the chain rule,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{\\partial {\\tilde{\\varphi}}_r}{\\partial X} =\n",
    "\\frac{\\partial {\\varphi}_i}{\\partial X} =\n",
    "\\frac{\\partial {\\varphi}_i}{\\partial x}\\frac{\\partial x}{\\partial X} +\n",
    "\\frac{\\partial {\\varphi}_i}{\\partial y}\\frac{\\partial y}{\\partial X},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{\\partial {\\tilde{\\varphi}}_r}{\\partial Y} =\n",
    "\\frac{\\partial {\\varphi}_i}{\\partial Y} =\n",
    "\\frac{\\partial {\\varphi}_i}{\\partial x}\\frac{\\partial x}{\\partial Y} +\n",
    "\\frac{\\partial {\\varphi}_i}{\\partial y}\\frac{\\partial y}{\\partial Y}\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write these two equations as a vector equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left[\\begin{array}{c}\n",
    "\\frac{\\partial {\\tilde{\\varphi}}_r}{\\partial X}\\\\ \n",
    "\\frac{\\partial {\\tilde{\\varphi}}_r}{\\partial Y}\n",
    "\\end{array}\\right]\n",
    "=\n",
    "\\left[\\begin{array}{cc}\n",
    "\\frac{\\partial x}{\\partial X} & \\frac{\\partial y}{\\partial X}\\\\ \n",
    "\\frac{\\partial x}{\\partial Y} & \\frac{\\partial y}{\\partial Y}\n",
    "\\end{array}\\right]\n",
    "\\left[\\begin{array}{c}\n",
    "\\frac{\\partial {\\varphi}_i}{\\partial x}\\\\ \n",
    "\\frac{\\partial {\\varphi}_i}{\\partial y}\n",
    "\\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Identifying"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\nabla_{\\boldsymbol{X}}{\\tilde{\\varphi}}_r = \\left[\\begin{array}{c}\n",
    "\\frac{\\partial {\\tilde{\\varphi}}_r}{\\partial X}\\\\ \n",
    "\\frac{\\partial {\\tilde{\\varphi}}_r}{\\partial Y}\n",
    "\\end{array}\\right],\n",
    "\\quad\n",
    "J =\n",
    "\\left[\\begin{array}{cc}\n",
    "\\frac{\\partial x}{\\partial X} & \\frac{\\partial y}{\\partial X}\\\\ \n",
    "\\frac{\\partial x}{\\partial Y} & \\frac{\\partial y}{\\partial Y}\n",
    "\\end{array}\\right],\n",
    "\\quad\n",
    "\\nabla_{\\boldsymbol{x}}{\\varphi}_r =\n",
    "\\left[\\begin{array}{c}\n",
    "\\frac{\\partial {\\varphi}_i}{\\partial x}\\\\ \n",
    "\\frac{\\partial {\\varphi}_i}{\\partial y}\n",
    "\\end{array}\\right],\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we have the relation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\nabla_{\\boldsymbol{X}}{\\tilde{\\varphi}}_r = J\\cdot\\nabla_{\\boldsymbol{x}}{\\varphi}_i,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which we can solve with respect to $\\nabla_{\\boldsymbol{x}}{\\varphi}_i$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto40\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\nabla_{\\boldsymbol{x}}{\\varphi}_i = J^{-1}\\cdot\\nabla_{\\boldsymbol{X}}{\\tilde{\\varphi}}_r{\\thinspace .}\n",
    "\\label{_auto40} \\tag{96}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the reference cell, ${\\varphi}_i(\\boldsymbol{x}) = {\\tilde{\\varphi}}_r(\\boldsymbol{X})$, so"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto41\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\nabla_{\\boldsymbol{x}}{\\tilde{\\varphi}}_r(\\boldsymbol{X}) = J^{-1}(\\boldsymbol{X})\\cdot\\nabla_{\\boldsymbol{X}}{\\tilde{\\varphi}}_r(\\boldsymbol{X}){\\thinspace .}\n",
    "\\label{_auto41} \\tag{97}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means that we have the following transformation of the\n",
    "integral in the physical domain to its counterpart over the reference cell:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto42\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\int_{\\Omega^{(e)}} {\\alpha}(\\boldsymbol{x})\\nabla_{\\boldsymbol{x}}{\\varphi}_i\\cdot\\nabla_{\\boldsymbol{x}}{\\varphi}_j{\\, \\mathrm{d}x} =\n",
    "\\int_{\\tilde\\Omega^r} {\\alpha}(\\boldsymbol{x}(\\boldsymbol{X}))(J^{-1}\\cdot\\nabla_{\\boldsymbol{X}}{\\tilde{\\varphi}}_r)\\cdot\n",
    "(J^{-1}\\cdot\\nabla{\\tilde{\\varphi}}_s)\\det J{\\, \\mathrm{d}X}\n",
    "\\label{_auto42} \\tag{98}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical integration\n",
    "\n",
    "Integrals are normally computed by numerical integration rules.\n",
    "For multi-dimensional cells, various families of rules exist.\n",
    "All of them are similar to what is shown in 1D:\n",
    "$\\int f {\\, \\mathrm{d}x}\\approx \\sum_jw_if(\\boldsymbol{x}_j)$, where $w_j$ are weights and\n",
    "$\\boldsymbol{x}_j$ are corresponding points.\n",
    "\n",
    "The file [numint.py](src/numint.py) contains the functions\n",
    "`quadrature_for_triangles(n)` and `quadrature_for_tetrahedra(n)`,\n",
    "which returns lists of points and weights corresponding to integration\n",
    "rules with `n` points over the reference triangle\n",
    "with vertices $(0,0)$, $(1,0)$, $(0,1)$, and the reference tetrahedron\n",
    "with vertices $(0,0,0)$, $(1,0,0)$, $(0,1,0)$, $(0,0,1)$,\n",
    "respectively. For example, the first two rules for integration over\n",
    "a triangle have 1 and 3 points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numint\n",
    "x, w = numint.quadrature_for_triangles(num_points=1)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, w = numint.quadrature_for_triangles(num_points=3)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rules with 1, 3, 4, and 7 points over the triangle will exactly integrate\n",
    "polynomials of degree 1, 2, 3, and 4, respectively.\n",
    "In 3D, rules with 1, 4, 5, and 11 points over the tetrahedron will\n",
    "exactly integrate polynomials of degree 1, 2, 3, and 4, respectively.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "## Convenient formulas for P1 elements in 2D\n",
    "\n",
    "We shall now provide some formulas for piecewise linear ${\\varphi}_i$ functions\n",
    "and their integrals *in the physical coordinate system*.\n",
    "These formulas make it convenient to compute with P1 elements without\n",
    "the need to work in the reference coordinate system and deal with mappings\n",
    "and Jacobians.\n",
    "A lot of computational and algorithmic details are hidden by this approach.\n",
    "\n",
    "Let $\\Omega^{(e)}$ be cell number $e$, and let the three vertices\n",
    "have global vertex numbers $I$, $J$, and $K$.\n",
    "The corresponding coordinates are\n",
    "$(x_{I},y_{I})$, $(x_{J},y_{J})$, and $(x_{K},y_{K})$.\n",
    "The basis function ${\\varphi}_I$ over $\\Omega^{(e)}$ have the explicit\n",
    "formula"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:approx:fe:2D:phi:I\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "{\\varphi}_I (x,y) = \\frac{1}{2}\\Delta \\left( \\alpha_I + \\beta_Ix\n",
    "+ \\gamma_Iy\\right),\n",
    "\\label{fem:approx:fe:2D:phi:I} \\tag{99}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where\n",
    "<!-- must split align in two because we need an array with & and \\\\ -->\n",
    "<!-- (sphinx, ipynb, pandoc requires splitting of align and & in the -->\n",
    "<!-- array confuses the splitting) -->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:approx:fe:2D:phi:alpha:I\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\alpha_I = x_{J}y_{K} - x_{K}y_{J},\n",
    "\\label{fem:approx:fe:2D:phi:alpha:I} \\tag{100}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:approx:fe:2D:phi:beta:I\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "\\beta_I = y_{J} - y_{K},\n",
    "\\label{fem:approx:fe:2D:phi:beta:I} \\tag{101}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:approx:fe:2D:phi:gamma:I\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "\\gamma_I = x_{K} - x_{J},\n",
    "\\label{fem:approx:fe:2D:phi:gamma:I} \\tag{102},\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:approx:fe:2D:phi:Delta\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "2\\Delta = \\det\\left(\\begin{array}{rrr}\n",
    "1 & x_{I} & y_{I} \\\\ \n",
    "1 & x_{J} & y_{J} \\\\ \n",
    "1 & x_{K} & y_{K} \\end{array}\\right)\n",
    "{\\thinspace .}\n",
    "\\label{fem:approx:fe:2D:phi:Delta} \\tag{103}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quantity $\\Delta$ is the area of the cell.\n",
    "\n",
    "The following formula is often convenient when computing element matrices\n",
    "and vectors:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:approx:fe:2D:phi:integral\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\int_{\\Omega^{(e)}} {\\varphi}_I^{p}{\\varphi}_J^{q}{\\varphi}_K^{r} dx dy =\n",
    "{p!q!r!\\over (p+q+r+2)!}2\\Delta\n",
    "\\label{fem:approx:fe:2D:phi:integral} \\tag{104}\n",
    "{\\thinspace .}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Note that the $q$ in this formula is not to be mixed with the $q(e,r)$\n",
    "mapping of degrees of freedom.)\n",
    "\n",
    "As an example, the element matrix entry\n",
    "$\\int_{\\Omega^{(e)}} {\\varphi}_I{\\varphi}_J{\\, \\mathrm{d}x}$\n",
    "can be computed by setting\n",
    "$p=q=1$ and $r=0$, when $I\\neq J$, yielding $\\Delta/12$, and\n",
    "$p=2$ and $q=r=0$, when $I=J$, resulting in $\\Delta/6$.\n",
    "We collect these numbers in a local element matrix:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{\\Delta}{12}\n",
    "\\left[\\begin{array}{ccc}\n",
    "2 & 1 & 1\\\\ \n",
    "1 & 2 & 1\\\\ \n",
    "1 & 1 & 2\n",
    "\\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The common element matrix entry $\\int_{\\Omega^{(e)}} \\nabla{\\varphi}_I\\cdot\\nabla{\\varphi}_J{\\, \\mathrm{d}x}$, arising from a Laplace term $\\nabla^2u$, can also easily be\n",
    "computed by the formulas above. We have"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\nabla{\\varphi}_I\\cdot\\nabla{\\varphi}_J =\n",
    "\\frac{\\Delta^2}{4}(\\beta_I\\beta_J + \\gamma_I\\gamma_J) = \\hbox{const},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "so that the element matrix entry becomes\n",
    "$\\frac{1}{4}\\Delta^3(\\beta_I\\beta_J + \\gamma_I\\gamma_J)$.\n",
    "\n",
    "From an implementational point of view, one will work with local vertex\n",
    "numbers $r=0,1,2$, parameterize the coefficients in the basis\n",
    "functions by $r$, and look up vertex coordinates through $q(e,r)$.\n",
    "\n",
    "Similar formulas exist for integration of P1 elements in 3D."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementation in 2D and 3D via FEniCS\n",
    "<div id=\"fem:varform:fenics\"></div>\n",
    "\n",
    "From a principle of view, we have seen that variational forms of the\n",
    "type: find $a(u,v)=L\\ \\forall v\\in V$ (and even general nonlinear problems\n",
    "$F(u;v)=0$), can apply the computational machinery of introduced for\n",
    "the approximation problem $u=f$. We actually need two extensions only:\n",
    "\n",
    "1. specify Dirichlet boundary conditions as part of $V$\n",
    "\n",
    "2. incorporate Neumann flux boundary conditions in the variational form\n",
    "\n",
    "The algorithms are all the same in any space dimension, we only need to\n",
    "choose the element type and associated integration rule. Once we know\n",
    "how to compute things in 1D, and made the computer code sufficiently\n",
    "flexible, the method and code should work for any variational form in\n",
    "any number of space dimensions! This fact is exactly the idea behind\n",
    "the [FEniCS](http://fenicsproject.org) finite element software.\n",
    "\n",
    "Therefore, if we know how to set up an approximation problem in any\n",
    "dimension in FEniCS, and know how to derive variational forms in higher\n",
    "dimensions, we are (in principle!) very close to solving a\n",
    "PDE problem in FEniCS. We shall now solve a quite general 1D/2D/3D Poisson problem in FEniCS.\n",
    "There is much more to FEniCS than what is shown in this example, but it\n",
    "illustrates the fact that when we go beyond 1D, there exists\n",
    "software which leverage the full power of the finite element method as\n",
    "a method for solving any problem on any mesh in any number of\n",
    "space dimensions.\n",
    "\n",
    "## Mathematical problem\n",
    "<div id=\"fem:varform:fenics:problem\"></div>\n",
    "\n",
    "The following model describes the pressure $u$ in the flow around a\n",
    "bore hole of radius $a$ in a porous medium. If the hole is long in the\n",
    "vertical direction (z-direction) then it is natural to assume that the \n",
    "vertical changes are small and $u_z\\approx\\mbox{constant}$. \n",
    "Therefore, we can model it by a 2D domain in the cross section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:varform:fenics:problem:pde\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\nabla\\cdot \\left( {\\alpha}\\nabla u\\right) = 0,  \\quad a < ||\\boldsymbol{x}|| < b,\n",
    "\\label{fem:varform:fenics:problem:pde} \\tag{105}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:varform:fenics:problem:ua\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "u(\\boldsymbol{x}) = U_a, \\quad   ||\\boldsymbol{x}|| = a,\n",
    "\\label{fem:varform:fenics:problem:ua} \\tag{106}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:varform:fenics:problem:ub\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "u(\\boldsymbol{x}) = U_b  \\quad   ||\\boldsymbol{x}|| = b{\\thinspace .}\n",
    "\\label{fem:varform:fenics:problem:ub} \\tag{107}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is, we have a hollow circular 2D domain with inner radius $a$ and\n",
    "outer radius $b$. The pressure is known on these two boundaries, so\n",
    "this is a pure Dirichlet problem.\n",
    "\n",
    "### Symmetry\n",
    "\n",
    "The first thing we should observe is that the problem is radially\n",
    "symmetric, so we can change to polar coordinates and obtain a 1D\n",
    "problem in the radial direction:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(r{\\alpha} u')' = 0,\\quad u(a)=U_a, u(b)=U_b{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is not very exciting beyond being able to find an analytical solution\n",
    "and compute the true error of a finite element approximation.\n",
    "\n",
    "However, many software packages solve problems in Cartesian coordinates, and\n",
    "FEniCS basically do this, so we want to take advantage of symmetry in\n",
    "Cartesian coordinates and reformulate the problem in a smaller\n",
    "domain.\n",
    "\n",
    "Looking at the domain as a cake with a hole, any piece of the\n",
    "cake will be a candidate for a reduced-size domain.  The solution is\n",
    "symmetric about any line $\\theta = \\hbox{const}$ in polar coordinates,\n",
    "so at such lines we have the symmetry boundary condition $\\partial\n",
    "u/\\partial n=0$, i.e., a homogeneous Neumann condition.  In [Figure](#fem:varform:fenics:problem:meshfig) we have plotted a possible\n",
    "mesh of cells as triangles, here with dense refinement toward the bore\n",
    "hole, because we know the solution will decay most rapidly toward the\n",
    "origin.  This mesh is a piece of the cake with four sides: Dirichlet\n",
    "conditions on the inner and outer boundary, named $\\Gamma_{D_a}$ and\n",
    "$\\Gamma_{D_b}$, and $\\partial u/\\partial n=0$\n",
    "on the two other sides, named $\\Gamma_N$.\n",
    "In this particular example, the arc of the piece\n",
    "of the cake is 45 degrees, but any value of the arc will work.\n",
    "\n",
    "<!-- dom:FIGURE: [fig/borehole_mesh1.png, width=700 frac=1.2] Mesh of a hollow cylinder, with refinement and utilizing symmetry. <div id=\"fem:varform:fenics:problem:meshfig\"></div> -->\n",
    "<!-- begin figure -->\n",
    "<div id=\"fem:varform:fenics:problem:meshfig\"></div>\n",
    "\n",
    "<p>Mesh of a hollow cylinder, with refinement and utilizing symmetry.</p>\n",
    "<img src=\"fig/borehole_mesh1.png\" width=700>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "The boundary problem can then be expressed as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:varform:fenics:problem:pde2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\nabla\\cdot \\left( {\\alpha}\\nabla u\\right) = 0,  \\quad \\boldsymbol{x}\\in\\Omega,\n",
    "\\label{fem:varform:fenics:problem:pde2} \\tag{108}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:varform:fenics:problem:Ga\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "u(\\boldsymbol{x}) = U_a,  \\quad   \\boldsymbol{x}\\in\\Gamma_{D_a},\n",
    "\\label{fem:varform:fenics:problem:Ga} \\tag{109}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:varform:fenics:problem:Gb\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "u(\\boldsymbol{x}) = U_b,  \\quad   \\boldsymbol{x}\\in\\Gamma_{D_b},\n",
    "\\label{fem:varform:fenics:problem:Gb} \\tag{110}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:varform:fenics:problem:GN\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "\\frac{\\partial u}{\\partial n} =0,\\quad  \\boldsymbol{x}\\in\\Gamma_N{\\thinspace .}\n",
    "\\label{fem:varform:fenics:problem:GN} \\tag{111}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variational formulation\n",
    "<div id=\"fem:varform:fenics:varform\"></div>\n",
    "\n",
    "To obtain the variational formulation, we multiply the PDE by a test\n",
    "function $v$ and integrate the second-order derivatives by part:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\int_\\Omega \\nabla\\cdot ({\\alpha}\\nabla u) v {\\, \\mathrm{d}x} \n",
    "&= -\\int_\\Omega {\\alpha}\\nabla u\\cdot\\nabla v{\\, \\mathrm{d}x} + \\int_{\\Gamma_N}{\\alpha}\n",
    "\\frac{\\partial u}{\\partial n}v{\\, \\mathrm{d}s}\\\\ \n",
    "&= -\\int_\\Omega {\\alpha}\\nabla u\\cdot\\nabla v{\\, \\mathrm{d}x},\\quad \\forall v\\in V\\\\ \n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are left with a problem of the form: find $u$ such that\n",
    "$a(u,v)=L(v)\\ \\forall v\\in V$, with"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto43\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "a(u,v) = \\int_\\Omega {\\alpha}\\nabla u\\cdot\\nabla v{\\, \\mathrm{d}x},\n",
    "\\label{_auto43} \\tag{112}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto44\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "L(v) = \\int_\\Omega 0v{\\, \\mathrm{d}x} {\\thinspace .}\n",
    "\\label{_auto44} \\tag{113}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We write the integrand as $0v{\\, \\mathrm{d}x}$ even though $L=0$, because it is necessary\n",
    "in FEniCS to specify $L$ as a linear form (i.e., a\n",
    "test function and some form of integration need to be present) and not the number zero.\n",
    "The Dirichlet conditions make a nonzero solution.\n",
    "\n",
    "## The FEniCS solver\n",
    "\n",
    "We suppose that we have a function `make_mesh` that can make the mesh for us. More\n",
    "details about this function will be provided later. \n",
    "A next step is then to define proper Dirichlet conditions.\n",
    "This might seem a bit complicated, but we introduce *markers* at the\n",
    "boundary for marking the Dirichlet boundaries. The inner boundary has\n",
    "marker 1 and the outer has marker 2. In this way, we can recognize the\n",
    "nodes that are on the boundary. It is usually a part of the mesh making\n",
    "process to compute both the mesh and its markers, so `make_mesh` returns\n",
    "a `Mesh` object as `mesh` and a `MeshFunction` object `markers`.\n",
    "Setting Dirichlet conditions in the solver is then a matter of\n",
    "introducing `DirichletBC` objects, one for each part of the boundary\n",
    "marked by `markers`, and then we collect all such Dirichlet conditions\n",
    "in a list that is used by the assembly process to incorporate the\n",
    "Dirichlet  conditions in the linear system.\n",
    "The code goes like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = FunctionSpace(mesh, 'P', degree)\n",
    "bc_inner = DirichletBC(V, u_a, markers, 1)\n",
    "bc_outer = DirichletBC(V, u_b, markers, 2)\n",
    "bcs = [bc_inner, bc_outer]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, `u_a` and `u_b` are constants (floats) set by the user. In general, \n",
    "anything that can be evaluated pointwise can be used, such as `Expression`, \n",
    "`Function`, and `Constant`. \n",
    "The next step is to define the variational problem and solve it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define variational problem\n",
    "u = TrialFunction(V)\n",
    "v = TestFunction(V)\n",
    "a = alpha*dot(grad(u), grad(v))*dx\n",
    "L = Constant(0)*v*dx  # L = 0*v*dx = 0 does not work...\n",
    "\n",
    "# Compute solution\n",
    "u = Function(V)\n",
    "solve(a == L, u, bcs)\n",
    "\n",
    "f = File(\"mesh.xml\")\n",
    "f << mesh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to avoid `L=0` (`L` equal to the float zero), we have to\n",
    "tell FEniCS that is a linear form, so zero must be specified as `Constant(0)`.\n",
    "\n",
    "Note that everything is the same as for the approximation problem in the section [fe:approx:fenics](#fe:approx:fenics),\n",
    "except for the Dirichlet conditions and the formulas for `a` and\n",
    "`L`. FEniCS has, of course, access to very efficient solution methods,\n",
    "so we could add arguments to the `solve` call to apply\n",
    "state-of-the-art iterative methods and preconditioners for large-scale\n",
    "problems. However, for this little 2D case a standard sparse Gaussian\n",
    "elimination, as implied by `solve(a = L, u, bcs)` is a sufficient approach.\n",
    "\n",
    "Finally, we can save the solution to file for using professional\n",
    "visualization software and, if desired, add a quick plotting using the\n",
    "built-in FEniCS tool `plot`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save solution to file in VTK format\n",
    "vtkfile = File(filename + '.pvd')\n",
    "vtkfile << u\n",
    "\n",
    "u.rename('u', 'u'); plot(u); plot(mesh)\n",
    "import matplotlib.pyplot as plt\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(The `u.rename` call is just for getting a more readable title in the plot.)\n",
    "\n",
    "The above statements are collected in a function `solver` in the\n",
    "file [`borehole_fenics.py`](${fem_src}/borehole_fenics.py):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solver(\n",
    "    mesh,\n",
    "    markers,  # MeshFunctions for Dirichlet conditions\n",
    "    alpha,    # Diffusion coefficient\n",
    "    u_a,      # Inner pressure\n",
    "    u_b,      # Outer pressure\n",
    "    degree,   # Element polynomial degree\n",
    "    filename, # Name of VTK file\n",
    "    ):\n",
    "    V = FunctionSpace(mesh, 'P', degree)\n",
    "    bc_inner = DirichletBC(V, u_a, markers, 1)\n",
    "    bc_outer = DirichletBC(V, u_b, markers, 2)\n",
    "    bcs = [bc_inner, bc_outer]\n",
    "\n",
    "    # Define variational problem\n",
    "    u = TrialFunction(V)\n",
    "    v = TestFunction(V)\n",
    "    a = alpha*dot(grad(u), grad(v))*dx\n",
    "    L = Constant(0)*v*dx  # L = 0*v*dx = 0 does not work...\n",
    "\n",
    "    # Compute solution\n",
    "    u = Function(V)\n",
    "    solve(a == L, u, bcs)\n",
    "\n",
    "    f = File(\"mesh.xml\")\n",
    "    f << mesh\n",
    "\n",
    "    # Save solution to file in VTK format\n",
    "    vtkfile = File(filename + '.pvd')\n",
    "    vtkfile << u\n",
    "\n",
    "    u.rename('u', 'u'); plot(u); plot(mesh)\n",
    "    import matplotlib.pyplot as plt\n",
    "    plt.show()\n",
    "    return u\n",
    "\n",
    "def problem():\n",
    "    mesh, markers = make_mesh(Theta=25*pi/180, a=1, b=2,\n",
    "                              nr=20, nt=20, s=1.9)\n",
    "    beta = 5\n",
    "    solver(mesh, markers, alpha=1, u_a=1, u_b=0, degree=1, filename='tmp')\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    problem()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Be careful with name clashes!**\n",
    "\n",
    "It is easy when coding mathematics to use variable names that correspond\n",
    "to one-letter names in the mathematics. For example, in the mathematics\n",
    "of this problem there are two $a$ variables: the radius of the inner\n",
    "boundary and the bilinear form in the variational formulation.\n",
    "Using `a` for the inner boundary in `solver` does not work: it is\n",
    "quickly overwritten by the bilinear form. We therefore have to introduce\n",
    "`x_a`. Long variable names are to be preferred for safe programming,\n",
    "though short names corresponding to the mathematics are often nicer.\n",
    "\n",
    "\n",
    "\n",
    "## Making the mesh\n",
    "\n",
    "The hardest part of a finite element problem is very often to make the\n",
    "mesh.  This is particularly the case in large industrial projects, but\n",
    "also often academic projects quickly lead to time-consuming work with\n",
    "constructing finite element meshes. In the present example we create\n",
    "the mesh for the symmetric problem by deforming an originally\n",
    "rectangular mesh.  The rectangular mesh is made by the FEniCS object\n",
    "`RectangleMesh` on $[a,b]\\times [0,1]$.  Therefore, we stretch the\n",
    "mesh towards the left before we bend the rectangle onto to \"a piece\n",
    "of cake\". [Figure](#fem:varform:fenics:mesh1:fig1) shows an example\n",
    "on the resulting mesh. The stretching gives us refinement in the\n",
    "radial direction because we expect the variations to be quite large in\n",
    "this direction, but uniform in $\\theta$ direction.\n",
    "\n",
    "We first make the rectangle and set boundary markers here for the inner\n",
    "and outer boundaries (since these are particularly simple: $x=a$ and $x=b$).\n",
    "Here is how we make\n",
    "the rectangular mesh from lower left corner $(a,0)$ to upper left corner\n",
    "$(b,1)$ with `nr` quadrilateral cells in $x$ direction (later to become\n",
    "the radial direction) and `nt` quadrilateral cells in the $y$ direction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = RectangleMesh(Point(a, 0), Point(b, 1), nr, nt, 'crossed')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each quadrilateral cell is divided into two triangles with right or left\n",
    "going diagonals, or four triangles using both diagonals. These choices\n",
    "of producing triangles from rectangles are named `right`, `left`, and\n",
    "`crossed`. Recall that FEniCS can only work with cells of triangular\n",
    "shape only and where the sides are straight. This means that we need a\n",
    "good resolution in $\\theta$ direction to represent a circular boundary.\n",
    "With isoparametric elements, it is easier to get a higher-order polynomial\n",
    "approximation of curved boundaries.\n",
    "\n",
    "We must then mark the boundaries for boundary conditions. Since we do not need\n",
    "to do anything with the homogeneous Neumann conditions, we can just mark\n",
    "the inner and outer boundary of the hole cylinder. This is very easy to do\n",
    "as long as the mesh is a rectangle since then the specifications of the\n",
    "boundaries are $x=a$ and $x=b$. The relevant FEniCS code requires the\n",
    "user to define a subclass of `SubDomain` and implement a function `inside`\n",
    "for indicating whether a given point `x` is on the desired boundary or not:\n",
    "\n",
    "<!-- dom:FIGURE: [fig/borehole_mesh1.png, width=600 frac=0.8] Finite element mesh for a porous medium outside a bore hole (hollow cylinder). <div id=\"fem:varform:fenics:mesh1:fig1\"></div> -->\n",
    "<!-- begin figure -->\n",
    "<div id=\"fem:varform:fenics:mesh1:fig1\"></div>\n",
    "\n",
    "<p>Finite element mesh for a porous medium outside a bore hole (hollow cylinder).</p>\n",
    "<img src=\"fig/borehole_mesh1.png\" width=600>\n",
    "\n",
    "<!-- end figure -->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# x=a becomes the inner borehole boundary\n",
    "class Inner(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        return on_boundary and abs(x[0] - a) < tol\n",
    "\n",
    "# x=b becomes the outer borehole boundary\n",
    "class Outer(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        return on_boundary and abs(x[0] - b) < tol\n",
    "\n",
    "inner = Inner(); outer = Outer();\n",
    "markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)\n",
    "markers.set_all(0)\n",
    "inner.mark(markers, 1)\n",
    "outer.mark(markers, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the instances `inner` and `outer` we fill a marker object, called\n",
    "`MeshFunction` in FEniCS. For this purpose we must introduce our own\n",
    "convention of numbering boundaries: here we use 1 for all points on\n",
    "the inner boundary and 2 for the outer boundary, while all other points\n",
    "are marked by 0. The solver applies the `markers` object to set\n",
    "the right Dirichlet boundary conditions.\n",
    "\n",
    "The next step is to deform the mesh. Given coordinates $x$, we can map\n",
    "these onto a stretched coordinate $\\bar x$ by"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto45\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\bar x = a + (b-a)\\left({x-a\\over b-a}\\right)^s,\n",
    "\\label{_auto45} \\tag{114}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $s$ is a parameter that controls the amount of stretching.\n",
    "The formula above gives a stretching towards $x=a$, while the next one\n",
    "stretches the coordinates towards $x=b$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto46\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\bar x = a + (b-a)\\left({x-a\\over b-a}\\right)^{1/s}{\\thinspace .}\n",
    "\\label{_auto46} \\tag{115}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code shown later\n",
    "shows the details of mapping coordinates in a FEniCS mesh object.\n",
    "\n",
    "The final step is to map the rectangle to a part of a hollow cylinder.\n",
    "Mathematically, a point $(x,y)$ in the rectangle is mapped onto $\\bar x,\n",
    "\\bar y)$ in our final geometry:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\hat x = \\bar x\\cos (\\Theta \\bar y),\\quad \\hat y = \\bar x\\sin (\\Theta \\bar y){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The relevant FEniCS code becomes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --- Deform mesh ---\n",
    "\n",
    "# First make a denser mesh towards r=a\n",
    "x = mesh.coordinates()[:,0]\n",
    "y = mesh.coordinates()[:,1]\n",
    "\n",
    "def denser(x, y):\n",
    "    return [a + (b-a)*((x-a)/(b-a))**s, y]\n",
    "\n",
    "x_bar, y_bar = denser(x, y)\n",
    "xy_bar_coor = np.array([x_bar, y_bar]).transpose()\n",
    "mesh.coordinates()[:] = xy_bar_coor\n",
    "\n",
    "# Then map onto to a \"piece of cake\"\n",
    "\n",
    "def cylinder(r, s):\n",
    "    return [r*np.cos(Theta*s), r*np.sin(Theta*s)]\n",
    "\n",
    "x_hat, y_hat = cylinder(x_bar, y_bar)\n",
    "xy_hat_coor = np.array([x_hat, y_hat]).transpose()\n",
    "mesh.coordinates()[:] = xy_hat_coor\n",
    "return mesh, markers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fortunately, the solver is independent of all the details of the mesh making.\n",
    "We could also have used the mesh tool `mshr` in FEniCS, but with our\n",
    "approach here we have full control of the refinement towards the hole.\n",
    "\n",
    "\n",
    "## Solving a problem\n",
    "\n",
    "We assume that ${\\alpha}$ is constant.  Before solving such a specific\n",
    "problem, it can be wise to scale the problem since it often reduces\n",
    "the amount of input data in the model. Here, the variation in $u$ is\n",
    "typically $|u_a-u_b|$ so we use that as characteristic pressure. The\n",
    "coordinates may be naturally scaled by the bore hole radius, so we\n",
    "have new, scaled variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\bar u = \\frac{u-u_a}{u_a-u_b},\\quad \\bar x = \\frac{x}{a},\\quad\n",
    "\\bar y = \\frac{y}{a}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we expect $\\bar u\\in [0,1]$, which is a goal of scaling.\n",
    "Inserting this in the problem gives the PDE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\nabla^2 \\bar u = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "in a domain with inner radius 1 and $\\bar u=0$, and outer radius"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\beta = \\frac{a}{b},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with $\\bar u = 1$. Our solver can solve this problem by setting\n",
    "`alpha=1`, `u_a=1`, and `u_b=0`. We see that the dimensionless\n",
    "parameter $\\beta$ goes to the mesh and not to the solver.\n",
    "[Figure](#fem:varform:fenics:sol1:fig2) shows a solution for $\\beta =2$\n",
    "on a mesh with $4\\cdot 20\\cdot 20 = 1600$ triangles, 25 degree opening,\n",
    "and P1 elements.\n",
    "Switching to higher-order, say P3, is a matter of changing the `degree`\n",
    "parameter that goes to the function `V` in the solver:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh, markers = make_mesh(Theta=25*pi/180, a=1, b=2,\n",
    "                          nr=20, nt=20, s=1.9)\n",
    "beta = 2\n",
    "solver(mesh, markers,\n",
    "       alpha=1, u_a=1, u_b=0, degree=3, filename='borehole1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The complete code is found in [`borehole_fenics.py`](${fem_src}/borehole_fenics.py). All fluids flow in the same way as long as the geometry is the same!\n",
    "\n",
    "<!-- dom:FIGURE: [fig/borehole1.png, width=800 frac=1] Solution for (scaled) fluid pressure around a bore hole in a porous medium. <div id=\"fem:varform:fenics:sol1:fig2\"></div> -->\n",
    "<!-- begin figure -->\n",
    "<div id=\"fem:varform:fenics:sol1:fig2\"></div>\n",
    "\n",
    "<p>Solution for (scaled) fluid pressure around a bore hole in a porous medium.</p>\n",
    "<img src=\"fig/borehole1.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "How can we solve a 3D version of this problem? Then we would make a long\n",
    "cylinder. The assumption is that nothing changes in the third direction, so\n",
    "$\\partial/\\partial z=0$. This means that the cross sections at the end of\n",
    "the cylinder have homogeneous Neumann conditions $\\partial u/\\partial n=0$.\n",
    "Therefore, nothing changes in the variational form.\n",
    "Actually, all we have to do is to a generate a 3D box and use the same\n",
    "stretching and mapping to make the cylinder, and run the solver without\n",
    "changes!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary\n",
    "\n",
    " * When approximating $f$ by $u = \\sum_j c_j{\\varphi}_j$, the least squares\n",
    "   method and the Galerkin/projection method give the same result.\n",
    "   The interpolation/collocation method is simpler and yields different\n",
    "   (mostly inferior) results.\n",
    "\n",
    " * Fourier series expansion can be viewed as a least squares or Galerkin\n",
    "   approximation procedure with sine and cosine functions.\n",
    "\n",
    " * Basis functions should optimally be orthogonal or almost orthogonal,\n",
    "   because this makes the coefficient matrix become diagonal or sparse and \n",
    "   results in little round-off errors when solving the linear\n",
    "   system.    One way to obtain almost orthogonal basis functions is to localize the\n",
    "   basis by making the basis functions that have local support in only a few\n",
    "   cells of a mesh. This is utilized in the finite element method.\n",
    "\n",
    " * Finite element basis functions are *piecewise* polynomials, normally\n",
    "   with discontinuous derivatives at the cell boundaries. The basis\n",
    "   functions overlap very little, leading to stable and efficient numerics involving sparse\n",
    "   matrices.\n",
    "\n",
    " * To use the finite element method for differential equations, we use\n",
    "   the Galerkin method or the method of weighted residuals\n",
    "   to arrive at a variational form. Technically, the differential equation\n",
    "   is multiplied by a test function and integrated over the domain.\n",
    "   Second-order derivatives are integrated by parts to allow for typical finite\n",
    "   element basis functions that have discontinuous derivatives.\n",
    "\n",
    " * The least squares method is not much used for finite element solution\n",
    "   of differential equations of second order, because\n",
    "   it then involves second-order derivatives which cause trouble for\n",
    "   basis functions with discontinuous derivatives.\n",
    "\n",
    " * We have worked with two common finite element terminologies and\n",
    "   associated data structures\n",
    "   (both are much used, especially the first one, while the other is more\n",
    "   general):\n",
    "\n",
    "a. *elements*, *nodes*, and *mapping between local and global\n",
    "       node numbers*\n",
    "\n",
    "b. an extended element concept consisting of *cell*, *vertices*,\n",
    "       *degrees of freedom*, *local basis functions*,\n",
    "       *geometry mapping*, and *mapping between\n",
    "       local and global degrees of freedom*\n",
    "\n",
    "\n",
    " * The meaning of the word \"element\" is multi-fold: the geometry of a finite\n",
    "   element (also known as a cell), the geometry and its basis functions,\n",
    "   or all information listed under point 2 above.\n",
    "\n",
    " * One normally computes integrals in the finite element method element\n",
    "   by element (cell by cell), either in a local reference coordinate\n",
    "   system or directly in the physical domain.\n",
    "\n",
    " * The advantage of working in the reference coordinate system is that\n",
    "   the mathematical expressions for the basis functions depend on the\n",
    "   element type only, not the geometry of that element in the physical\n",
    "   domain.  The disadvantage is that a mapping must be used, and\n",
    "   derivatives must be transformed from reference to physical\n",
    "   coordinates.\n",
    "\n",
    " * Element contributions to the global linear system are collected in\n",
    "   an element matrix and vector, which must be assembled into the\n",
    "   global system using the degree of freedom mapping (`dof_map`) or\n",
    "   the node numbering mapping (`elements`), depending on which terminology\n",
    "   that is used.\n",
    "\n",
    " * Dirichlet conditions, involving prescribed values of $u$ at the\n",
    "   boundary, are mathematically taken care of via a boundary function that takes\n",
    "   on the right Dirichlet values, while the basis functions vanish at\n",
    "   such boundaries. The finite element method features a general\n",
    "   expression for the boundary function. In software implementations,\n",
    "   it is easier to drop the boundary function and the requirement that\n",
    "   the basis functions must vanish on Dirichlet boundaries and instead\n",
    "   manipulate the global matrix system (or the element matrix and vector)\n",
    "   such that the Dirichlet values are imposed on the unknown parameters.\n",
    "\n",
    " * Neumann conditions, involving prescribed values of the derivative\n",
    "   (or flux) of $u$, are incorporated in boundary terms arising from\n",
    "   integrating terms with second-order derivatives by part.\n",
    "   Forgetting to account for the boundary terms implies the\n",
    "   condition $\\partial u/\\partial n=0$ at parts of the boundary where\n",
    "   no Dirichlet condition is set."
   ]
  }
 ],
 "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
}