How to compute the strain & stress matrices?

A projection is simply solving
\int_\Omega u \cdot v ~\mathrm{d}x = \int_\Omega f \cdot v ~\mathrm{d} x
which is easy to implement in DOLFINx:

A more scalable approach than the route taken in DOLFINy is to create a projector class, that caches the mass matrix, so that you do not have to re-assemble (and refactorize in the case of LU) if you use it multiple times).

For the quadrature points (in physical space), consider the first bit of: Numerical values from ufl.SpatialCoordinate - #2 by dokken