%
%  my_delsqdemo.m
%
%  solving the poisson equation   del^2 u = 1  with u=0 on bndry
%    - a hacked version of "delsqdemo.m"
%

n =  32;
%n =  64;
%n = 128;

R = 'S'

clf reset 

    % Generate and display the grid.
    G = numgrid(R,n);
    spy(G)
    title('A finite difference grid')
    g = numgrid(R,12)
    disp('pause'), disp(' '), pause

    % Generate and display the discrete Laplacian.
    D = delsq(G);
    spy(D)
    title('The 5-point Laplacian')
    disp('pause'), disp(' '), pause
   
    % Number of interior points
    N = sum(G(:)>0);

    % Solve the Dirichlet boundary value problem:
    %    delsq(u) = 1 in the interior,
    %    u = 0 on the boundary.
    disp('Solving the sparse linear system ...')
    rhs = ones(N,1);
    flops(0)
    tic
    if R == 'N'
        % For nested dissection, turn off minimum degree ordering.
        spparms('autommd',0)
        u = D\rhs;
        spparms('autommd',1)
    else
        u = D\rhs;
    end
    flops
    toc

    % Map the solution onto the grid and display it.
    U = G;
    U(G>0) = full(u(G(G>0)));
    clabel(contour(U));
    prism
    axis('square'), axis('ij')
    disp('pause'), disp(' '), pause

    colormap((cool+1)/2);
    mesh(U)
    axis([0 n 0 n 0 max(max(U))])
    axis('square'), axis('ij')
    disp('pause'), disp(' '), pause