It is quite trivial to port these demos to DOLFINx, as you have demos for mixed elements: dolfinx/demo_stokes.py at main · FEniCS/dolfinx · GitHub

```
import dolfinx
from mpi4py import MPI
import ufl
import numpy as np
# Create mesh
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 25, 25)
# Define function spaces and mixed (product) space
BDM = ufl.FiniteElement("BDM", mesh.ufl_cell(), 1)
DG = ufl.FiniteElement("DG", mesh.ufl_cell(), 0)
el = ufl.MixedElement([BDM ,DG])
W = dolfinx.fem.FunctionSpace(mesh, el)
# Define trial and test functions
(sigma, u) = ufl.TrialFunctions(W)
(tau, v) = ufl.TestFunctions(W)
# Define source function
x = ufl.SpatialCoordinate(mesh)
f = 10*ufl.exp(-((x[0] - 0.5)**2 + (x[1] - 0.5)**2) / 0.02)
# Define variational form
a = (ufl.dot(sigma, tau) + ufl.div(tau)*u + ufl.div(sigma)*v)*ufl.dx
L = - f*v*ufl.dx
# Define function G such that G \cdot n = g
def bottom_boundary(x):
return np.isclose(x[1], 0)
def top_boundary(x):
return np.isclose(x[1], 1)
def G(x):
return np.sin(5*x[0])
def g(x, n):
return (np.sin(5*x[0])*n[0], np.sin(5*x[0])*n[1])
W0, _ = W.sub(0).collapse()
w_bc = dolfinx.fem.Function(W0)
top_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim-1, top_boundary)
top_cells = dolfinx.mesh.compute_incident_entities(mesh, top_facets, mesh.topology.dim-1, mesh.topology.dim)
w_bc.interpolate(lambda x: g(x, (0,1)), top_cells)
bottom_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim-1, bottom_boundary)
bottom_cells = dolfinx.mesh.compute_incident_entities(mesh, bottom_facets, mesh.topology.dim-1, mesh.topology.dim)
w_bc.interpolate(lambda x: g(x, (0, -1)), bottom_cells)
w_bc.x.scatter_forward()
bc_dofs = dolfinx.fem.locate_dofs_topological((W.sub(0), W0), mesh.topology.dim-1, np.hstack([top_facets, bottom_facets]))
bc = dolfinx.fem.dirichletbc(w_bc, bc_dofs, W.sub(0))
# Compute solution
w = dolfinx.fem.Function(W)
solver = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc], u=w)
solver.solve()
(sigma, u) = w.split()
BDM_out = dolfinx.fem.VectorFunctionSpace(mesh, ("DG", 1))
sigma_out = dolfinx.fem.Function(BDM_out)
sigma_out.interpolate(sigma)
with dolfinx.io.VTXWriter(mesh.comm, "sigma.bp", [sigma_out]) as vtx:
vtx.write(0)
with dolfinx.io.XDMFFile(mesh.comm, "u.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(u)
```

yielding the solutions