function val = interp2d(x,y,t,xref,yref)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Computes interpolation error 
%% Call:  interp2d(x,y,t,xref,yref) 
%%%%%% where x,y,t are nodes and element info
%%%% example: (No plotting): comment out the plotting
%% x = [0 1 0 1];
%% y = [0 0 1 1];
%% t = [ 4 2 3; 1 2 3];
%%%%% another example with a Delanuay tesselation
%% nx =10; ny = 10;
%% xx = linspace(0,1,nx+1);
%% xx = exp(20*xx)/exp(20);
%% yy=linspace(0,1,ny+1);
%% [x,y] = meshgrid(xx,yy);
%% t = delaunay(x,y);
%% xref = linspace(0,1,10*nx+1);yref = xref;
%% [xr,yr]=meshgrid(xref,yref);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

trimesh(t,x,y);title('Triangulation');pause
surfc(x,y,myfun(x,y));title('Plot of the function that is interpolated');pause

%% compute the interpolation error in L2 and H1 norm 
sval = 0; svald = 0;
for numel=1:size(t,1)
    [v,vd] = myinterp2d (numel,t,x,y);
    %% fprintf('numel=%d err =%g %g\n',numel,v,vd);
    sval = sval +v.^2;
    svald = svald + vd.^2;
end
sval = sqrt(sval);
svald = sqrt(svald);

fprintf('L2err = %g H01 err =%g\n',sval,svald);    

if nargin > 3
    surfc(xref,yref,myfun(xref,yref)-interp2(x,y,myfun(x,y),xref,yref));
    title('Plot of interpolation error');
end
function [val,vald]=myinterp2d(numel,t,x,y)    
    coords = mycoord(numel,t,x,y);
    xc = x(t(numel,:)); yc = y(t(numel,:));
    %%%%%%%%%%%%% get the coefficients of the quadrature
    [wq,xq,yq] = tri_quadcofs(coords,13);
    %%%%%%%%%%%%%% find coefficients of the basis functions
    vander = [ones(3,1),coords];
    psi = vander \ eye(3,3); %% solve vandermonde system 
    psi1 = psi(:,1);
    psi2 = psi(:,2);
    psi3 = psi(:,3);

    uval = myfun(xc,yc); %% evaluate u at the nodes of triangle
    basisfun1 = @(xx,yy)(psi1(1) + psi1(2)*xx + psi1(3)*yy);
    basisfun2 = @(xx,yy)(psi2(1) + psi2(2)*xx + psi2(3)*yy);
    basisfun3 = @(xx,yy)(psi3(1) + psi3(2)*xx + psi3(3)*yy);
    interpolant = @(xx,yy) uval(1)*basisfun1(xx,yy) ...
        + uval(2)*basisfun2(xx,yy)+uval(3)*basisfun3(xx,yy);
    dxinterpolant = @(xx,yy) uval(1)*psi1(2)+uval(2)*psi2(2)+uval(3)*psi3(2);
    dyinterpolant = @(xx,yy) uval(1)*psi1(3)+uval(2)*psi2(3)+uval(3)*psi3(3);

    s = 0; sd = 0;
    for r=1:length(wq)
        [val,dx,dy] = myfun(xq(r),yq(r));
        myval = interpolant(xq(r),yq(r)) - val;
        mydxval = dxinterpolant(xq(r),yq(r)) - dx;
        mydyval = dyinterpolant(xq(r),yq(r)) - dy;
        s = s + wq(r)*myval*myval;
        sd = sd + wq(r)*(mydxval*mydxval + mydyval*mydyval);
    end
    val = sqrt(s);
    vald = sqrt(sd);

function carr = mycoord(num,t,x,y)
n1 = t(num,1); n2 = t(num,2) ; n3 = t(num,3);
carr = [x(n1) y(n1); x(n2) y(n2); x(n3) y(n3)];


function [z,dx,dy]=myfun(x,y)
    z = x.^2 + y.^2;
    dx = 2*x;
    dy = 2*y;
    