Approximation with continuous and discontinuous finite elements#
We introduced the notion of a projection in The UFL forms where we want to find the best approximation of an expression in a finite element space.
The goal of this section is to approximate
where \(\alpha\) is a pre-defined constant.
We will use ufl.conditional
as explained in the the previous section.
import sys
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx.fem.petsc
import ufl
def h(alpha, mesh: dolfinx.mesh.Mesh):
x = ufl.SpatialCoordinate(mesh)
return ufl.conditional(ufl.lt(x[0], alpha), ufl.cos(ufl.pi * x[0]), -ufl.sin(x[0]))
Reusable projector#
Imagine we want to solve a sequence of such post-processing steps for functions \(f\) and \(g\).
If the mesh is not changing between each projection, the left hand side is constant.
Thus, it would make sense to assemble the matrix once.
Following this, we create a general projector class (class Projector
)
With this class, we can send in any expression written in UFL to the projector, and then generate code for the right hand side assembly, and then solve the linear system. If we use LU factorization, most of the cost will be in the first projection, when the matrix is factorized. This is then cached, so that the solution algorithm is reduced to solving to linear problems; one upper diagonal matrix and one lower diagonal matrix.
Non-aligning discontinuity#
We will start considering the case where \(\alpha\) is not aligned with the mesh. We choose \(\alpha = \frac{\pi}{10}\) and get the following \(h\):
The partial-function
For the next cases, we would like to fix alpha in our function h
, but we want to keep the mesh as a variable.
We can use the partial
function from the functools
module to fix the first argument of a function.
It creates a new object that behaves as a function where mesh
is the only input parameter.
i.e.
from functools import partial
def f(x, y, z):
print(f"{x=}, {y=}, {z=}")
h = partial(f, 3, 4)
h(7) # prints x=3, y=4, z=7
from functools import partial
h_nonaligned = partial(h, np.pi / 10)
Let us now try to use the re-usable projector to approximate this function with a continuous Lagrange space of order 1
Nx = 20
mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, Nx)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
petsc_options = {"ksp_type": "preonly", "pc_type": "lu"}
V_projector = Projector(V, petsc_options=petsc_options)
uh = V_projector.project(h_nonaligned(V.mesh))
We can now repeat the study for a discontinuous Lagrange space of order 1
W = dolfinx.fem.functionspace(mesh, ("DG", 1))
W_projector = Projector(W, petsc_options=petsc_options)
wh = W_projector.project(h_nonaligned(W.mesh))
We compare the two solutions side by side

We observe that both solutions overshoot and undershoot around the discontinuity Let us refine the mesh several times to see if the solution converges



We still see overshoots with either space. This is known as Gibbs phenomenon and is discussed in detail in [Zha22].
Grid-aligned discontinuity#
Next, we choose \(\alpha = 0.2\) and choose grid sizes such that the discontinuity is aligned with a cell boundary.
h_aligned = partial(h, 0.2)



Interpolation of functions and UFL-expressions#
Above we have defined a function that is dependent on the spatial coordinates of ufl
, and it is a purely symbolic expression.
If we want to evaluate this expression, either at a given point or interpolate it into a function space, we need to compile code
similar to the code generated with dolfinx.fem.form
or calling FFCx.
The main difference is that for an expression, there is no summation over quadrature.
To perform this compilation for a given point in the reference cell, we call
mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 7)
compiled_h = dolfinx.fem.Expression(h_aligned(mesh), np.array([0.5]))
We can now evaluate the expression at the point 0.5 in the reference element for any cell (this coordinate is then pushed forward to the given input cell). For instance, we can evaluate this expression in the cell with index 0 with
compiled_h.eval(mesh, np.array([0], dtype=np.int32))
array([[0.97492791]])
Interpolate expressions#
We can also use expressions for post-processing, by interpolating into an appropriate
finite element function space (Q
). To do so, we compile the UFL-expression to be evaluated
at the interpolation points of Q
.
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
u = dolfinx.fem.Function(V)
Q = dolfinx.fem.functionspace(mesh, ("DG", 1))
dudx = ufl.grad(u)[0]
compile_dudx = dolfinx.fem.Expression(dudx, Q.element.interpolation_points())
We populate u
with some data on some part of the domain
left_cells = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, lambda x: x[0] >= 0.3 + 1e-14)
u.interpolate(lambda x: x[0] ** 2, cells0=left_cells)
We can then interpolate dudx
into Q
with
q = dolfinx.fem.Function(Q)
q.interpolate(compile_dudx)
and plot the result


Bibliography#
Shun Zhang. Battling Gibbs phenomenon: On finite element approximations of discontinuous solutions of PDEs. Computers & Mathematics with Applications, 122:35–47, 2022. doi:10.1016/j.camwa.2022.07.014.