Computing grad(u) over each triangle

General Discussion
1

Hi there!

I am trying to find the components of grad(w) in each triangle where w is a P1 element. Here is my code

for (int i=0;i<nbtriangles;i++){ 
	// For each triangle T_i in the mesh, we have w(x)=a_ix+b_iy+c_i
	// and thus grad w = (a_i,b_i). 	
	real xp1= Th[i][0].x;
	real xp2= Th[i][1].x;
	real xp3= Th[i][2].x;
	real yp1= Th[i][0].y;
	real yp2= Th[i][1].y;
	real yp3= Th[i][2].y;
	real wp1=w(xp1, yp1); // Interpolation  
	real wp2=w(xp2, yp2); // is happening  
	real wp3=w(xp3, yp3); // here. 
	M = [[xp1,yp1,1],[xp2,yp2,1],[xp3,yp3,1]];					
	real[int] bb = [wp1,wp2,wp3];
	real[int] xx = [0,0,0];	
	set(M,solver=sparsesolver);
	xx=M^-1 * bb; 
}

Is there a faster way to do this? Interpolation is making the code slow!

Also, if I remove the line set(M,solver=sparsesolver) then I get the error"

sizestack + 1024 =16896 ( 15872 )

– Square mesh : nb vertices =1089 , nb triangles = 2048 , nb boundary edges 128
– Solve :
min -0.137402 max 0.137394
– Solve :
min -0.0755225 max 0.0755066

2

Hello,
Did you try using “dx(w)” and “dy(w)” ?
I did not check if it needs to be “transferred” from function type to vector type…
HTH
Denis

3

I tried “dx(w)” and it outputs a number not a vector.

On another note, with regards to the code in the question, is there a way to get w(x,y) if (x,y) is a mesh point, without doing interpolation?

4

dx(w) gives the partial derivative w.r.t. x and dy(w) gives the partial derivative w.r.t. y.

Yes, if u is a FE function then u is a real[int] which gives the values at the degrees of freedom. So something like u[Th[i][0]] should give you the value of u at the point Th[i][0]. But you should check the FreeFem manual to be sure how this works.

5

But w is a P1 element, not the constant function. Therefore dx(w) should be a function of x and y and not just a number. More specifically, grad(w) should be a different constant on each triangle in the mesh Th.
So what is dx(w) actually outputting?

6

dx(w) is not a function, but a Finite element function.

In order to get a different constant on every triangle in the mesh use a P0 finite element to capture it. In fact, if you try to store dx(w) as a P1 function you get inaccurate results (as shown in the code below). If you store it as a P0 function then I guess you’ll recover what you want.

int NN = 20;
mesh Th=square(NN,NN,[-1+2x,-1+2y]);

fespace Vh(Th,P1);
fespace Vh0(Th,P0);

Vh u = sqrt(x^2+y^2);

Vh gx = dx(u);
plot(gx,fill=1,wait=1,cmm=“grad P1”);

Vh0 gx0 = dx(u);
plot(gx0,fill=1,wait=1,cmm=“grad P0”);