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

using namespace std;

# include "poisson_1d.hpp"

int main ( );

void test01 ( );
void test02 ( );
double exact1 ( double x );
double force1 ( double x );
double exact2 ( double x );
double force2 ( double x );
void timestamp ( );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    poisson_1d_test() tests poisson_1d().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 October 2024
//
//  Author:
//
//    John Burkardt
//
{
  timestamp ( );
  cout << "\n";
  cout << "poisson_1d_test():\n";
  cout << "  C++ version\n";
  cout << "  Test poisson_1d().\n";

  test01 ( );
  test02 ( );
//
//  Terminate.
//
  cout << "\n";
  cout << "poisson_1d_test():\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

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

void test01 ( ) 

//****************************************************************************80
//
//  Purpose:
//
//    test01() tests poisson_1d() on test case #1.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    07 December 2011
//
//  Author:
//
//    John Burkardt
//
{
  double a;
  double b;
  double difmax;
  int i;
  int it_num;
  int k;
  int n;
  double *u;
  double ua;
  double ub;
  double *x;

  cout << "\n";
  cout << "test01()\n";
  cout << "  poisson_1d() solves a 1D Poisson BVP\n";
  cout << "  using the Gauss-Seidel method.\n";

  a = 0.0;
  b = 1.0;
  ua = 0.0;
  ub = 0.0;

  cout << "\n";
  cout << "  -u''(x) = 1, for 0 < x < 1\n";
  cout << "  u(0) = u(1) = 0.\n";
  cout << "  Solution is u(x) = ( -x^2 + x ) / 2\n";

  for ( k = 5; k <= 5; k++ )
  {
    n = pow ( 2, k );

    u = new double[n+1];
    x = new double[n+1];
    for ( i = 0; i <= n; i++ )
    {
      x[i] = ( ( n - i ) * a + i * b ) / n;
    }

    cout << "\n";
    cout << "  Mesh index K = " << k << "\n";
    cout << "  Number of intervals N=2^K = " << n << "\n";
    cout << "  Number of nodes = 2^K+1 =   " << n + 1 << "\n";

    poisson_1d ( n, a, b, ua, ub, force1, exact1, it_num, u );

    cout << "\n";
    cout << "     I        X(I)      U(I)         U Exact(X(I))\n";
    cout << "\n";
    for ( i = 0; i < n + 1; i++ )
    {
      cout << "  " << setw(4) << i
           << "  " << setw(10) << x[i]
           << "  " << setw(14) << u[i]
           << "  " << setw(14) << exact1 ( x[i] ) << "\n";
    }

    cout << "\n";

    difmax = 0.0;
    for ( i = 0; i < n + 1; i++ )
    {
      difmax = fmax ( difmax, fabs ( u[i] - exact1 ( x[i] ) ) );
    } 
    cout << "  Maximum error = " << difmax << "\n";
    cout << "  Number of Gauss-Seidel iterations = " << it_num << "\n";

    delete [] u;
    delete [] x;
  }
  return;
}
//****************************************************************************80

double exact1 ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    exact1() evaluates the exact solution #1.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    07 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    William Hager,
//    Applied Numerical Linear Algebra,
//    Prentice-Hall, 1988,
//    ISBN13: 978-0130412942,
//    LC: QA184.H33.
//
//  Input:
//
//    double X, the evaluation point.
//
//  Output:
//
//    double EXACT1, the value of the exact solution at X.
//
{
  double value;

  value = 0.5 * ( - x * x + x );

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

double force1 ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    force1() evaluates the forcing function for test case #1.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    07 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    William Hager,
//    Applied Numerical Linear Algebra,
//    Prentice-Hall, 1988,
//    ISBN13: 978-0130412942,
//    LC: QA184.H33.
//
//  Input:
//
//    double X, the evaluation point.
//
//  Output:
//
//    double FORCE1, the value of the forcing function at X.
//
{
  double value;

  value = 1.0;

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

void test02 ( ) 

//****************************************************************************80
//
//  Purpose:
//
//    test02() tests poisson_1d() on test case 2.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 July 2014
//
//  Author:
//
//    John Burkardt
//
{
  double a;
  double b;
  double difmax;
  int i;
  int it_num;
  int k;
  int n;
  double *u;
  double ua;
  double ub;
  double *x;

  cout << "\n";
  cout << "test02():\n";
  cout << "  poisson_1d() solves a 1D Poisson BVP\n";
  cout << "  using the Gauss-Seidel method.\n";

  a = 0.0;
  b = 1.0;
  ua = 0.0;
  ub = 0.0;

  cout << "\n";
  cout << "  -u''(x) = - x * (x+3) * exp(x), for 0 < x < 1\n";
  cout << "  u(0) = u(1) = 0.\n";
  cout << "  Solution is u(x) = x * (x-1) * exp(x)\n";

  for ( k = 5; k <= 5; k++ )
  {
    n = pow ( 2, k );

    u = new double[n+1];
    x = new double[n+1];
    for ( i = 0; i <= n; i++ )
    {
      x[i] = ( ( n - i ) * a + i * b ) / n;
    }

    cout << "\n";
    cout << "  Mesh index K = " << k << "\n";
    cout << "  Number of intervals N=2^K = " << n << "\n";
    cout << "  Number of nodes = 2^K+1 =   " << n + 1 << "\n";

    poisson_1d ( n, a, b, ua, ub, force2, exact2, it_num, u );

    cout << "\n";
    cout << "     I        X(I)      U(I)         U Exact(X(I))\n";
    cout << "\n";
    for ( i = 0; i < n + 1; i++ )
    {
      cout << "  " << setw(4) << i
           << "  " << setw(10) << x[i]
           << "  " << setw(14) << u[i]
           << "  " << setw(14) << exact2 ( x[i] ) << "\n";
    }

    cout << "\n";

    difmax = 0.0;
    for ( i = 0; i < n + 1; i++ )
    {
      difmax = fmax ( difmax, fabs ( u[i] - exact2 ( x[i] ) ) );
    } 
    cout << "  Maximum error = " << difmax << "\n";
    cout << "  Number of Gauss-Seidel iterations = " << it_num << "\n";

    delete [] u;
    delete [] x;
  }
  return;
}
//****************************************************************************80

double exact2 ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    exact2() evaluates the exact solution for test case #2.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    07 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    William Hager,
//    Applied Numerical Linear Algebra,
//    Prentice-Hall, 1988,
//    ISBN13: 978-0130412942,
//    LC: QA184.H33.
//
//  Input:
//
//    double X, the evaluation point.
//
//  Output:
//
//    double EXACT2, the value of the exact solution at X.
//
{
  double value;

  value = x * ( x - 1.0 ) * exp ( x );

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

double force2 ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    force2() evaluates the forcing function for test case #2.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    07 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    William Hager,
//    Applied Numerical Linear Algebra,
//    Prentice-Hall, 1988,
//    ISBN13: 978-0130412942,
//    LC: QA184.H33.
//
//  Input:
//
//    double X, the evaluation point.
//
//  Output:
//
//    double FORCE2, the value of the forcing function at X.
//
{
  double value;

  value = - x * ( x + 3.0 ) * exp ( x );

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