Computing strains and stresses: per mesh vertice or element?

General Discussion
1

Hi,

So I am not sure how FreeFem++ works, but I know that Abaqus for example has:

  • Displacements computed in the vertices
  • While strains and stresses are computed per finite element in Gaussian points.

However, I am not sure how FreeFem++ does that. Assuming we use

fespace Vh(Th, P1);
Vh u1, u2, v1, v2;

The displacements will be known in the mesh vertices. But computing strains as

macro e(u1, u2)
	[
		dx(u1),
		0.5 * (dy(u1) + dx(u2)),
		0.5 * (dx(u2) + dy(u1)), 
		dy(u2)
	] // EOM

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?

2

Hi Mitja,

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.

1 Like
3

Looking into the old documentation, when you use the P1 elements for the integrals it suffices to use:

int2d(Th, qft=qf1pT) ( ... ) \\ or
int2d(Th, qforder=2) ( ... )

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.

qf5pT

2 Likes