Accessing discretized gradient matrix (B-matrix)

dolfinx
1

In dolfinx, is it somehow possible to access the discretized gradient matrices (“B-matrices”)? With B-matrix I refer to the matrix which relates the nodal displacements to local strains strains=B*displacements. Either on the element level or in assembled form.

Or is this concept not used in dolfinx?

Thank you in advance.

2

See, e.g., https://github.com/FEniCS/dolfinx/blob/main/python/demo/demo_elasticity.py.

3

Thanks for the answer. However, I was unfortunately not able see how this example provides the answer I was looking for. I see that a ‘near-nullspace’ consisting of translations and rotations is computed, but not explicitly the ‘B-matrix’ which relates the local strains to the displacement field in linear elasticity. But maybe I am missing something.

4

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

5

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