% This function provides a starting point for solving homework problem 5.2.
%
% DRIVER_dragrid: Simply a function to draw the grid.
%
% DRIVER_computeschur: Builds the Schur complement for the basic 5-point
% stencil associated with the Laplace equation.
%
% DRIVER_randomcond: Builds the Schur complement associated with a grid
% conduction problem. The conductivities of the bars are drawn
% independently, so there is no "smoothness" at all. 

function HW05p2

%DRIVER_drawgrid
DRIVER_computeschur
%DRIVER_randomcond

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_drawgrid

n  = 5;
n1 = 2*n+3;
n2 = n+2;
h  = 1/(n+1);
[XX1,XX2] = meshgrid(linspace(-1,1,n1),linspace(0,1,n2));
xx = [reshape(XX1,1,numel(XX1));...
      reshape(XX2,1,numel(XX2))];

I1 = find((abs(xx(1,:)) < (0.5*h)) & ...
          (abs(xx(2,:) - 0.5) < (0.5 - 0.5*h)));
I2 = find((abs(xx(1,:) + 0.5) < (0.5 - 0.5*h)) & ...
          (abs(xx(2,:) - 0.5) < (0.5 - 0.5*h)));
I3 = find((abs(xx(1,:) - 0.5) < (0.5 - 0.5*h)) & ...
          (abs(xx(2,:) - 0.5) < (0.5 - 0.5*h)));
han = plot([-1,1,1,-1,-1],[0,0,1,1,0],'k',...
            xx(1,:),xx(2,:),'g.',...
            xx(1,I1),xx(2,I1),'k.',...
            xx(1,I2),xx(2,I2),'b.',...
            xx(1,I3),xx(2,I3),'r.',...
            'MarkerSize',40,'LineWidth',2);
set(han(2),'MarkerEdgeColor',0.7*[1,1,1])
axis equal
axis off
print -dpng fig_twobox.png
keyboard


return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_computeschur

%%% Set the problem size and build the basic grid.
n  = 100;
if (mod(n,2) == 1)
  fprintf(1,'ERROR: n should be even.\n')
  keyboard
end
n1 = 2*n+3;
n2 = n+2;
h  = 1/(n+1);
[XX1,XX2] = meshgrid(linspace(-1,1,n1),linspace(0,1,n2));
xx = [reshape(XX1,1,numel(XX1));...
      reshape(XX2,1,numel(XX2))];

%%% Partition the grid.
I1 = find((abs(xx(1,:)) < (0.5*h)) & ...
          (abs(xx(2,:) - 0.5) < (0.5 - 0.5*h)));
I2 = find((abs(xx(1,:) + 0.5) < (0.5 - 0.5*h)) & ...
          (abs(xx(2,:) - 0.5) < (0.5 - 0.5*h)));
I3 = find((abs(xx(1,:) - 0.5) < (0.5 - 0.5*h)) & ...
          (abs(xx(2,:) - 0.5) < (0.5 - 0.5*h)));

%%% Build the finite difference matrix.
D1 = spdiags((1/(h*h))*[-ones(n1,1),2*ones(n1,1),-ones(n1,1)],[-1,0,1],n1,n1);
D2 = spdiags((1/(h*h))*[-ones(n2,1),2*ones(n2,1),-ones(n2,1)],[-1,0,1],n2,n2);
A  = kron(speye(n1),D2) + kron(D1,speye(n2));

%%% Form the Schur complement S and compute the svds of the quadrant.
S      = full(A(I1,I1) - A(I1,I2)*(A(I2,I2)\A(I2,I1)) ...
                       - A(I1,I3)*(A(I3,I3)\A(I3,I1)));
nhalf  = round(n/2);
Sblock = S(1:nhalf,nhalf + (1:nhalf));
ss     = svd(Sblock);

%%% Plot the results.
semilogy(1:nhalf,ss/ss(1),'r.-')
axis tight
              
keyboard
return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_randomcond

%%% Set the problem size and build the basic grid.
n  = 100;
if (mod(n,2) == 1)
  fprintf(1,'ERROR: n should be even.\n')
  keyboard
end
n1 = 2*n+3;
n2 = n+2;
h  = 1/(n+1);
[XX1,XX2] = meshgrid(linspace(-1,1,n1),linspace(0,1,n2));
xx = [reshape(XX1,1,numel(XX1));...
      reshape(XX2,1,numel(XX2))];

%%% Partition the grid.
I1 = find((abs(xx(1,:)) < (0.5*h)) & ...
          (abs(xx(2,:) - 0.5) < (0.5 - 0.5*h)));
I2 = find((abs(xx(1,:) + 0.5) < (0.5 - 0.5*h)) & ...
          (abs(xx(2,:) - 0.5) < (0.5 - 0.5*h)));
I3 = find((abs(xx(1,:) - 0.5) < (0.5 - 0.5*h)) & ...
          (abs(xx(2,:) - 0.5) < (0.5 - 0.5*h)));

%%% Build the finite difference matrix.
% The object "BARS" holds a list of all links in the lattice.
%   BARS(1,i) = index of starting point of link i.
%   BARS(2,i) = index of end point of link i.
%   BARS(3,i) = conductivityof link i.
BARS      = LOCAL_get_links(n1,n2);
nbars     = size(BARS,2);
BARS(3,:) = 1 + rand(1,nbars); % Let each conductivity be a random number.
A         = sparse([BARS(1,:), BARS(1,:), BARS(2,:), BARS(2,:)],...
                   [BARS(1,:), BARS(2,:), BARS(1,:), BARS(2,:)],...
                   [BARS(3,:),-BARS(3,:),-BARS(3,:), BARS(3,:)]);

%%% Form the Schur complement S and compute the svds of the quadrant.
S      = full(A(I1,I1) - A(I1,I2)*(A(I2,I2)\A(I2,I1)) ...
                       - A(I1,I3)*(A(I3,I3)\A(I3,I1)));
nhalf  = round(n/2);
Sblock = S(1:nhalf,nhalf + (1:nhalf));
ss     = svd(Sblock);

%%% Plot the results.
semilogy(1:nhalf,ss/ss(1),'k.-')
axis tight

keyboard
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This supplementary functions returns a list of all links in a
% rectangular grid of size n1 x n2.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function BARS = LOCAL_get_links(n1,n2)

if (nargin == 1)
  n2 = n1;
end

BARS  = zeros(3,2*n1*n2 - n1 - n2);

ndone = 0;
for i1 = 1:(n1-1)
  for i2 = 1:(n2-1)
    j_c = n2*(i1-1) + i2;
    j_n = n2*(i1-1) + i2 + 1;
    j_e = n2*i1 + i2;
    BARS(1:2,ndone + (1:2)) = [j_c, j_c; j_n, j_e];
    ndone = ndone + 2;
  end
  j_c = n2*i1;
  j_e = n2*(i1+1);
  BARS(1:2,ndone + 1) = [j_c; j_e];
  ndone = ndone + 1;
end
for i2 = 1:(n2-1)
  j_c = n2*(n1-1) + i2;
  j_n = n2*(n1-1) + i2 + 1;
  BARS(1:2,ndone + 1) = [j_c; j_n];
  ndone = ndone + 1;
end

BARS(3,:) = ones(1,ndone);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

