function biharmonic_fd2d_test ()

%*****************************************************************************80
%
%% biharmonic_fd2d_test() tests biharmonic_fd2d().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 June 2021
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../biharmonic_fd2d' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'biharmonic_fd2d_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Test biharmonic_fd2d()\n' );

  biharmonic_fd2d_solution_test ( );
% biharmonic_fd2d_solve_test ( );
  biharmonic_fd2d_stencil_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'biharmonic_fd2d_test:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../biharmonic_fd2d' );

  return
end
function biharmonic_fd2d_solution_test ( )

%*****************************************************************************80
%
%% biharmonic_fd2d_solution_test() tests biharmonic_fd2d_solution().
%
%  Discussion:
%
%    For each example, evaluate and plot the solution function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    31 May 2021
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'biharmonic_fd2d_solution_test():\n' );
  fprintf ( 1, '  biharmonic_fd2d_solution() evaluates an example\n' );
  fprintf ( 1, '  solution of the 2D biharmonic equation.\n' );
%
%  Get the default parameter values, so you can reset them.
%
  [ k, xmin, xmax, ymin, ymax ] = biharmonic_fd2d_parameters ( );
%
%  Set the test case.
%
  for k = 1 : 4

    [ k, xmin, xmax, ymin, ymax ] = biharmonic_fd2d_parameters ( ...
      k, xmin, xmax, ymin, ymax );

    n = 17;
    x = linspace ( xmin, xmax, n );
    y = linspace ( ymin, ymax, n );

    [ X, Y ] = ndgrid ( x, y );
    U = biharmonic_fd2d_solution ( X, Y );
    surf ( X, Y, U );
    title_string = sprintf ( 'biharmonic fd2d case %d\n', k );
    title ( title_string );
    filename = sprintf ( 'biharmonic_fd2d_case%d.png', k );
    print ( '-dpng', filename );
    fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  end

  return
end
function biharmonic_fd2d_solve_test ( )

%*****************************************************************************80
%
%% biharmonic_fd2d_solve_test() tests biharmonic_fd2d_solve().
%
%  Discussion:
%
%    For each example, evaluate and plot the solution function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 June 2021
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'biharmonic_fd2d_solve_test():\n' );
  fprintf ( 1, '  biharmonic_fd2d_solve() uses a 13 point finite difference\n' );
  fprintf ( 1, '  stencil to approximate the solution of a 2D biharmonic equation.\n' );
%
%  Get the default parameter values, so you can reset them.
%
  [ k, xmin, xmax, ymin, ymax ] = biharmonic_fd2d_parameters ( );
%
%  Set the test case.
%
  for k = 1 : 1

    [ k, xmin, xmax, ymin, ymax ] = biharmonic_fd2d_parameters ( ...
      k, xmin, xmax, ymin, ymax );

    xn = 17;
    x = linspace ( xmin, xmax, xn );
    yn = 17
    y = linspace ( ymin, ymax, yn );

    [ X, Y ] = ndgrid ( x, y );
    F = biharmonic_fd2d_rhs ( X, Y );

%   A = ?

    U1 = A \ F;

    U2 = biharmonic_fd2d_solution ( X, Y );
%
%  I want a single number, not a vector of results!
%
    err = rms ( U1 - U2 );

  end

  return
end
function biharmonic_fd2d_stencil_test ( )

%*****************************************************************************80
%
%% biharmonic_fd2d_stencil_test() tests biharmonic_fd2d_stencil().
%
%  Discussion:
%
%    For each example, evaluate and plot the solution function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    31 May 2021
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'biharmonic_fd2d_stencil_test():\n' );
  fprintf ( 1, '  biharmonic_fd2d_stencil() defines a finite\n' );
  fprintf ( 1, '  difference stencil approximation the 2D biharmonic\n' );
  fprintf ( 1, '  operator.\n' );
%
%  Get the default parameter values, so you can reset them.
%
  [ k, xmin, xmax, ymin, ymax ] = biharmonic_fd2d_parameters ( );
%
%  Set the test case.
%
  A = biharmonic_fd2d_stencil ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Stencil matrix:\n' );
  fprintf ( 1, '\n' );
  disp ( A )
%
%  Compare A*U and RHS.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compare A*U and RHS\n' );
  fprintf ( 1, '\n' );

  for k = 1 : 4

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Case = %d\n', k );
    fprintf ( 1, '\n' );

    [ k, xmin, xmax, ymin, ymax ] = biharmonic_fd2d_parameters ( ...
      k, xmin, xmax, ymin, ymax );

    n = 17;
    x = linspace ( xmin, xmax, n );
    y = linspace ( ymin, ymax, n );

    [ X, Y ] = ndgrid ( x, y );
    h = ( xmax - xmin ) / ( n - 1 );

    U = biharmonic_fd2d_solution ( X, Y );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  U(1:5,1:5)\n' );
    fprintf ( 1, '\n' );
    U(1:5,1:5)
    F = biharmonic_fd2d_rhs ( X, Y );
    S = zeros ( n, n );

    for i = 3 : n - 2
      for j = 3 : n - 2
        S(i,j) = sum ( sum ( A .* U(i-2:i+2,j-2:j+2) ) ) / h^4;        
        fprintf ( 1, '  (%2d,%2d)  %14.6g  %14.6g\n', i, j, S(i,j), F(i,j) );
      end
    end

  end

  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

