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

# include "poisson_1d_multigrid.h"

int main ( );

void test01_mono ( );
void test01_multi ( );
void test02_mono ( );
void test02_multi ( );
double exact1 ( double x );
double force1 ( double x );
double exact2 ( double x );
double force2 ( double x );
void timestamp ( );

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

int main ( )

/******************************************************************************/
/*
  Purpose:

    poisson_1d_multigrid_test() tests poisson_1d_multigrid().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    07 December 2011

  Author:

    John Burkardt
*/
{
  timestamp ( );

  printf ( "\n" );
  printf ( "poisson_1d_multigrid_test():\n" );
  printf ( "  C version\n" );
  printf ( "  Test poisson_1d_multigrid().\n" );

  test01_mono ( );
  test01_multi ( );
  test02_mono ( );
  test02_multi ( );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "poisson_1d_multigrid_test():\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

void test01_mono ( ) 

/******************************************************************************/
/*
  Purpose:

    test01_mono() tests poisson_1d_monogrid() on test case 1.

  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;

  printf ( "\n" );
  printf ( "test01_mono()\n" );
  printf ( "  poisson_1d_monogrid() solves a 1D Poisson BVP\n" );
  printf ( "  using the Gauss-Seidel method.\n" );

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

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

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

    u = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) );
    x = r8vec_linspace_new ( n + 1, a, b );

    printf ( "\n" );
    printf ( "  Mesh index K = %d\n", k );
    printf ( "  Number of intervals N=2^K = %d\n", n );
    printf ( "  Number of nodes = 2^K+1 =   %d\n", n + 1 );

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

    printf ( "\n" );
    printf ( "     I        X(I)      U(I)         U Exact(X(I))\n" );
    printf ( "\n" );
    for ( i = 0; i < n + 1; i++ )
    {
      printf ( "  %4d  %10f  %14g  %14g\n", i, x[i], u[i], exact1 ( x[i] ) );
    }

    printf ( "\n" );

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

    free ( u );
    free ( x );
  }
  return;
}
/******************************************************************************/

void test01_multi ( ) 

/******************************************************************************/
/*
  Purpose:

    test01_multi() tests poisson_1d_multigrid() on test case 1.

  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;

  printf ( "\n" );
  printf ( "test01_multi():\n" );
  printf ( "  poisson_1d_multigrid() solves a 1D Poisson BVP\n" );
  printf ( "  using the multigrid method.\n" );

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

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

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

    u = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) );
    x = r8vec_linspace_new ( n + 1, a, b );

    printf ( "\n" );
    printf ( "  Mesh index K = %d\n", k );
    printf ( "  Number of intervals N=2^K = %d\n", n );
    printf ( "  Number of nodes = 2^K+1 =   %d\n", n + 1 );

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

    printf ( "\n" );
    printf ( "     I        X(I)      U(I)         U Exact(X(I))\n" );
    printf ( "\n" );
    for ( i = 0; i < n + 1; i++ )
    {
      printf ( "  %4d  %10f  %14g  %14g\n", i, x[i], u[i], exact1 ( x[i] ) );
    }

    printf ( "\n" );

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

    free ( u );
    free ( x );
  }
  return;
}
/******************************************************************************/

double exact1 ( double x )

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

double force1 ( double x )

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

void test02_mono ( ) 

/******************************************************************************/
/*
  Purpose:

    test02_mono() tests poisson_1d_monogrid() 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;

  printf ( "\n" );
  printf ( "test02_mono()\n" );
  printf ( "  poisson_1d_monogrid() solves a 1D Poisson BVP\n" );
  printf ( "  using the Gauss-Seidel method.\n" );

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

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

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

    u = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) );
    x = r8vec_linspace_new ( n + 1, a, b );

    printf ( "\n" );
    printf ( "  Mesh index K = %d\n", k );
    printf ( "  Number of intervals N=2^K = %d\n", n );
    printf ( "  Number of nodes = 2^K+1 =   %d\n", n + 1 );

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

    printf ( "\n" );
    printf ( "     I        X(I)      U(I)         U Exact(X(I))\n" );
    printf ( "\n" );
    for ( i = 0; i < n + 1; i++ )
    {
      printf ( "  %4d  %10f  %14g  %14g\n", i, x[i], u[i], exact2 ( x[i] ) );
    }

    printf ( "\n" );

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

    free ( u );
    free ( x );
  }
  return;
}
/******************************************************************************/

void test02_multi ( ) 

/******************************************************************************/
/*
  Purpose:

    test02_multi() tests poisson_1d_multigrid() 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;

  printf ( "\n" );
  printf ( "test02_multi():\n" );
  printf ( "  poisson_1d_multigrid() solves a 1D Poisson BVP\n" );
  printf ( "  using the multigrid method.\n" );

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

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

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

    u = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) );
    x = r8vec_linspace_new ( n + 1, a, b );

    printf ( "\n" );
    printf ( "  Mesh index K = %d\n", k );
    printf ( "  Number of intervals N=2^K = %d\n", n );
    printf ( "  Number of nodes = 2^K+1 =   %d\n", n + 1 );

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

    printf ( "\n" );
    printf ( "     I        X(I)      U(I)         U Exact(X(I))\n" );
    printf ( "\n" );
    for ( i = 0; i < n + 1; i++ )
    {
      printf ( "  %4d  %10f  %14g  %14g\n", i, x[i], u[i], exact2 ( x[i] ) );
    }

    printf ( "\n" );

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

    free ( u );
    free ( x );
  }
  return;
}
/******************************************************************************/

double exact2 ( double x )

/******************************************************************************/
/*
  Purpose:

    exact2() evaluates the exact solution #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;
}
/******************************************************************************/

double force2 ( double x )

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

void timestamp ( )

/******************************************************************************/
/*
  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:

    24 September 2003

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

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}