# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <time.h>

# include "poisson_2d.h"

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 ( );

/******************************************************************************/

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

/******************************************************************************/
/*
  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 ( );
  printf ( "\n" );
  printf ( "poisson_2d_test():\n" );
  printf ( "  C version\n" );
  printf ( "  Test poisson_2d().\n" );
/*
  Problem 1.
*/
  t_start = wtime ( );
  nx = 161;
  ny = 161;
  poisson_2d ( nx, ny, rhs1, u_exact1 );
  t_stop = wtime ( );
  printf ( "\n" );
  printf ( "  Wall clock time in seconds is %g\n", t_stop - t_start );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "poisson_2d_test():\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

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

/******************************************************************************/
/*
  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.
*/
{
  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 ( i = 0; i < ny; i++ )
  {
    y = ( double ) ( i ) / ( double ) ( ny - 1 );
    for ( j = 0; j < nx; j++ )
    {
      x = ( double ) ( j ) / ( double ) ( nx - 1 );
      if ( i == 0 || i == ny - 1 || j == 0 || j == nx - 1 )
      {
        f[i][j] = u_exact1 ( x, y );
      }
      else
      {
        f[i][j] = - uxxyy_exact1 ( x, y );
      }
    }
  }

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

void timestamp ( )

/******************************************************************************/
/*
  Purpose:
 
    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    17 June 2014 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    01 May 2021

  Author:

    John Burkardt
*/
{
# define TIME_SIZE 40

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

  now = time ( NULL );
  tm = localtime ( &now );

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

  printf ( "%s\n", time_buffer );

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

double u_exact1 ( double x, double y )

/******************************************************************************/
/*
  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;
}
/******************************************************************************/

double uxxyy_exact1 ( double x, double y )

/******************************************************************************/
/*
  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;
}
/******************************************************************************/

double wtime ( )

/******************************************************************************/
/*
  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;
}