Is there a more efficient way to many integrals over subdomains?

dolfinx
1

For postprocessing, I am evaluating the average strain/stress over many subdomains, e.g. 100. The code snippet below actually works, but it seems it could be more efficient. The code snippet evaluates 1200 integrals (2 variables X 100 subdomains X 6 components), each using its own form. I am thinking there is probably a way to use a one or two forms with the ufl.action() function.

This is not urgent, as the snippet works and is reallly not that slow. Just curious.

        ng = 100 # for example
        sig_avg = np.zeros((ng, 3, 3))
        eps_avg = np.zeros((ng, 3, 3))
        dx = ufl.Measure("dx", subdomain_data=ldr.cell_tags)
        for g in range(ng):
            for i in range(3):
                for j in range(i, 3):
                    eps_avg[g, i, j] = ldr.mesh.comm.allreduce(assemble(
                        form(strain[i, j] * dx(g))), op=MPI.SUM
                    )
                    sig_avg[g, i, j] = ldr.mesh.comm.allreduce(assemble(
                        form(stress[i, j] * dx(g))), op=MPI.SUM
                    )
2

This might be faster;

  1. Add an indicator function d from DG 0 to your form.
  2. Precompile the forms
  3. Change the values in d to be 0 on all cells except the one you want to integrate over.
    Following is a mwe
import ufl
import dolfinx
from mpi4py import MPI
import numpy as np
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
left_cells = dolfinx.mesh.locate_entities(
    mesh, mesh.topology.dim, lambda x: x[0] <= 0.5+1e-10)
markers = np.arange(len(left_cells), dtype=np.int32)
ct = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, left_cells, markers)

Q = dolfinx.fem.TensorFunctionSpace(mesh, ("DG", 1))
strain = dolfinx.fem.Function(Q)


def s(x):
    return (x[0], 3*x[1], 5*x[2], 0*x[0], x[2], x[1], x[2], x[1], x[0])


strain.interpolate(s)
D = dolfinx.fem.FunctionSpace(mesh, ("DG", 0))
d = dolfinx.fem.Function(D)


forms = []
dx = ufl.Measure("dx", domain=mesh)
for i in range(3):
    forms.append([])
    for j in range(3):
        forms[i].append(dolfinx.fem.form(d*strain[i, j]*dx))


ng = len(markers)
sig_avg = np.zeros((ng, 3, 3))
for g in range(ng):
    for i in range(3):
        for j in range(i, 3):
            d.x.set(0)
            d.x.array[ct.find(g)] = 1

            sig_avg[g, i, j] = dolfinx.fem.assemble_scalar(forms[i][j])
3

Thanks. This was a big improvement.