Accessing low-level FEM details in FEniCSx for teaching purposes

General
1

Hi all,

I am considering using FEniCSx in a nonlinear FEM course. For didactic reasons, I would like to first show the “classical” FEM structure:

  • explicit loop over cells,
  • explicit loop over quadrature points,
  • evaluation of shape functions and their gradients at a given quadrature point in a given cell,
  • computation of the cell Jacobian,
  • manual assembly of local element tensors into a global PETSc matrix.

My questions:

  1. Is it straightforward in FEniCSx (Python API) to access these low-level ingredients (cell geometry, quadrature points, basis tabulation via Basix, manual insertion into global matrices)?
  2. Are there references or minimal examples that showcase this more “textbook-style” FEM implementation in FEniCSx?

The goal is purely pedagogical, to first expose the algorithmic structure, and then show how FEniCSx automates it.

Thanks in advance!

2

I’m sure @dokken is going to give you a much nicer answer and better context, but just to help you on your way already, here is a workshop he wrote that also starts with the basics of FEM, basis functions, element mappings, integration points, etc: Motivation — FEniCS Workshop

1 Like
3

As @Stein already pointed out, this is possible with the components described in the workshop notes.
Below is an example where I combine all the features described in the workshop notes to do assembly of a stiffness matrix over a manifold (line embedded in 2D)
by manual assembly into a petsc matrix.
Note that this code also runs in parallel:

"""
Handcoded assembly of a stiffness matrix using basix for tabulaton of basis functions and
DOLFINx for dofmap handling.

Copyright (2026) Jørgen S. Dokken
SPDX-License-Identifier:    MIT
"""
from mpi4py import MPI
from petsc4py import PETSc

import dolfinx.fem.petsc
import basix
import ufl
import numpy as np


# Start by creating a mesh consistig of nodes and the cells, where `cells[i]` indicates which nodes in `nodes`
# is in the i-th cell. For example if `cells[i] = [j, k]`, then the i-th cells consists of the points `points[j], points[k]`
if MPI.COMM_WORLD.rank == 0:
    cells = np.array([[0,1], [1,2], [2, 3], [3, 4]], dtype=np.int64)
    nodes = np.array([[0.1, 0.0], [0.2, 0.1], [0.4, 0.2], [0.5,0.4],[0.6, 0.2]], dtype=np.float64)
else:
    # On any other process than the rank 0 process, we start with no cells, and let
    # DOLFINx distribute the cells from rank 0 across all processes
    cells = np.zeros((0, 2), dtype=np.int64)
    nodes = np.zeros((0, 2), dtype=np.float64)

# Create the distributed DOLFINx mesh
c_el = basix.ufl.element("Lagrange", "interval", 1, shape=(nodes.shape[1], ))
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells=cells, e=ufl.Mesh(c_el), x=nodes)

# Note that we have chosen to use intervals embedded in R^2, i.e. a manifold surface, thus
# the geometric dimesion of the mesh is 2, while the topological dimension is 1 (intervals)

gdim = mesh.geometry.dim
tdim = mesh.topology.dim

# The number of cells owned by this process can be extracted from the index map
num_cells = mesh.topology.index_map(tdim).size_local

# We create the function space for our unknown, pick whatever degree you'd like
degree = 3
el = basix.ufl.element("Lagrange", mesh.basix_cell(), degree)
V = dolfinx.fem.functionspace(mesh, el)

# To get exact integration of straight intervals for a stiffness matrix, we pick a
# sufficiently high quadrature degree

q_deg = 2*(degree - 1)
q_pts, q_weights = basix.make_quadrature(mesh.basix_cell(), q_deg)

# Next, we will tabulate the basis functions required to compute the jacobian
# and gradientes on the space chosen for the unknown on the reference element.

# Compute Jacobians on reference elements
# The first entry to tabulate is the number of derivatives (we skip the basis values phi_i(x)...
# by using [1:])
coordinate_basis_derivatives = c_el.sub_elements[0].tabulate(1, q_pts)[1:]

# Compute basis functions and derivatives on reference element
# Similar procedure to what we saw above
el_basis_derivatives = el.tabulate(1, q_pts)[1:] # shape: (derivative direction, num_pts, num_basis_funcs)
num_basis_funcs = el_basis_derivatives.shape[2]

# Create local matrix to populate in assembly
A_loc = np.zeros((num_basis_funcs, num_basis_funcs), dtype=np.float64)

# Create a petsc matrix where we allocate the potential non-zero entries by looping through the dofmap and inserting them 
# into a sparsity pattern
pattern = dolfinx.fem.SparsityPattern(mesh.comm, [V.dofmap.index_map, V.dofmap.index_map], [V.dofmap.index_map_bs, V.dofmap.index_map_bs])
for i in range(num_cells):
    local_dofs = V.dofmap.cell_dofs(i)
    pattern.insert(local_dofs, local_dofs)
pattern.finalize()

# Create matrix
A = dolfinx.cpp.la.petsc.create_matrix(mesh.comm, pattern)
A.zeroEntries()

# Assemble over cells local to process
for i in range(num_cells):
    # Zero local tensor
    A_loc[:,:] = 0
    
    # Extract cell nodal coordinates
    cell_node_indices = mesh.geometry.dofmap[i]
    local_mesh_geometry = mesh.geometry.x[cell_node_indices, :mesh.geometry.dim]

    # Compute jacobian for cell
    dphi_dxi = []
    for j in range(tdim):
        dphi_dxi.append(coordinate_basis_derivatives[j] @ local_mesh_geometry)
    jacobian = np.transpose(np.hstack(dphi_dxi).reshape(q_pts.shape[0], tdim, gdim), (0, 2, 1))

    # Loop over quadrature points
    for j in range(q_pts.shape[0]): 
        # Since we are working on a manifold we need to compute the Moore-Penrose pseudo-inverse and
        # Gram determinant rather than the standard inverse and determinant

        # Moore-Penrose pseudo-inverse   
        jacobian_T = jacobian[j].T
        jacobian_pinv = np.linalg.pinv(jacobian_T)
        # Gram determinant
        jac_det = np.sqrt(np.linalg.det(jacobian_T@jacobian[j]))

        # Insert into local tensor
        for k in range(num_basis_funcs):
            for l in range(num_basis_funcs):
                dphi_k = el_basis_derivatives[:, j, k]
                dphi_l = el_basis_derivatives[:, j, l]
                A_loc[k, l] += q_weights[j] * np.dot(jacobian_pinv @  dphi_k, jacobian_pinv @ dphi_l) * jac_det
    
    # Insert local matrix into global matrix
    cell_dofs = V.dofmap.cell_dofs(i)
    A.setValuesLocal(cell_dofs, cell_dofs, A_loc, PETSc.InsertMode.ADD_VALUES)
# Accumulate contributions from all processors
A.assemble()


# Compare with solution made with UFL
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
A_ref = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(ufl.inner(ufl.grad(u), ufl.grad(v))*ufl.dx))
A_ref.assemble()

if MPI.COMM_WORLD.size == 0:
    np.testing.assert_allclose(A_ref.to_dense(), A[:,:])

D = A - A_ref
PETSc.Sys.Print(A.norm(2), A_ref.norm(2), D.norm(2))
assert np.isclose(D.norm(2), 0)
3 Likes