# Computing strains and stresses: per mesh vertice or element?

**URL:** <https://community.freefem.org/t/computing-strains-and-stresses-per-mesh-vertice-or-element/1286>\
**Category:** General Discussion\
**Created:** [October 27, 2021, 7:15am UTC](https://community.freefem.org/t/computing-strains-and-stresses-per-mesh-vertice-or-element/1286 "2021-10-27T07:15:18Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![distractor](https://avatars.discourse-cdn.com/v4/letter/d/57b2e6/32.png) [@distractor](https://community.freefem.org/u/distractor)\
**Post date:** [October 27, 2021, 7:15am UTC](https://community.freefem.org/t/computing-strains-and-stresses-per-mesh-vertice-or-element/1286/1 "2021-10-27T07:15:18Z")

</div>

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

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

```

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

```auto
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?

---

<div class="post-metadata">

**Author:** ![ivan](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ivan/32/823_2.png) [@ivan](https://community.freefem.org/u/ivan)\
**Post date:** [November 4, 2021, 6:09pm UTC](https://community.freefem.org/t/computing-strains-and-stresses-per-mesh-vertice-or-element/1286/2 "2021-11-04T18:09:34Z")

</div>

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:

```cpp
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.

---

<div class="post-metadata">

**Author:** ![ivan](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ivan/32/823_2.png) [@ivan](https://community.freefem.org/u/ivan)\
**Post date:** [November 4, 2021, 6:37pm UTC](https://community.freefem.org/t/computing-strains-and-stresses-per-mesh-vertice-or-element/1286/3 "2021-11-04T18:37:54Z")

</div>

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

```auto
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](https://en.wikipedia.org/wiki/Triangle#Using_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](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/f3534b02907aedc69ea2dd9507337887996a6ea2.png)