% Solve - Lap u = 1 on a unit circle with zero bc
% approx solution uh is a column vector
%  xi, yi,  ex, ey, area are row vectors
% This code requires the files circleg.m circleb1.m
% containing the geometry and boundary conditions.
% These can be formed using pdetool
% graderror computes the gradient error between u_h and u_I
format long
[p,e,t] = initmesh('circleg', 'Hmax', 1);
maxerror = [] ; maxerr =1;
fcnerror = [];
graderror = [];
maxexerror = [];
maxeyerror = [];
while maxerr > 0.0005
     [p,e,t] = refinemesh('circleg', p,e,t);
     uh = assempde('circleb1',p,e,t,1,0,1);
% p(1,:) denotes the first row of the matrix p
% that contains the x coordinates of the vertices of the mesh
% p(2,:) denotes the second row of the matrix p
% that contains the y coordinates of the vertices of the mesh
     exact = -(p(1,:).^2 + p(2,:).^2-1)/4;
% The expression p(1,:).^2 gives a row vector each of whose
% components is the square of the corresponding component of
% the row vector p(1,:)
     maxerr = norm(uh-exact',inf);
     maxerror = [maxerror maxerr];
%  Compute the gradient error directly using pdegrad
     [ex,ey] = pdegrad(p,t,exact'-uh);
% Compute the area of the triangles directly using pdetrg (area)
     [area,a1,a2,a3]=pdetrg(p,t);
     graderr = sqrt(sum(area.*(ex.^2 + ey.^2)));
     graderror = [graderror, graderr];
% Compute the max gradient error for this mesh
     maxexerr = norm(ex',inf);
     maxexerror = [maxexerror maxexerr];
     maxeyerr = norm(ey',inf);
     maxeyerror = [maxeyerror maxeyerr];
% Compute the L2 function error
% Compute the error at the midpoints of the edges of the triangles
% Note: result is a row vector
     em1 = (exact(t(2,:))+ exact(t(3,:)) - uh(t(2,:))'- uh(t(3,:))')/2;
     em2 = (exact(t(3,:))+ exact(t(1,:)) - uh(t(3,:))'- uh(t(1,:))')/2;
     em3 = (exact(t(1,:))+ exact(t(2,:)) - uh(t(1,:))'- uh(t(2,:))')/2;
     fcnerr = sqrt((em1*(area.*em1)' + em2*(area.*em2)' + em3*(area.*em3)')/3);
     fcnerror = [fcnerror, fcnerr];
end