% An example using adaptmesh to resolve a solution that has
% a singularity. 
% This is a (slightly) annotated version of the example at
% http://www.mathworks.com/access/helpdesk/help/toolbox/pde/adaptmesh.html
% DPO 03/2004

% Solve the Laplace equation over a circle sector, with 
% Dirichlet boundary conditions u = cos(2/3atan2(y,x)) along the arc, 
% and u = 0 along the straight lines, and compare to the exact 
% solution. We refine the triangles using the worst error criterion 
% until we obtain a mesh with at least 500 triangles:

[ua,pa,ea,ta]=adaptmesh('cirsg','cirsb',1,0,0,'maxt',500,... 
                       'tripick','pdeadworst','ngen',inf); 
xa=pa(1,:); ya=pa(2,:); 
exacta=((xa.^2+ya.^2).^(1/3).*cos(2/3*atan2(ya,xa)))'; 
disp('The true solution is u(x,y) = ((x^2+y^2)^(1/3)*cos(2/3*atan2(y,x)))')
disp('Near the origin, the solution is approx u = r^{1/3}.')
disp(sprintf('Adaptive mesh with %d triangles: max error at the mesh points is %e',size(ta,2),max(abs(ua-exacta)) ))
% ans = 0.0058 
% ans = 534 
figure(1)
pdemesh(pa,ea,ta)
title('Mesh used by adaptmesh')
figure(2)
pdeplot(pa,ea,ta,'xydata',ua)
title('Solution computed by adaptmesh')

[p,e,t]=initmesh('cirsg'); 
[p,e,t]=refinemesh('cirsg',p,e,t); 
u=assempde('cirsb',p,e,t,1,0,0); 
x=p(1,:); y=p(2,:); 
exact=((x.^2+y.^2).^(1/3).*cos(2/3*atan2(y,x)))'; 
disp(sprintf('Uniform mesh with %d triangles: max error at the mesh points is %e',size(t,2),max(abs(u-exact)) ))
%  ans = 0.0085 
% ans = 1640 

% We test how many refinements we have to use with a uniform mesh
% to get error comparable to that obtained by adaptmesh.

[p,e,t]=refinemesh('cirsg',p,e,t); 
u=assempde('cirsb',p,e,t,1,0,0); 
x=p(1,:); y=p(2,:); 
exact=((x.^2+y.^2).^(1/3).*cos(2/3*atan2(y,x)))'; 
disp(sprintf('Uniform mesh with %d triangles, max error at the mesh points is %e',size(t,2),max(abs(u-exact)) ))
% ans = 0.0054 
% ans = 6560 
figure(3)
pdemesh(p,e,t)
title('Uniformly refined mesh')
%figure(4)
%pdeplot(p,e,t,'xydata',u)
%title('Solution for uniformly refined mesh')