function [uest] = esterr(elements3,coordinates,u)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% esterr for P1 on triangular elements
%% works with ACF code (Alberty/Carstensen/Funken)
%% MPeszynska, for MTH 655-9, Winter 2017
%% Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Remarks: estimates the H1 errror given nodal values <u>
%% This function should be called from fem2d.m or fem2d*.m 
%% Functions tri_quadcofs and basisfun must be available
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     nels = size(elements3,1);
     uel = zeros(nels,2);

     for j = 1:nels

  %%%%%%%%%%%%%%%%%%%%%% compute the actual integral of uex,unum

       midp = (sum(coordinates(elements3(j,:),:))/3)';

       [dum,dnum] = basisfun (midp(1),midp(2), ...
           coordinates(elements3(j,:),:),u([elements3(j,:)]));

       uel(j,:)= dnum;
       ff(j,:) = f([midp(1),midp(2)]);

     end
     %%%%%%%%%%%%%%%%%%%%%%%%
     for j = 1:nels

        [n,l] = edge(coordinates(elements3(j,:),:));
        v = (n*uel(j,:)');
        uerr(j) = l.^2*v.^2/2;

        [wq,xq,yq]= tri_quadcofs(coordinates(elements3(j,:),:),1);

        uerr(j) = uerr(j) + sum(wq)*ff(j,:)^2;

     end

uest = uerr;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [normal, length] = edge(vertices)

e(1,:) = vertices(1,:)-vertices(2,:);
e(2,:) = vertices(2,:)-vertices(3,:);
e(3,:) = vertices(3,:)-vertices(1,:);

for i=1:3 length(i)=norm(e(i,:)); end;

for i=1:3 un(i,:) = [e(i,2),-e(i,1)]; end;

normal=un;
end
