The standard way of compiling code with DOLFINx#
In this section, we will focus on the approach most users use to interact with UFL, FFCx and basix. Here we will start by creating the domain we want to solve a problem on. In this case, we will use a unit square
from mpi4py import MPI
import dolfinx
import ufl
import numpy as np
import scipy
N = 10
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N)
tdim = mesh.topology.dim
Problem specification: Non-linear Poisson.#
Next, let’s consider the problem
where \(p(u)=1+u^2\). We choose to use a first order Lagrange space for the unknown.
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
Alternative declaration of finite elements
When working in DOLFINx, we can also supply the element information directly to the function space
with a tuple (family, degree)
, or a triplet (family, degree, shape)
, i.e. if we want a M dimensional
Lagrange 5th order vector space, one could specify this as ("Lagrange", 5, (M, ))
.
We do as in the previous section, and define a manufactured solution
def u_exact(module, x):
return module.sin(2 * module.pi * x[1]) + x[0] * x[1] ** 2
x = ufl.SpatialCoordinate(mesh)
u_ex = u_exact(ufl, x)
As the problem is non-linear we need to prepare for solving this with a Newton-solver.
uh = dolfinx.fem.Function(V)
We define the residual
def p(u):
return 1 + u**2
f = -ufl.div(p(u_ex) * ufl.grad(u_ex))
v = ufl.TestFunction(V)
F = ufl.inner(p(uh) * ufl.grad(uh), ufl.grad(v)) * ufl.dx - ufl.inner(f, v) * ufl.dx
gh = dolfinx.fem.Function(V)
gh.interpolate(lambda x: u_exact(np, x))
mesh.topology.create_entities(tdim - 1)
mesh.topology.create_connectivity(tdim - 1, tdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, tdim - 1, boundary_facets)
bc = dolfinx.fem.dirichletbc(gh, boundary_dofs)
bcs = [bc]
We compute the Jacobian of the system using UFL.
J = ufl.derivative(F, uh)
Now that we have associated the forms with the discrete problem already, we use
dolfinx.fem.form
residual = dolfinx.fem.form(F)
jacobian = dolfinx.fem.form(J)
A = dolfinx.fem.create_matrix(jacobian)
b = dolfinx.fem.create_vector(residual)
We are now ready to use these forms to solve the variational problem.
As we want \(u_{k+1}=g\), but we do not know if \(u_k=g\), we need to take this into account when assembling the right hand side for the Jacobian equation.
We will go through how we apply this boundary condition using lifting
Lifting#
Let us split the degrees of freedom into two disjoint sets, \(u_d\), and \(u_{bc}\), and set up the corresponding linear system
In the identity row approach, we set the rows corresponding to the Dirichlet conditions to the identity row and set the appropriate dofs on the right hand side to contain the Dirichlet values:
where \(g\) is the vector satisfying the various Dirichlet conditions. We can now reduce this to a smaller system
which is symmetric if \(A_{d,d}\) is symmetric. However, we do not want to remove the degrees of freedom from the sparse matrix, as it makes the matrix less adaptable for varying boundary conditions, so we set up the system