% This program constructs a stiffness matrix for a five-point-stencil
% on a simple rectangular lattice of size
%    nside1 x nside2
% 
% The figure illustrates the ordering of points in the grid
% for the case nside1=5 and nside2=8.
%
%        5-----10----15----20----25----30----35----40
%        |     |     |     |     |     |     |     |
%   n    4-----9-----14----19----24----29----34----39
%   s    |     |     |     |     |     |     |     |
%   i    3-----8-----13----18----23----28----33----38
%   d    |     |     |     |     |     |     |     |
%   e    2-----7-----12----17----22----27----32----37
%   1    |     |     |     |     |     |     |     |
%        1-----6-----11----16----21----26----31----36
%
%                        nside2
%
% To demonstrate how things work, the program solves a Dirichlet and a
% Neumann problem on the grid.
%
% Note that the function that constructs the stiffness matrix has several
% different calling options:
%
%   flag_matrix='orig'   The usual 5-point discrete Laplacian 4/-1/-1/-1/-1
%
%   flag_matrix='rcon'   A grid Laplacian where the conductivities of all
%                        bars are drawn at random.
%
%   flag_matrix='rcut'   A grid Laplacian for a grid where some bars have
%                        been cut at random.
%
%   flag_matrix='peri'   Conductivities vary periodically.
%
%   flag_matrix='crac'   Two long cracks in the grid.

function main_stiffness

%LOCAL_plot_lattice
%keyboard
%LOCAL_solve_dirichlet
%keyboard
%LOCAL_solve_neumann
%keyboard
%LOCAL_compute_evals
%keyboard
LOCAL_block_tessellate
keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% An auxiliary function that plots all bars in the lattice.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_plot_lattice()

rng('default')
rng(0)

nside1      = 25;               %
nside2      = 12;
N           = nside1*nside2;
flag_matrix = 'rcon';

fprintf(1,'=== Plot lattce - domain size: %d x %d = %d\n',nside1,nside2,N)

%%% Construct the geometry.
[XX1,XX2] = meshgrid((1:nside1)-0.5,(1:nside2)-0.5);
xx        = [reshape(XX1,1,N);reshape(XX2,1,N)];

%%% Construct the stiffness matrix (this is a sparse matrix).
flag_anchor = 0;
A           = LOCAL_get_A(nside1,nside2,flag_matrix,flag_anchor);

%%% Identify boundary nodes and interior nodes.
I1 = find( (xx(1,:) < 1) | ...
           (xx(1,:) > (nside1-1)) | ...
           (xx(2,:) < 1) | ...
           (xx(2,:) > (nside2-1)) );
I2 = LOCAL_index_complement(I1,N);

%%% Plot the lattice
hold off
LOCAL_draw_lattice(xx,A)
hold on
han = plot(xx(1,I1),xx(2,I1),'b.',...
           xx(1,I2),xx(2,I2),'r.');
set(han,'MarkerSize',25)
legend(han,'Boundary nodes','Interior nodes')
axis equal
axis([0,nside1,0,nside2])

return

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

function LOCAL_solve_dirichlet

rng('default')
rng(0)

nside1      = 25;               %
nside2      = 12;
N           = nside1*nside2;
flag_matrix = 'rcon';

fprintf(1,'=== Domain size: %d x %d = %d\n',nside1,nside2,N)

%%% Construct the geometry.
[XX1,XX2] = meshgrid((1:nside1)-0.5,(1:nside2)-0.5);
xx        = [reshape(XX1,1,N);reshape(XX2,1,N)];

%%% Construct the stiffness matrix (this is a sparse matrix).
flag_anchor = 0;
A           = LOCAL_get_A(nside1,nside2,flag_matrix,flag_anchor);

%%% Identify boundary nodes and interior nodes.
I1 = find( (xx(1,:) < 1) | ...
           (xx(1,:) > (nside1-1)) | ...
           (xx(2,:) < 1) | ...
           (xx(2,:) > (nside2-1)) );
I2 = LOCAL_index_complement(I1,N);

%%% Solve a Dirichlet problem
fprintf(1,'Solving a Dirichlet problem.\n')
func  = inline('x1.*x1 + x2.*x1','x1','x2');
u1 = func(xx(1,I1),xx(2,I1))';
u2 = A(I2,I2)\(-A(I2,I1)*u1);
uu = zeros(N,1);
uu(I1) = u1;
uu(I2) = u2;
fprintf(1,'Residual after Dirichlet solve = %12.5e\n',...
        max(abs(A(I2,:)*uu)))
