Isn’t that just the domain integration (surface integration)?
\bar\sigma = \frac{\int \sigma dS}{\int dS}
i.e. (not exact syntax, but similar to)
area = scalar_assemble(Constant(1)*dx) # note this is not the correct syntax for Constant
stress = assemble(sigma*dx)
average_stress = stress / area