function [ nodeInfo, neighborInfo, cInteriorNodes ] = createFDGrid( gamma, box, n )
%CREATEFDGRID Summary of this function goes here
% Input:
%   gamma   Function handle describing the boundary of domain.
%           It is gamma<0 outside the domain, gamma>0 inside and gamma=0 on
%           the boundary.
%   box     A vector that describes a box containing the domain.
%           The box in 2D is given by (box[1,1], box[1,2]) x (box[2,1],
%           box[2,2]).
%   n       Number of edges in one direction. Can be a vector or scalar.
%
% Output:
%   nodeInfo        A matrix which contains in its columns all the nodes of
%                   the grid. The first cInteriorNodes nodes are interior
%                   nodes, all others are boundary nodes.
%   neighborInfo    A matrix containing in each column a list of all
%                   neighbor indices for each interior node.
%   cInteriorNodes  Used for interpreting nodeInfo.
%
% Remarks:
%   - There will be an error if the box does not contain all of the domain.
%     This case would have to be treated separately.
%   - Uniqueness of boundary nodes is not necessarily given.
%   - Boundary nodes can be outside the box.
%   - Behavior for multiple zeros of gamma on one edge is not defined.
%   - No serious error processing.
%   - There may be unstable situations when interior points are very close
%     to the boundary.
%

    % Do basic setup.
    [~,dim] = size(box);
    if( length(n) == 1 )
        n   = repmat(n,dim,1);
    end
    h               = (box(:,2)-box(:,1)) ./ n;
    n               = n+1; % n now describes the nodes in one direction.
    cNodes          = prod(n);
    nodeInfo        = zeros( dim+1, 2*cNodes );
    neighborInfo    = zeros( 2^dim, cNodes );
    cInteriorNodes  = 0;
    cBoundaryNodes  = 0;
    
    % Iterate over the uniform grid.
    for i = 1:cNodes
        
        % Convert i to multi-index.
        I   = toMultiIndex( i-1, n );
        
        % Is the current point inside the grid?
        x   = box(:,1)  + I.*h;
        if( gamma(x) > 0 )
            cInteriorNodes  = cInteriorNodes+1;
            nodeInfo(:,i)   = [x; cInteriorNodes];
            
            % Insert neighbor information.
            J   = I;
            y   = x;
            for k = 1:dim
                for s = 1:2
                    J(k)    = I(k)+(-1)^s;
                    y(k)    = x(k)+h(k)*(-1)^s;
                    if( gamma(y) > 0 )
                        neighborInfo(2*(k-1)+s,i) = fromMultiIndex( J, n ) + 1;
                    else
                        alpha   = fzero( @(t) gamma( x + t*(y-x) ), [0,1] );
                        z       = x;
                        z(k)    = x(k) + alpha*h(k)*(-1)^s;
                        cBoundaryNodes = cBoundaryNodes+1;
                        nodeInfo(:,cNodes+cBoundaryNodes) = [z; -1];
                        neighborInfo(2*(k-1)+s,i) = cNodes+cBoundaryNodes;
                    end
                end
                y(k)    = x(k);
                J(k)    = I(k);
            end
        end
    end
    
    % Update indices of boundary nodes.
    for i = 1:cBoundaryNodes
        nodeInfo(dim+1,cNodes+i)    = cInteriorNodes+i;
    end

    % Generate set of reduced indices.
    globalIdxActive                 = (nodeInfo(dim+1,:)~=0);
    reducedIdx                      = zeros(2*cNodes,1);
    reducedIdx(globalIdxActive)     = 1:(cInteriorNodes+cBoundaryNodes);
    
    % Update indices from global to reduced ones and drop unused parts.
    nodeInfo            = nodeInfo(1:dim, globalIdxActive);
	neighborInfo        = neighborInfo(:,globalIdxActive(1:cNodes));
    neighborInfo(:)     = reducedIdx( neighborInfo(:) );
end

function i = fromMultiIndex( I, n )
%FROMMULTIINDEX Converts a multi-index to a number
% To-do: Add explanation.

    dim         = length(n);
    i           = 0;
    multiplier  = 1;
    for k = 1:dim
        i           = i + I(k)*multiplier;
        multiplier  = multiplier * n(k);
    end

end

function I = toMultiIndex( i, n )
%TOMULTIINDEX Converts a number i to a multi-index given a vector n.
% To-do: Add explanation.

    dim = length(n);
    I   = zeros( dim, 1 );
    for k = 1:dim
        I(k)    = mod(i,n(k));
        i       = floor(i/n(k));
    end
    
end

