%
% this is a test of the quality of fish2d. We make a 2d problem, pass it
% to the solver, and measure the increase in work and the decrease in
% error as the mesh is refined. 
%
% This is similar to what we measured for the 1D mg codes.
%
lmax=8; lmin=3; ld=lmax-lmin+1; errdata=zeros(ld,5);
        for l=lmin:lmax
	nt=2^l+1; h=2^(-l); x=0:nt-1; x=x'*h;
	errdata(l-lmin+1,1)=h;
%
% Fix the right side so that the solution is 
%
% u(x,y) = sin(pi x) sin(pi y) e^x 
%
	wx=sin(pi*x); wxp=pi*cos(pi*x); wxpp=-pi*pi*wx;
	wy=wx; wypp=wxpp;
	u=(wx.*exp(x))*wy';
	wxx=(wx + 2*wxp + wxpp).*exp(x); uxx=wxx*wy';
	wyy=wypp; uyy=(wx.*exp(x))*wyy';
%
% rhs2d is a 2d grid function that includes the boundary conditions
%
	rhs2d=-(uxx+uyy); 
%
% rhs1d is a linear array that does not include the boundary conditions
%
	tmp=rhs2d(2:nt-1,2:nt-1);
	rhs1d=tmp(:);
%
% send rhs1d to the Poisson solver, get a 1d array back, map to a 2d
% array that includes the boundary data
%
	fin=flops;
	v1d=fish2d(rhs1d);
        fout=flops-fin;
	tmp(:)=v1d; v2d=zeros(nt,nt); v2d(2:nt-1,2:nt-1)=tmp;
        errdata(l-lmin+1,4)=fout;
%
%       Don't use norm(v2d-u,inf) on a 2d array. MATLAB will think it's
%       a matrix!
%
        errdata(l-lmin+1,2)=max(max(v2d-u));
end
errdata(2:ld,3)=errdata(2:ld,2)./errdata(1:ld-1,2);
errdata(2:ld,5)=errdata(2:ld,4)./errdata(1:ld-1,4)