UU = reshape(uu,nside2,nside1);
hold off
mesh(XX1,XX2,UU)
title('Solution to Dirichlet problem')

return

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

function LOCAL_solve_neumann

rng('default')
rng(0)

nside1      = 25;               %
nside2      = 12;
N           = nside1*nside2;
flag_matrix = 'rcon';

fprintf(1,'=== Domain size: %d x %d = %d\n',nside1,nside2,N)

%%% Construct the geometry.
[XX1,XX2] = meshgrid((1:nside1)-0.5,(1:nside2)-0.5);
xx        = [reshape(XX1,1,N);reshape(XX2,1,N)];

%%% Construct the stiffness matrix (this is a sparse matrix).
flag_anchor = 0;
A           = LOCAL_get_A(nside1,nside2,flag_matrix,flag_anchor);

%%% Identify boundary nodes and interior nodes.
I1 = find( (xx(1,:) < 1) | ...
           (xx(1,:) > (nside1-1)) | ...
           (xx(2,:) < 1) | ...
           (xx(2,:) > (nside2-1)) );
I2 = LOCAL_index_complement(I1,N);

%%% Solve a Neumann problem.
%%% The load consists of a unit source and a unit sink placed randomly on
%%% the boundary.
hold off
fprintf(1,'Solving a Neumann problem.\n')
A(1,1) = A(1,1) + 1; % "Anchor" a corner node to avoid null-space issue.
ff = zeros(N,1);                  
i_in  = ceil(length(I1)*rand(1));
i_out = ceil(length(I1)*rand(1));
ff(I1(i_in))  =  1;
ff(I1(i_out)) = -1;
uu = A\ff;
fprintf(1,'Residual after Neumann solve = %12.5e\n',...
        max(abs(A*uu - ff)))
UU = reshape(uu,nside2,nside1);
mesh(XX1,XX2,UU)
title('Solution to Neumann problem')
return

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

function LOCAL_compute_evals()
rng('default')
rng(0)

nside1      = 25;               %
nside2      = 12;
N           = nside1*nside2;
flag_matrix = 'orig';

fprintf(1,'=== Domain size: %d x %d = %d\n',nside1,nside2,N)

%%% Construct the geometry.
[XX1,XX2] = meshgrid((1:nside1)-0.5,(1:nside2)-0.5);
xx        = [reshape(XX1,1,N);reshape(XX2,1,N)];

%%% Construct the stiffness matrix (this is a sparse matrix).
flag_anchor = 0;
A           = LOCAL_get_A(nside1,nside2,flag_matrix,flag_anchor);

%%% Identify boundary nodes and interior nodes.
I1 = find( (xx(1,:) < 1) | ...
           (xx(1,:) > (nside1-1)) | ...
           (xx(2,:) < 1) | ...
           (xx(2,:) > (nside2-1)) );
I2 = LOCAL_index_complement(I1,N);

%%% Compute the eigenvalues of A.
hold off
[V,D] = eig(full(A(I2,I2)));
ee = diag(D);
figure(1)
plot(1:length(I2),real(ee),'x-')
title('Eigenvalues of A (real part)')
for imode = 1:10
  u  = zeros(N,1);
  u(I2) = real(V(:,imode));
  UU = reshape(u,nside2,nside1);
  figure(2)
  mesh(XX1,XX2,UU)
  title(sprintf('Eigenmode %d',imode))
  keyboard
end

return

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

function LOCAL_block_tessellate()

rng('default')
rng(10)

nleaf  = 10; % Each patch has nleaf x nleaf nodes.
m1     = 6; % The grid consists of m1 x m2 patches.
m2     = 4; % The grid consists of m1 x m2 patches.
nside1 = nleaf*m1;
nside2 = nleaf*m2;
N      = nside1*nside2;
fprintf(1,'\n=== Parameters: m1=%2d   m2=%2d   nleaf=%3d   ntot=%d\n',...
        m1,m2,nleaf,N)

% Construct the stiffness matrix (this is a sparse matrix).
flag_matrix = 'rcon';
flag_anchor = 1;
A = LOCAL_get_A(nside1,nside2,flag_matrix,flag_anchor);

%%% Construct the geometry.
[XX1,XX2] = meshgrid((1:nside1)-0.5,(1:nside2)-0.5);
xx        = [reshape(XX1,1,N);reshape(XX2,1,N)];

