\chapter{Examples}\label{sec:examples}\index{examples}The following examples illustrate basic usage of the form language forthe definition of a collection of standard multilinear forms. Weassume that \texttt{dx} has been declared as an integral over theinterior of $\Omega$ and that both \texttt{i} and \texttt{j} havebeen declared as a free \texttt{Index}.The examples presented below can all be found in the subdirectory\texttt{src/demo} of the \ffc{} source tree together with numerousother examples.%------------------------------------------------------------------------------\section{The mass matrix}\index{mass matrix}As a first example, consider the bilinear form corresponding to amass matrix,\begin{equation}  a(v, u) = \int_{\Omega} v \, u \dx,\end{equation}which can be implemented in \ffc{} as follows:\begin{code}element = FiniteElement("Lagrange", "triangle", 1)v = TestFunction(element)u = TrialFunction(element)    a = v*u*dx\end{code}This example is implemented in the file \texttt{Mass.form} in thecollection of demonstration forms included with the \ffc{} sourcedistribution.%------------------------------------------------------------------------------\section{Poisson's equation}\index{Poisson's equation}The bilinear and linear forms form for Poisson's equation,\begin{eqnarray}  a(v, u) &=& \int_{\Omega} \nabla v \cdot \nabla u \dx, \\  L(v)    &=& \int_{\Omega} v \, f \dx,\end{eqnarray}can be implemented as follows:\begin{code}element = FiniteElement("Lagrange", "triangle", 1)v = TestFunction(element)u = TrialFunction(element)f = Function(element)  a = dot(grad(v), grad(u))*dxL = v*f*dx\end{code}Alternatively, index notation can be used to express the scalarproduct:\begin{code}a = D(v, i)*D(u, i)*dx\end{code}This example is implemented in the file \texttt{Poisson.form} in thecollection of demonstration forms included with the \ffc{} sourcedistribution.%------------------------------------------------------------------------------\section{Vector-valued Poisson}\index{vector-valued Poisson}The bilinear and linear forms for a system of (independent) Poissonequations,\begin{eqnarray}  a(v, u) &=& \int_{\Omega} \nabla v : \nabla u \dx, \\  L(v)    &=& \int_{\Omega} v \cdot f \dx,\end{eqnarray}with $v$, $u$ and $f$ vector-valued can be implemented as follows:\begin{code}element = VectorElement("Lagrange", "triangle", 1)  v = TestFunction(element)u = TrialFunction(element)f = Function(element) a = dot(grad(v), grad(u))*dxL = dot(v, f)*dx\end{code}Alternatively, index notation may be used:\begin{code}a = D(v[i], j)*D(u[i], j)*dxL = v[i]*f[i]*dx\end{code}This example is implemented in the file \texttt{PoissonSystem.form} inthe collection of demonstration forms included with the \ffc{} sourcedistribution.%------------------------------------------------------------------------------\section{The strain-strain term of linear elasticity}\index{elasticity}\index{linear elasticity}\index{strain}The strain-strain term of linear elasticity,\begin{equation}  a(v, u) = \int_{\Omega} \epsilon(v) : \epsilon(u) \dx,\end{equation}where\begin{equation}  \epsilon(v) = \frac{1}{2}(\nabla v + (\nabla v)^{\top})\end{equation}can be implemented as follows:\begin{code}element = VectorElement("Lagrange", "tetrahedron", 1)  v = TestFunction(element)u = TrialFunction(element)  def epsilon(v):    return 0.5*(grad(v) + transp(grad(v)))a = dot(epsilon(v), epsilon(u))*dx\end{code}Alternatively, index notation can be used to define the form:\begin{code}a = 0.25*(D(v[i], j) + D(v[j], i))* \         (D(u[i], j) + D(u[j], i))*dx\end{code}This example is implemented in the file \texttt{Elasticity.form} inthe collection of demonstration forms included with the \ffc{} sourcedistribution.%------------------------------------------------------------------------------\section{The nonlinear term of Navier--Stokes}\index{Navier-Stokes}\index{fixed-point iteration}The bilinear form for fixed-point iteration on the nonlinear term ofthe incompressible Navier--Stokes equations,\begin{equation}  a(v, u) = \int_{\Omega} v \cdot ((w \cdot \nabla) u) \dx,\end{equation}with $w$ the frozen velocity from a previous iteration, can beconveniently implemented using index notation as follows:\begin{code}element = FiniteElement("Vector Lagrange", "tetrahedron", 1)  v = TestFunction(element)u = TrialFunction(element)w = Function(element)a = v[i]*w[j]*D(u[i], j)*dx\end{code}This example is implemented in the file \texttt{NavierStokes.form} inthe collection of demonstration forms included with the \ffc{} sourcedistribution.%------------------------------------------------------------------------------\section{The heat equation}\index{heat equation}\index{time-stepping}\index{backward Euler}Discretizing the heat equation,\begin{equation}  \dot{u} - \nabla \cdot (c \nabla u) = f,\end{equation}in time using the $\mathrm{dG}(0)$ method (backward Euler), weobtain the following variational problem for the discrete solution $u_h= u_h(x, t)$: Find $u_h^n = u_h(\cdot, t_n)$ with$u_h^{n-1} = u_h(\cdot, t_{n-1})$ given such that\begin{equation}  \frac{1}{k_n} \int_{\Omega} v \, (u_h^n - u_h^{n-1}) \dx +  \int_{\Omega} c \, \nabla v \cdot \nabla u_h^n \dx =  \int_{\Omega} v \, f^n \dx\end{equation}for all test functions $v$, where $k = t_n - t_{n-1}$ denotes the timestep . In the example below, we implement this variational problemwith piecewise linear test and trial functions, but other choices arepossible (just choose another finite element).Rewriting the variational problem in the standard form $a(v, u_h) =L(v)$ for all $v$, we obtain the following pair of bilinear and linearforms:\begin{eqnarray}  a(v, u_h^n) &=& \int_{\Omega} v \, u_h^n \dx +  k_n \int_{\Omega} c \, \nabla v \cdot \nabla u_h^n \dx, \\  L(v) &=& \int_{\Omega} v \, u_h^{n-1} \dx + k_n \int_{\Omega} v \, f^n \dx,\end{eqnarray}which can be implemented as follows:\begin{code}element = FiniteElement("Lagrange", "triangle", 1)v  = TestFunction(element) # Test functionu1 = TrialFunction(element) # Value at t_nu0 = Function(element)      # Value at t_n-1c  = Function(element)      # Heat conductivityf  = Function(element)      # Heat sourcek  = Constant()             # Time stepa = v*u1*dx + k*c*dot(grad(v), grad(u1))*dxL = v*u0*dx + k*v*f*dx\end{code}%------------------------------------------------------------------------------\section{Mixed formulation of Stokes}\index{Stokes' equations}\index{mixed formulation}\index{Taylor-Hood element}To solve Stokes' equations,\begin{eqnarray}  - \Delta u + \nabla p &=& f, \\  \nabla \cdot u &=& 0,\end{eqnarray}we write the variational problem in standard form$a(v, u) = L(v)$ for all $v$ to obtain the following pair of bilinearand linear forms:\begin{eqnarray}  a((v, q), (u, p)) &=& \int_{\Omega} \nabla v : \nabla u - (\nabla \cdot v) \, p +  q \, (\nabla \cdot u) \dx, \\  L((v, q)) &=& \int_{\Omega} v \cdot f \dx.\end{eqnarray}Using a mixed formulation with Taylor-Hood elements, this can beimplemented as follows:\begin{code}P2 = FiniteElement("Vector Lagrange", "triangle", 2)P1 = FiniteElement("Lagrange", "triangle", 1)TH = P2 + P1(v, q) = TestFunctions(TH)(u, p) = TrialFunctions(TH) f = Function(P2) a = (dot(grad(v), grad(u)) - div(v)*P + q*div(u))*dxL = dot(v, f)*dx\end{code}This example is implemented in the file \texttt{Heat.form} in thecollection of demonstration forms included with the \ffc{} sourcedistribution.%------------------------------------------------------------------------------\section{Mixed formulation of Poisson}\index{mixed Poisson}\index{BDM elements}\index{Brezzi--Douglas--Marini elements}We next consider the following formulation of Poisson's equation as apair of first order equations for $\sigma \in H(\mathrm{div})$and $u \in L_2$:\begin{eqnarray}    \sigma + \nabla u &= 0, \\    \nabla \cdot \sigma &= f.\end{eqnarray}We multiply the two equations by a pair of test functions $\tau$ and$w$ and integrate by parts to obtain the following variationalproblem: Find $(\sigma, u) \in V = H(\mathrm{div}) \times L_2$ such that\begin{equation}  a((\tau, w), (\sigma, u)) = L((\tau, w)) \quad \forall \, (\tau, w) \in V,\end{equation}where\begin{eqnarray}  a((\tau, w), (\sigma, u)) &=& \int_{\Omega} \tau \cdot \sigma - \nabla \cdot \tau \, u + w \nabla \cdot \sigma \dx, \\  L((\tau, w)) &=& \int_{\Omega} w \cdot f \dx.\end{eqnarray}We may implement the corresponding forms in the \ffc{} form languageusing first order BDM $H(\mathrm{div})$-conformingelements for $\sigma$ and piecewise constant $L_2$-conforming elementsfor $u$ as follows:\begin{code}BDM1 = FiniteElement("Brezzi-Douglas-Marini", "triangle", 1)DG0  = FiniteElement("Discontinuous Lagrange", "triangle", 0)element = BDM1 + DG0(tau, w) = TestFunctions(element)(sigma, u) = TrialFunctions(element)f = Function(DG0)a = (dot(tau, sigma) - div(tau)*u + w*div(sigma))*dxL = w*f*dx\end{code}This example is implemented in the file \texttt{MixedPoisson.form} inthe collection of demonstration forms included with the \ffc{} sourcedistribution.%------------------------------------------------------------------------------\section{Poisson's equation with DG elements}\index{Discontinuous Galerkin}We consider again Poisson's equation, but now in an (interior penalty)discontinuous Galerkin formulation: Find $u \in V = L_2$ such that\begin{displaymath}  a(v, u) = L(v) \quad \forall v \in V,\end{displaymath}where\begin{equation}  \begin{split}    a(v, u) &= \int_{\Omega} \nabla v \cdot \nabla u \dx \\    &+ \sum_S \int_S     - \langle \nabla v \rangle \cdot \llbracket u \rrbracket_n    - \llbracket v \rrbracket_n \cdot \langle \nabla u \rangle    + (\alpha/h) \llbracket v \rrbracket_n \cdot \llbracket u \rrbracket_n \dS \\    &+ \int_{\partial\Omega}    - \nabla v \cdot \llbracket u \rrbracket_n - \llbracket v \rrbracket_n \cdot \nabla u    + (\gamma/h) v u \ds \\    L(v) &= \int_{\Omega} v f \dx + \int_{\partial\Omega} v g \ds.  \end{split}\end{equation}The corresponding finite element variational problem for discontinuousfirst order elements may be implemented as follows:\begin{code}DG1 = FiniteElement("Discontinuous Lagrange", "triangle", 1)v = TestFunction(DG1)u = TrialFunction(DG1)f = Function(DG1)g = Function(DG1)n = FacetNormal("triangle")h = MeshSize("triangle")a = dot(grad(v), grad(u))*dx \  - dot(avg(grad(v)), jump(u, n))*dS \  - dot(jump(v, n), avg(grad(u)))*dS \  + alpha/h('+')*dot(jump(v, n), jump(u, n))*dS \  - dot(grad(v), jump(u, n))*ds \  - dot(jump(v, n),\ grad(u))*ds \  + gamma/h*v*u*dsL = v*f*dx + v*g*ds\end{code}This example is implemented in the file \texttt{PoissonDG.form} in thecollection of demonstration forms included with the \ffc{} sourcedistribution.