# Computing grad(u) over each triangle

**URL:** https://community.freefem.org/t/computing-grad-u-over-each-triangle/287
**Category:** General Discussion
**Created:** [March 3, 2020, 2:17pm UTC](https://community.freefem.org/t/computing-grad-u-over-each-triangle/287 "2020-03-03T14:17:27Z")
**Posts on this page:** 6
**Page:** 1

<div class="post-metadata">

### Author: ![Neutrino](https://avatars.discourse-cdn.com/v4/letter/n/bb73d2/32.png) [@Neutrino](https://community.freefem.org/u/Neutrino)
#### Post date: [March 3, 2020, 2:17pm UTC](https://community.freefem.org/t/computing-grad-u-over-each-triangle/287/1 "2020-03-03T14:17:27Z")

</div>

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

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

---

<div class="post-metadata">

### Author: ![denis.crete](https://avatars.discourse-cdn.com/v4/letter/d/9d8465/32.png) [@denis.crete](https://community.freefem.org/u/denis.crete)
#### Post date: [March 4, 2020, 8:29pm UTC](https://community.freefem.org/t/computing-grad-u-over-each-triangle/287/2 "2020-03-04T20:29:20Z")

</div>

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

---

<div class="post-metadata">

### Author: ![Neutrino](https://avatars.discourse-cdn.com/v4/letter/n/bb73d2/32.png) [@Neutrino](https://community.freefem.org/u/Neutrino)
#### Post date: [March 5, 2020, 12:38pm UTC](https://community.freefem.org/t/computing-grad-u-over-each-triangle/287/3 "2020-03-05T12:38:24Z")

</div>

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?

---

<div class="post-metadata">

### Author: ![BeniaminBogosel](https://avatars.discourse-cdn.com/v4/letter/b/3ec8ea/32.png) [@BeniaminBogosel](https://community.freefem.org/u/BeniaminBogosel)
#### Post date: [March 5, 2020, 1:48pm UTC](https://community.freefem.org/t/computing-grad-u-over-each-triangle/287/4 "2020-03-05T13:48:28Z")

</div>

> [@Neutrino](#):
>
> I tried “dx(w)” and it outputs a number not a vector.

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

> [@Neutrino](#):
>
> is there a way to get w(x,y) if (x,y) is a mesh point, without doing interpolation?

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.

---

<div class="post-metadata">

### Author: ![Neutrino](https://avatars.discourse-cdn.com/v4/letter/n/bb73d2/32.png) [@Neutrino](https://community.freefem.org/u/Neutrino)
#### Post date: [March 5, 2020, 8:01pm UTC](https://community.freefem.org/t/computing-grad-u-over-each-triangle/287/5 "2020-03-05T20:01:30Z")

</div>

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?

---

<div class="post-metadata">

### Author: ![BeniaminBogosel](https://avatars.discourse-cdn.com/v4/letter/b/3ec8ea/32.png) [@BeniaminBogosel](https://community.freefem.org/u/BeniaminBogosel)
#### Post date: [March 5, 2020, 9:47pm UTC](https://community.freefem.org/t/computing-grad-u-over-each-triangle/287/6 "2020-03-05T21:47:24Z")

</div>

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+2_x,-1+2_y]);

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”);