Hello. I’m starting with FreeFem, so sorry if my question is too simple, or if this is not the right place to ask something like this. But I haven’t found the answer in the tutorials.

I’m implementing a code to solve the Stokes Problem with velocity and pressure in P1 (space of continuous functions in the domain. Each element of the space is a linear polynomial in each element of the mesh). The mesh is the typical 2D mesh composed of triangles.

I short, the numerical scheme to solve the Stokes problem is:

```
solve Stokes( [u1, u2, p], [v1, v2, q] ) //, solver=LU)
= int2d(Th)(
something + hK*something
)
- int2d(Th)(
somethign + hK*something
)
+ on(1,2,3,4, u1=g1, u2=g2);
```

where

```
hK=4*AreaT/(3*r)
```

with AreaT is the area of the current triangle and

```
r = sqrt((x1-xc)^2+(x2-xc)^2+(x3-xc)^2+(y1-yc)^2+(y2-yc)^2+(y3-yc)^2)
```

where (x1,y1), (x2,y2), and (x3,y3) are the vertices of the current element (triangle) of the mesh, and

```
(xc,yc) = (x1+x2+x3,y1+y2+y3)/3.
```

Since each hK is different on each element, I don’t know how I could implement this stabilization parameter in the integral. In particular, how can I calculate “r”? Maybe it is not possible, and that is why I did not find a way to implement it.

Thanks in advance for your help and advice.