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

# include "ppc_mesh.h"
# include "ppc_twb_quad.h"

int main ( );
void demo1 ( );
void demo2 ( );
double f ( double x, double y );
void timestamp ( );

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

int main ( )

/******************************************************************************/
/*
  Purpose:
 
    ppc_twb_quad_minimal() tests ppc_twb_quad().

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 June 2024

  Author:

    John Burkardt
*/
{
  timestamp ( );

  printf ( "\n" );
  printf ( "ppc_twb_quad_minimal():\n" );
  printf ( "  C version.\n" );
  printf ( "  Test ppc_twb_quad()\n" );

  demo1 ( );
  demo2 ( );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "ppc_twb_quad_minimal():\n" );
  printf ( "  Normal end of execution.\n" );

  timestamp ( );

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

void demo1 ()

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

    demo1() retrieves and prints a quadrature rule.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    20 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt.
*/
{
  int d;
  int i;
  int n;
  struct TWB_qdat *qdat;

  printf ( "\n" );
  printf ( "demo1():\n" );
  printf ( "  Retrieve and print a single quadrature rule.\n" );	
  d = 6;
	
  qdat = twb_qdat ( &d, &n );

  printf ( "  strength = %d,  number of points = %d\n", d, n );
  printf ( "\n" );
  printf ( "            x                    y                    z                    w\n" );
  printf ( "\n" );

  for ( i = 0; i < n; i++ )
  {
    printf("%4d %20.13f %20.13f %20.13f %20.13f\n",
      i, qdat[i].lambda1,
      qdat[i].lambda2,
      qdat[i].lambda3,
      qdat[i].weight );
  }

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

void demo2 ()

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

    demo2() uses a quadrature rule to estimate an integral over a triangle.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    20 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt.
*/
{
  double area = 0.5;
  int d = 2;
  int i;
  int n;
  struct TWB_qdat *qdat;
  double sum;
  double v1[] = { 0.0, 0.0 };
  double v2[] = { 1.0, 0.0 };
  double v3[] = { 0.0, 1.0 };

  printf ( "\n" );
  printf ( "demo2():\n" );
  printf ( "  Estimate integral of x^2+y^2 over reference triangle.\n" );	
	
  qdat = twb_qdat ( &d, &n );

  sum = 0.0;
  for ( i = 0; i < n; i++ )
  {
    double lambda1 = qdat[i].lambda1;
    double lambda2 = qdat[i].lambda2;
    double lambda3 = qdat[i].lambda3;

    double w = qdat[i].weight;

    double x = lambda1 * v1[0] + lambda2 * v2[0] + lambda3 * v3[0];
    double y = lambda1 * v1[1] + lambda2 * v2[1] + lambda3 * v3[1];

    sum = sum + w * f ( x, y );
  }

  sum = sum * area / TWB_STANDARD_AREA;

  printf ( "  Integral estimate = %g\n", sum );

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

double f ( double x, double y )

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

    f() evaluates an integrand function in a triangle.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    20 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt.

  Input:

    double x, y: the position of a point.

  Output:

    double value: the value of the integrand function there.
*/
{
  double value;
  value = x * x + y * y;
  return value;
}
/******************************************************************************/

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
}