function [K1,rhs1,h1,K2,rhs2,h2] = slit2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [K1,rhs1,h1,K2,rhs2,h2] = slit2
% This program generates two linear systems:
%    K1 u1 = rhs1        K1: 1208 x 1208  h1 = mesh parameter
%    K2 u2 = rhs2        K2: 4931 x 4931  h2 = mesh parameter
% related to solving Laplace's equation on a circle sector.
% They are generated by adaptive mesh generation on a problem
% whose solution has a singularity.
% This is a modified and annotated version of the example at
% http://www.mathworks.com/access/helpdesk/help/toolbox/pde/adaptmesh.html
% A different version, slit.m, was used to demonstrate adaptive 
% mesh refinement.  This version is used in Homework 3.
%
% You may comment out the plotting if that is more convenient.
%
% Dianne O'Leary 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(1)
pdeplot(pa,ea,ta,'xydata',ua)
% [K,F1,B1,UD]=ASSEMPDE(B,P,E,T,C,A,F) assembles the PDE problem by
%         eliminating the Dirichlet boundary conditions from the system of
%         linear equations. UN=K\F1 returns the solution on the non-Dirichlet
%         points.  The solution to the full PDE problem can be obtained by
%         the MATLAB command U=B1*UN+UD.
[K,rhs,B,ud]=assempde('cirsb',pa,ea,ta,1,0,0);
un = K \ rhs;
uc = B*un+ud;
pdeplot(pa,ea,ta,'xydata',uc)

% Refine the grid once (uniformly) to make a bigger matrix.

figure(2)
[p1,e1,t1]=refinemesh('cirsg',pa,ea,ta); 
[K1,rhs1,B1,ud1]=assempde('cirsb',p1,e1,t1,1,0,0);
un1 = K1 \ rhs1;
uc1 = B1*un1+ud1;
pdeplot(p1,e1,t1,'xydata',uc1)
h1 = 1/sqrt(size(K1,1));

figure(1)
pdemesh(p1,e1,t1)

% Refine the grid once more.

figure(3)
[p2,e2,t2]=refinemesh('cirsg',p1,e1,t1); 
[K2,rhs2,B2,ud2]=assempde('cirsb',p2,e2,t2,1,0,0);
un2 = K2 \ rhs2;
uc2 = B2*un2+ud2;
pdeplot(p2,e2,t2,'xydata',uc2)
h2 = 1/sqrt(size(K2,1));