from mpi4py import MPI
import numpy as np
import os
= os.getcwd()
CMM_DIR "PYTHONPATH"] = f"{CMM_DIR}:{os.environ.get('PYTHONPATH', '')}"
os.environ[
import library.plot
import library.helpers
= library.helpers.make_unique_dir(os.path.join(CMM_DIR, 'ex01')) output_dir
My first FenicsX program
This notebook gives a small introduction to FenicsX, an open-source FEM library. We use the Python API of FenicsX to solve an example PDE problem.
1 Learning goals
- Understand the basic structure of a FenicsX program
- Learn the basic building blocks of modern FEM software: mesh, function space, boundary conditions, weak form, and solvers
2 Poisson model
In this example, we solve the classical Poisson equation for a scalar \(u\),
\[ -\Delta u = f \text{ in } \Omega = \left[0, \,1\right] \times \left[0,\,1\right], \quad u = g \text{ on } \partial\Omega \quad , \]
where \(\Omega\) is the domain of interest, \(f\) is a given scalar source term, and \(g\) is a Dirichlet boundary condition.
3 Structure
Below, we learn how to
- create a mesh for a simple domain
- pick discrete FEM function spaces (which functions do we use to approximate the PDE solution?)
- define the weak form of the PDE
- apply boundary conditions
- assemble and solve the linear system resulting from the discretization
- visualize the solution
4 Let’s go!
First, we import some of the necessary python libraries
Make sure to download and unpack the library.tar.gz file in the same directory as this notebook.
5 Mesh and computational domain
The basis for our discretization is a discrete approximation \(\Omega_h\) of the computational domain.
Here, he have the choice between different discretization types, e.g. quadrilateral, triangular etc.
from dolfinx import mesh
= 5
number_elements_x = 1
number_elements_y = mesh.create_unit_square(
domain =mesh.CellType.quadrilateral
MPI.COMM_WORLD, number_elements_x, number_elements_y, cell_type )
6 The variational problem
A key step in building finite element methods is the weak form of the PDE, upon which the discrete variational problem is based.
In all cases, the weak form is obtained by multiplying the PDE with a test function \(v\) and integrating over the domain \(\Omega\).
Further operations on the resulting integral equations, such as integration by parts, can then be applied. This is often done improve the properties of the method, for example to - reduce the required regularity (we need less derivatives of the solution to exist) and hence allow for lower-order (=cheaper) finite element spaces - introduce symmetry to the (bi-)linear forms, which is favored by many linear solvers - introduce certain natural boundary conditions, e.g. the terms that appear as boundary integrals in the weak form
A first key realization is that the weak form for a given PDE is not unique.
7 A weak form for the Poisson equation
A typical variational problem for the Poisson equation is given by
Find \(u \in V\) such that
\[ \int_\Omega \nabla u_h \cdot \nabla v_h \, dx - \cancel{\int_{\partial\Omega} \mathbf{\nabla} u_h \cdot \mathbf{n} \, v_h \, ds} = \int_\Omega f v_h \, dx \]
for all \(v \in V_{0}\). Here \(\mathbf{n}\) is the outward normal vector on the boundary \(\partial\Omega\) and \(V_{0}\) is the space of functions that vanish on the boundary, hence the boundary integral vanishes.
The discrete variational problem is obtained by replacing the continous domain \(\Omega\) with its discrete approximation \(\Omega_h\) (our mesh) and the continuous trial and test spaces \(V\) and \(V_{0}\) with the discrete finite element spaces \(V_h\) and \(V_{h,0}\), respectively.
Technical note: By default, FenicsX test functions are zero on parts of the boundary where Dirichlet conditions are applied. Since we only use Dirichlet condiitons in this example, we don’t have to worry about the boundary integral term and we omit the \(V_{h,0}\) notation.
8 Our first finite element
The term “finite element” refers to the basis functions defined on our mesh. The idea of the finite element method is to find the best approximation of the solution \(u\) in a finite-dimensional function space \(V_h\).
We thus have to pick a finite element function space \(V_h\) on our mesh.
Typically, the choice of the function space is crucial for the accuracy and convergence of the method (lego block analogy)
Let’s start with a simple linear polynomial function space on the quadrilateral mesh that we created (check out what it looks like here)
# used for plotting
from dolfinx.fem import functionspace, Function, Constant, form
from basix.ufl import element
= domain.topology.cell_name()
mesh_element_name """ type of mesh element, e.g. "quadrilateral" """
= 1
basis_functions_degree """ degree of the finite basis functions """
= element(
finite_element ="CG", cell=domain.topology.cell_name(), degree=basis_functions_degree
family
)= functionspace(mesh=domain, element=finite_element) discrete_fem_space
Lets have a look at our beautiful mesh:
=discrete_fem_space, title="Mesh", figurepath=output_dir, figurename='mesh', show_edges=True) library.plot.mesh(functionspace
image saved in /home/is086873/CMM/ex01/mesh.png
2025-04-23 16:49:19.813 ( 0.802s) [ 14E0A30E5740]vtkXOpenGLRenderWindow.:1416 WARN| bad X server connection. DISPLAY=
=os.path.join(output_dir, 'mesh.png')) library.plot.display_image(image_path
9 From Math to Code: the weak form and UFL
The main advantage of software like FenicsX is that it provides a way to work in a syntax that is very close to the mathematical notation of the weak form. The underlying syntax is called UFL (Unified Form Language) and is a domain-specific language for finite element variational forms.
Let’s see how we can write the weak form in UFL.
from ufl import TestFunction, TrialFunction, dot, ds, dx, grad
= discrete_fem_space
V_h """ Discrete finite element space """
= TrialFunction(V_h)
u_h """ discrete trial function """
= TestFunction(V_h)
v_h """ discrete test function """
= Constant(domain, 0.0)
f """ right-hand side source term """
# Define the weak form lhs == rhs
= dot(grad(u_h), grad(v_h)) * dx
weak_form_lhs = f * v_h * dx weak_form_rhs
10 Boundary conditions
There are a couple of ways to apply boundary conditions in FenicsX. We cover the simplest way to set Dirichlet boundary conditions.
The basic idea is to tell FenicsX where we want to apply the Dirichlet boundary condition in terms of the mesh and then let the software figure out which nodes and facets of the mesh coincide with that geometrical location.
We therefore start by defining functions to tell FenicsX what we consider to be the left and right boundary of the domain (\(x=0\) and \(x=1\)).
# return 1 on the left boundary of the unit square, e.g when x = (x,y)[0] close to 0
def left_boundary_marker(x):
return np.isclose(x[0], 0.0)
# return 1 on the right boundary of the unit square, e.g when x = (x,y)[0] close to 1
def right_boundary_marker(x):
return np.isclose(x[0], 1.0)
We now create two boundary conditions for a constant value of \(g=1\) on the left boundary and \(g=0\) on the right boundary.
from dolfinx.fem import dirichletbc, locate_dofs_topological
from dolfinx.mesh import locate_entities_boundary
= domain.topology.dim
mesh_topology_dim = mesh_topology_dim - 1
facet_geometrical_dimension
= locate_entities_boundary(
facets_on_left_boundary
domain, facet_geometrical_dimension, left_boundary_marker
)= locate_dofs_topological(
dofs_on_left_boundary
V_h, facet_geometrical_dimension, facets_on_left_boundary
)
= Constant(domain, 1.0)
value_left """ value on the left boundary """
= dirichletbc(value_left, dofs_on_left_boundary, V_h)
bc_left
= locate_entities_boundary(
facets_on_right_boundary
domain, facet_geometrical_dimension, right_boundary_marker
)= locate_dofs_topological(
dofs_on_right_boundary
V_h, facet_geometrical_dimension, facets_on_right_boundary
)= Constant(domain, 0.0)
value_right """ value on the right boundary """
= dirichletbc(value_right, dofs_on_right_boundary, V_h)
bc_right
= [bc_left, bc_right]
dirichlet_boundary_conditions """ list of Dirichlet boundary conditions """
' list of Dirichlet boundary conditions '
11 Assemble and solve the linear system
We start by creating a discrete function to store our solution vector. The function lives in the same approximation space as our trial function \(u_h\) and is initialized to zero:
= Function(V_h)
discrete_solution """ discrete solution, a function in the finite element space. Our solution vector will be stored here """
' discrete solution, a function in the finite element space. Our solution vector will be stored here '
Then, we create a writer for the VTK file to store our results on disk.
from dolfinx.io import VTKFile
from library.helpers import make_unique_dir
= os.path.join(output_dir, "simulation.pvd")
vtk_file_abs_path_name = VTKFile(
vtk_writer "w+"
domain.comm, vtk_file_abs_path_name,
)=0.0)
vtk_writer.write_function(discrete_solution, t
print("Writing solution to file " + vtk_file_abs_path_name)
Writing solution to file /home/is086873/CMM/ex01/simulation.pvd
12 This is where the magic happens
We now reap the true benefit of modern FEM backends: we can assemble the linear system and solve it in a single line of code.
from dolfinx.fem.petsc import LinearProblem
= {"ksp_type": "preonly", "pc_type": "lu"}
linear_solver_options
= LinearProblem(
problem
weak_form_lhs,
weak_form_rhs,=dirichlet_boundary_conditions,
bcs= linear_solver_options,
petsc_options
)
= problem.solve()
discrete_solution
# Write the solution to disk
=0.5)
vtk_writer.write_function(discrete_solution, t
=1.0)
vtk_writer.write_function(discrete_solution, t vtk_writer.close()
Ok… something happened … I guess? Let’s see if we can visualize the solution.
13 Visualization
There are different ways to visualize the solution. For a check of the solution from within our python script, we can visualize the solution using pyvista
.
The result is not super pretty but we can check that the solution indeed “looks” like a linear connection between the boundary conditions
library.plot.scalar_field(=discrete_solution,
function="Solution",
title=output_dir,
figurepath='solution',
figurename=True,
show_edges )
image saved in /home/is086873/CMM/ex01/solution.png
=os.path.join(output_dir, 'solution.png')) library.plot.display_image(image_path
14 Post-processing with Paraview
We can also visualize the solution using Paraview. To do that, we open the files written to the following location
print(vtk_file_abs_path_name)
/home/is086873/CMM/ex01/simulation.pvd
15 Verification: Comparing with the exact solution
For the simple 1D Poisson equation, we can easily compare with the analytical solution, which is just a linear connection of the boundary values, e.g.
\[ u_{\text{exact}}(x,y) = 1 - x \quad \text{ for } x \in [0,1], y \in [0,1] \]
from dolfinx.fem import assemble_scalar, form
from petsc4py import PETSc
def u_exact(x):
= np.zeros((1 , x.shape[1]), dtype=PETSc.ScalarType)
values 0] = 1.0 - x[0]
values[return values
= Function(V_h)
u_ex
u_ex.interpolate(u_exact)
= form(dot(discrete_solution - u_ex, discrete_solution - u_ex) * dx)
L2_error
= np.sqrt(domain.comm.allreduce(assemble_scalar(L2_error), op=MPI.SUM))
error_L2 = domain.comm.allreduce(
error_max max(discrete_solution.x.petsc_vec.array - u_ex.x.petsc_vec.array), op=MPI.MAX
np.
)
print("L2 error with respect to the analytical solution: " + str(error_L2))
print("Maximum error at the degrees of freedom: " + str(error_max))
L2 error with respect to the analytical solution: 1.4256072914824157e-16
Maximum error at the degrees of freedom: 4.440892098500626e-16