<jupyter_start><jupyter_text>Variational formulations with finite elementsWe shall now take the ideas from previous chapters and put together such thatwe can solve PDEs using the flexible finite element basis functions.This is quite a machinery with many details, but the chapter is mostly anassembly of concepts and details we have already met. Computing with finite elementsThe purpose of this section is to demonstrate in detail howthe finite element method can then be applied to the model problem $$-u''(x) = 2,\quad x\in (0,L),\ u(0)=u(L)=0,$$ with variational formulation $$(u',v') = (2,v)\quad\forall v\in V{\thinspace .}$$ Any $v\in V$ must obey $v(0)=v(L)=0$ because of the Dirichlet conditionson $u$. Finite element mesh and basis functionsWe introduce a finite element mesh with $N_e$ cells, allwith length $h$, and numberthe cells from left to right.Choosing P1 elements, there are twonodes per cell, and the coordinates of the nodes become $$x_{i} = i h,\quad h=L/N_e,\quad i=0,\ldots,N_n-1=N_e{\thinspace .}$$ Any node $i$ is associated with a finite element basis function${\varphi}_i(x)$.  When approximating a given function $f$ by a finiteelement function $u$, we expand $u$ using finite element basisfunctions associated with *all* nodes in the mesh. The parameter$N$, which counts the unknowns from 0 to $N$, is then equal to$N_n-1$ such that the total number of unknowns, $N+1$, is thetotal number of nodes.However, when solving differential equations we will often have$N<N_n-1$ because of Dirichlet boundary conditions. The reason issimple: we know what $u$ are at some (here two) nodes, and the number ofunknown parameters is naturally reduced.In our case with homogeneous Dirichlet boundary conditions, we do notneed any boundary function $B(x)$, so we can work with the expansion $$\begin{equation}u(x) = \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j(x){\thinspace .}\label{fem:deq:1D:fem1:ex:u} \tag{64}\end{equation}$$ Because of the boundary conditions, we must demand${\psi}_i(0)={\psi}_i(L)=0$, $i\in{\mathcal{I}_s}$. When ${\psi}_i$ for all$i=0,\ldots,N$ is to be selected among the finite element basisfunctions ${\varphi}_j$, $j=0,\ldots,N_n-1$, we have to avoid using${\varphi}_j$ functions that do not vanish at $x_{0}=0$ and$x_{N_n-1}=L$. However, all ${\varphi}_j$ vanish at these two nodes for$j=1,\ldots,N_n-2$.  Only basis functions associated with the end nodes,${\varphi}_0$ and ${\varphi}_{N_n-1}$, violate the boundary conditions ofour differential equation. Therefore, we select the basis functions${\varphi}_i$ to be the set of finite element basis functions associatedwith all the interior nodes in the mesh: $${\psi}_i={\varphi}_{i+1},\quad i=0,\ldots,N{\thinspace .}$$ The $i$ index runs over all the unknowns $c_i$ in the expansionfor $u$, and in this case $N=N_n-3$.In the general case, and in particular on domains in higher dimensions,the nodes are not necessarily numbered from leftto right, so we introduce a mapping from the node numbering, or moreprecisely the degree of freedom numbering, to the numbering ofthe unknowns in the final equation system. These unknowns take onthe numbers $0,\ldots,N$. Unknown number $j$ in the linear systemcorresponds to degree of freedom number $\nu (j)$, $j\in{\mathcal{I}_s}$.We can then write $${\psi}_i={\varphi}_{\nu(i)},\quad i=0,\ldots,N{\thinspace .}$$ With a regular numbering as in the present example,$\nu(j) = j+1$, $j=0,\ldots,N=N_n-3$. Computation in the global physical domainWe shall first perform a computation in the $x$coordinate system because the integrals can be easily computedhere by simple, visual,geometric considerations. This is called a global approachsince we work in the $x$ coordinate system and compute integrals onthe global domain $[0,L]$.The entries in the coefficient matrix and right-hand side are $$A_{i,j}=\int_0^L{\psi}_i'(x){\psi}_j'(x) {\, \mathrm{d}x},\quadb_i=\int_0^L2{\psi}_i(x) {\, \mathrm{d}x}, \quad i,j\in{\mathcal{I}_s}{\thinspace .}$$ Expressed in terms of finite element basis functions ${\varphi}_i$ weget the alternative expressions $$A_{i,j}=\int_0^L{\varphi}_{i+1}'(x){\varphi}_{j+1}'(x) {\, \mathrm{d}x},\quadb_i=\int_0^L2{\varphi}_{i+1}(x) {\, \mathrm{d}x},\quad i,j\in{\mathcal{I}_s}{\thinspace .}$$ For the following calculations the subscripts on the finiteelement basis functions are more conveniently written as$i$ and $j$ instead of $i+1$ and $j+1$, so our notation becomes $$A_{i-1,j-1}=\int_0^L{\varphi}_{i}'(x){\varphi}_{j}'(x) {\, \mathrm{d}x},\quadb_{i-1}=\int_0^L2{\varphi}_{i}(x) {\, \mathrm{d}x},$$ where the $i$ and $j$ indices run as $i,j=1,\ldots,N+1=N_n-2$.The ${\varphi}_i(x)$ function is a hat function with peak at $x=x_{i}$and a linear variation in $[x_{i-1},x_{i}]$ and$[x_{i},x_{i+1}]$.The derivative is $1/h$ to the left of $x_{i}$ and $-1/h$ tothe right, or more formally, $$\begin{equation}{\varphi}_i'(x) = \left\lbrace\begin{array}{ll}0, & x < x_{i-1},\\ h^{-1},& x_{i-1} \leq x < x_{i},\\ -h^{-1},& x_{i} \leq x < x_{i+1},\\ 0, & x\geq x_{i+1}\end{array}\right.\label{fem:approx:fe:Dphi:1:formula2} \tag{65}\end{equation}$$ [Figure](fem:approx:fe:fig:dP1:2) shows ${\varphi}_2'(x)$ and ${\varphi}_3'(x)$. -->Illustration of the derivative of piecewise linear basis functions associated with nodes in cell 2.We realize that ${\varphi}_i'$ and ${\varphi}_j'$ has no overlap, and hence theirproduct vanishes, unless $i$ and $j$ are nodes belonging to the samecell. The only nonzero contributions to the coefficient matrix aretherefore $$\begin{align*}A_{i-1,i-2} &=\int_0^L{\varphi}_i'(x) {\varphi}_{i-1}'(x) {\, \mathrm{d}x},\\ A_{i-1,i-1}&=\int_0^L{\varphi}_{i}'(x)^2 {\, \mathrm{d}x}, \\ A_{i-1,i}&=\int_0^L{\varphi}_{i}'(x){\varphi}_{i+1}'(x) {\, \mathrm{d}x},\end{align*}$$ for $i=1,\ldots,N+1$, but for $i=1$, $A_{i-1,i-2}$ is not defined,and for $i=N+1$, $A_{i-1,i}$ is not defined.From [Figure](fem:approx:fe:fig:dP1:2),we see that ${\varphi}_{i-1}'(x)$ and ${\varphi}_i'(x)$ have overlap of onecell $\Omega^{(i-1)}=[x_{i-1},x_{i}]$ and that their productthen is $-1/h^{2}$. The integrand is constant and therefore$A_{i-1,i-2}=-h^{-2}h=-h^{-1}$.A similar reasoning can be applied to$A_{i-1,i}$, which also becomes $-h^{-1}$. The integral of${\varphi}_i'(x)^2$ gets contributions from two cells,$\Omega^{(i-1)}=[x_{i-1},x_{i}]$ and$\Omega^{(i)}=[x_{i},x_{i+1}]$, but ${\varphi}_i'(x)^2=h^{-2}$ inboth cells, and the length of the integration interval is $2h$ sowe get$A_{i-1,i-1}=2h^{-1}$.The right-hand side involves an integral of $2{\varphi}_i(x)$,$i=1,\ldots,N_n-2$,which is just the area under a hat function of height 1 and width$2h$, i.e., equal to $h$. Hence, $b_{i-1}=2h$.To summarize the linear system, we switch from $i$ to $i+1$ such thatwe can write $$A_{i,i-1}=A_{i,i+1}=-h^{-1},\quad A_{i,i}=2h^{-1},\quadb_i = 2h{\thinspace .}$$ The equation system to be solved only involves the unknowns$c_i$ for $i\in{\mathcal{I}_s}$. With our numbering of unknowns andnodes, we have that $c_i$ equals $u(x_{i+1})$.The complete matrix system then takes the following form: $$\begin{equation}\frac{1}{h}\left(\begin{array}{ccccccccc}1 & -1 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -1 & 2 & -1 & \ddots &   & &  & &  \vdots \\ 0 & -1 & 2 & -1 &\ddots & &  &  & \vdots \\ \vdots & \ddots &  & \ddots & \ddots & 0 &  & & \vdots \\ \vdots &  & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & &  & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & &  & \ddots & \ddots & \ddots &\ddots  & 0 \\ \vdots & & & &  &\ddots  & \ddots &\ddots  & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots  & 0 & -1 & 1\end{array}\right)\left(\begin{array}{c}c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N}\end{array}\right)=\left(\begin{array}{c}2h \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ 2h\end{array}\right)\label{fem:deq:1D:ex1:Ab:glob} \tag{66}\end{equation}$$  Comparison with a finite difference discretizationA typical row in the matrix system ([66](fem:deq:1D:ex1:Ab:glob))can be written as $$\begin{equation}-\frac{1}{h}c_{i-1} + \frac{2}{h}c_{i} - \frac{1}{h}c_{i+1} = 2h{\thinspace .}\label{fem:deq:1D:fem:ex1:c} \tag{67}\end{equation}$$ Let us introduce the notation $u_j$ for the value of $u$ at node $j$:$u_j=u(x_{j})$, since we have the interpretation$u(x_{j})=\sum_jc_j{\varphi}(x_{j})=\sum_j c_j\delta_{ij}=c_j$.The unknowns $c_0,\ldots,c_N$ are $u_1,\ldots,u_{N_n-2}$.Shifting $i$ with $i+1$ in ([67](fem:deq:1D:fem:ex1:c)) and inserting$u_i = c_{i-1}$, we get $$\begin{equation}-\frac{1}{h}u_{i-1} + \frac{2}{h}u_{i} - \frac{1}{h}u_{i+1} = 2h,\label{fem:deq:1D:fem:ex1} \tag{68}\end{equation}$$ A finite difference discretization of $-u''(x)=2$ by a centered,second-order finite difference approximation $u''(x_i)\approx [D_x D_x u]_i$with $\Delta x = h$yields $$\begin{equation}-\frac{u_{i-1} - 2u_{i} + u_{i+1}}{h^2} = 2,\label{_auto35} \tag{69}\end{equation}$$ which is, in fact, equal to ([68](fem:deq:1D:fem:ex1)) if([68](fem:deq:1D:fem:ex1)) is divided by $h$.Therefore, the finite difference and the finite element method areequivalent in this simple test problem.Sometimes a finite element method generates the finite differenceequations on a uniform mesh, and sometimes the finite element methodgenerates equations that are different.  The differences are modest,but may influence the numerical quality of the solution significantly,especially in time-dependent problems.It depends on the problem at handwhether a finite element discretization is more or less accurate thana corresponding finite difference discretization. Cellwise computationsSoftware for finite element computations normally employsthe cell by cell computational procedure wherean element matrix and vector are calculated for each cell andassembled in the global linear system.Let us go through the details of this type of algorithm.All integrals are mapped to the local reference coordinate system$X\in [-1,1]$.In the present case, the matrix entries contain derivativeswith respect to $x$, $$A_{i-1,j-1}^{(e)}=\int_{\Omega^{(e)}} {\varphi}_i'(x){\varphi}_j'(x) {\, \mathrm{d}x}= \int_{-1}^1 \frac{d}{dx}{\tilde{\varphi}}_r(X)\frac{d}{dx}{\tilde{\varphi}}_s(X)\frac{h}{2} {\, \mathrm{d}X},$$ where the global degree of freedom $i$ is related to the localdegree of freedom $r$ through $i=q(e,r)$. Similarly,$j=q(e,s)$. The local degrees of freedom run as $r,s=0,1$ for a P1element. The integral for the element matrixThere are simple formulas for the basis functions ${\tilde{\varphi}}_r(X)$ asfunctions of $X$.However, we nowneed to find the derivative of ${\tilde{\varphi}}_r(X)$ with respect to $x$.Given $${\tilde{\varphi}}_0(X)=\frac{1}{2}(1-X),\quad{\tilde{\varphi}}_1(X)=\frac{1}{2}(1+X),$$ we can easily compute $d{\tilde{\varphi}}_r/ dX$: $$\frac{d{\tilde{\varphi}}_0}{dX} = -\frac{1}{2},\quad  \frac{d{\tilde{\varphi}}_1}{dX} = \frac{1}{2}{\thinspace .}$$ From the chain rule, $$\begin{equation}\frac{d{\tilde{\varphi}}_r}{dx} = \frac{d{\tilde{\varphi}}_r}{dX}\frac{dX}{dx}= \frac{2}{h}\frac{d{\tilde{\varphi}}_r}{dX}{\thinspace .}   \label{_auto36} \tag{70}\end{equation}$$ The transformed integral is then $$A_{i-1,j-1}^{(e)}=\int_{\Omega^{(e)}} {\varphi}_i'(x){\varphi}_j'(x) {\, \mathrm{d}x}= \int_{-1}^1 \frac{2}{h}\frac{d{\tilde{\varphi}}_r}{dX}\frac{2}{h}\frac{d{\tilde{\varphi}}_s}{dX}\frac{h}{2} {\, \mathrm{d}X}{\thinspace .}$$  The integral for the element vectorThe right-hand side is transformed according to $$b_{i-1}^{(e)} = \int_{\Omega^{(e)}} 2{\varphi}_i(x) {\, \mathrm{d}x} =\int_{-1}^12{\tilde{\varphi}}_r(X)\frac{h}{2} {\, \mathrm{d}X},\quad i=q(e,r),\ r=0,1{\thinspace .}$$  Detailed calculations of the element matrix and vectorSpecifically for P1 elements we arrive at the following calculations forthe element matrix entries: $$\begin{align*}\tilde A_{0,0}^{(e)} &= \int_{-1}^1\frac{2}{h}\left(-\frac{1}{2}\right)\frac{2}{h}\left(-\frac{1}{2}\right)\frac{h}{2} {\, \mathrm{d}X} = \frac{1}{h}\\ \tilde A_{0,1}^{(e)} &= \int_{-1}^1\frac{2}{h}\left(-\frac{1}{2}\right)\frac{2}{h}\left(\frac{1}{2}\right)\frac{h}{2} {\, \mathrm{d}X} = -\frac{1}{h}\\ \tilde A_{1,0}^{(e)} &= \int_{-1}^1\frac{2}{h}\left(\frac{1}{2}\right)\frac{2}{h}\left(-\frac{1}{2}\right)\frac{h}{2} {\, \mathrm{d}X} = -\frac{1}{h}\\ \tilde A_{1,1}^{(e)} &= \int_{-1}^1\frac{2}{h}\left(\frac{1}{2}\right)\frac{2}{h}\left(\frac{1}{2}\right)\frac{h}{2} {\, \mathrm{d}X} = \frac{1}{h}\end{align*}$$ The element vector entries become $$\begin{align*}\tilde b_0^{(e)} &= \int_{-1}^12\frac{1}{2}(1-X)\frac{h}{2} {\, \mathrm{d}X} = h\\ \tilde b_1^{(e)} &= \int_{-1}^12\frac{1}{2}(1+X)\frac{h}{2} {\, \mathrm{d}X} = h{\thinspace .}\end{align*}$$ Expressing these entries in matrix and vector notation, we have $$\begin{equation}\tilde A^{(e)} =\frac{1}{h}\left(\begin{array}{rr}1 & -1\\ -1 & 1\end{array}\right),\quad\tilde b^{(e)} = h\left(\begin{array}{c}1\\ 1\end{array}\right){\thinspace .}\label{fem:deq:1D:ex1:Ab:elm} \tag{71}\end{equation}$$  Contributions from the first and last cellThe first and last cell involve only one unknown and one basis functionbecause of the Dirichlet boundary conditions at the first and lastnode.The element matrix therefore becomes a $1\times 1$ matrix and thereis only one entry in the element vector. On cell 0, only ${\psi}_0={\varphi}_1$is involved, corresponding to integration with ${\tilde{\varphi}}_1$. On cell $N_e$,only ${\psi}_N={\varphi}_{N_n-2}$ is involved, corresponding tointegration with ${\tilde{\varphi}}_0$.We then get the special end-cell contributions $$\begin{equation}\tilde A^{(e)} =\frac{1}{h}\left(\begin{array}{r}1\end{array}\right),\quad\tilde b^{(e)} = h\left(\begin{array}{c}1\end{array}\right),\label{fem:deq:1D:ex1:Ab:elm:ends} \tag{72}\end{equation}$$ for $e=0$ and $e=N_e$. In these cells, we have only one degree offreedom, not two as in the interior cells. AssemblyThe next step is to assemble the contributions from the various cells.The assembly of an element matrix and vector into the global matrixand right-hand side can be expressed as $$A_{q(e,r),q(e,s)} = A_{q(e,r),q(e,s)} + \tilde A^{(e)}_{r,s},\quadb_{q(e,r)} = b_{q(e,r)} + \tilde b^{(e)}_{r},\quad$$ for $r$ and $s$ running over all local degrees of freedom in cell $e$.To make the assembly algorithm more precise, it is convenient to set upPython data structures and a code snippet for carrying out all detailsof the algorithm.For a mesh of four equal-sized P1 elements and $L=2$ we have<jupyter_code>vertices = [0, 0.5, 1, 1.5, 2]
cells = [[0, 1], [1, 2], [2, 3], [3, 4]]
dof_map = [[0], [0, 1], [1, 2], [2]]<jupyter_output><empty_output><jupyter_text>The total number of degrees of freedom is 3, being the functionvalues at the internal 3 nodes where $u$ is unknown.In cell 0 we have global degree of freedom 0, the nextcell has $u$ unknown at its two nodes, which becomeglobal degrees of freedom 0 and 1, and so forth according tothe `dof_map` list. The mathematical $q(e,r)$ quantity is nothingbut the `dof_map` list.Assume all element matrices are stored in a list `Ae` such that`Ae[e][i,j]` is $\tilde A_{i,j}^{(e)}$. A corresponding listfor the element vectors is named `be`, where `be[e][r]` is$\tilde b_r^{(e)}$.A Python code snippetillustrates all details of the assembly algorithm: ```python A[i,j]: coefficient matrix, b[i]: right-hand sidefor e in range(len(Ae)):    for r in range(Ae[e].shape[0]):        for s in range(Ae[e].shape[1]):            A[dof_map[e,r],dof_map[e,s]] += Ae[e][i,j]        b[dof_map[e,r]] += be[e][i,j]``` The general case with `N_e` P1 elements of length `h` has ```pythonN_n = N_e + 1vertices = [i*h for i in range(N_n)]cells = [[e, e+1] for e in range(N_e)]dof_map = [[0]] + [[e-1, e] for i in range(1, N_e)] + [[N_n-2]]``` Carrying out the assembly results in a linear system that is identicalto ([66](fem:deq:1D:ex1:Ab:glob)), which is not surprising, sincethe procedures is mathematically equivalent to the calculationsin the physical domain.So far, our technique for computing the matrix system have assumedthat $u(0)=u(L)=0$. The next section deals with the extension tononzero Dirichlet conditions. Boundary conditions: specified nonzero valueWe have to take special actions to incorporate nonzeroDirichlet conditions,such as $u(L)=D$, into the computational procedures. The presentsection outlines alternative, yet mathematically equivalent, methods. General construction of a boundary functionIn previous sections, we introduced a boundary function $B(x)$to deal with nonzero Dirichlet boundary conditions for $u$. Theconstruction of such a function is not always trivial, especially notin multiple dimensions. However, a simple and general constructive ideaexists when thebasis functions have the property $${\varphi}_i(x_{j}) = \delta_{ij},\quad\delta_{ij} = \left\lbrace\begin{array}{ll}1, & i=j,\\ 0, & i\neq j,\end{array}\right.$$ where $x_{j}$ is a boundary point. Examples on suchfunctions are the Lagrange interpolating polynomials and finiteelement functions.Suppose now that $u$ has Dirichlet boundary conditions at nodeswith numbers $i\in{I_b}$. For example, ${I_b} = \{0,N_n-1\}$ in a 1Dmesh with node numbering from left to right and Dirichlet conditionsat the end nodes $i=0$ and $i=N_n-1$.Let $U_i$ be the corresponding prescribed values of $u(x_{i})$.We can then, in general, use $$\begin{equation}B(x) = \sum_{j\in{I_b}} U_j{\varphi}_j(x){\thinspace .}\label{_auto37} \tag{73}\end{equation}$$ It is easy to verify that$B(x_{i})= \sum_{j\in{I_b}} U_j{\varphi}_j(x_{i})=U_i$.The unknown function can then be written as $$\begin{equation}u(x) = \sum_{j\in{I_b}} U_j{\varphi}_j(x) + \sum_{j\in{\mathcal{I}_s}}c_j{\varphi}_{\nu(j)},\label{_auto38} \tag{74}\end{equation}$$ where $\nu(j)$ maps unknown number $j$ in the equation system tonode $\nu(j)$, ${I_b}$ is the set of indices corresponding to basis functionsassociated with nodes where Dirichlet conditions apply, and ${\mathcal{I}_s}$ is theset of indices used to number the unknowns from zero to $N$.We can easily show that with this $u$, a Dirichletcondition $u(x_{k})=U_k$ is fulfilled: $$u(x_{k}) = \sum_{j\in{I_b}} U_j\underbrace{{\varphi}_j(x_k)}_{\neq 0\,\Leftrightarrow\,j=k} +\sum_{j\in{\mathcal{I}_s}} c_j\underbrace{{\varphi}_{\nu(j)}(x_{k})}_{=0,\ k\not\in{\mathcal{I}_s}}= U_k$$ Some examples will further clarify the notation. With a regularleft-to-right numbering of nodes in a mesh with P1 elements,and Dirichlet conditions at $x=0$, we use finite element basisfunctions associated with the nodes $1, 2, \ldots, N_n-1$, implyingthat  $\nu(j)=j+1$, $j=0,\ldots,N$, where $N=N_n-2$. Considera particular mesh:The expansion associated with this mesh becomes $$u(x) = U_0{\varphi}_0(x) + c_0{\varphi}_1(x) +c_1{\varphi}_2(x) + \cdots + c_4{\varphi}_5(x){\thinspace .}$$ Switching to the more standard case of left-to-right numbering andboundary conditions $u(0)=C$, $u(L)=D$, we have $N=N_n-3$ and $$\begin{align*}u(x) &= C{\varphi}_0 + D{\varphi}_{N_n-1} + \sum_{j\in{\mathcal{I}_s}} c_j{\varphi}_{j+1}\\ &= C{\varphi}_0 + D{\varphi}_{N_n} + c_0{\varphi}_1 + c_1{\varphi}_2 +\cdots+ c_N{\varphi}_{N_n-2}{\thinspace .}\end{align*}$$ Finite element meshes in non-trivial 2D and 3D geometries usually leadsto an irregular cell and node numbering. Let us therefore take a lookat an irregular numbering in 1D:Say we in this mesh have Dirichlet conditions on the left-most andright-most node, with numbers 3 and 1, respectively.  We can numberthe unknowns at the interior nodes as we want, e.g., from left toright, resulting in $\nu(0)=0$, $\nu(1)=4$, $\nu(2)=5$, $\nu(3)=2$.This gives $$B(x) = U_3{\varphi}_3(x) + U_1{\varphi}_1(x),$$ and $$u(x) = B(x) + \sum_{j=0}^3 c_j{\varphi}_{\nu(j)}= U_3{\varphi}_3 + U_1{\varphi}_1 + c_0{\varphi}_0 + c_1{\varphi}_4+ c_2{\varphi}_5 + c_3{\varphi}_2{\thinspace .}$$ The idea of constructing $B$, described here, generalizes almosttrivially to 2D and 3D problems: $B=\sum_{j\in{I_b}}U_j{\varphi}_j$,where ${I_b}$ is the index set containing the numbers of all thenodes on the boundaries where Dirichlet values are prescribed. Example on computing with a finite element-based boundary functionLet us see how the model problem $-u''=2$, $u(0)=C$, $u(L)=D$,is affected by a $B(x)$ to incorporate boundary values.Inserting the expression $$u(x) = B(x) + \sum_{j\in{\mathcal{I}_s}}c_j{\psi}_j(x)$$ in $-(u'',{\psi}_i)=(f,{\psi}_i)$ andintegrating by parts results in a linear system with $$A_{i,j} = \int_0^L {\psi}_i'(x){\psi}_j'(x) {\, \mathrm{d}x},\quadb_i = \int_0^L (f(x){\psi}_i(x) - B'(x){\psi}_i'(x)) {\, \mathrm{d}x}{\thinspace .}$$ We choose ${\psi}_i={\varphi}_{i+1}$, $i=0,\ldots,N=N_n-3$if the node numbering is from leftto right. (Later we also need the assumption that cells tooare numbered from left to right.)The boundary function becomes $$B(x) = C{\varphi}_0(x) + D{\varphi}_{N_n-1}(x){\thinspace .}$$ The expansion for $u(x)$ is $$u(x)  = B(x) + \sum_{j\in{\mathcal{I}_s}} c_j{\varphi}_{j+1}(x){\thinspace .}$$ We can write the matrix and right-hand side entries as $$\begin{align*}A_{i-1,j-1} &= \int_0^L {\varphi}_i'(x){\varphi}_j'(x) {\, \mathrm{d}x},\\ b_{i-1} &=\int_0^L (f(x){\varphi}_i'(x) - (C{\varphi}_{0}'(x) + D{\varphi}_{N_n-1}'(x)){\varphi}_i'(x) ){\, \mathrm{d}x},\end{align*}$$ for $i,j = 1,\ldots,N+1=N_n-2$. Note that we have here used$B'=C{\varphi}_0' + D{\varphi}_{N_n-1}'$. Computations in physical coordinatesMost of the terms in the linear system have already been computedso we concentrate on the new contribution from the boundary function.The integral $C\int_0^L {\varphi}_{0}'(x)){\varphi}_i'(x) {\, \mathrm{d}x}$, associatedwith the Dirichlet condition in $x=0$,  can only geta nonzero contribution from the first cell,$\Omega^{(0)}=[x_{0},x_{1}]$since ${\varphi}_{0}'(x)=0$ on all other cells. Moreover,${\varphi}_{0}'(x){\varphi}_i'(x) {\, \mathrm{d}x} \neq 0$ only for $i=0$ and $i=1$(but node $i=0$ is excluded from the formulation),since ${\varphi}_{i}=0$ on the first cell if $i>1$.With a similar reasoning we realize that$D\int_0^L {\varphi}_{N_n-1}'(x)){\varphi}_i'(x) {\, \mathrm{d}x}$ can only geta nonzero contribution from the last cell. $$\begin{align*}\int_0^L {\varphi}_{0}'(x){\varphi}_{1}'(x) {\, \mathrm{d}x} &=(-\frac{1}{h})\cdot\frac{1}{h}\cdot h = -\frac{1}{h},\\ \int_0^L {\varphi}_{N_n-1}'(x){\varphi}_{N_n-2}'(x) {\, \mathrm{d}x} &=\frac{1}{h}\cdot(-\frac{1}{h})\cdot h = -\frac{1}{h}{\thinspace .}\end{align*}$$ With these expressions we get $$b_0 = \int_0^Lf(x){\varphi}_1{\, \mathrm{d}x} - C(-\frac{1}{h}),\quadb_N = \int_0^L f(x){\varphi}_{N_n-2}{\, \mathrm{d}x} - D(-\frac{1}{h}){\thinspace .}$$  Cellwise computations on the reference elementAs an equivalent alternative, we now turn to cellwise computations.The element matrices and vectors are calculated as in the section [Cellwise computations](fem:deq:1D:comp:elmwise), so we concentrate on the impact ofthe new term involving $B(x)$. This new term,$B'=C{\varphi}_0' + D{\varphi}_{N_n-1}'$, vanishes on all cells exceptfor $e=0$ and $e=N_e$. Over the first cell ($e=0$) the $B'(x)$ functionin local coordinates reads $$\frac{dB}{dx} = C\frac{2}{h}\frac{d{\tilde{\varphi}}_0}{dX},$$ while over the last cell ($e=N_e$) it looks like $$\frac{dB}{dx} = D\frac{2}{h}\frac{d{\tilde{\varphi}}_1}{dX}{\thinspace .}$$ For an arbitrary interior cell, we have the formula $$\tilde b_r^{(e)} = \int_{-1}^1 f(x(X)){\tilde{\varphi}}_r(X)\frac{h}{2}{\, \mathrm{d}X},$$ for an entry in the local element vector.In the first cell, the value at local node 0 is known so only the valueat local node 1 is unknown. The associated element vector entry becomes $$\begin{align*}\tilde b_0^{(1)} &= \int_{-1}^1 \left(f{\tilde{\varphi}}_1 -C\frac{2}{h}\frac{d{\tilde{\varphi}}_0}{dX}\frac{2}{h}\frac{d{\tilde{\varphi}}_1}{dX}\right)\frac{h}{2} {\, \mathrm{d}X} \\ & = \frac{h}{2} 2\int_{-1}^1 {\tilde{\varphi}}_1  {\, \mathrm{d}X}- C\frac{2}{h}(-\frac{1}{2})\frac{2}{h}\frac{1}{2}\frac{h}{2}\cdot 2= h + C\frac{1}{h}{\thinspace .}\end{align*}$$ The value at local node 1 in the last cell is known so theelement vector here is $$\begin{align*}\tilde b_0^{N_e} &= \int_{-1}^1 \left(f{\tilde{\varphi}}_0 -D\frac{2}{h}\frac{d{\tilde{\varphi}}_1}{dX}\frac{2}{h}\frac{d{\tilde{\varphi}}_0}{dX}\right)\frac{h}{2} {\, \mathrm{d}X}\\ & = \frac{h}{2} 2\int_{-1}^1 {\tilde{\varphi}}_0  {\, \mathrm{d}X}- D\frac{2}{h}\frac{1}{2}\frac{2}{h}(-\frac{1}{2})\frac{h}{2}\cdot 2= h + D\frac{1}{h}{\thinspace .}\end{align*}$$ The contributions from the $B(x)$ function to the global right-hand sidevector becomes $C/h$ for $b_0$ and $D/h$ for $b_N$, exactly as wecomputed in the physical domain. Modification of the linear systemFrom an implementational point of view, there is a convenient alternativeto adding the $B(x)$ function and using only the basis functions associatedwith nodes where $u$ is truly unknown.Instead of seeking $$\begin{equation}u(x) = \sum_{j\in{I_b}} U_j{\varphi}_j(x)+ \sum_{j\in{\mathcal{I}_s}}c_j{\varphi}_{\nu(j)}(x),\label{fem:deq:1D:fem:essBC:Bfunc:modsys:utrad} \tag{75}\end{equation}$$ we use the sum over all degrees of freedom, including the known boundaryvalues: $$\begin{equation}u(x) = \sum_{j\in{\mathcal{I}_s}}c_j{\varphi}_j(x){\thinspace .}\label{fem:deq:1D:fem:essBC:Bfunc:modsys:uall} \tag{76}\end{equation}$$ Note that the collections of unknowns$\left\{ {c}_i \right\}_{i\in{\mathcal{I}_s}}$ in ([75](fem:deq:1D:fem:essBC:Bfunc:modsys:utrad))and ([76](fem:deq:1D:fem:essBC:Bfunc:modsys:uall)) are different.The index set ${\mathcal{I}_s}=\{0,\ldots,N\}$ always goes to $N$, and thenumber of unknowns is $N+1$, butin ([75](fem:deq:1D:fem:essBC:Bfunc:modsys:utrad)) the unknownscorrespond to nodes where $u$ is not known, whilein ([76](fem:deq:1D:fem:essBC:Bfunc:modsys:uall)) the unknownscover $u$ values at all the nodes. So, if the index set${I_b}$ contains $N_b$ node numbers where $u$ is prescribed,we have that $N=N_n-N_b$ in([75](fem:deq:1D:fem:essBC:Bfunc:modsys:utrad)) and$N=N_n$ in ([76](fem:deq:1D:fem:essBC:Bfunc:modsys:uall)).The idea is to compute the entries in the linear system as if noDirichlet values are prescribed. Afterwards, we modify the linear systemto ensure that the known $c_j$ values are incorporated.A potential problem arises for the boundary term $[u'v]_0^L$ from theintegration by parts: imagining no Dirichlet conditions means that weno longer require $v=0$ at Dirichlet points, and the boundary term isthen nonzero at these points. However, when we modify the linearsystem, we will erase whatever the contribution from $[u'v]_0^L$should be at the Dirichlet points in the right-hand side of the linearsystem. We can therefore safely forget $[u'v]_0^L$ at any point wherea Dirichlet condition applies. Computations in the physical systemLet us redo the computations in the example inthe section [General construction of a boundary function](fem:deq:1D:fem:essBC:Bfunc). We solve $-u''=2$ with$u(0)=0$ and $u(L)=D$. The expressions for $A_{i,j}$ and $b_i$are the same, but the numbering is different as the numbering ofunknowns and nodes now coincide: $$A_{i,j} = \int_0^L {\varphi}_i'(x){\varphi}_j'(x) {\, \mathrm{d}x},\quadb_{i} = \int_0^L f(x){\varphi}_i(x) {\, \mathrm{d}x},$$ for $i,j = 0,\ldots,N=N_n-1$.The integrals involving basis functionscorresponding to interior mesh nodes, $i,j=1,\ldots,N_n-2$, areobviously the same as before. We concentrate on the contributionsfrom ${\varphi}_0$ and ${\varphi}_{N_n-1}$: $$\begin{align*}A_{0,0} &= \int_0^L ({\varphi}_0')^2{\, \mathrm{d}x} = \int_{0}^{x_{1}}= ({\varphi}_0')^2{\, \mathrm{d}x} \frac{1}{h},\\ A_{0,1} &= \int_0^L {\varphi}_0'{\varphi}_1'{\, \mathrm{d}x}= \int_{0}^{x_{1}} {\varphi}_0'{\varphi}_1'{\, \mathrm{d}x} = -\frac{1}{h},\\ A_{N,N} &= \int_0^L ({\varphi}_N')^2{\, \mathrm{d}x}= \int_{x_{N_n-2}}^{x_{N_n-1}} ({\varphi}_N')^2{\, \mathrm{d}x} = \frac{1}{h},\\ A_{N,N-1} &= \int_0^L {\varphi}_N'{\varphi}_{N-1}'{\, \mathrm{d}x}=\int_{x_{N_n-2}}^{x_{N_n-1}} {\varphi}_N'{\varphi}_{N-1}'{\, \mathrm{d}x} = -\frac{1}{h}{\thinspace .}\end{align*}$$ The new terms on the right-hand side are also those involving${\varphi}_0$ and ${\varphi}_{N_n-1}$: $$\begin{align*}b_0 &= \int_0^L 2{\varphi}_0(x) {\, \mathrm{d}x} = \int_0^{x_{1}} 2{\varphi}_0(x){\, \mathrm{d}x} = h,\\ b_N &=  \int_0^L 2{\varphi}_{N_n-1}{\, \mathrm{d}x} =\int_{x_{N_n-2}}^{x_{N_n-1}} 2{\varphi}_{N_n-1}{\, \mathrm{d}x} = h{\thinspace .}\end{align*}$$ The complete matrix system, involving all degrees of freedom, takes the form $$\begin{equation}\frac{1}{h}\left(\begin{array}{ccccccccc}1 & -1 & 0&\cdots &\cdots & \cdots & \cdots &\cdots & 0 \\ -1 & 2 & -1 & \ddots &   & &  & &  \vdots \\ 0 & -1 & 2 & -1 &\ddots & &  &  & \vdots \\ \vdots & \ddots &  & \ddots & \ddots & 0 &  & & \vdots \\ \vdots &  & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & &  & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & &  & \ddots & \ddots & \ddots &\ddots  & 0 \\ \vdots & & & &  &\ddots  & \ddots &\ddots  & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots  & 0 & -1 & 1\end{array}\right)\left(\begin{array}{c}c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N}\end{array}\right)=\left(\begin{array}{c}h \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h\\ h\end{array}\right)\label{fem:deq:1D:ex1:Ab:glob2} \tag{77}\end{equation}$$ Incorporation of Dirichlet values can now be done by replacingthe first and last equation bythe very simple equations $c_0=0$ and $c_N=D$, respectively.Note that the factor $1/h$ in front of the matrix then requiresa factor $h$ to be introduce appropriately on the diagonal in the first and lastrow of the matrix. $$\begin{equation}\frac{1}{h}\left(\begin{array}{ccccccccc}h & 0 & 0&\cdots &\cdots & \cdots & \cdots &\cdots & 0 \\ -1 & 2 & -1 & \ddots &   & &  & &  \vdots \\ 0 & -1 & 2 & -1 &\ddots & &  &  & \vdots \\ \vdots & \ddots &  & \ddots & \ddots & 0 &  & & \vdots \\ \vdots &  & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & &  & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & &  & \ddots & \ddots & \ddots &\ddots  & 0 \\ \vdots & & & &  &\ddots  & \ddots &\ddots  & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots  & 0 & 0 & h\end{array}\right)\left(\begin{array}{c}c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N}\end{array}\right)=\left(\begin{array}{c}0 \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h\\ D\end{array}\right)\label{fem:deq:1D:ex1:Ab:glob3} \tag{78}\end{equation}$$ Note that because we do not require ${\varphi}_i(0)=0$ and${\varphi}_i(L)=0$, $i\in{\mathcal{I}_s}$, the boundary term $[u'v]_0^L$,in principle, givescontributions $u'(0){\varphi}_0(0)$ to $b_0$ and$u'(L){\varphi}_N(L)$ to $b_N$ ($u'{\varphi}_i$ vanishes for $x=0$ or$x=L$ for $i=1,\ldots,N-1$).  Nevertheless, we erase thesecontributions in $b_0$ and $b_N$ and insert boundary values instead. Thisargument shows why we can drop computing $[u'v]_0^L$ at Dirichletnodes when we implement the Dirichlet values by modifying the linearsystem. Symmetric modification of the linear systemThe original matrix system ([66](fem:deq:1D:ex1:Ab:glob)) is symmetric,but the modifications in ([78](fem:deq:1D:ex1:Ab:glob3)) destroy thissymmetry. Our described modification will in general destroy aninitial symmetry in the matrix system. This is not a particularcomputational disadvantage for tridiagonal systems arising in 1Dproblems, but may be more serious in 2D and 3D problems when thesystems are large and exploiting symmetry can be important for halvingthe storage demands and speeding up computations. Methods for solvingsymmetric matrix are also usually more stable and efficient than thosefor non-symmetric systems.  Therefore, an alternative modificationwhich preserves symmetry is attractive.One can formulate a general algorithm for incorporating a Dirichletcondition in a symmetric way.Let $c_k$ be a coefficient corresponding to a known value$u(x_{k}) = U_k$.We want to replace equation $k$ in the system by $c_k=U_k$, i.e.,insert zeroes in row number $k$ in the coefficient matrix,set 1 on the diagonal, and replace $b_k$ by $U_k$.A symmetry-preserving modification consists in firstsubtracting column number $k$ in the coefficient matrix, i.e., $A_{i,k}$for $i\in{\mathcal{I}_s}$, times the boundary value $U_k$, from theright-hand side: $b_i \leftarrow b_i - A_{i,k}U_k$,$i=0,\ldots,N$. Then we putzeroes in both row number $k$ *and* column number $k$ in the coefficient matrix,and finally set $b_k=U_k$. The steps in algorithmic form becomes1. $b_i \leftarrow b_i - A_{i,k}U_k$ for $i\in{\mathcal{I}_s}$2. $A_{i,k} = A_{k,i} = 0$ for $i\in{\mathcal{I}_s}$3. $A_{k,k}=1$4. $b_i = U_k$This modification goes as follows for the specific linear systemwritten out in ([77](fem:deq:1D:ex1:Ab:glob2)) inthe section [Modification of the linear system](fem:deq:1D:fem:essBC:Bfunc:modsys). First wesubtract the first column in the coefficient matrix, times the boundaryvalue, from the right-hand side. Because $c_0=0$, this subtractionhas no effect. Then we subtract the last column, times the boundary value $D$,from the right-hand side. This action results in $b_{N-1}=2h+D/h$ and$b_N=h-2D/h$. Thereafter, we place zeros in the first and last row andcolumn in the coefficient matrix and 1 on the two corresponding diagonalentries. Finally, we set $b_0=0$ and $b_N=D$. The result becomes $$\begin{equation}\frac{1}{h}\left(\begin{array}{ccccccccc}h & 0 & 0&\cdots &\cdots & \cdots & \cdots &\cdots & 0 \\ 0 & 2 & -1 & \ddots &   & &  & &  \vdots \\ 0 & -1 & 2 & -1 &\ddots & &  &  & \vdots \\ \vdots & \ddots &  & \ddots & \ddots & 0 &  & & \vdots \\ \vdots &  & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & &  & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & &  & \ddots & \ddots & \ddots &\ddots  & 0 \\ \vdots & & & &  &\ddots  & \ddots &\ddots  & 0 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots  & 0 & 0 & h\end{array}\right)\left(\begin{array}{c}c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N}\end{array}\right)=\left(\begin{array}{c}0 \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h +D/h\\ D\end{array}\right)\label{fem:deq:1D:ex1:Ab:glob3:symm} \tag{79}\end{equation}$$  Modification of the element matrix and vectorThe modifications of the global linear system can alternativelybe done for the element matrix and vector. Let us perform theassociated calculations in the computational example wherethe element matrix and vector is given by([71](fem:deq:1D:ex1:Ab:elm)). The modifications are needed incells where one of the degrees of freedom is known. In thepresent example, this meansthe first and last cell. We compute the element matrixand vector as if there were no Dirichlet conditions. The boundary term$[u'v]_0^L$ is simply forgotten at nodes that have Dirichlet conditionsbecause the modification of the element vector will anyway erase thecontribution from the boundary term. In the first cell,local degree of freedom number 0is known and the modification becomes $$\begin{equation}\tilde A^{(0)} =A = \frac{1}{h}\left(\begin{array}{rr}h & 0\\ -1 & 1\end{array}\right),\quad\tilde b^{(0)} = \left(\begin{array}{c}0\\ h\end{array}\right){\thinspace .}\label{fem:deq:1D:ex1:Ab:elm:bc:0} \tag{80}\end{equation}$$ In the last cell we set $$\begin{equation}\tilde A^{(N_e)} =A = \frac{1}{h}\left(\begin{array}{rr}1 & -1\\ 0 & h\end{array}\right),\quad\tilde b^{(N_e)} = \left(\begin{array}{c}h\\ D\end{array}\right){\thinspace .}\label{fem:deq:1D:ex1:Ab:elm:bc:N} \tag{81}\end{equation}$$ We can also perform the symmetric modification. This operation affectsonly the last cell with a nonzero Dirichlet condition. The algorithmis the same as for the global linear system, resulting in $$\begin{equation}\tilde A^{(N_e)} =A = \frac{1}{h}\left(\begin{array}{rr}1 & 0\\ 0 & h\end{array}\right),\quad\tilde b^{(N_e)} = \left(\begin{array}{c}h + D/h\\ D\end{array}\right){\thinspace .}\label{fem:deq:1D:ex1:Ab:elm:bc:N:symm} \tag{82}\end{equation}$$ The reader is encouraged to assemble the element matrices and vectors andcheck that the result coincides with the system([79](fem:deq:1D:ex1:Ab:glob3:symm)). Boundary conditions: specified derivativeSuppose our model problem $-u''(x)=f(x)$ featuresthe boundary conditions $u'(0)=C$ and $u(L)=D$.As already indicated,the former condition can be incorporated through the boundary termthat arises from integration by parts. The details of this method will now beillustrated in the context of finite element basis functions. The variational formulationStarting with the Galerkin method, $$\int_0^L(u''(x)+f(x)){\psi}_i(x) {\, \mathrm{d}x} = 0,\quad i\in{\mathcal{I}_s},$$ integrating $u''{\psi}_i$ by parts results in $$\int_0^Lu'(x)'{\psi}_i'(x) {\, \mathrm{d}x} -(u'(L){\psi}_i(L) - u'(0){\psi}_i(0)) =\int_0^L f(x){\psi}_i(x) {\, \mathrm{d}x}, \quad i\in{\mathcal{I}_s}{\thinspace .}$$ The first boundary term, $u'(L){\psi}_i(L)$,vanishes because $u(L)=D$.The second boundaryterm, $u'(0){\psi}_i(0)$, can be used to implement the condition $u'(0)=C$,provided ${\psi}_i(0)\neq 0$ for some $i$ (but with finite elementswe fortunately have ${\psi}_0(0)=1$).The variational form of the differential equation then becomes $$\int_0^Lu'(x){\varphi}_i'(x) {\, \mathrm{d}x} + C{\varphi}_i(0) =\int_0^L f(x){\varphi}_i(x) {\, \mathrm{d}x},\quad i\in{\mathcal{I}_s}{\thinspace .}\label{fem:deq:1D:BC:nat:varform} \tag{83}$$  Boundary term vanishes because of the test functionsAt points where $u$ is known we may require ${\psi}_i$ to vanish.Here, $u(L)=D$ and then ${\psi}_i(L)=0$, $i\in{\mathcal{I}_s}$. Obviously,the boundary term $u'(L){\psi}_i(L)$ then vanishes.The set of basis functions $\left\{ {{\psi}}_i \right\}_{i\in{\mathcal{I}_s}}$ contains, in thiscase, all the finite element basis functions on the mesh, exceptthe one that is 1 at $x=L$. The basis function that is left out isused in a boundary function $B(x)$ instead.With a left-to-right numbering,${\psi}_i = {\varphi}_i$, $i=0,\ldots,N_n-2$, and $B(x)=D{\varphi}_{N_n-1}$: $$u(x) = D{\varphi}_{N_n-1}(x) + \sum_{j=0}^{N=N_n-2} c_j{\varphi}_j(x){\thinspace .}$$ Inserting this expansion for $u$ in the variational form([83](fem:deq:1D:BC:nat:varform)) leads to the linear system $$\begin{equation}\sum_{j=0}^{N}\left(\int_0^L {\varphi}_i'(x){\varphi}_j'(x) {\, \mathrm{d}x} \right)c_j =\int_0^L\left(f(x){\varphi}_i(x) -D{\varphi}_{N_n-1}'(x){\varphi}_i'(x)\right) {\, \mathrm{d}x} - C{\varphi}_i(0),\label{fem:deq:1D:natBC} \tag{84}\end{equation}$$ for $i=0,\ldots,N=N_n-2$. Boundary term vanishes because of linear system modificationsWe may, as an alternative to the approach in the previous section, usea basis $\left\{ {{\psi}}_i \right\}_{i\in{\mathcal{I}_s}}$ which contains all the finite elementfunctions on the mesh: ${\psi}_i={\varphi}_i$, $i=0,\ldots,N_n-1=N$.  Inthis case, $u'(L){\psi}_i(L)=u'(L){\varphi}_i(L)\neq 0$ for the $i$corresponding to the boundary node at $x=L$ (where ${\varphi}_i=1$).The number of this node is $i=N_n-1=N$ if a left-to-right numbering ofnodes is utilized.However, even though $u'(L){\varphi}_{N_n-1}(L)\neq 0$, we do not need tocompute this term.  For $i<N_n-1$ we realize that ${\varphi}_i(L)=0$.  Theonly nonzero contribution to the right-hand side comes from $i=N$($b_N$). Without a boundary function we must implement thecondition $u(L)=D$ by the equivalent statement $c_N=D$ and modify thelinear system accordingly. This modification will erase the lastrow and replace $b_N$ by another value. Any attempt to computethe boundary term $u'(L){\varphi}_{N_n-1}(L)$ and store it in $b_N$ will belost. Therefore, we can safely forget about boundary termscorresponding to Dirichlet boundary conditions also when we usethe methods from the section [Modification of the linear system](fem:deq:1D:fem:essBC:Bfunc:modsys)or the section [Symmetric modification of the linear system](fem:deq:1D:fem:essBC:Bfunc:modsys:symm).The expansion for $u$ reads $$u(x) = \sum_{j\in{\mathcal{I}_s}} c_j{\varphi}_j(x){\thinspace .}$$ Insertion in the variational form([83](fem:deq:1D:BC:nat:varform)) leads tothe linear system $$\begin{equation}\sum_{j\in{\mathcal{I}_s}}\left(\int_0^L {\varphi}_i'(x){\varphi}_j'(x) {\, \mathrm{d}x} \right)c_j =\int_0^L\left(f(x){\varphi}_i(x)\right) {\, \mathrm{d}x} - C{\varphi}_i(0),\quad i\in{\mathcal{I}_s}{\thinspace .}\label{fem:deq:1D:natBC2} \tag{85}\end{equation}$$ After having computed the system, we replace the last row by$c_N=D$, either straightforwardly as inthe section [Modification of the linear system](fem:deq:1D:fem:essBC:Bfunc:modsys) or in a symmetricfashion as in the section [Symmetric modification of the linear system](fem:deq:1D:fem:essBC:Bfunc:modsys:symm).These modifications can also be performed in the element matrix andvector for the right-most cell. Direct computation of the global linear systemWe now turn to actual computations with P1 finite elements.The focus is on how the linear system andthe element matrices and vectors are modified by thecondition $u'(0)=C$.Consider first the approach where Dirichlet conditions are incorporatedby a $B(x)$ function and the known degree of freedom$C_{N_n-1}$ is left out of the linear system(see the section [Boundary term vanishes because of the test functions](fem:deq:1D:BC:nat:uLtest)).The relevant formula for the linearsystem is given by ([84](fem:deq:1D:natBC)).There are three differences compared to the extensivelycomputed 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).First, because we do not have a Dirichletcondition at the left boundary, we need to extend the linear system([66](fem:deq:1D:ex1:Ab:glob)) with an equation associated with the node$x_{0}=0$.According to the section [Modification of the linear system](fem:deq:1D:fem:essBC:Bfunc:modsys), thisextension consists of including $A_{0,0}=1/h$, $A_{0,1}=-1/h$, and $b_0=h$.For $i>0$ we have $A_{i,i}=2/h$, $A_{i-1,i}=A_{i,i+1}=-1/h$.Second, we need to includethe extra term$-C{\varphi}_i(0)$ on the right-hand side. Since all ${\varphi}_i(0)=0$for $i=1,\ldots,N$, this term reduces to $-C{\varphi}_0(0)=-C$ andaffects only the first equation ($i=0$). We simply add $-C$ to $b_0$such that $b_0=h - C$.Third, the boundary term $-\int_0^L D{\varphi}_{N_n-1}(x){\varphi}_i{\, \mathrm{d}x}$must be computed. Since $i=0,\ldots,N=N_n-2$, this integral can onlyget a nonzero contribution with $i=N_n-2$ over the last cell.The result becomes $-Dh/6$.The resulting linear system can be summarized in the form $$\begin{equation}\frac{1}{h}\left(\begin{array}{ccccccccc}1 & -1 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -1 & 2 & -1 & \ddots &   & &  & &  \vdots \\ 0 & -1 & 2 & -1 &\ddots & &  &  & \vdots \\ \vdots & \ddots &  & \ddots & \ddots & 0 &  & & \vdots \\ \vdots &  & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & &  & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & &  & \ddots & \ddots & \ddots &\ddots  & 0 \\ \vdots & & & &  &\ddots  & \ddots &\ddots  & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots  & 0 & -1 & 2\end{array}\right)\left(\begin{array}{c}c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N}\end{array}\right)=\left(\begin{array}{c}h - C \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ 2h - Dh/6\end{array}\right){\thinspace .}\label{fem:deq:1D:BC:nat:Aub:system} \tag{86}\end{equation}$$ Next we consider the technique where we modify the linear system toincorporate Dirichlet conditions(see the section [Boundary term vanishes because of linear system modifications](fem:deq:1D:BC:nat:uLmod)). Now $N=N_n-1$.The two differences from thecase above is that the $-\int_0^LD{\varphi}_{N_n-1}{\varphi}_i{\, \mathrm{d}x}$ term isleft out of the right-hand side and an extra last row associatedwith the node $x_{N_n-1}=L$ where the Dirichlet condition appliesis appended to the system.This last row is anyway replaced by the condition $c_N=D$ or thiscondition can be incorporated in a symmetric fashion. Using the simplest,former approach gives $$\begin{equation}\frac{1}{h}\left(\begin{array}{ccccccccc}1 & -1 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -1 & 2 & -1 & \ddots &   & &  & &  \vdots \\ 0 & -1 & 2 & -1 &\ddots & &  &  & \vdots \\ \vdots & \ddots &  & \ddots & \ddots & 0 &  & & \vdots \\ \vdots &  & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & &  & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & &  & \ddots & \ddots & \ddots &\ddots  & 0 \\ \vdots & & & &  &\ddots  & -1 & 2  & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots  & 0 & 0 & h\end{array}\right)\left(\begin{array}{c}c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N}\end{array}\right)=\left(\begin{array}{c}h - C \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h \\ D\end{array}\right){\thinspace .}\label{fem:deq:1D:BC:nat:Aub:system:mod} \tag{87}\end{equation}$$  Cellwise computationsNow we compute with one element at a time, working in the referencecoordinate system $X\in [-1,1]$.We need to see how the$u'(0)=C$ condition affects the element matrix and vector.The extra term $-C{\varphi}_i(0)$ in the variational formulationonly affects the element vector in the first cell.On the reference cell, $-C{\varphi}_i(0)$ is transformed to$-C{\tilde{\varphi}}_r(-1)$, where $r$ counts local degrees of freedom.We have ${\tilde{\varphi}}_0(-1)=1$ and ${\tilde{\varphi}}_1(-1)=0$ sowe are left with the contribution$-C{\tilde{\varphi}}_0(-1)=-C$ to $\tilde b^{(0)}_0$: $$\begin{equation}\tilde A^{(0)} =A = \frac{1}{h}\left(\begin{array}{rr}1 & 1\\ -1 & 1\end{array}\right),\quad\tilde b^{(0)} = \left(\begin{array}{c}h - C\\ h\end{array}\right){\thinspace .}\label{fem:deq:1D:ex1:Ab:elm:bc:nat} \tag{88}\end{equation}$$ No other element matrices or vectors are affected by the $-C{\varphi}_i(0)$boundary term.There are two alternative ways of incorporating the Dirichlet condition.Following the section [Boundary term vanishes because of the test functions](fem:deq:1D:BC:nat:uLtest), we geta $1\times 1$ element matrix in the last cell andan element vector with an extra term containing $D$: $$\begin{equation}\tilde A^{(e)} =\frac{1}{h}\left(\begin{array}{r}1\end{array}\right),\quad\tilde b^{(e)} = h\left(\begin{array}{c}1 - D/6\end{array}\right),\label{fem:deq:1D:ex1:Ab:elm:ends2} \tag{89}\end{equation}$$ Alternatively, we include the degree of freedom at the node with$u$ specified. The element matrix and vector must then be modifiedto constrain the $\tilde c_1 = c_N$ value at local node $r=1$: $$\begin{equation}\tilde A^{(N_e)} =A = \frac{1}{h}\left(\begin{array}{rr}1 & 1\\ 0 & h\end{array}\right),\quad\tilde b^{(N_e)} = \left(\begin{array}{c}h\\ D\end{array}\right){\thinspace .}\label{fem:deq:1D:ex1:Ab:elm:bc:nat:mod} \tag{90}\end{equation}$$  Implementation of finite element algorithmsAt this point, it is sensible to create a program with symboliccalculations to perform all the steps in the computational machinery,both for automating the work and for documenting the completealgorithms.  As we have seen, there are quite many details involvedwith finite element computations and incorporation of boundaryconditions.  An implementation will also act as a structured summaryof all these details. Extensions of the code for approximationImplementation of the finite element algorithms for differentialequations follows closely the algorithm for approximation of functions.The new additional ingredients are1. other types of integrands (as implied by the variational formulation)2. additional boundary terms in the variational formulation for   Neumann boundary conditions3. modification of element matrices and vectors due to Dirichlet   boundary conditionsPoint 1 and 2 can be taken care of by letting the user supplyfunctions defining the integrands and boundary terms on theleft- and right-hand side of the equation system: * Integrand on the left-hand side: `ilhs(e, phi, r, s, X, x, h)` * Integrand on the right-hand side: `irhs(e, phi, r, X, x, h)` * Boundary term on the left-hand side: `blhs (e, phi, r, s, X, x, h)` * Boundary term on the right-hand side: `brhs (e, phi, r, s, X, x, h)`Here, `phi` is a dictionary where `phi[q]` holds a list of thederivatives of order `q` of the basis functions with respect to thephysical coordinate $x$.  The derivatives are available as Pythonfunctions of `X`.  For example, `phi[0][r](X)` means ${\tilde{\varphi}}_r(X)$,and `phi[1][s](X, h)` means $d{\tilde{\varphi}}_s (X)/dx$ (we refer tothe file [fe1D.py](src/fe1D.py) for detailsregarding the function `basis` that computes the `phi`dictionary).  The `r` and `s`arguments in the above functions correspond to the index in theintegrand contribution from an integration point to $\tildeA^{(e)}_{r,s}$ and $\tilde b^{(e)}_r$. The variables `e` and `h` arethe current element number and the length of the cell, respectively.Specific examples below will make it clear how to construct thesePython functions.Given a mesh represented by `vertices`, `cells`, and `dof_map` asexplained before, we can write a pseudo Python code to list allthe steps in the computational algorithm for finite element solutionof a differential equation. ```pythonfor e in range(len(cells)):     Compute element matrix and vector    n = len(dof_map[e])   no of dofs in this element    h = vertices[cells[e][1]] - vertices[cells[e][1]]         Integrate over the reference cell    points, weights =     for X, w in zip(points, weights):        phi =         detJ = h/2        dX = detJ*w        x =         for r in range(n):            for s in range(n):                A_e[r,s] += ilhs(e, phi, r, s, X, x, h)*dX            b_e[r] += irhs(e, phi, r, X, x, h)*dX     Add boundary terms    for r in range(n):        for s in range(n):            A_e[r,s] += blhs(e, phi, r, s, X, x)*dX        b_e[r] += brhs(e, phi, r, X, x, h)*dX     Incorporate essential boundary conditions    for r in range(n):        global_dof = dof_map[e][r]        if global_dof in essbc:             local dof r is subject to an essential condition            value = essbc[global_dof]             Symmetric modification            b_e -= value*A_e[:,r]            A_e[r,:] = 0            A_e[:,r] = 0            A_e[r,r] = 1            b_e[r] = value     Assemble    for r in range(n):        for s in range(n):            A[dof_map[e][r], dof_map[e][s]] += A_e[r,s]        b[dof_map[e][r] += b_e[r]```           Having implemented this function, the user only has supply the `ilhs`,`irhs`, `blhs`, and `brhs` functions that specify the PDE and boundaryconditions. The rest of the implementation forms a genericcomputational engine. The big finite element packages Diffpack and FEniCSare structured exactly the same way. A sample implementation of `ilhs`for the 1D Poisson problem is:<jupyter_code>def integrand_lhs(phi, i, j):
        return phi[1][i]*phi[1][j]<jupyter_output><empty_output><jupyter_text>which returns  $d{\tilde{\varphi}}_i (X)/dx d{\tilde{\varphi}}_j (X)/dx$.  Reducing theamount of code the user has to supply and making the code as close aspossible to the mathematical formulation makes the software user-friendly andeasy to debug.A complete function `finite_element1D_naive` for the 1D finitealgorithm above, is found in the file [fe1D.py](src/fe1D.py). The term "naive" refers to a version of thealgorithm where we use a standard dense square matrix as global matrix`A`. The implementation also has a verbose mode for printing out theelement matrices and vectors as they are computed.  Below is thecomplete function without the print statements.  You should study indetail since it contains all the steps in the finite elementalgorithm.<jupyter_code>def finite_element1D_naive(
    vertices, cells, dof_map,     # mesh
    essbc,                        # essbc[globdof]=value
    ilhs,                         # integrand left-hand side
    irhs,                         # integrand right-hand side
    blhs=lambda e, phi, r, s, X, x, h: 0,
    brhs=lambda e, phi, r, X, x, h: 0,
    intrule='GaussLegendre',      # integration rule class
    verbose=False,                # print intermediate results?
    ):
    N_e = len(cells)
    N_n = np.array(dof_map).max() + 1

    A = np.zeros((N_n, N_n))
    b = np.zeros(N_n)

    for e in range(N_e):
        Omega_e = [vertices[cells[e][0]], vertices[cells[e][1]]]
        h = Omega_e[1] - Omega_e[0]

        d = len(dof_map[e]) - 1  # Polynomial degree
        # Compute all element basis functions and their derivatives
        phi = basis(d)

        # Element matrix and vector
        n = d+1  # No of dofs per element
        A_e = np.zeros((n, n))
        b_e = np.zeros(n)

        # Integrate over the reference cell
        if intrule == 'GaussLegendre':
            points, weights = GaussLegendre(d+1)
        elif intrule == 'NewtonCotes':
            points, weights = NewtonCotes(d+1)

        for X, w in zip(points, weights):
            detJ = h/2
            x = affine_mapping(X, Omega_e)
            dX = detJ*w

            # Compute contribution to element matrix and vector
            for r in range(n):
                for s in range(n):
                    A_e[r,s] += ilhs(phi, r, s, X, x, h)*dX
                b_e[r] += irhs(phi, r, X, x, h)*dX

        # Add boundary terms
        for r in range(n):
            for s in range(n):
                A_e[r,s] += blhs(phi, r, s, X, x, h)
            b_e[r] += brhs(phi, r, X, x, h)

        # Incorporate essential boundary conditions
        modified = False
        for r in range(n):
            global_dof = dof_map[e][r]
            if global_dof in essbc:
                # dof r is subject to an essential condition
                value = essbc[global_dof]
                # Symmetric modification
                b_e -= value*A_e[:,r]
                A_e[r,:] = 0
                A_e[:,r] = 0
                A_e[r,r] = 1
                b_e[r] = value
                modified = True

        # Assemble
        for r in range(n):
            for s in range(n):
                A[dof_map[e][r], dof_map[e][s]] += A_e[r,s]
            b[dof_map[e][r]] += b_e[r]

    c = np.linalg.solve(A, b)
    return c, A, b, timing<jupyter_output><empty_output><jupyter_text>The `timing` object is a dictionary holding the CPU spent on computing`A` and the CPU time spent on solving the linear system. (We haveleft out the timing statements.) Utilizing a sparse matrixA potential efficiency problem with the `finite_element1D_naive` functionis that it uses dense $(N+1)\times (N+1)$ matrices, while we know thatonly $2d+1$ diagonals around the main diagonal are different from zero.Switching to a sparse matrix is very easy. Using the DOK (dictionary ofkeys) format, we declare `A` as ```pythonimport scipy.sparseA = scipy.sparse.dok_matrix((N_n, N_n))``` Assignments or in-place arithmetics are done as for a dense matrix,         A[i,j] += term        A[i,j]  = term but only the index pairs `(i,j)` we have used in assignments orin-place arithmetics are actually stored.A tailored solution algorithm is needed. The most reliable issparse Gaussian elimination. SciPy gives access to the[UMFPACK](https://en.wikipedia.org/wiki/UMFPACK) algorithmfor this purpose: ```pythonimport scipy.sparse.linalgc = scipy.sparse.linalg.spsolve(A.tocsr(), b, use_umfpack=True)``` The declaration of `A` and the solve statement are the onlychanges needed in the `finite_element1D_naive` to utilizesparse matrices. The resulting modification is found in thefunction `finite_element1D`. Application to our model problemLet us demonstrate the finite element software on $$-u''(x)=f(x),\quad x\in (0,L),\quad u'(0)=C,\ u(L)=D{\thinspace .}$$ This problem can be analytically solved by the`model2` function from the section "Simple model problems and their solutions" in previous chapter.Let $f(x)=x^2$. Calling `model2(x**2, L, C, D)` gives $$u(x) = D + C(x-L) + \frac{1}{12}(L^4 - x^4)$$ The variational formulation reads $$(u', v) = (x^2, v) - Cv(0){\thinspace .}$$ The entries in the element matrix and vector,which we need to set up the `ilhs`, `irhs`,`blhs`, and `brhs` functions, becomes $$\begin{align*}A^{(e)}_{r,s} &= \int_{-1}^1 \frac{d{\tilde{\varphi}}_r}{dx}\frac{{\tilde{\varphi}}_s}{dx}(\det J{\, \mathrm{d}X}),\\ b^{(e)} &= \int_{-1}^1 x^2{\tilde{\varphi}}_r\det J{\, \mathrm{d}X} - C{\tilde{\varphi}}_r(-1)I(e,0),\end{align*}$$ where $I(e)$ is an indicator function: $I(e,q)=1$ if $e=q$, otherwise $I(e)=0$.We use this indicator function to formulate that the boundary term$Cv(0)$, which in the local element coordinate system becomes $C{\tilde{\varphi}}_r(-1)$,is only included for the element $e=0$.The functions for specifying the element matrix and vector entriesmust contain the integrand, but without the $\det J{\, \mathrm{d}X}$ termas this term is taken care of by the quadrature loop, andthe derivatives $d{\tilde{\varphi}}_r(X)/dx$with respect to the physical $x$ coordinates arecontained in `phi[1][r](X)`, computed by the function `basis`.<jupyter_code>def ilhs(e, phi, r, s, X, x, h):
    return phi[1][r](X, h)*phi[1][s](X, h)

def irhs(e, phi, r, X, x, h):
    return x**2*phi[0][r](X)

def blhs(e, phi, r, s, X, x, h):
    return 0

def brhs(e, phi, r, X, x, h):
    return -C*phi[0][r](-1) if e == 0 else 0<jupyter_output><empty_output><jupyter_text>We can then make the call to `finite_element1D_naive` or `finite_element1D`to solve the problem with two P1 elements:<jupyter_code>from src.fe1D import finite_element1D_naive, mesh_uniform
C = 5;  D = 2;  L = 4
d = 1

vertices, cells, dof_map = mesh_uniform(
    N_e=2, d=d, Omega=[0,L], symbolic=False)
essbc = {}
essbc[dof_map[-1][-1]] = D

c, A, b, timing = finite_element1D_naive(
    vertices, cells, dof_map, essbc,
    ilhs=ilhs, irhs=irhs, blhs=blhs, brhs=brhs,
    intrule='GaussLegendre')<jupyter_output><empty_output><jupyter_text>It remains to plot the solution (with high resolution in each element).To this end, we use the `u_glob` function imported from`fe1D`, which imports it from `fe_approx1D_numit` (the`u_glob` function in `fe_approx1D.py`works with `elements` and `nodes`, while `u_glob` in`fe_approx1D_numint` works with `cells`, `vertices`,and `dof_map`):<jupyter_code>u_exact = lambda x: D + C*(x-L) + (1./6)*(L**3 - x**3)
from src.fe1D import u_glob
x, u, nodes = u_glob(c, cells, vertices, dof_map)
u_e = u_exact(x)
# print((u_exact(nodes) - c))  # difference at the nodes

# import matplotlib.pyplot as plt
# plt.plot(x, u, 'b-', x, u_e, 'r--')
# plt.legend(['finite elements, d=%d' %d, 'exact'], loc='upper left')
# plt.show()<jupyter_output><empty_output><jupyter_text>The result is shown in [Figure](fem:deq:1D:code:fe:fig1). We seethat the solution using P1 elements is exact at the nodes, butfeature considerable discrepancy between the nodes.[Exercise 10: Investigate exact finite element solutions](fem:deq:exer:1D:exact_numerics) asks you to explorethis problem further using other $m$ and $d$ values. -->Finite element and exact solution using two cells. Variational formulations in 2D and 3DThe major difference between deriving variational formulations in 2Dand 3D compared to 1D is the rule for integrating by parts.  The cellshave shapes different from an interval, so basis functions look a bitdifferent, and there is a technical difference in actually calculatingthe integrals over cells. Otherwise, going to 2D and 3D is not a bigstep from 1D. All the fundamental ideas still apply. Integration by partsA typical second-order term in a PDE may be written in dimension-independentnotation as $$\nabla^2 u \quad\hbox{or}\quad \nabla\cdot\left( {\alpha}(\boldsymbol{x})\nabla u\right){\thinspace .}$$ The explicit forms in a 2D problem become $$\nabla^2 u = \nabla\cdot\nabla u =\frac{\partial^2 u}{\partial x^2} +\frac{\partial^2 u}{\partial y^2},$$ and $$\nabla\cdot\left( a(\boldsymbol{x})\nabla u\right) =\frac{\partial}{\partial x}\left( {\alpha}(x,y)\frac{\partial u}{\partial x}\right) +\frac{\partial}{\partial y}\left( {\alpha}(x,y)\frac{\partial u}{\partial y}\right){\thinspace .}$$ We shall continue with the latter operator as the former arises fromjust setting ${\alpha} =1$.**The integration by parts formula for $\int\nabla\cdot({\alpha}\nabla)$.**The general rule for integrating by parts is often referred to as[Green's first identity](http://en.wikipedia.org/wiki/Green's_identities): $$\begin{equation}-\int_{\Omega} \nabla\cdot ({\alpha}(\boldsymbol{x})\nabla u) v{\, \mathrm{d}x} =\int_{\Omega} {\alpha}(\boldsymbol{x})\nabla u\cdot\nabla v {\, \mathrm{d}x} -\int_{\partial\Omega} a\frac{\partial u}{\partial n} v {\, \mathrm{d}s},\label{fem:deq:2D:int:by:parts} \tag{91}\end{equation}$$ where $\partial\Omega$ is the boundary of $\Omega$ and$\partial u/\partial n = \boldsymbol{n}\cdot\nabla u$ is the derivativeof $u$ in the outward normal direction, $\boldsymbol{n}$ being an outwardunit normal to $\partial\Omega$. The integrals $\int_\Omega (){\, \mathrm{d}x}$ arearea integrals in 2D and volume integrals in 3D, while$\int_{\partial\Omega} (){\, \mathrm{d}s}$ is a line integral in 2D and a surfaceintegral in 3D.It will be convenient to divide the boundary into two parts: * $\partial\Omega_N$, where we have Neumann conditions   $-a\frac{\partial u}{\partial n} = g$, and * $\partial\Omega_D$, where we have Dirichlet conditions   $u = u_0$.The test functions $v$ are (as usual) required to vanish on$\partial\Omega_D$. Example on a multi-dimensional variational problemHere is a quite general, stationary, linear PDE arising in many problems: $$\begin{equation}\label{varform:conv:diff:pde:pre} \tag{92}\boldsymbol{v}\cdot\nabla u + \beta u = \nabla\cdot\left( {\alpha}\nabla u\right) + f,\quad\boldsymbol{x}\in\Omega,\end{equation}$$ $$\begin{equation}  \label{varform:conv:diff:bc1:pre} \tag{93}u = u_0,\quad\boldsymbol{x}\in\partial\Omega_D,\end{equation}$$ $$\begin{equation}  \label{varform:conv:diff:bc2:pre} \tag{94}-{\alpha}\frac{\partial u}{\partial n} = g,\quad\boldsymbol{x}\in\partial\Omega_N{\thinspace .}\end{equation}$$ The vector field $\boldsymbol{v}$ and the scalar functions $a$, $\alpha$, $f$, $u_0$, and$g$ may vary with the spatial coordinate $\boldsymbol{x}$ and must be known.Such a second-order PDE needs exactly one boundary condition at eachpoint of the boundary, so $\partial\Omega_N\cup\partial\Omega_D$must be the complete boundary $\partial\Omega$.Assume that the boundary function $u_0(\boldsymbol{x})$ is defined for all $\boldsymbol{x}\in\Omega$.The unknown function can then be expanded as $$u = B + \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j,\quad B = u_0 {\thinspace .}$$ As long as any ${\psi}_j=0$ on $\partial\Omega_D$, we realize that $u=u_0$on $\partial\Omega_D$.The variational formula is obtained from Galerkin's method, whichtechnically means multiplying the PDE by a testfunction $v$ and integrating over $\Omega$: $$\int_{\Omega} (\boldsymbol{v}\cdot\nabla u + \beta u)v{\, \mathrm{d}x} =\int_{\Omega} \nabla\cdot\left( {\alpha}\nabla u\right){\, \mathrm{d}x} + \int_{\Omega}fv {\, \mathrm{d}x}{\thinspace .}$$ The second-order term is integrated by parts, according to the formula([91](fem:deq:2D:int:by:parts)): $$\int_{\Omega} \nabla\cdot\left( {\alpha}\nabla u\right)v {\, \mathrm{d}x} =-\int_{\Omega} {\alpha}\nabla u\cdot\nabla v{\, \mathrm{d}x}+ \int_{\partial\Omega} {\alpha}\frac{\partial u}{\partial n} v{\, \mathrm{d}s}{\thinspace .}$$ Galerkin's method therefore leads to $$\int_{\Omega} (\boldsymbol{v}\cdot\nabla u + \beta u)v{\, \mathrm{d}x} =-\int_{\Omega} {\alpha}\nabla u\cdot\nabla v{\, \mathrm{d}x}+ \int_{\partial\Omega} {\alpha}\frac{\partial u}{\partial n} v{\, \mathrm{d}s}+ \int_{\Omega} fv {\, \mathrm{d}x}{\thinspace .}$$ The boundary term can be developed further by noticing that $v\neq 0$only on $\partial\Omega_N$, $$\int_{\partial\Omega} {\alpha}\frac{\partial u}{\partial n} v{\, \mathrm{d}s}= \int_{\partial\Omega_N} {\alpha}\frac{\partial u}{\partial n} v{\, \mathrm{d}s},$$ and that on $\partial\Omega_N$, we have the condition$a\frac{\partial u}{\partial n}=-g$, so the term becomes $$-\int_{\partial\Omega_N} gv{\, \mathrm{d}s}{\thinspace .}$$ The final variational form is then $$\int_{\Omega} (\boldsymbol{v}\cdot\nabla u + \beta u)v{\, \mathrm{d}x} =-\int_{\Omega} {\alpha}\nabla u\cdot\nabla v {\, \mathrm{d}x}- \int_{\partial\Omega_N} g v{\, \mathrm{d}s}+ \int_{\Omega} fv {\, \mathrm{d}x}{\thinspace .}$$ Instead of using the integral signs, we may use the inner productnotation: $$(\boldsymbol{v}\cdot\nabla u, v) + (\beta u,v) =- ({\alpha}\nabla u,\nabla v) - (g,v)_{N} + (f,v){\thinspace .}$$ The subscript $\,{}_N$ in $(g,v)_{N}$ is a notation for a line or surfaceintegral over $\partial\Omega_N$, while $(\cdot,\cdot)$ is the area/volumeintegral over $\Omega$.We can derive explicit expressions for the linear system for $\left\{ {c}_j \right\}_{j\in{\mathcal{I}_s}}$that arises from the variational formulation.Inserting the $u$ expansion results in $$\begin{align*}\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 = \\ & (g,{\psi}_i)_{N} + (f,{\psi}_i) -(\boldsymbol{v}\cdot\nabla u_0, {\psi}_i) + (\beta u_0 ,{\psi}_i) +({\alpha}\nabla u_0,\nabla {\psi}_i){\thinspace .}\end{align*}$$ This is a linear system with matrix entries $$A_{i,j} = (\boldsymbol{v}\cdot\nabla {\psi}_j, {\psi}_i) + (\beta {\psi}_j ,{\psi}_i) + ({\alpha}\nabla {\psi}_j,\nabla {\psi}_i)$$ and right-hand side entries $$b_i = (g,{\psi}_i)_{N} + (f,{\psi}_i) -(\boldsymbol{v}\cdot\nabla u_0, {\psi}_i) + (\beta u_0 ,{\psi}_i) +({\alpha}\nabla u_0,\nabla {\psi}_i),$$ for $i,j\in{\mathcal{I}_s}$.In the finite element method, we usually express $u_0$ in terms ofbasis functions and restrict $i$ and $j$ to run over the degrees offreedom that are not prescribed as Dirichlet conditions.However, we can also keep all the $\left\{ {c}_j \right\}_{j\in{\mathcal{I}_s}}$ as unknowns,drop the $u_0$ in the expansion for $u$, and incorporate all theknown $c_j$ values in the linear system. This has been explainedin detail in the 1D case, and the technique is the same for 2D and3D problems. Transformation to a reference cell in 2D and 3DThe real power of the finite element method first becomes evident whenwe want to solve partial differential equations posed on two- andthree-dimensional domains of non-trivial geometric shape.  As in 1D,the domain $\Omega$ is divided into $N_e$ non-overlapping cells. Theelements have simple shapes: triangles and quadrilaterals are popularin 2D, while tetrahedra and box-shapes elements dominate in 3D.  Thefinite element basis functions ${\varphi}_i$ are, as in 1D, polynomialsover each cell.  The integrals in the variational formulation are, asin 1D, split into contributions from each cell, and thesecontributions are calculated by mapping a physical cell, expressed inphysical coordinates $\boldsymbol{x}$, to a reference cell in a local coordinatesystem $\boldsymbol{X}$. This mapping will now be explained in detail.We consider an integral of the type $$\begin{equation}\int_{{\Omega}^{(e)}} {\alpha}(\boldsymbol{x})\nabla{\varphi}_i\cdot\nabla{\varphi}_j{\, \mathrm{d}x},\label{_auto39} \tag{95}\end{equation}$$ where the ${\varphi}_i$ functions are finite element basis functions in2D or 3D, defined in the physical domain.Suppose we want to calculate this integral over a reference cell,denoted by $\tilde\Omega^r$, in a coordinate system with coordinates$\boldsymbol{X} = (X_0, X_1)$ (2D) or $\boldsymbol{X} = (X_0, X_1, X_2)$ (3D).The mapping between a point $\boldsymbol{X}$ in the reference coordinate system  andthe corresponding point $\boldsymbol{x}$ in the physical coordinate system isgiven by a vector relation $\boldsymbol{x}(\boldsymbol{X})$.The corresponding Jacobian, $J$, of this mapping has entries $$J_{i,j}=\frac{\partial x_j}{\partial X_i}{\thinspace .}$$ The change of variables requires ${\, \mathrm{d}x}$ to be replaced by $\det J{\, \mathrm{d}X}$.The derivatives in the $\nabla$ operator in the variational form arewith respect to $\boldsymbol{x}$, which we may denote by $\nabla_{\boldsymbol{x}}$.The ${\varphi}_i(\boldsymbol{x})$ functions in the integralare replaced by local basis functions ${\tilde{\varphi}}_r(\boldsymbol{X})$ sothe integral features $\nabla_{\boldsymbol{x}}{\tilde{\varphi}}_r(\boldsymbol{X})$. We readily have$\nabla_{\boldsymbol{X}}{\tilde{\varphi}}_r(\boldsymbol{X})$ from formulas for the basis functions inthe reference cell, butthe desired quantity $\nabla_{\boldsymbol{x}}{\tilde{\varphi}}_r(\boldsymbol{X})$ requires some effortsto compute. All the details are provided below.Let $i=q(e,r)$ and consider two space dimensions. By the chain rule, $$\frac{\partial {\tilde{\varphi}}_r}{\partial X} =\frac{\partial {\varphi}_i}{\partial X} =\frac{\partial {\varphi}_i}{\partial x}\frac{\partial x}{\partial X} +\frac{\partial {\varphi}_i}{\partial y}\frac{\partial y}{\partial X},$$ and $$\frac{\partial {\tilde{\varphi}}_r}{\partial Y} =\frac{\partial {\varphi}_i}{\partial Y} =\frac{\partial {\varphi}_i}{\partial x}\frac{\partial x}{\partial Y} +\frac{\partial {\varphi}_i}{\partial y}\frac{\partial y}{\partial Y}{\thinspace .}$$ We can write these two equations as a vector equation $$\left[\begin{array}{c}\frac{\partial {\tilde{\varphi}}_r}{\partial X}\\ \frac{\partial {\tilde{\varphi}}_r}{\partial Y}\end{array}\right]=\left[\begin{array}{cc}\frac{\partial x}{\partial X} & \frac{\partial y}{\partial X}\\ \frac{\partial x}{\partial Y} & \frac{\partial y}{\partial Y}\end{array}\right]\left[\begin{array}{c}\frac{\partial {\varphi}_i}{\partial x}\\ \frac{\partial {\varphi}_i}{\partial y}\end{array}\right]$$ Identifying $$\nabla_{\boldsymbol{X}}{\tilde{\varphi}}_r = \left[\begin{array}{c}\frac{\partial {\tilde{\varphi}}_r}{\partial X}\\ \frac{\partial {\tilde{\varphi}}_r}{\partial Y}\end{array}\right],\quadJ =\left[\begin{array}{cc}\frac{\partial x}{\partial X} & \frac{\partial y}{\partial X}\\ \frac{\partial x}{\partial Y} & \frac{\partial y}{\partial Y}\end{array}\right],\quad\nabla_{\boldsymbol{x}}{\varphi}_r =\left[\begin{array}{c}\frac{\partial {\varphi}_i}{\partial x}\\ \frac{\partial {\varphi}_i}{\partial y}\end{array}\right],$$ we have the relation $$\nabla_{\boldsymbol{X}}{\tilde{\varphi}}_r = J\cdot\nabla_{\boldsymbol{x}}{\varphi}_i,$$ which we can solve with respect to $\nabla_{\boldsymbol{x}}{\varphi}_i$: $$\begin{equation}\nabla_{\boldsymbol{x}}{\varphi}_i = J^{-1}\cdot\nabla_{\boldsymbol{X}}{\tilde{\varphi}}_r{\thinspace .}\label{_auto40} \tag{96}\end{equation}$$ On the reference cell, ${\varphi}_i(\boldsymbol{x}) = {\tilde{\varphi}}_r(\boldsymbol{X})$, so $$\begin{equation}\nabla_{\boldsymbol{x}}{\tilde{\varphi}}_r(\boldsymbol{X}) = J^{-1}(\boldsymbol{X})\cdot\nabla_{\boldsymbol{X}}{\tilde{\varphi}}_r(\boldsymbol{X}){\thinspace .}\label{_auto41} \tag{97}\end{equation}$$ This means that we have the following transformation of theintegral in the physical domain to its counterpart over the reference cell: $$\begin{equation}\int_{\Omega^{(e)}} {\alpha}(\boldsymbol{x})\nabla_{\boldsymbol{x}}{\varphi}_i\cdot\nabla_{\boldsymbol{x}}{\varphi}_j{\, \mathrm{d}x} =\int_{\tilde\Omega^r} {\alpha}(\boldsymbol{x}(\boldsymbol{X}))(J^{-1}\cdot\nabla_{\boldsymbol{X}}{\tilde{\varphi}}_r)\cdot(J^{-1}\cdot\nabla{\tilde{\varphi}}_s)\det J{\, \mathrm{d}X}\label{_auto42} \tag{98}\end{equation}$$  Numerical integrationIntegrals are normally computed by numerical integration rules.For multi-dimensional cells, various families of rules exist.All of them are similar to what is shown in 1D:$\int f {\, \mathrm{d}x}\approx \sum_jw_if(\boldsymbol{x}_j)$, where $w_j$ are weights and$\boldsymbol{x}_j$ are corresponding points.The file [numint.py](src/numint.py) contains the functions`quadrature_for_triangles(n)` and `quadrature_for_tetrahedra(n)`,which returns lists of points and weights corresponding to integrationrules with `n` points over the reference trianglewith vertices $(0,0)$, $(1,0)$, $(0,1)$, and the reference tetrahedronwith vertices $(0,0,0)$, $(1,0,0)$, $(0,1,0)$, $(0,0,1)$,respectively. For example, the first two rules for integration overa triangle have 1 and 3 points:<jupyter_code>import numint
x, w = numint.quadrature_for_triangles(num_points=1)
x
w
x, w = numint.quadrature_for_triangles(num_points=3)
x
w<jupyter_output><empty_output><jupyter_text>Rules with 1, 3, 4, and 7 points over the triangle will exactly integratepolynomials of degree 1, 2, 3, and 4, respectively.In 3D, rules with 1, 4, 5, and 11 points over the tetrahedron willexactly integrate polynomials of degree 1, 2, 3, and 4, respectively. Convenient formulas for P1 elements in 2DWe shall now provide some formulas for piecewise linear ${\varphi}_i$ functionsand their integrals *in the physical coordinate system*.These formulas make it convenient to compute with P1 elements withoutthe need to work in the reference coordinate system and deal with mappingsand Jacobians.A lot of computational and algorithmic details are hidden by this approach.Let $\Omega^{(e)}$ be cell number $e$, and let the three verticeshave global vertex numbers $I$, $J$, and $K$.The corresponding coordinates are$(x_{I},y_{I})$, $(x_{J},y_{J})$, and $(x_{K},y_{K})$.The basis function ${\varphi}_I$ over $\Omega^{(e)}$ have the explicitformula $$\begin{equation}{\varphi}_I (x,y) = \frac{1}{2}\Delta \left( \alpha_I + \beta_Ix+ \gamma_Iy\right),\label{fem:approx:fe:2D:phi:I} \tag{99}\end{equation}$$ where $$\begin{equation}\alpha_I = x_{J}y_{K} - x_{K}y_{J},\label{fem:approx:fe:2D:phi:alpha:I} \tag{100}\end{equation}$$ $$\begin{equation}  \beta_I = y_{J} - y_{K},\label{fem:approx:fe:2D:phi:beta:I} \tag{101}\end{equation}$$ $$\begin{equation}  \gamma_I = x_{K} - x_{J},\label{fem:approx:fe:2D:phi:gamma:I} \tag{102},\end{equation}$$ and $$\begin{equation}2\Delta = \det\left(\begin{array}{rrr}1 & x_{I} & y_{I} \\ 1 & x_{J} & y_{J} \\ 1 & x_{K} & y_{K} \end{array}\right){\thinspace .}\label{fem:approx:fe:2D:phi:Delta} \tag{103}\end{equation}$$ The quantity $\Delta$ is the area of the cell.The following formula is often convenient when computing element matricesand vectors: $$\begin{equation}\int_{\Omega^{(e)}} {\varphi}_I^{p}{\varphi}_J^{q}{\varphi}_K^{r} dx dy ={p!q!r!\over (p+q+r+2)!}2\Delta\label{fem:approx:fe:2D:phi:integral} \tag{104}{\thinspace .}\end{equation}$$ (Note that the $q$ in this formula is not to be mixed with the $q(e,r)$mapping of degrees of freedom.)As an example, the element matrix entry$\int_{\Omega^{(e)}} {\varphi}_I{\varphi}_J{\, \mathrm{d}x}$can be computed by setting$p=q=1$ and $r=0$, when $I\neq J$, yielding $\Delta/12$, and$p=2$ and $q=r=0$, when $I=J$, resulting in $\Delta/6$.We collect these numbers in a local element matrix: $$\frac{\Delta}{12}\left[\begin{array}{ccc}2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2\end{array}\right]$$ 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 becomputed by the formulas above. We have $$\nabla{\varphi}_I\cdot\nabla{\varphi}_J =\frac{\Delta^2}{4}(\beta_I\beta_J + \gamma_I\gamma_J) = \hbox{const},$$ so that the element matrix entry becomes$\frac{1}{4}\Delta^3(\beta_I\beta_J + \gamma_I\gamma_J)$.From an implementational point of view, one will work with local vertexnumbers $r=0,1,2$, parameterize the coefficients in the basisfunctions by $r$, and look up vertex coordinates through $q(e,r)$.Similar formulas exist for integration of P1 elements in 3D.  Implementation in 2D and 3D via FEniCSFrom a principle of view, we have seen that variational forms of thetype: find $a(u,v)=L\ \forall v\in V$ (and even general nonlinear problems$F(u;v)=0$), can apply the computational machinery of introduced forthe approximation problem $u=f$. We actually need two extensions only:1. specify Dirichlet boundary conditions as part of $V$2. incorporate Neumann flux boundary conditions in the variational formThe algorithms are all the same in any space dimension, we only need tochoose the element type and associated integration rule. Once we knowhow to compute things in 1D, and made the computer code sufficientlyflexible, the method and code should work for any variational form inany number of space dimensions! This fact is exactly the idea behindthe [FEniCS](http://fenicsproject.org) finite element software.Therefore, if we know how to set up an approximation problem in anydimension in FEniCS, and know how to derive variational forms in higherdimensions, we are (in principle!) very close to solving aPDE problem in FEniCS. We shall now solve a quite general 1D/2D/3D Poisson problem in FEniCS.There is much more to FEniCS than what is shown in this example, but itillustrates the fact that when we go beyond 1D, there existssoftware which leverage the full power of the finite element method asa method for solving any problem on any mesh in any number ofspace dimensions. Mathematical problemThe following model describes the pressure $u$ in the flow around abore hole of radius $a$ in a porous medium. If the hole is long in thevertical direction (z-direction) then it is natural to assume that the vertical changes are small and $u_z\approx\mbox{constant}$. Therefore, we can model it by a 2D domain in the cross section. $$\begin{equation}\nabla\cdot \left( {\alpha}\nabla u\right) = 0,  \quad a < ||\boldsymbol{x}|| < b,\label{fem:varform:fenics:problem:pde} \tag{105}\end{equation}$$ $$\begin{equation}  u(\boldsymbol{x}) = U_a, \quad   ||\boldsymbol{x}|| = a,\label{fem:varform:fenics:problem:ua} \tag{106}\end{equation}$$ $$\begin{equation}  u(\boldsymbol{x}) = U_b  \quad   ||\boldsymbol{x}|| = b{\thinspace .}\label{fem:varform:fenics:problem:ub} \tag{107}\end{equation}$$ That is, we have a hollow circular 2D domain with inner radius $a$ andouter radius $b$. The pressure is known on these two boundaries, sothis is a pure Dirichlet problem. SymmetryThe first thing we should observe is that the problem is radiallysymmetric, so we can change to polar coordinates and obtain a 1Dproblem in the radial direction: $$(r{\alpha} u')' = 0,\quad u(a)=U_a, u(b)=U_b{\thinspace .}$$ This is not very exciting beyond being able to find an analytical solutionand compute the true error of a finite element approximation.However, many software packages solve problems in Cartesian coordinates, andFEniCS basically do this, so we want to take advantage of symmetry inCartesian coordinates and reformulate the problem in a smallerdomain.Looking at the domain as a cake with a hole, any piece of thecake will be a candidate for a reduced-size domain.  The solution issymmetric about any line $\theta = \hbox{const}$ in polar coordinates,so at such lines we have the symmetry boundary condition $\partialu/\partial n=0$, i.e., a homogeneous Neumann condition.  In [Figure](fem:varform:fenics:problem:meshfig) we have plotted a possiblemesh of cells as triangles, here with dense refinement toward the borehole, because we know the solution will decay most rapidly toward theorigin.  This mesh is a piece of the cake with four sides: Dirichletconditions on the inner and outer boundary, named $\Gamma_{D_a}$ and$\Gamma_{D_b}$, and $\partial u/\partial n=0$on the two other sides, named $\Gamma_N$.In this particular example, the arc of the pieceof the cake is 45 degrees, but any value of the arc will work. -->Mesh of a hollow cylinder, with refinement and utilizing symmetry.The boundary problem can then be expressed as $$\begin{equation}\nabla\cdot \left( {\alpha}\nabla u\right) = 0,  \quad \boldsymbol{x}\in\Omega,\label{fem:varform:fenics:problem:pde2} \tag{108}\end{equation}$$ $$\begin{equation}  u(\boldsymbol{x}) = U_a,  \quad   \boldsymbol{x}\in\Gamma_{D_a},\label{fem:varform:fenics:problem:Ga} \tag{109}\end{equation}$$ $$\begin{equation}  u(\boldsymbol{x}) = U_b,  \quad   \boldsymbol{x}\in\Gamma_{D_b},\label{fem:varform:fenics:problem:Gb} \tag{110}\end{equation}$$ $$\begin{equation}  \frac{\partial u}{\partial n} =0,\quad  \boldsymbol{x}\in\Gamma_N{\thinspace .}\label{fem:varform:fenics:problem:GN} \tag{111}\end{equation}$$  Variational formulationTo obtain the variational formulation, we multiply the PDE by a testfunction $v$ and integrate the second-order derivatives by part: $$\begin{align*}\int_\Omega \nabla\cdot ({\alpha}\nabla u) v {\, \mathrm{d}x} &= -\int_\Omega {\alpha}\nabla u\cdot\nabla v{\, \mathrm{d}x} + \int_{\Gamma_N}{\alpha}\frac{\partial u}{\partial n}v{\, \mathrm{d}s}\\ &= -\int_\Omega {\alpha}\nabla u\cdot\nabla v{\, \mathrm{d}x},\quad \forall v\in V\\ \end{align*}$$ We are left with a problem of the form: find $u$ such that$a(u,v)=L(v)\ \forall v\in V$, with $$\begin{equation}a(u,v) = \int_\Omega {\alpha}\nabla u\cdot\nabla v{\, \mathrm{d}x},\label{_auto43} \tag{112}\end{equation}$$ $$\begin{equation}  L(v) = \int_\Omega 0v{\, \mathrm{d}x} {\thinspace .}\label{_auto44} \tag{113}\end{equation}$$ We write the integrand as $0v{\, \mathrm{d}x}$ even though $L=0$, because it is necessaryin FEniCS to specify $L$ as a linear form (i.e., atest function and some form of integration need to be present) and not the number zero.The Dirichlet conditions make a nonzero solution. The FEniCS solverWe suppose that we have a function `make_mesh` that can make the mesh for us. Moredetails about this function will be provided later. A next step is then to define proper Dirichlet conditions.This might seem a bit complicated, but we introduce *markers* at theboundary for marking the Dirichlet boundaries. The inner boundary hasmarker 1 and the outer has marker 2. In this way, we can recognize thenodes that are on the boundary. It is usually a part of the mesh makingprocess to compute both the mesh and its markers, so `make_mesh` returnsa `Mesh` object as `mesh` and a `MeshFunction` object `markers`.Setting Dirichlet conditions in the solver is then a matter ofintroducing `DirichletBC` objects, one for each part of the boundarymarked by `markers`, and then we collect all such Dirichlet conditionsin a list that is used by the assembly process to incorporate theDirichlet  conditions in the linear system.The code goes like this:<jupyter_code>V = FunctionSpace(mesh, 'P', degree)
bc_inner = DirichletBC(V, u_a, markers, 1)
bc_outer = DirichletBC(V, u_b, markers, 2)
bcs = [bc_inner, bc_outer]<jupyter_output><empty_output><jupyter_text>Here, `u_a` and `u_b` are constants (floats) set by the user. In general, anything that can be evaluated pointwise can be used, such as `Expression`, `Function`, and `Constant`. The next step is to define the variational problem and solve it:<jupyter_code># Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = alpha*dot(grad(u), grad(v))*dx
L = Constant(0)*v*dx  # L = 0*v*dx = 0 does not work...

# Compute solution
u = Function(V)
solve(a == L, u, bcs)

f = File("mesh.xml")
f << mesh<jupyter_output><empty_output><jupyter_text>In order to avoid `L=0` (`L` equal to the float zero), we have totell FEniCS that is a linear form, so zero must be specified as `Constant(0)`.Note that everything is the same as for the approximation problem in the section [fe:approx:fenics](fe:approx:fenics),except for the Dirichlet conditions and the formulas for `a` and`L`. FEniCS has, of course, access to very efficient solution methods,so we could add arguments to the `solve` call to applystate-of-the-art iterative methods and preconditioners for large-scaleproblems. However, for this little 2D case a standard sparse Gaussianelimination, as implied by `solve(a = L, u, bcs)` is a sufficient approach.Finally, we can save the solution to file for using professionalvisualization software and, if desired, add a quick plotting using thebuilt-in FEniCS tool `plot`:<jupyter_code># Save solution to file in VTK format
vtkfile = File(filename + '.pvd')
vtkfile << u

u.rename('u', 'u'); plot(u); plot(mesh)
import matplotlib.pyplot as plt
plt.show()<jupyter_output><empty_output><jupyter_text>(The `u.rename` call is just for getting a more readable title in the plot.)The above statements are collected in a function `solver` in thefile [`borehole_fenics.py`](${fem_src}/borehole_fenics.py):<jupyter_code>def solver(
    mesh,
    markers,  # MeshFunctions for Dirichlet conditions
    alpha,    # Diffusion coefficient
    u_a,      # Inner pressure
    u_b,      # Outer pressure
    degree,   # Element polynomial degree
    filename, # Name of VTK file
    ):
    V = FunctionSpace(mesh, 'P', degree)
    bc_inner = DirichletBC(V, u_a, markers, 1)
    bc_outer = DirichletBC(V, u_b, markers, 2)
    bcs = [bc_inner, bc_outer]

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    a = alpha*dot(grad(u), grad(v))*dx
    L = Constant(0)*v*dx  # L = 0*v*dx = 0 does not work...

    # Compute solution
    u = Function(V)
    solve(a == L, u, bcs)

    f = File("mesh.xml")
    f << mesh

    # Save solution to file in VTK format
    vtkfile = File(filename + '.pvd')
    vtkfile << u

    u.rename('u', 'u'); plot(u); plot(mesh)
    import matplotlib.pyplot as plt
    plt.show()
    return u

def problem():
    mesh, markers = make_mesh(Theta=25*pi/180, a=1, b=2,
                              nr=20, nt=20, s=1.9)
    beta = 5
    solver(mesh, markers, alpha=1, u_a=1, u_b=0, degree=1, filename='tmp')

if __name__ == '__main__':
    problem()<jupyter_output><empty_output><jupyter_text>**Be careful with name clashes!**It is easy when coding mathematics to use variable names that correspondto one-letter names in the mathematics. For example, in the mathematicsof this problem there are two $a$ variables: the radius of the innerboundary and the bilinear form in the variational formulation.Using `a` for the inner boundary in `solver` does not work: it isquickly overwritten by the bilinear form. We therefore have to introduce`x_a`. Long variable names are to be preferred for safe programming,though short names corresponding to the mathematics are often nicer. Making the meshThe hardest part of a finite element problem is very often to make themesh.  This is particularly the case in large industrial projects, butalso often academic projects quickly lead to time-consuming work withconstructing finite element meshes. In the present example we createthe mesh for the symmetric problem by deforming an originallyrectangular mesh.  The rectangular mesh is made by the FEniCS object`RectangleMesh` on $[a,b]\times [0,1]$.  Therefore, we stretch themesh towards the left before we bend the rectangle onto to "a pieceof cake". [Figure](fem:varform:fenics:mesh1:fig1) shows an exampleon the resulting mesh. The stretching gives us refinement in theradial direction because we expect the variations to be quite large inthis direction, but uniform in $\theta$ direction.We first make the rectangle and set boundary markers here for the innerand outer boundaries (since these are particularly simple: $x=a$ and $x=b$).Here is how we makethe rectangular mesh from lower left corner $(a,0)$ to upper left corner$(b,1)$ with `nr` quadrilateral cells in $x$ direction (later to becomethe radial direction) and `nt` quadrilateral cells in the $y$ direction:<jupyter_code>mesh = RectangleMesh(Point(a, 0), Point(b, 1), nr, nt, 'crossed')<jupyter_output><empty_output><jupyter_text>Each quadrilateral cell is divided into two triangles with right or leftgoing diagonals, or four triangles using both diagonals. These choicesof producing triangles from rectangles are named `right`, `left`, and`crossed`. Recall that FEniCS can only work with cells of triangularshape only and where the sides are straight. This means that we need agood resolution in $\theta$ direction to represent a circular boundary.With isoparametric elements, it is easier to get a higher-order polynomialapproximation of curved boundaries.We must then mark the boundaries for boundary conditions. Since we do not needto do anything with the homogeneous Neumann conditions, we can just markthe inner and outer boundary of the hole cylinder. This is very easy to doas long as the mesh is a rectangle since then the specifications of theboundaries are $x=a$ and $x=b$. The relevant FEniCS code requires theuser to define a subclass of `SubDomain` and implement a function `inside`for indicating whether a given point `x` is on the desired boundary or not: -->Finite element mesh for a porous medium outside a bore hole (hollow cylinder).<jupyter_code># x=a becomes the inner borehole boundary
class Inner(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - a) < tol

# x=b becomes the outer borehole boundary
class Outer(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - b) < tol

inner = Inner(); outer = Outer();
markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
markers.set_all(0)
inner.mark(markers, 1)
outer.mark(markers, 2)<jupyter_output><empty_output><jupyter_text>With the instances `inner` and `outer` we fill a marker object, called`MeshFunction` in FEniCS. For this purpose we must introduce our ownconvention of numbering boundaries: here we use 1 for all points onthe inner boundary and 2 for the outer boundary, while all other pointsare marked by 0. The solver applies the `markers` object to setthe right Dirichlet boundary conditions.The next step is to deform the mesh. Given coordinates $x$, we can mapthese onto a stretched coordinate $\bar x$ by $$\begin{equation}\bar x = a + (b-a)\left({x-a\over b-a}\right)^s,\label{_auto45} \tag{114}\end{equation}$$ where $s$ is a parameter that controls the amount of stretching.The formula above gives a stretching towards $x=a$, while the next onestretches the coordinates towards $x=b$: $$\begin{equation}\bar x = a + (b-a)\left({x-a\over b-a}\right)^{1/s}{\thinspace .}\label{_auto46} \tag{115}\end{equation}$$ The code shown latershows the details of mapping coordinates in a FEniCS mesh object.The final step is to map the rectangle to a part of a hollow cylinder.Mathematically, a point $(x,y)$ in the rectangle is mapped onto $\bar x,\bar y)$ in our final geometry: $$\hat x = \bar x\cos (\Theta \bar y),\quad \hat y = \bar x\sin (\Theta \bar y){\thinspace .}$$ The relevant FEniCS code becomes<jupyter_code># --- Deform mesh ---

# First make a denser mesh towards r=a
x = mesh.coordinates()[:,0]
y = mesh.coordinates()[:,1]

def denser(x, y):
    return [a + (b-a)*((x-a)/(b-a))**s, y]

x_bar, y_bar = denser(x, y)
xy_bar_coor = np.array([x_bar, y_bar]).transpose()
mesh.coordinates()[:] = xy_bar_coor

# Then map onto to a "piece of cake"

def cylinder(r, s):
    return [r*np.cos(Theta*s), r*np.sin(Theta*s)]

x_hat, y_hat = cylinder(x_bar, y_bar)
xy_hat_coor = np.array([x_hat, y_hat]).transpose()
mesh.coordinates()[:] = xy_hat_coor
return mesh, markers<jupyter_output><empty_output><jupyter_text>Fortunately, the solver is independent of all the details of the mesh making.We could also have used the mesh tool `mshr` in FEniCS, but with ourapproach here we have full control of the refinement towards the hole. Solving a problemWe assume that ${\alpha}$ is constant.  Before solving such a specificproblem, it can be wise to scale the problem since it often reducesthe amount of input data in the model. Here, the variation in $u$ istypically $|u_a-u_b|$ so we use that as characteristic pressure. Thecoordinates may be naturally scaled by the bore hole radius, so wehave new, scaled variables $$\bar u = \frac{u-u_a}{u_a-u_b},\quad \bar x = \frac{x}{a},\quad\bar y = \frac{y}{a}{\thinspace .}$$ Now, we expect $\bar u\in [0,1]$, which is a goal of scaling.Inserting this in the problem gives the PDE $$\nabla^2 \bar u = 0$$ in a domain with inner radius 1 and $\bar u=0$, and outer radius $$\beta = \frac{a}{b},$$ with $\bar u = 1$. Our solver can solve this problem by setting`alpha=1`, `u_a=1`, and `u_b=0`. We see that the dimensionlessparameter $\beta$ goes to the mesh and not to the solver.[Figure](fem:varform:fenics:sol1:fig2) shows a solution for $\beta =2$on a mesh with $4\cdot 20\cdot 20 = 1600$ triangles, 25 degree opening,and P1 elements.Switching to higher-order, say P3, is a matter of changing the `degree`parameter that goes to the function `V` in the solver:<jupyter_code>mesh, markers = make_mesh(Theta=25*pi/180, a=1, b=2,
                          nr=20, nt=20, s=1.9)
beta = 2
solver(mesh, markers,
       alpha=1, u_a=1, u_b=0, degree=3, filename='borehole1')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16464040633273896, 'fraction_numerical': 0.01685848597827849, 'mean_word_length': 4.561391910571548, 'pattern_counts': {'":': 0, '<': 78, '<?xml version=': 0, '>': 67, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 2}, 'pii_count': 0, 'substrings_counts': 6, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '530353', 'license': 'MIT', 'n_tokens_mistral': 31287, 'n_tokens_neox': 29946, 'n_words': 10164, 'path': '5- varform-fe.ipynb', 'repo_name': 'mbarzegary/finite-element-intro', 'size': '92535'}