# include <cmath>
# include <cstdlib>
# include <ctime>
# include <iomanip>
# include <iostream>

using namespace std;

# include "poisson_2d.hpp"

int main ( int argc, char *argv[] );
void rhs1 ( int nx, int ny, double **f );
void timestamp ( );
double u_exact1 ( double x, double y );
double uxxyy_exact1 ( double x, double y );
double wtime ( );

//****************************************************************************80

int main ( int argc, char *argv[] )

//****************************************************************************80
//
//  Purpose:
//
//    poisson_2d_test() tests poisson_2d().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    03 October 2024
//
//  Author:
//
//    John Burkardt
//
{
  int nx;
  int ny;
  double t_start;
  double t_stop;
 
  timestamp ( );
  cout << "\n";
  cout << "poisson_2d_test():\n";
  cout << "  C++ version\n";
  cout << "  Test poisson_2d().\n";
//
//  Problem 1.
//
  t_start = wtime ( );
  nx = 161;
  ny = 161;
  poisson_2d ( nx, ny, rhs1, u_exact1 );
  t_stop = wtime ( );
  cout << "\n";
  cout << "  Wall clock time in seconds is " << t_stop - t_start << "\n";
//
//  Terminate.
//
  cout << "\n";
  cout << "poisson_2d_test():\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

  return 0;
}
//****************************************************************************80

void rhs1 ( int nx, int ny, double **f )

//****************************************************************************80
//
//  Purpose:
//
//    rhs1() initializes the right hand side "vector".
//
//  Discussion:
//
//    It is convenient for us to set up RHS as a 2D array.  However, each
//    entry of RHS is really the right hand side of a linear system of the
//    form
//
//      A * U = F
//
//    In cases where U(I,J) is a boundary value, then the equation is simply
//
//      U(I,J) = F(i,j)
//
//    and F(I,J) holds the boundary data.
//
//    Otherwise, the equation has the form
//
//      (1/DX^2) * ( U(I+1,J)+U(I-1,J)+U(I,J-1)+U(I,J+1)-4*U(I,J) ) = F(I,J)
//
//    where DX is the spacing and F(I,J) is the value at X(I), Y(J) of
//
//      pi^2 * ( x^2 + y^2 ) * sin ( pi * x * y )
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    28 October 2011
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int NX, NY, the X and Y grid dimensions.
//
//  Output:
//
//    double **F, the initialized right hand side data.
//
{
  double fnorm;
  int i;
  int j;
  double x;
  double y;
//
//  The "boundary" entries of F store the boundary values of the solution.
//  The "interior" entries of F store the right hand sides of the Poisson equation.
//
  for ( j = 0; j < ny; j++ )
  {
    y = ( double ) ( j ) / ( double ) ( ny - 1 );
    for ( i = 0; i < nx; i++ )
    {
      x = ( double ) ( i ) / ( double ) ( nx - 1 );
      if ( i == 0 || i == nx - 1 || j == 0 || j == ny - 1 )
      {
        f[i][j] = u_exact1 ( x, y );
      }
      else
      {
        f[i][j] = - uxxyy_exact1 ( x, y );
      }
    }
  }

  fnorm = rms_norm ( nx, ny, f );

  cout << "  RMS of F = " << fnorm << "\n";

  return;
}
//****************************************************************************80

void timestamp ( )

//****************************************************************************80
//
//  Purpose:
//
//    timestamp() prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    31 May 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    08 July 2009
//
//  Author:
//
//    John Burkardt
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct std::tm *tm_ptr;
  std::time_t now;

  now = std::time ( NULL );
  tm_ptr = std::localtime ( &now );

  std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );

  std::cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}
//****************************************************************************80

double u_exact1 ( double x, double y )

//****************************************************************************80
//
//  Purpose:
//
//    u_exact1() evaluates the exact solution.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 October 2011
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double X, Y, the coordinates of a point.
//
//  Output:
//
//    double U_EXACT1, the value of the exact solution at (X,Y).
//
{
  double r8_pi = 3.141592653589793;
  double value;

  value = sin ( r8_pi * x * y );

  return value;
}
//****************************************************************************80

double uxxyy_exact1 ( double x, double y )

//****************************************************************************80
//
//  Purpose:
//
//    uxxyy_exact1() evaluates ( d/dx d/dx + d/dy d/dy ) of the exact solution.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 October 2011
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double X, Y, the coordinates of a point.
//
//  Output:
//
//    double UXXYY_EXACT1, the value of 
//    ( d/dx d/dx + d/dy d/dy ) of the exact solution at (X,Y).
//
{
  double r8_pi = 3.141592653589793;
  double value;

  value = - r8_pi * r8_pi * ( x * x + y * y ) * sin ( r8_pi * x * y );

  return value;
}
//****************************************************************************80

double wtime ( )

//****************************************************************************80
//
//  Purpose:
// 
//    wtime() reports the elapsed wallclock time.
//
//  Discussion:
//
//    The reliability of this function depends in part on the value of
//    CLOCKS_PER_SEC.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    27 April 2009
//
//  Author:
//
//    John Burkardt
//
//  Output:
//
//    double VALUE, the reading of the wall clock, in seconds.
//
{
  double value;

  value = ( double ) clock ( ) 
        / ( double ) CLOCKS_PER_SEC;

  return value;
}