will again result in some values computed in mesh vertices, while, according to my understanding the strains should be given per finite element instead.

Question

So how does this work in FreeFem++? The displacements are known in the vertices. Does the same somehow apply to the strains and stresses too? Are is there any additional postprocessing required to get the correct values, e.g. averaging nodal stress/strain values per finite element?

With P1 elements I am not sure if there is a need to use Gauss quadrature points. The integrals over an element can be computed exactly just by taking a weighted average of the vertex values. In each element the displacement derivatives will be constant.

The way I used to compute strains and stresses was to use a lower order fespace, for example something like:

macro u [ux,uy] // displacement
macro v [vx,vy] // test function
macro stress [Sxx,Syy,Sxy] // stress tensor
macro e(u) [dx(u[0]), dy(u[1]), (dx(u[1]) + dy(u[0]))/2] // strain tensor
macro C [[2*mu+lambda,lambda,0],[lambda,2*mu+lambda,0],[0,0,2*mu]] // plane strain stress-strain matrix
fespace Sh(Th,P1); // scalars
fespace Vh(Th,[P2,P2]); // vectors
fespace Zh(Th,[P1,P1,P1]); // tensors
Vh u,v;
// compute u
Zh stress;
stress = C*e(u);

But I am not sure what exactly happens upon assigment when the left-hand side and right-hand side use differing finite element spaces. This seemed the easiest way to do it.

In both cases FreeFem will use a numerical integration formula with a single node (1/3,1/3) (presumably in barycentric triangle coordinates) and weight |T_k| equal to the triangle area that can be computed from the vertex coordinates.

Edit: By default FreeFem uses the formula qf5pt which is exact for polynomials of degree 5 on triangles or edges (in 3d). This formula uses 7 quadrature nodes.