%%% Construct the index vectors.
I_STORE = cell(1,m1*m2);  % These are the lists of BOUNDARY NODES.
J_STORE = cell(1,m1*m2);  % These are the lists of INTERIOR NODES.
for i1 = 1:m1
  for i2 = 1:m2
    icell = m2*(i1-1) + i2;
    x1min = (i1-1)*nleaf;
    x1max = i1*nleaf;
    x2min = (i2-1)*nleaf;
    x2max = i2*nleaf;
    % Find the index vector "indbox" marking the nodes on the rim.
    indbox = find( (x1min < xx(1,:)) .* ...
                  (xx(1,:) < x1max) .* ...
                  (x2min < xx(2,:)) .* ...
                  (xx(2,:) < x2max) );
    % Find the indices within indbox that are interior.
    Jrel = find( ((x1min+1) < xx(1,indbox)) .* ...
                 (xx(1,indbox) < (x1max-1)) .* ...
                 ((x2min+1) < xx(2,indbox)) .* ...
                 (xx(2,indbox) < (x2max-1)) );
    Jrel_comp = LOCAL_index_complement(Jrel,length(indbox));
    J_STORE{icell} = indbox(Jrel);      % These are the interior nodes.
    I_STORE{icell} = indbox(Jrel_comp); % These are the exterior nodes.
  end
end

%%% Assemble the compressed stiffness matrix.
k = length(I_STORE{1});     % This is the block size of the compressed matrix.
ACOMP = sparse(k*m1*m2,k*m1*m2);     % This is the compressed matrix.
indcomp = zeros(1,k*m1*m2); % This is an index vector representing the "surviving nodes".
for ibox = 1:(m1*m2);
  indi = k*(ibox-1) + (1:k);  % This is an index vector into the REDUCED matrix.
  I    = I_STORE{ibox};
  J    = J_STORE{ibox};
  % Inser the Schur complement on the diagonal block.
  ACOMP(indi,indi) = full(A(I,I) - A(I,J)*(A(J,J)\A(J,I)));  
  % Insert the "I" vector
  indcomp(indi)    = I_STORE{ibox};                          
  for jbox = 1:(m1*m2)
    if ~(ibox == jbox)
      indj = k*(jbox-1) + (1:k);
      ACOMP(indi,indj) = A(I,I_STORE{jbox});
    end
  end
end

figure(1)
plot(xx(1,indcomp),xx(2,indcomp),'.')
axis equal
title('The surviving nodes')

%%% Set up and solve a global problem.
f      = zeros(N,1);
indrim = find( (xx(1,:) < 1) | ...
               ((nside1-1) < xx(1,:)) | ...
               (xx(2,:) < 1) | ...
               ((nside2-1) < xx(2,:)) );
i_in  = ceil(length(indrim)*rand(1)); % We place a point source here.
i_out = ceil(length(indrim)*rand(1)); % We place a point sink here.
f(indrim(i_in))  =  1;
f(indrim(i_out)) = -1;
tic
u = A\f;
fprintf(1,'Time for solving original   problem = %0.3f\n',toc)

%%% Solve the compressed problem.
fcomp = f(indcomp);
tic
ucomp = ACOMP\fcomp;
fprintf(1,'Time for solving compressed problem = %0.3f\n',toc)
fprintf(1,'max|u(indcomp) - ucomp| = %12.5e\n',max(abs(u(indcomp) - ucomp)))

%%% Reconstruct the full solution from ucomp.
u_recon = zeros(N,1);
for ibox = 1:(m1*m2)
  ind   = k*(ibox-1) + (1:k);
  I     = I_STORE{ibox};
  J     = J_STORE{ibox};
  u_rim = ucomp(ind);
  u_int = A(J,J)\(-A(J,I)*u_rim);
  u_recon(I) = u_rim;
  u_recon(J) = u_int;
end
fprintf(1,'max|u - u_recon|        = %12.5e\n',max(abs(u - u_recon)))

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is the key function in this program. 
% For documention on the call see the top of the file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_get_A(nside1,nside2,flag_matrix,flag_anchor)

if (nargin < 4)
  flag_anchor = 0;
end

if (nargin == 1)
  nside2 = nside1;
end

BARS  = LOCAL_get_links(nside1,nside2);
nbars = size(BARS,2);

