My approach to get the B-matrix which relates the strain field to the displacement field by $$\varepsilon_i = B_{ij} u_j $$ is to compute $$B_{ij} =\partial \varepsilon_i/\partial u_j$$.

Here is my code:

```
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import mesh, fem
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([1.0, 1.0, 1.0])],
[2,2,2], cell_type=mesh.CellType.hexahedron)
element = ufl.VectorElement("CG", domain.ufl_cell(), 1, 3)
V = fem.FunctionSpace(domain, element)
element_gp = ufl.VectorElement("CG", domain.ufl_cell(), 1, 3)
# element_gp = ufl.VectorElement("DG", domain.ufl_cell(), 0, 3)
V_gp = fem.FunctionSpace(domain, element_gp)
def epsilon(u):
return u[0].dx(0)
u = fem.Function(V)
# u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
uv = ufl.variable(u)
b = ufl.diff(epsilon(uv),uv)
# b_apply = ufl.algorithms.apply_derivatives.apply_derivatives(b)
b_expr = fem.Expression(b,V_gp.element.interpolation_points())
B_fun = fem.Function(V_gp)
B_fun.interpolate(b_expr)
print(B_fun.x.array)
```

For simplicity, I have considered only a single, scalar strain component.

However, the output for B_fun is just an array of zeros.

I note that most “classical” FEM programs define these B-matrices explicitly, while dolfinx operates using automatic differentiation, as I understand it. I am still wondering if it is nevertheless possible to recover the B-matrix.

Any help would be greatly appreciated. Thanks!

Best

Erik