# The UFL forms#

In this section we will show how to set up the following problem

for some known \(g\). This problem is known as the \(L^2\)-projection of \(g\) into \(u\). We can observe this by looking at what solving this problem implies: To find the minimum, we need to find where \(\frac{\mathrm{d}G}{\mathrm{d}u}\) is \(0\). We define the residual \(F\) as

We could compute this by hand and get: Find \(u_h\in V\) such that

Thus we can define the bi-linear form \(a(u, v)\) and linear form \(L(v)\)

We start by defining and abstract domain and a finite element

```
import basix.ufl
import ufl
cell = "triangle"
c_el = basix.ufl.element("Lagrange", cell, 1, shape=(2,))
domain = ufl.Mesh(c_el)
el = basix.ufl.element("Lagrange", cell, 2)
```

## The functions space#

We define the function space \(V\), which will depend on the domain and the finite element

```
V = ufl.FunctionSpace(domain, el)
```

Now we are ready to create the formulations above.
We can start with defining our unknown \(u_h\).
A function in UFL is described as a `ufl.Coefficient`

, this means that it can be written as
\(u(x)=\sum_{i=0}^N u_i \phi_i(x)\).

```
uh = ufl.Coefficient(V)
```

## Spatial derivatives#

If we want to compute the derivative in the \(i\)th spatial direction, we can call `ufl.Coefficient.dx`

```
duh_dx = uh.dx(0)
duh_dy = uh.dx(1)
```

If we want the gradient, we can call `ufl.grad(operator)`

```
grad_uh = ufl.grad(uh)
```

We can also use `ufl.as_vector`

to create our own vectors, for instance

```
alt_grad_uh = ufl.as_vector((duh_dx, duh_dy))
```

## Representing \(g\)#

Next, we define our \(g\). We will consider four different cases for \(g\):

\(g\) is a constant over the whole domain

\(g\) is an expression depending on the spatial coordinates of the domain, i.e. \(g=g(x,y,z)\)

\(g\) is a function form another finite element function space

\(g\) is a mixture of the above, combined with spatial derivatives

For case 1., we create a symbolic representation of a constant

```
g = ufl.Constant(domain)
```

For case 2., we use `ufl.SpatialCoordinate`

to create a symbolic representation of the coordinates
in the domain \(\Omega\)

```
x = ufl.SpatialCoordinate(domain)
```

`x`

is here a tuple of 2 components, representing the `x`

and `y`

coordinate in 2D, i.e. `x[0]`

represents the x-coordinate, while `x[1]`

represent the y-coordinate.
The unified form language consist of many
geometrical expressions, for instance `sin`

`pi`

and `cos`

.

```
g = ufl.sin(x[0]) + ufl.cos(ufl.pi * x[1])
```

As the third case is easy to represent, we move directly to the fourth, where we choose \(g=\nabla \cdot f where \)f$ in a vector DG-3 space

```
el_f = basix.ufl.element("DG", cell, 3, shape=(2,))
Q = ufl.FunctionSpace(domain, el_f)
f = ufl.Coefficient(Q)
g = ufl.div(f)
```

## The functional#

Now that we have defined \(u_h\) and \(g\) we can define the functional \(G\)

```
G = (uh - g) ** 2 * ufl.dx
```

The integration measure

In the code above, we have used `ufl.dx`

to indicate that we want to integrate over the cells in our domain.

An alternative definition is

```
dx = ufl.Measure("dx", domain=domain)
```

We will come back to why this can be smart later.

## Symbolic differentiation#

Now that we have our \(G\), we can use `ufl.derivative`

to compute the derivative of \(G\) with respect to a coefficient.
As we want to differentiate with respect to all functions \(\delta u \in V\), we define a `ufl.TestFunction`

```
du = ufl.TestFunction(V)
F = ufl.derivative(G, uh, du)
```

### Non-linear functional#

What if our functional had not been a projection, but something non-linear? How would we set up a symbolic representation of the operators needed to solve the problem? We could then compute the derivative of the residual, known as the Jacobian (not to be confused with the Jacobian of the mapping between the reference and physical cell).

We could use a Newton method to solve the problem in an iterative fashion: Given \(u_{k}\) solve

