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