Accessing discretized gradient matrix (B-matrix)

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()
1 Like