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