where \( H= \frac{\mathrm{d}^2G}{\mathrm{d}u^2}[\delta u, \delta v]\) We obtain \(J\) with

```
dv = ufl.TrialFunction(V)
J = ufl.derivative(F, uh, dv)
```

Special case for quadratic functional

For a quadratic functonal, the residual is linear, and thus the Newton method converges in one iteration

```
forms = [G, J, F]
```

## Further analysis of the variational form#

We next consider the steps we would have to implement if we wanted to solve this problem by hand. For illustrative purposes, we choose \(g=\frac{f}{h}\) where \(f\) and \(h\) are two known functions in respectve finite element spaces \(Q\), \(T\), where both \(Q\) and \(T\) uses the scalar valued elements.

We next use the map \(F_K:K_{ref}\mapsto K\) to map the integrals over each cell in the domain back to the reference cell.

where \(K\) is each element in the physical space, \(J_K\) the Jacobian of the mapping.

Next, we can insert the expansion of \(u, v, f, h\) into the formulation: \(u=\sum_{i=0}^{\mathcal{N}}u_i\phi_i(x)\qquad v=\sum_{i=0}^{\mathcal{N}}v_i\phi_i(x)\qquad f=\sum_{k=0}^{\mathcal{M}}f_k\psi_k(x)\qquad g=\sum_{l=0}^{\mathcal{T}}g_l\varphi_l(x)\)

and identify the matrix system \(Au=b\), where

Warning

Next, one can choose an appropriate quadrature rule with points and weights, include the correct mapping/restrictions of degrees of freedom for each cell. All of this becomes quite tedious and error prone work, and has to be repeated for every variational form!

```
Q = ufl.FunctionSpace(domain, basix.ufl.element("Lagrange", cell, 3))
f = ufl.Coefficient(Q)
T = ufl.FunctionSpace(domain, basix.ufl.element("DG", cell, 1))
h = ufl.Coefficient(T)
g = f / g
v = ufl.TestFunction(V)
L = f / g * v * ufl.dx
```

We can use UFL to analyse the linear form above

```
pulled_back_L = ufl.algorithms.compute_form_data(
L,
do_apply_function_pullbacks=True,
do_apply_integral_scaling=True,
do_apply_geometry_lowering=True,
preserve_geometry_types=(ufl.classes.Jacobian,),
)
print(pulled_back_L.integral_data[0])
```

```
IntegralData over domain(cell, ('otherwise',)) with integrals:
{ weight * |(J[0, 0] * J[1, 1] + -1 * J[0, 1] * J[1, 0])| * (reference_value(w_2)) / (reference_value(w_2)) / (sum_{i_0} ({ A | A_{i_{10}, i_{11}} = sum_{i_{12}} ({ A | A_{i_{13}, i_{14}} = ([
[J[1, 1], -1 * J[0, 1]],
[-1 * J[1, 0], J[0, 0]]
])[i_{13}, i_{14}] / (J[0, 0] * J[1, 1] + -1 * J[0, 1] * J[1, 0]) })[i_{12}, i_{11}] * (reference_grad(reference_value(w_1)))[i_{10}, i_{12}] })[i_0, i_0] ) * (reference_value(v_0)) } * dx(<Mesh #0>[('otherwise',)], {'estimated_polynomial_degree': 10})
and metadata:
{metadata}
```

## Inner products and dot products#

When working with vectors and tensors, we often want to compute the inner product or dot product
between two values `u`

and `v`

.

```
V = ufl.FunctionSpace(domain, basix.ufl.element("N1curl", cell, 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
```

The dot product of two vectors \(u\) and \(v\) is defined as `dot(u,v)=u[i]v[i]`

using Einstein summation notation

```
i = ufl.Index()
f = u[i] * v[i] * dx
f_equivalent = ufl.dot(u, v) * dx
```

For two tensors `r`

, `s`

we have that

The inner product is a contraction over all axes

where \(s_{ij}^*\) is the complex conjugate of \(s_{ij}\).

How would the above change if we choose “N1curl” elements as the receiving space?

One would have to add in the covariant Piola mapping to map from the physical to reference basis function.

Can you get UFL to do the computations for you? (Expand for hint)

