function delsq_test ( )

%*****************************************************************************80
%
%% delsq_test() tests delsq().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    30 August 2023
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../delsq' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'delsq_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Test delsq()\n' );
%
%  Start with small sample.
%
  R = 'L';
  n = 10;
  G_10 = numgrid ( R, n );
  fprintf ( 1, '\n' )
  fprintf ( 1, '  numgrid region "%c" with n = %d\n', R, n );
  spy ( G_10 );
  title ( 'L region: 10x10 finite difference grid' );
  filename = 'grid_10x10.png';
  print ( '-dpng', filename );
  fprintf ( '  Graphics saved as "%s"\n', filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Node indices:\n' );
  disp ( G_10 )
%
%  Now go to an interesting size.
%
  R = 'L';
  n = 32;
  G = numgrid ( R, n );
  fprintf ( 1, '\n' )
  fprintf ( 1, '  numgrid region "%c" with n = %d\n', R, n );
  spy ( G );
  title ( 'L region: 32x32 finite difference grid' );
  filename = 'grid_32x32.png';
  print ( '-dpng', filename );
  fprintf ( '  Graphics saved as "%s"\n', filename );
%
%  Generate the Laplacian operator.
%
  D = delsq ( G );
  fprintf ( 1, '\n' )
  fprintf ( 1, '  Sparse Laplacian operator D\n')
  spy ( D );
  title ( 'L region: 32x32 Laplacian operator' );
  filename = 'delsq_10x10.png';
  print ( '-dpng', filename );
  fprintf ( '  Graphics saved as "%s"\n', filename );
%
%  Count the grid points.
%
  grid_num = sum ( G(:) > 0 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This grid uses %d interior nodes.\n', grid_num );
%
%  Define a boundary value problem for grid function u in the region:
%    delsq(u) = 1 in interior
%          u  = 0 on boundary.
%
  rhs = ones ( grid_num, 1 );
  u = D \ rhs;
%
%  Map the solution back onto the grid.
%  
  U = G;
  U ( G > 0 ) = full ( u ( G ( G > 0 ) ) );
%  
%  Create a contour plot, labeling the contour lines.
%  Use the prism colormap.
%
  clabel ( contour ( U ) );
  prism ( );
  axis ( 'square', 'ij' );
  filename = 'contour_plot.png';
  print ( '-dpng', filename );
  fprintf ( '  Graphics saved as "%s"\n', filename );
%
%  Now show solution as a mesh plot.
%
  mesh ( U );
  axis ( [ 0, n, 0, n, 0, max(max(U)) ] );
  axis ( 'square', 'ij' );
  filename = 'mesh_plot.png';
  print ( '-dpng', filename );
  fprintf ( '  Graphics saved as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'delsq_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../delsq' );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

