I have worked out a solution, just in case somebody else is interested in the future:
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import mesh, fem
from dolfinx.fem import petsc#this seems to be required since dolfinx 0.8
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([1,1,1])],
[2,1,1], cell_type=mesh.CellType.hexahedron)
element = ufl.VectorElement("CG", domain.ufl_cell(), 1, 3)
V = fem.FunctionSpace(domain, element)
element_strain = ufl.VectorElement("DG", domain.ufl_cell(), 0, 6)
V_strain = fem.FunctionSpace(domain, element_strain)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V_strain)
#volumes of integration points
V_vol = fem.FunctionSpace(domain, ("DG", 0))
vol = fem.Function(V_vol)
fem.petsc.assemble_vector(vol.vector, fem.form(
1*ufl.TestFunction(V_vol)*ufl.dx(metadata={"quadrature_degree": 1})))
vol.x.scatter_forward()
def epsilon(u):
return ufl.as_vector([u[0].dx(0),
u[1].dx(1),
u[2].dx(2),
(u[0].dx(1) + u[1].dx(0)),
(u[0].dx(2) + u[2].dx(0)),
(u[1].dx(2) + u[2].dx(1))])
b = fem.form(ufl.inner(epsilon(u),v)/vol*ufl.dx)
B = fem.petsc.assemble_matrix(b)
B.assemble()