Try setting up the bi-linear form `a`

for your problem and use `ufl.algorithms.compute_form_data`

to see what expression you get out.
Try turning various options on and off and see how the result looks like

Expand the below to investigate the solution

## Show code cell source

```
V = ufl.FunctionSpace(domain, basix.ufl.element("N1curl", cell, 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(u, v) * ufl.dx
pulled_back_a = ufl.algorithms.compute_form_data(
a,
do_apply_function_pullbacks=True,
do_apply_integral_scaling=False,
do_apply_geometry_lowering=False,
preserve_geometry_types=(ufl.classes.Jacobian,),
)
print(pulled_back_a.integral_data[0])
```

## Show code cell output

```
IntegralData over domain(cell, ('otherwise',)) with integrals:
{ sum_{i_0} ({ A | A_{i_{19}} = sum_{i_{20}} K[i_{20}, i_{19}] * (reference_value(v_0))[i_{20}] })[i_0] * ({ A | A_{i_{17}} = sum_{i_{18}} K[i_{18}, i_{17}] * (reference_value(v_1))[i_{18}] })[i_0] } * dx(<Mesh #0>[('otherwise',)], {'estimated_polynomial_degree': 2})
and metadata:
{metadata}
```

## Quadrature rule#

The next step in the analysis of (2) is to choose a suitable quadrature rule, and expand the integral as a sum of evaluations of the basis functions at points.

As you might have spotted in the previous exercise, UFL can estimate want degree you should use for your quadrature rule. We can get this number explicitly with the following code.

```
pulled_back_a = ufl.algorithms.compute_form_data(
a,
do_apply_function_pullbacks=True,
do_apply_integral_scaling=True,
do_apply_geometry_lowering=True,
preserve_geometry_types=(ufl.classes.Jacobian,),
)
print(list(itg.metadata() for itg in pulled_back_a.integral_data[0].integrals))
```

```
[{'estimated_polynomial_degree': 2}]
```

Would the estimate change if we use a different cell for the domain?

Yes, choosing a quadrilateral cell will increase the quadrature degree, as the map from the reference element is non-affine.

Expand to see example

## Show code cell source

```
cell = "quadrilateral"
domain = ufl.Mesh(basix.ufl.element("Lagrange", cell, 1, shape=(2,)))
V = ufl.FunctionSpace(domain, basix.ufl.element("N1curl", cell, 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(u, v) * ufl.dx
pulled_back_a = ufl.algorithms.compute_form_data(
a,
do_apply_function_pullbacks=True,
do_apply_integral_scaling=True,
do_apply_geometry_lowering=True,
preserve_geometry_types=(ufl.classes.Jacobian,),
)
print(list(itg.metadata() for itg in pulled_back_a.integral_data[0].integrals))
```

## Show code cell output

```
[{'estimated_polynomial_degree': 4}]
```

Would the estimate change if use a higher order triangle?

Yes, as above, the Jacobian becomes non-constant as the map is non-affine.

## Show code cell source

```
cell = "triangle"
domain = ufl.Mesh(basix.ufl.element("Lagrange", cell, 2, shape=(2,)))
V = ufl.FunctionSpace(domain, basix.ufl.element("N1curl", cell, 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(u, v) * ufl.dx
pulled_back_a = ufl.algorithms.compute_form_data(
a,
do_apply_function_pullbacks=True,
do_apply_integral_scaling=True,
do_apply_geometry_lowering=True,
preserve_geometry_types=(ufl.classes.Jacobian,),
)
print(list(itg.metadata() for itg in pulled_back_a.integral_data[0].integrals))
```

## Show code cell output

```
[{'estimated_polynomial_degree': 4}]
```

We can fix the number of quadraure points by setting the `quadrature_degree`

in the `Measure`

```
dx = ufl.Measure("dx", domain=domain, metadata={"quadrature_degree": 4})
```

Warning

Variational crimes
Reducing the accuracy of the integration by lowering the quadrature rule is considered to be a
**variational crime** [Suli12] (Chapter 3.4) and should be done with caution.

## References#

Endre Süli. Lecture notes on finite element methods for partial differential equations. 2012. URL: https://people.maths.ox.ac.uk/suli/fem.pdf.