Implementing a non-constant stabilized parameter

General Discussion
1

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.

2

Isn’t that a good question? Or maybe it’s not the place to ask it?

Maybe I wasn’t clear enough. Because I doubt FreeFem doesn’t have the option.

3

No problem , AreaT is defined in FreeFEM it is area.

Now computation of r

fespce Ph(Th,P0); // function constant par element
Ph xc=x,yc=x; // Barycenter of all element 
varf vr(unused,v) = int2d(Th,qft=qf1pTlump) ( 3*((x-xc)^2+(y-yc)^2)/area); 
Ph r;
r []= vr(0,Ph); 
r = sqrt(r);

I am sorry the technics is tricking : qf1pTlump quadrature sum on all vertices of triangle /3
and xc and yc is the barycenter of the triangle because the dof of Ph are all barycenter of triangle.

Remark: not tester so a bug is possible.

4

Many thanks prof. Hecht!!