Stabilized FEM Method

General Discussion
1

Hello everyone,
I’m new in FreeFem and i’m trying to implement a stabilized FEM method, in particular i have the following linear form:
image
My problem is in that summation over the triangles of the triangularization of the domain, how can i integrate over every one of them in a individual way in order to get the vector corresponding to that linear form?
The parameter \tau_K its a stabilization parameter that depends on the size of each triangle.

Thanks in advance!

Juan.

2

Hello Juan,
I understand that you want to compute the piecewise gradient and Laplacian of v, obtained by taking the derivatives inside each triangle, and ignoring the jumps between triangles.
This is what does Freefem when you write dx() for a first-order derivative, or dxx() for a second-order derivative. For example dxx(v) is zero if v is P1, and dx(v) is zero if v is P0.
Therefore you can do for example

int nn=30;//number of intervals on boundary
border circle(t=0.,1.){x=cos(t*2.*pi);y=sin(t*2.*pi);label=1;};
mesh Th=buildmesh(circle(nn));//unit disc

fespace Vh(Th,P2);//space for the unknown or test function (choose P1,P2,...)
fespace Rh(Th,P0);//space for functions constant by triangle

Vh v;//this declaration is optional
Rh tau;

func f=x;//choose an appropriate value
func ax=x-y;//choose an appropriate value
func ay=-1.;//choose an appropriate value
tau=hTriangle;// example of size of triangle
real sigma=1.;
real cnu=1.;

varf linform(unused,v)=
       int2d(Th)(f*v)
       -int2d(Th)(f*tau*(sigma*v-ax*dx(v)-ay*dy(v)
                         -cnu*(dxx(v)+dyy(v))))
       ;
real[int] linf(Vh.ndof);
linf=linform(0,Vh);//the vector of the linear form, obtained by evaluating the variational formulation for v in the basis of Vh
1 Like
3

Hello François,
Thanks for the help, your explanation was clear. Just another question:
When you write “Rh tau”, is it possible to assign to the tau function a value that depends on each element? For example

\tau_K = \frac{h_K ^2}{\sigma h_K ^2 + \nu}

Where “h_K ^2” is the size of the triangle.

Thanks!

4

Rh tau just declares tau to be piecewise constant, since Rh is P0

For your required formula, instead of tau=hTriangle; just set
tau=hTriangle^2/(sigma*hTriangle^2+cnu);
Of course sigma and cnu need to be declared before.
hTriangle is a piecewise constant predefined function in FreeFem taking the value that you guess.

1 Like
5

Thank you, François!
You were very helpful.

Best regards, Juan.