function [ graduTri, gradu ]  = gradientReconstruction( mesh, u )

[ graduTri, ~ ]  = gradientPerElement( mesh, u );
gradu = zeros(2,mesh.NuDoF);

for i=1:mesh.NuDoF
  index    = mesh.Point2Element(i).Elements;
  gradTri  = graduTri(:,index);
  gradu(1,i) = sum(gradTri(1,:))/size(gradTri,2);
  gradu(2,i) = sum(gradTri(2,:))/size(gradTri,2);
end