%
%  solving the laplace equation   del^2 u = 0
%    - with u given on bndry as in the assignment!
%    - a hacked version of "delsqdemo.m"
%

n =  32;  h = 1/n;

A = 1;  B = exp(i*2*pi/3);   C = exp(i*4*pi/3); 

clf reset

% generate & plot the vectorizing indices
G = numgrid('S',n);

image(G,'cdatamapping','scaled');  
colormap(jet);  colorbar;  axis square, axis image
title('\bf colorized ordering of the vectorizing indices')
xlabel('\bf x-axis matrix indices')
ylabel('\bf y-axis matrix indices')

disp('pause');  pause

% generate & display the discrete, vectorized Laplacian

lapl = -delsq(G)/(h^2);
spy(lapl)
title('\bf banded matrix laplacian')
xlabel('\bf vector indices')
ylabel('\bf vector indices')

disp('pause');  pause
   
% Number of interior points (count when G_index > 0)
   
N = sum(G(:)>0);

% solve the Laplace boundary value problem:
%    delsq(u) = 0 in the interior
%    u = A,B,C on the boundary

rhs = zeros(N,1);
rhs(round(G(  2,2:n-1))) = -B/(h^2);
rhs(round(G(n-1,2:n-1))) = -C/(h^2);
rhs(round(G(2:n-1,  2))) = -C/(h^2);
rhs(round(G(2:n-1,n-1))) = -A/(h^2);

% map "rhs" onto a grid & plot abs(rhs)

R = G;  R(G>0) = (full(rhs(G(G>0))) ~=0);

image(R,'cdatamapping','scaled');  
colormap(jet);  colorbar;  axis square;  axis image
title('\bf boundary data location check')
xlabel('\bf x-axis matrix indices')
ylabel('\bf y-axis matrix indices')

%  do the solve

flops(0);  tic; 

u = lapl\rhs;

disp(' ')
disp(['  # of flops = ' num2str(flops) '   ;   cpu time = ' num2str(toc)])
disp(' ')

disp('pause');  pause

%  plot the solution (with BCs & corners = 0)

U = complex(G);  U(G>0) = full(u(G(G>0)));
U(1,2:end-1) = B;  U(2:end-1,1) = C;
U(n,2:end-1) = C;  U(2:end-1,n) = A;

%image(real(U),'cdatamapping','scaled');  colormap(jet);  
%image(imag(U),'cdatamapping','scaled');  colormap(jet);  
%image( abs(U),'cdatamapping','scaled');  colormap(jet);  
image(angle(U),'cdatamapping','scaled');  colormap(hsv);

colorbar;  axis square;  axis image
title('\bf complex laplace solution')
xlabel('\bf x-axis matrix indices')
ylabel('\bf y-axis matrix indices')