if strcmp(flag_matrix,'orig')
    
  BARS(3,:) = ones(1,nbars);
  
elseif strcmp(flag_matrix,'rcon')
    
  BARS(3,:) = 1 + rand(1,nbars);
  
elseif strcmp(flag_matrix,'rcut')
    
  p = 0.04;
  BARS(3,:) = (rand(1,nbars) > p);
  
elseif strcmp(flag_matrix,'peri')
  
  [XX1,XX2] = meshgrid(0:(nside1-1),0:(nside2-1));
  xx        = [reshape(XX1,1,nside1*nside2);reshape(XX2,1,nside1*nside2)];
  h         = 1/(max([nside1,nside2])-1);
  Ncell     = round(min(nside1,nside2)/25);
  XX1sc     = h*xx(1,BARS(1,:));
  XX2sc     = h*xx(2,BARS(1,:));
  YY1sc     = h*xx(1,BARS(2,:));
  YY2sc     = h*xx(2,BARS(2,:));
  ZZ1sc     = 0.5*(XX1sc+YY1sc);
  ZZ2sc     = 0.5*(XX2sc+YY2sc);
  BARS(3,:) = 1 - 0.9*((cos(Ncell*pi*ZZ1sc).*cos(Ncell*pi*ZZ2sc)).^2);
    
elseif strcmp(flag_matrix,'crac')

  % We create a list of bars to be removed:
  n1a = round(0.4*nside1);
  n1b = round(0.6*nside1);
  n2a = round(0.4*nside2);
  n2b = round(0.7*nside2);
  gap = 4;
  
  ii1 = n1a*nside2 + (1:n2a);
  jj1 = ii1 + nside2;
  ii2 = n1a*nside2 + n2a + 1 + nside2*(1:(n1b-n1a));
  jj2 = ii2 - 1;
  ii3 = n1b*nside2 + 1 + (n2a:n2b);
  jj3 = ii3 + nside2;
  ii4 = (n1b-gap)*nside2 + (nside2:(-1):(n2a+gap+2));
  jj4 = ii4-nside2;
  ii5 = (n1a-1)*nside2 + (n2a+gap+1) + nside2*(1:(n1b-n1a-gap));
  jj5 = ii5+1;
  
  ii = [ii1,ii2,ii3,ii4,ii5];
  jj = [jj1,jj2,jj3,jj4,jj5];
  
  ni   = length(ii);
  BARS(3,:) = ones(1,nbars);
  BARS = [BARS,[ii;jj;-ones(1,ni)]];
      
end

if (flag_anchor == 0)
    
  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,:)]);
         
else
    
  ii_rim = [1:nside1,...
            (nside1+1):nside1:(nside1*(nside2-2)+1),...
            (2*nside1):nside1:(nside1*(nside2-1)),...
             nside1*(nside2-1) + (1:nside1)];
     
  A = sparse([BARS(1,:), BARS(1,:), BARS(2,:), BARS(2,:),ii_rim],...
             [BARS(1,:), BARS(2,:), BARS(1,:), BARS(2,:),ii_rim],...
             [BARS(3,:),-BARS(3,:),-BARS(3,:), BARS(3,:),ones(1,2*(nside1+nside2-2))]);
    
end

if strcmp(flag_matrix,'crac')
  [ii,jj,aa] = find(A);
  ind = find(abs(aa) < 1e-12);
  aa(ind) = zeros(size(ind));
  A = sparse(ii,jj,aa);
end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function draws the grid corresponding to the sparse matrix A.
% Node "i" has coordinates xx(:,i) in the plane.
% Two nodes i and j are connected iff A(i,j) is non-zero.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_draw_lattice(xx,A)

[ii,jj,aa] = find(A);

hold off
plot([xx(1,ii);xx(1,jj)],[xx(2,ii);xx(2,jj)],'k')
hold on
plot(xx(1,:),xx(2,:),'b.')
hold off

bord = 2;
axis([min(xx(1,:))-bord,...
      max(xx(1,:))+bord,...
      min(xx(2,:))-bord,...
      max(xx(2,:))+bord])
axis equal

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function constructs the complement of the index vector 
% ind1 in the set [1,2,3,...,ntot].
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ind2 = LOCAL_index_complement(ind1,ntot)

ind       = 1:ntot;
ind(ind1) = 2*ntot*ones(1,length(ind1));
ind       = sort(ind);
ind2      = ind(1:(ntot-length(ind1)));

return

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

