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

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
)

``````

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

``````

Thanks. This was a big improvement.