Example forms

The following examples illustrate basic usage of the form language for the definition of a collection of standard multilinear forms. We assume that dx has been declared as an integral over the interior of $\mathrm{\Omega }$ and that both i and j have been declared as a free Index.

The examples presented below can all be found in the subdirectory demo/ of the UFL source tree together with numerous other examples.

The mass matrix

As a first example, consider the bilinear form corresponding to a mass matrix,

$$a(v,u)={\int }_{\mathrm{\Omega }}vudx,$$

which can be implemented in UFL as follows:

element = FiniteElement("Lagrange", triangle, 1)

v = TestFunction(element)
u = TrialFunction(element)

a = v*u*dx

This example is implemented in the file Mass.py in the collection of demonstration forms included with the UFL source distribution.

Poisson equation

The bilinear and linear forms form for Poisson’s equation,

$$\begin{array}{rl}a(v,u)& ={\int }_{\mathrm{\Omega }}\mathrm{\nabla }v\cdot \mathrm{\nabla }udx,\\ L(v;f)& ={\int }_{\mathrm{\Omega }}vfdx,\end{array}$$

can be implemented as follows:

element = FiniteElement("Lagrange", triangle, 1)

v = TestFunction(element)
u = TrialFunction(element)
f = Coefficient(element)

a = dot(grad(v), grad(u))*dx
L = v*f*dx

Alternatively, index notation can be used to express the scalar product like this:

a = Dx(v, i)*Dx(u, i)*dx

or like this:

a = v.dx(i)*u.dx(i)*dx

This example is implemented in the file Poisson.py in the collection of demonstration forms included with the UFL source distribution.

Vector-valued Poisson

The bilinear and linear forms for a system of (independent) Poisson equations,

$$\begin{array}{rl}a(v,u)& ={\int }_{\mathrm{\Omega }}\mathrm{\nabla }v:\mathrm{\nabla }udx,\\ L(v;f)& ={\int }_{\mathrm{\Omega }}v\cdot fdx,\end{array}$$

with $v$, $u$ and $f$ vector-valued can be implemented as follows:

element = VectorElement("Lagrange", triangle, 1)

v = TestFunction(element)
u = TrialFunction(element)
f = Coefficient(element)

a = inner(grad(v), grad(u))*dx
L = dot(v, f)*dx

Alternatively, index notation may be used like this:

a = Dx(v[i], j)*Dx(u[i], j)*dx
L = v[i]*f[i]*dx

or like this:

a = v[i].dx(j)*u[i].dx(j)*dx
L = v[i]*f[i]*dx

This example is implemented in the file PoissonSystem.py in the collection of demonstration forms included with the UFL source distribution.

The strain-strain term of linear elasticity

The strain-strain term of linear elasticity,

$$a(v,u)={\int }_{\mathrm{\Omega }}ϵ(v):ϵ(u)dx,$$

where

$$ϵ(v)=\frac{1}{2}(\mathrm{\nabla }v+(\mathrm{\nabla }v{)}^{\mathrm{\top }})$$

can be implemented as follows:

element = VectorElement("Lagrange", tetrahedron, 1)

v = TestFunction(element)
u = TrialFunction(element)

def epsilon(v):
    Dv = grad(v)
    return 0.5*(Dv + Dv.T)

a = inner(epsilon(v), epsilon(u))*dx

Alternatively, index notation can be used to define the form:

a = 0.25*(Dx(v[j], i) + Dx(v[i], j))* \
         (Dx(u[j], i) + Dx(u[i], j))*dx

or like this:

a = 0.25*(v[j].dx(i) + v[i].dx(j))* \
         (u[j].dx(i) + u[i].dx(j))*dx

This example is implemented in the file Elasticity.py in the collection of demonstration forms included with the UFL source distribution.

The nonlinear term of Navier–Stokes

The bilinear form for fixed-point iteration on the nonlinear term of the incompressible Navier–Stokes equations,

$$a(v,u;w)={\int }_{\mathrm{\Omega }}(w\cdot \mathrm{\nabla }u)\cdot vdx,$$

with $w$ the frozen velocity from a previous iteration, can be implemented as follows:

element = VectorElement("Lagrange", tetrahedron, 1)

v = TestFunction(element)
u = TrialFunction(element)
w = Coefficient(element)

a = dot(grad(u)*w, v)*dx

alternatively using index notation like this:

a = v[i]*w[j]*Dx(u[i], j)*dx

or like this:

a = v[i]*w[j]*u[i].dx(j)*dx

