function [l2error,h1error]=computeError(T,uh,u,ux,uy)

% [l2error,h1error]=computeError(T,uh,u,ux,uy)
%
% Input:
%      T    : basic triangulation
%      uh   : vector of nodal values of u_h
%      u    : function
%      ux   : partial_x u
%      uy   : partial_y u
% Output:
%      l2error : \| u-u_h\|_{L^2}
%      h1error : \|\grad u -\grad u_h\|_{L^2}

% Seven point quadrature formula with order 5

nodes = [0.333333333333333 0.333333333333333 0.333333333333333; ...
	   0.059715871789770 0.470142064105115 0.470142064105115; ...
	   0.470142064105115 0.059715871789770 0.470142064105115; ...
	   0.470142064105115 0.470142064105115 0.059715871789770; ...
	   0.797426985353087 0.101286507323456 0.101286507323456; ...
	   0.101286507323456 0.797426985353087 0.101286507323456; ...
	   0.101286507323456 0.101286507323456 0.797426985353087]';
weights = 0.5*[0.225000000000000 ...	   
	           0.132394152788506 ...
	           0.132394152788506 ...
	           0.132394152788506 ...
	           0.125939180544827 ...
	           0.125939180544827 ...
	           0.125939180544827]';

x=T.coordinates(:,1); x=x(T.elements)*nodes;      % x coordinates of quadrature points
y=T.coordinates(:,2); y=y(T.elements)*nodes;      % y coordinates of quadrature points
T=expandGrid(T);

% L2 error

uint=uh(T.elements)*nodes;    % Value of u_h at all the quadrature points
u   =u(x,y);                  % Value of u at all the quadrature points
l2error=sqrt(dot( T.detB, (uint-u).^2*weights ));

% H1 error

uhx = uh(T.elements(:,1)).*T.g11 + uh(T.elements(:,2)).*T.g12...
         +uh(T.elements(:,3)).*T.g13;    % One value per triangle
uhy = uh(T.elements(:,1)).*T.g21 + uh(T.elements(:,2)).*T.g22...
         +uh(T.elements(:,3)).*T.g23;
ux  = ux(x,y);
uy  = uy(x,y);
error=bsxfun(@minus,uhx,ux).^2+bsxfun(@minus,uhy,uy).^2;   % Errors in quad points
h1error=sqrt(dot( T.detB, error*weights ));

return