function [p] = poisson_3d(rho, delta)
% Solves poisson equation for a 3D cuboid domain
%   poisson_3d solves the poisson equation \Delta(p) = rho using a finite
%   differencing scheme as discussed in [1] p. 1024, with a direct matrix
%   method. Note that for large regions, this computation can become extremely
%   expensive.
%
%   Currently, the boundary condition applied is the simplest possible: a
%   dirichlet condition, with the value p = 0 everywhere on the boundary. It
%   would not be difficult to modify the Right Hand Side solution vector to
%   encompass a different boundary. Reference [2] is a useful link showing you
%   how to do that.
%
%   NB FOR THIS WEB VERSION I'VE REMOVED THE SUITESPARSE DEPENDANCY   
%        %   The SUITESPARSE toolbox, made available to the MATLAB community by Tim Davis
%        %   (posted on www.mathworks.com) is deployed, in order to speed up the 
%        %   initialisation of the sparse coefficients matrix [3].
%
%   There are also some notes in the code below, at the appropriate points.
%
% Syntax:  
%       [p] = poisson_3d(rho, delta)
%
% Inputs:
%       rho     [m x n x q] double  3D array containing the source distribution 
%                                   for which the equation is to be solved. This
%                                   is the RHS of the \Delta(p) = rho equation
%
%       delta   [1 x 1] double      The grid spacing, in whatever units are
%                                   required (typically m) of the rho array. The
%                                   grid spacing should be the same in all
%                                   directions.
%       
% Outputs:
%       p       [m x n x q] double  Solution to the equation. 
%
% References:
%
%   [1] Press, Teulosky, Vetterling, Flannery 'Numerical Recipes Third Edition'
%       Cambridge University Press 2007, ISBN 978-0-521-88068-8
%   
%   [2] http://en.wikipedia.org/wiki/Discrete_Poisson_equation
%
%   [3] http://www.mathworks.com/matlabcentral/fileexchange/12233
%
%
% Other m-files required:   SuiteSparse toolbox package [3]
% Subfunctions:             none
% Nested Functions:         none
% MAT-files required:       none
%
%
% Author:           T.H. Clark
% Work address:     Fluids Lab
%                   Cambridge University Engineering Department
%                   2 Trumpington Street
%                   Cambridge
%                   CB21PZ
% Email:            t.clark@cantab.net
% Website:          http://www.eng.cam.ac.uk/thc29
%
% Created:          10 December 2008 
% Last revised:     21 December 2008


%% CONTROL PARAMETERS

% Spymatrix: true = plot the sparsity pattern of the matrix, false = do not plot
opts.spymatrix = false;


%% CONSTRUCT THE RIGHT HAND SIDE

% We need to reshape the 'rho' input array to a column vector, to solve Ax=rho.
% Reshape works in a COLUMNAR direction - this is fine, since the coefficients
% matrix is symmetric
rhs = reshape(rho,numel(rho),1);
rhs = rhs*(delta.^2);


%% ADD THE BOUNDARY CONDITIONS

% Boundary conditions are applied in the form of corrections to the right hand
% side array. However, in this case, all corrections are zero since a dirichlet
% (value of p, rather than derivative) condition is applied; with a uniform
% boundary value equal to zero.
%
% Thus we do nothing here... see [2] for how to apply a non-zero boundary


%% CONSTRUCT THE COEFFICIENTS MATRIX

% The coefficients matrix is of dimension mn x mn:
[m,n,q] = size(rho);
mnq = m*n*q;

% Number of nonzero elements in the mxm D-Component matrix [3]
nnz_D = 2*(m-1) + m;

% Number of nonzero elements in the mxm Identity matrix [3] is simply m

% There are 2*q*(n-1) + 2*n*(q-1) identity matrices within the coefficients
% matrix, and there are nq D-component matrices within the coefficients matrix.

% Thus the total number of nonzeros in the coefficients matrix is:
nnz_total = nnz_D*n*q + 2*m*q*(n-1) + 2*m*n*(q-1);

% We can calculate the sparsity of the matrix.
sparsity = 1-(nnz_total/(mnq*mnq));
disp(['poisson_3d: Matrix A is ' num2str(sparsity*100) '% sparse'])
disp(['poisson_3d: Matrix A has ' num2str(nnz_total) ' nonzeros'])

% Construct the coefficients matrix.
%   Note that we use sign convention from ref [1] rather than ref [3], where
%   identity fringes are positive, and the main diagonal contains negative
%   coefficients.

% The pattern of the matrix is highly sparse, symmetric and tridiagonal,
% with fringes.
%
% While it would appear ideal, the spdiag function has irritating
% dimension-dependant functionality, not to mention a time-consuming and
% memory-heavy implementation. Instead, we'll directly index in.

% Indices of the subside outer identity fringe diagonal
i_inds = (m*n+1:mnq+1:mnq*m*n*(q-1))';
% i_inds = ( m+1 : mn+1 : m*m*n*(n-1)  )';

% Add the indices of the superside outer identity fringe diagonal
i_inds = [i_inds; i_inds + (m*n*mnq-m*n)];

% Add the indices of the subside inner identity fringe diagonal
tempvec = (m+1:mnq+1:(mnq^2 - mnq))';
deletervec = [];
for ctr = 1:1:m
    deletervec = [deletervec, m*(n-1)+ctr : (m*n) : numel(tempvec)]; %#ok<AGROW> NB FIX THIS LATER!!
end
tempvec(deletervec) = [];
i_inds = [i_inds; tempvec];

% Add the indices of the superside inner identity fringe diagonal
tempvec = tempvec -m + (n*mnq);
i_inds = [i_inds;tempvec];

% Add the indices of the subdiagonal
tempvec = ( 2 : mnq+1 : mnq*mnq )';
tempvec(m:m:end) = [];
i_inds = [i_inds; tempvec];

% Add the indices of the superdiagonal
i_inds = [i_inds; tempvec + mnq-1];
diag_start = numel(i_inds)+1;

% And the indices of the main diagonal
i_inds = [i_inds; ( 1 : mnq+1 : mnq*mnq )'];

% Define the coefficients:
coefficients = ones(size(i_inds));
coefficients(diag_start:end) = -6;

% We need the i,j indices, rather than a single-index form.
[I,J] = ind2sub([mnq mnq],i_inds);

% Clear up all unnecessary variables...
clearvars -except I J coefficients mnq m n q nnz_total rhs opts

% For the sake of web publishing, I've removed the dependancy on suitesparse,
% and am using ML's regular sparse() command
        % Assign the sparse matrix, with the required amount of memory. This uses the
        % suitesparse version of a sparse matrix initialiser, sparse2(), which behaves
        % in exactly the same manner as MATLAB's sparse() function, but is faster. This
        % is due to a different preliminary sorting algorithm implementation.
A = sparse(I,J,coefficients,mnq,mnq,nnz_total);

% Plot the matrix, if required
if opts.spymatrix
    figure()
    spy(A)
    title('Sparsity Pattern of 3D Poisson Matrix')
end


% Clear up again (absolutely maximising available memory!!)
clearvars -except A rhs m n q
disp('poisson_3d: Memory allocation immediately before solution:')
whos


%% SOLVE FOR p

% Matrix A is positive-definite and symmetric, with high sparsity, so we solve
% by simultaneous method. I've tried various methods in the suitesparse toolkit,
% but the matlab backslash operator prevails as the fastest for these matrices.
p = A\rhs;

% Reshape p into array form, to match the dimensions of the input rho
p = reshape(p,m,n,q);