% This program is modeled after Matlab's pdedemo1.m
% Dianne O'Leary  March 2005

%       Solve Poisson's equation
%        -div(grad(u))=1
%       on the unit disk with u=0 on the boundary.
%       Compare with exact solution.

%       Problem definition

g='circleg'; % The unit circle
b='circleb1'; % 0 on the boundary
c=1;
a=0;
f=1;

%       Take a coarse initial mesh
%       and perform regular refinement n-1 times

n = 6;
for i=1:n,
  if (i==1)
     [p,e,t]=initmesh(g,'hmax',1);
  else
     [p,e,t]=refinemesh(g,p,e,t);
  end
  u=assempde(b,p,e,t,c,a,f);
  exact=(1-p(1,:).^2-p(2,:).^2)'/4;
  er(i)=norm(u-exact,'inf');
  nt(i) = size(t,2);   % number of triangles
  fprintf('Error %d: %e. Number of triangles: %d\n',i,er(i),nt(i));...
end

% Plot the error and the mesh.

figure(1)
pdemesh(p,e,t)
disp('We solve the Poisson equation on the unit circle.')
disp('Figure 1 shows the final mesh.')
figure(2)
pdesurf(p,t,u-exact);
title('The error')

% Discuss the results.

disp('The true solution to this problem is (1-x^2-y^2)/4.')
disp('Figure 2 shows the error.  Notice the large error')
disp('at the origin, and also anomalies along the lines')
disp('where the sectors of the mesh meet.              ')
disp(' ')
disp('We estimate a mesh parameter h by assuming that')
disp('all triangles are equilateral right triangles with')
disp('equal area -- actually, h/4 works better, so that')
disp('is what we compute here (although I cannot explain why).')

h4 = sqrt(pi./nt)/sqrt(8);  % (h would be sqrt(2)*sqrt(pi./nt))

% Ratio of errors; asymptotically, should be 4

est = [0, er(1:n-1)./er(2:n)];
pred_h = [0,exp(-est(2:n)*log(2)./(4-est(2:n)))];

disp(' ')
disp('Now we display h/4 (computed from the number of trianges),')
disp('               the h predicted from the convergence ratio,')
disp('               and the convergence ratio of old error to new.')
disp(' ')
disp('grid  h/4    predicted_h  error ratio   ')
for i=1:n,
  disp(sprintf('%3d %7.2e   %7.2e   %7.2e   ', ...
           i,h4(i),pred_h(i),est(i)))
end
disp('  ')
disp('If your machine is fast enough to allow more refinements (n=9,')
disp('for example), you can see the predicted_h converging to  h/4.')