This example is implemented in the file NavierStokes.py in the collection of demonstration forms included with the UFL source distribution.

The heat equation

Discretizing the heat equation,

$$\dot{u}-\mathrm{\nabla }\cdot (c\mathrm{\nabla }u)=f,$$

in time using the $\mathrm{dG}(0)$ method (backward Euler), we obtain 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

$$\frac{1}{k_n}{\int }_{\mathrm{\Omega }}v(u_h^n-u_h^{n-1})dx+{\int }_{\mathrm{\Omega }}c\mathrm{\nabla }v\cdot \mathrm{\nabla }{u}_h^{n}dx={\int }_{\mathrm{\Omega }}v{f}^{n}dx$$

for all test functions $v$, where $k_n=t_n-t_{n-1}$ denotes the time step. In the example below, we implement this variational problem with piecewise linear test and trial functions, but other choices are possible (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 linear forms:

$$\begin{array}{rl}a(v,u_h^n;c,k)& ={\int }_{\mathrm{\Omega }}v{u}_h^{n}dx+k_n{\int }_{\mathrm{\Omega }}c\mathrm{\nabla }v\cdot \mathrm{\nabla }{u}_h^{n}dx,\\ L(v;u_h^{n-1},f,k)& ={\int }_{\mathrm{\Omega }}v{u}_h^{n-1}dx+k_n{\int }_{\mathrm{\Omega }}v{f}^{n}dx,\end{array}$$

which can be implemented as follows:

element = FiniteElement("Lagrange", triangle, 1)

v  = TestFunction(element)  # Test function
u1 = TrialFunction(element) # Value at t_n
u0 = Coefficient(element)   # Value at t_n-1
c  = Coefficient(element)   # Heat conductivity
f  = Coefficient(element)   # Heat source
k  = Constant("triangle")   # Time step

a = v*u1*dx + k*c*dot(grad(v), grad(u1))*dx
L = v*u0*dx + k*v*f*dx

This example is implemented in the file Heat.py in the collection of demonstration forms included with the UFL source distribution.

Mixed formulation of Stokes

To solve Stokes’ equations,

$$\begin{array}{rl}-\mathrm{\Delta }u+\mathrm{\nabla }p& =f,\\ \mathrm{\nabla }\cdot u& =0,\end{array}$$

we write the variational problem in standard form $a(v,u)=L(v)$ for all $v$ to obtain the following pair of bilinear and linear forms:

$$\begin{array}{rl}a((v,q),(u,p))& ={\int }_{\mathrm{\Omega }}\mathrm{\nabla }v:\mathrm{\nabla }u-(\mathrm{\nabla }\cdot v)p+q(\mathrm{\nabla }\cdot u)dx,\\ L((v,q);f)& ={\int }_{\mathrm{\Omega }}v\cdot fdx.\end{array}$$

Using a mixed formulation with Taylor-Hood elements, this can be implemented as follows:

cell = triangle
P2 = VectorElement("Lagrange", cell, 2)
P1 = FiniteElement("Lagrange", cell, 1)
TH = P2 * P1

(v, q) = TestFunctions(TH)
(u, p) = TrialFunctions(TH)

f = Coefficient(P2)

a = (inner(grad(v), grad(u)) - div(v)*p + q*div(u))*dx
L = dot(v, f)*dx

This example is implemented in the file Stokes.py in the collection of demonstration forms included with the UFL source distribution.

Mixed formulation of Poisson

We next consider the following formulation of Poisson’s equation as a pair of first order equations for $\sigma \in H(\mathrm{div})$ and $u\in L^2$:

$$\begin{array}{rl}\sigma +\mathrm{\nabla }u& =0,\\ \mathrm{\nabla }\cdot \sigma & =f.\end{array}$$

We multiply the two equations by a pair of test functions $\tau$ and $w$ and integrate by parts to obtain the following variational problem: Find $(\sigma ,u)\in V=H(\mathrm{div})\times L^2$ such that

$$a((\tau ,w),(\sigma ,u))=L((\tau ,w))\mathrm{\forall }(\tau ,w)\in V,$$

where

$$\begin{array}{rl}a((\tau ,w),(\sigma ,u))& ={\int }_{\mathrm{\Omega }}\tau \cdot \sigma -\mathrm{\nabla }\cdot \tau u+w\mathrm{\nabla }\cdot \sigma dx,\\ L((\tau ,w);f)& ={\int }_{\mathrm{\Omega }}w\cdot fdx.\end{array}$$

We may implement the corresponding forms in our form language using first order BDM H(div)-conforming elements for $\sigma$ and piecewise constant $L^2$-conforming elements for $u$ as follows:

cell = triangle
BDM1 = FiniteElement("Brezzi-Douglas-Marini", cell, 1)
DG0  = FiniteElement("Discontinuous Lagrange", cell, 0)

element = BDM1 * DG0

(tau, w) = TestFunctions(element)
(sigma, u) = TrialFunctions(element)

f = Coefficient(DG0)

a = (dot(tau, sigma) - div(tau)*u + w*div(sigma))*dx
L = w*f*dx

This example is implemented in the file MixedPoisson.py in the collection of demonstration forms included with the UFL source distribution.

Poisson equation with DG elements

We consider again Poisson’s equation, but now in an (interior penalty) discontinuous Galerkin formulation: Find $u\in V=L^2$ such that

$$a(v,u)=L(v)\mathrm{\forall }v\in V,$$

where

$$\begin{array}{rl}a(v,u;h)& ={\int }_{\mathrm{\Omega }}\mathrm{\nabla }v\cdot \mathrm{\nabla }udx\\ & +\sum _S{\int }_S-⟨\mathrm{\nabla }v⟩\cdot [[u]{]}_n-[[v]{]}_n\cdot ⟨\mathrm{\nabla }u⟩+(\alpha /h)[[v]{]}_n\cdot [[u]{]}_{n}dS\\ & +{\int }_{\partial \mathrm{\Omega }}-\mathrm{\nabla }v\cdot [[u]{]}_n-[[v]{]}_n\cdot \mathrm{\nabla }u+(\gamma /h)vuds\\ L(v;f,g)& ={\int }_{\mathrm{\Omega }}vfdx+{\int }_{\partial \mathrm{\Omega }}vgds.\end{array}$$

The corresponding finite element variational problem for discontinuous first order elements may be implemented as follows:

cell = triangle
DG1 = FiniteElement("Discontinuous Lagrange", cell, 1)

v = TestFunction(DG1)
u = TrialFunction(DG1)

f = Coefficient(DG1)
g = Coefficient(DG1)
h = 2.0*Circumradius(cell)
alpha = 1
gamma = 1

a = dot(grad(v), grad(u))*dx \
  - dot(avg(grad(v)), jump(u))*dS \
  - dot(jump(v), avg(grad(u)))*dS \
  + alpha/h('+')*dot(jump(v), jump(u))*dS \
  - dot(grad(v), jump(u))*ds \
  - dot(jump(v), grad(u))*ds \
  + gamma/h*v*u*ds
L = v*f*dx + v*g*ds

This example is implemented in the file PoissonDG.py in the collection of demonstration forms included with the UFL source distribution.

The Quadrature family

We consider here a nonlinear version of the Poisson’s equation to illustrate the main point of the Quadrature finite element family. The strong equation looks as follows:

$$-\mathrm{\nabla }\cdot (1+u^2)\mathrm{\nabla }u=f.$$

The linearised bilinear and linear forms for this equation,

$$\begin{array}{rl}a(v,u;u_0)& ={\int }_{\mathrm{\Omega }}(1+u_0^2)\mathrm{\nabla }v\cdot \mathrm{\nabla }udx+{\int }_{\mathrm{\Omega }}2u_{0}u\mathrm{\nabla }v\cdot \mathrm{\nabla }{u}_{0}dx,\\ L(v;u_0,f)& ={\int }_{\mathrm{\Omega }}vfdx-{\int }_{\mathrm{\Omega }}(1+u_0^2)\mathrm{\nabla }v\cdot \mathrm{\nabla }{u}_{0}dx,\end{array}$$

can be implemented in a single form file as follows:

element = FiniteElement("Lagrange", triangle, 1)

v = TestFunction(element)
u = TrialFunction(element)
u0 = Coefficient(element)
f = Coefficient(element)

a = (1+u0**2)*dot(grad(v), grad(u))*dx + 2*u0*u*dot(grad(v), grad(u0))*dx
L = v*f*dx - (1+u0**2)*dot(grad(v), grad(u0))*dx

Here, $u_0$ represents the solution from the previous Newton-Raphson iteration.

The above form will be denoted REF1 and serves as our reference implementation for linear elements. A similar form (REF2) using quadratic elements will serve as a reference for quadratic elements.

Now, assume that we want to treat the quantities $C=(1+u_0^2)$ and ${\sigma }_0=(1+u_0^2)\mathrm{\nabla }{u}_0$ as given functions (to be computed elsewhere). Substituting into the bilinear and linear forms, we obtain

$$\begin{array}{rl}a(v,u)& ={\int }_{\mathrm{\Omega }}\text{C}\mathrm{\nabla }v\cdot \mathrm{\nabla }udx+{\int }_{\mathrm{\Omega }}2u_{0}u\mathrm{\nabla }v\cdot \mathrm{\nabla }{u}_{0}dx,\\ L(v;{\sigma }_0,f)& ={\int }_{\mathrm{\Omega }}vfdx-{\int }_{\mathrm{\Omega }}\mathrm{\nabla }v\cdot {\sigma }_{0}dx.\end{array}$$

Then, two additional forms are created to compute the tangent C and the gradient of $u_0$. This situation shows up in plasticity and other problems where certain quantities need to be computed elsewhere (in user-defined functions). The three forms using the standard FiniteElement (linear elements) can then be implemented as

# NonlinearPoisson.py
from ufl import *

element = FiniteElement("Lagrange", triangle, 1)
DG = FiniteElement("Discontinuous Lagrange", triangle, 0)
sig = VectorElement("Discontinuous Lagrange", triangle, 0)

v    = TestFunction(element)
u    = TrialFunction(element)
u0   = Coefficient(element)
C    = Coefficient(DG)
sig0 = Coefficient(sig)
f    = Coefficient(element)

a = v.dx(i)*C*u.dx(i)*dx + v.dx(i)*2*u0*u*u0.dx(i)*dx
L = v*f*dx - dot(grad(v), sig0)*dx

and

# Tangent.py
from ufl import *

element = FiniteElement("Lagrange", triangle, 1)
DG = FiniteElement("Discontinuous Lagrange", triangle, 0)

v = TestFunction(DG)
u = TrialFunction(DG)
u0= Coefficient(element)

a = v*u*dx
L = v*(1.0 + u0**2)*dx

and

# Gradient.py
from ufl import *

element = FiniteElement("Lagrange", triangle, 1)
DG = VectorElement("Discontinuous Lagrange", triangle, 0)

v = TestFunction(DG)
u = TrialFunction(DG)
u0 = Coefficient(element)

a = dot(v, u)*dx
L = dot(v, (1.0 + u0**2)*grad(u0))*dx

The three forms can be implemented using the QuadratureElement in a similar fashion in which only the element declaration is different:

# QE1NonlinearPoisson.py
from ufl import *

element = FiniteElement("Lagrange", triangle, 1)
QE = FiniteElement("Quadrature", triangle, 2)
sig = VectorElement("Quadrature", triangle, 2)

and

# QE1Tangent.py
from ufl import *

element = FiniteElement("Lagrange", triangle, 1)
QE = FiniteElement("Quadrature", triangle, 2)

and

# QE1Gradient.py
from ufl import *

element = FiniteElement("Lagrange", triangle, 1)
QE = VectorElement("Quadrature", triangle, 2)

Note that we use two points when declaring the QuadratureElement. This is because the RHS of Tangent.form is second order and therefore we need two points for exact integration. Due to consistency issues, when passing functions around between the forms, we also need to use two points when declaring the QuadratureElement in the other forms.

Typical values of the relative residual for each Newton iteration for all three approaches are shown in the table below. It is to be noted that the convergence rate is quadratic as it should be for all three methods.

Relative residuals for each approach for linear elements:

Iteration REF1      FE1      QE1
========= ====      ===      ===
1         6.3e-02   6.3e-02  6.3e-02
2         5.3e-04   5.3e-04  5.3e-04
3         3.7e-08   3.7e-08  3.7e-08
4         2.9e-16   2.9e-16  2.5e-16

However, if quadratic elements are used to interpolate the unknown field $u$, the order of all elements in the above forms is increased by 1. This influences the convergence rate as seen in the table below. Clearly, using the standard FiniteElement leads to a poor convergence whereas the QuadratureElement still leads to quadratic convergence.

Relative residuals for each approach for quadratic elements:

Iteration REF2      FE2      QE2
========= ====      ===      ===
1         2.6e-01   3.9e-01  2.6e-01
2         1.1e-02   4.6e-02  1.1e-02
3         1.2e-05   1.1e-02  1.6e-05
4         1.1e-11   7.2e-04  9.1e-09

More examples

Feel free to send additional demo form files for your favourite PDE to the UFL mailing list.