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

# include "triangle_witherden_rule.h"

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

void comp_next ( int n, int k, int a[], int *more, int *h, int *t )

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

    comp_next() computes the compositions of the integer N into K parts.

  Discussion:

    A composition of the integer N into K parts is an ordered sequence
    of K nonnegative integers which sum to N.  The compositions (1,2,1)
    and (1,1,2) are considered to be distinct.

    The routine computes one composition on each call until there are no more.
    For instance, one composition of 6 into 3 parts is
    3+2+1, another would be 6+0+0.

    On the first call to this routine, set MORE = FALSE.  The routine
    will compute the first element in the sequence of compositions, and
    return it, as well as setting MORE = TRUE.  If more compositions
    are desired, call again, and again.  Each time, the routine will
    return with a new composition.

    However, when the LAST composition in the sequence is computed 
    and returned, the routine will reset MORE to FALSE, signaling that
    the end of the sequence has been reached.

    This routine originally used a STATIC statement to maintain the
    variables H and T.  I have decided (based on an wasting an
    entire morning trying to track down a problem) that it is safer
    to pass these variables as arguments, even though the user should
    never alter them.  This allows this routine to safely shuffle
    between several ongoing calculations.


    There are 28 compositions of 6 into three parts.  This routine will
    produce those compositions in the following order:

     I         A
     -     ---------
     1     6   0   0
     2     5   1   0
     3     4   2   0
     4     3   3   0
     5     2   4   0
     6     1   5   0
     7     0   6   0
     8     5   0   1
     9     4   1   1
    10     3   2   1
    11     2   3   1
    12     1   4   1
    13     0   5   1
    14     4   0   2
    15     3   1   2
    16     2   2   2
    17     1   3   2
    18     0   4   2
    19     3   0   3
    20     2   1   3
    21     1   2   3
    22     0   3   3
    23     2   0   4
    24     1   1   4
    25     0   2   4
    26     1   0   5
    27     0   1   5
    28     0   0   6

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    02 July 2008

  Author:

    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
    C version by John Burkardt.

  Reference:

    Albert Nijenhuis, Herbert Wilf,
    Combinatorial Algorithms for Computers and Calculators,
    Second Edition,
    Academic Press, 1978,
    ISBN: 0-12-519260-6,
    LC: QA164.N54.

  Parameters:

    Input, int N, the integer whose compositions are desired.

    Input, int K, the number of parts in the composition.

    Input/output, int A[K], the parts of the composition.

    Input/output, int *MORE.
    Set MORE = FALSE on first call.  It will be reset to TRUE on return
    with a new composition.  Each new call returns another composition until
    MORE is set to FALSE when the last composition has been computed
    and returned.

    Input/output, int *H, *T, two internal parameters needed for the
    computation.  The user should allocate space for these in the calling
    program, include them in the calling sequence, but never alter them!
*/
{
  int i;

  if ( !( *more ) )
  {
    *t = n;
    *h = 0;
    a[0] = n;
    for ( i = 1; i < k; i++ )
    {
       a[i] = 0;
    }
  }
  else
  {
    if ( 1 < *t )
    {
      *h = 0;
    }
    *h = *h + 1;
    *t = a[*h-1];
    a[*h-1] = 0;
    a[0] = *t - 1;
    a[*h] = a[*h] + 1;
  }

  *more = ( a[k-1] != n );

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

double *monomial_value ( int m, int n, int e[], double x[] )

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

    monomial_value() evaluates a monomial.

  Discussion:

    This routine evaluates a monomial of the form

      product ( 1 <= i <= m ) x(i)^e(i)

    The combination 0.0^0 is encountered is treated as 1.0.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 August 2014

  Author:

    John Burkardt

  Input:

    int M, the spatial dimension.

    int N, the number of evaluation points.

    int E[M], the exponents.

    double X[M*N], the point coordinates.

  Output:

    double MONOMIAL_VALUE[N], the monomial values.
*/
{
  int i;
  int j;
  double *v;

  v = ( double * ) malloc ( n * sizeof ( double ) );

  for ( j = 0; j < n; j++ )
  {
    v[j] = 1.0;
  }

  for ( i = 0; i < m; i++ )
  {
    if ( 0 != e[i] )
    {
      for ( j = 0; j < n; j++ )
      {
        v[j] = v[j] * pow ( x[i+j*m], e[i] );
      }
    }
  }

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

void r8vec_copy ( int n, double a1[], double a2[] )

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

    r8vec_copy() copies an R8VEC.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    03 July 2005

  Author:

    John Burkardt

  Input:

    int N, the number of entries in the vectors.

    double A1[N], the vector to be copied.

  Output:

    double A2[N], the copy of A1.
*/
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    a2[i] = a1[i];
  }
  return;
}
/******************************************************************************/

int rule_order ( int p )

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

    rule_order() returns the order of a quadrature rule of given precision.

  Discussion:

    The "order" of a quadrature rule is the number of points involved.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int p: the precision, 0 <= p <= 20.

  Output:

    int rule_order: the order of the rule.
*/
{
  int order;
  static int order_save[] = {
      1,
      1,   3,   6,   6,   7, 
     12,  15,  16,  19,  25, 
     28,  33,  37,  42,  49, 
     55,  60,  67,  73,  79 };

  if ( p < 0 ) 
  {
    printf ( "\n" );
    printf ( "rule_order(): Fatal error!\n" );
    printf ( "  Input p < 0.\n" );
    exit ( 1 );
  }

  if ( 20 < p )
  {
    printf ( "\n" );
    printf ( "rule_order(): Fatal error!\n" );
    printf ( "  Input 20 < p.\n" );
    exit ( 1 );
  }

  order = order_save[p];

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

void rule00 ( int n, double a[], double b[], double c[], double w[] )

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

    rule00() returns the rule of precision 0.

    The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.3333333333333334 };

  static double b_save[] = {
          0.3333333333333334 };

  static double c_save[] = {
          0.3333333333333333 };

  double w_save[] = {
          1.0000000000000000 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule01 ( int n, double a[], double b[], double c[], double w[] )

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

   rule01() returns the rule of precision 1.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.3333333333333334 };

  static double b_save[] = {
          0.3333333333333334 };

  static double c_save[] = {
          0.3333333333333333 };

  double w_save[] = {
          1.0000000000000000 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule02 ( int n, double a[], double b[], double c[], double w[] )

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

   rule02() returns the rule of precision 2.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.1666666666666667,
          0.6666666666666666,
          0.1666666666666667 };

  static double b_save[] = {
          0.6666666666666666,
          0.1666666666666667,
          0.1666666666666667 };

  static double c_save[] = {
          0.1666666666666666,
          0.1666666666666667,
          0.6666666666666666 };

  double w_save[] = {
          0.3333333333333333,
          0.3333333333333333,
          0.3333333333333333 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule03 ( int n, double a[], double b[], double c[], double w[] )

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

   rule03() returns the rule of precision 3.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.4459484909159649,
          0.1081030181680702,
          0.4459484909159649,
          0.0915762135097707,
          0.8168475729804585,
          0.0915762135097707 };

  static double b_save[] = {
          0.1081030181680702,
          0.4459484909159649,
          0.4459484909159649,
          0.8168475729804585,
          0.0915762135097707,
          0.0915762135097707 };

  static double c_save[] = {
          0.4459484909159649,
          0.4459484909159649,
          0.1081030181680702,
          0.0915762135097707,
          0.0915762135097707,
          0.8168475729804585 };

  double w_save[] = {
          0.2233815896780115,
          0.2233815896780115,
          0.2233815896780115,
          0.1099517436553219,
          0.1099517436553219,
          0.1099517436553219 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule04 ( int n, double a[], double b[], double c[], double w[] )

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

   rule04() returns the rule of precision 4.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.4459484909159649,
          0.1081030181680702,
          0.4459484909159649,
          0.0915762135097707,
          0.8168475729804585,
          0.0915762135097707 };

  static double b_save[] = {
          0.1081030181680702,
          0.4459484909159649,
          0.4459484909159649,
          0.8168475729804585,
          0.0915762135097707,
          0.0915762135097707 };

  static double c_save[] = {
          0.4459484909159649,
          0.4459484909159649,
          0.1081030181680702,
          0.0915762135097707,
          0.0915762135097707,
          0.8168475729804585 };

  double w_save[] = {
          0.2233815896780115,
          0.2233815896780115,
          0.2233815896780115,
          0.1099517436553219,
          0.1099517436553219,
          0.1099517436553219 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule05 ( int n, double a[], double b[], double c[], double w[] )

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

   rule05() returns the rule of precision 5.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.3333333333333334,
          0.1012865073234563,
          0.7974269853530873,
          0.1012865073234563,
          0.4701420641051151,
          0.0597158717897698,
          0.4701420641051151 };

  static double b_save[] = {
          0.3333333333333334,
          0.7974269853530873,
          0.1012865073234563,
          0.1012865073234563,
          0.0597158717897698,
          0.4701420641051151,
          0.4701420641051151 };

  static double c_save[] = {
          0.3333333333333333,
          0.1012865073234563,
          0.1012865073234563,
          0.7974269853530873,
          0.4701420641051151,
          0.4701420641051151,
          0.0597158717897698 };

  double w_save[] = {
          0.2250000000000000,
          0.1259391805448271,
          0.1259391805448271,
          0.1259391805448271,
          0.1323941527885062,
          0.1323941527885062,
          0.1323941527885062 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule06 ( int n, double a[], double b[], double c[], double w[] )

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

   rule06() returns the rule of precision 6.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.0630890144915022,
          0.8738219710169955,
          0.0630890144915022,
          0.2492867451709104,
          0.5014265096581791,
          0.2492867451709104,
          0.0531450498448169,
          0.6365024991213987,
          0.3103524510337844,
          0.6365024991213987,
          0.3103524510337844,
          0.0531450498448169 };

  static double b_save[] = {
          0.8738219710169955,
          0.0630890144915022,
          0.0630890144915022,
          0.5014265096581791,
          0.2492867451709104,
          0.2492867451709104,
          0.6365024991213987,
          0.0531450498448169,
          0.6365024991213987,
          0.3103524510337844,
          0.0531450498448169,
          0.3103524510337844 };

  static double c_save[] = {
          0.0630890144915022,
          0.0630890144915022,
          0.8738219710169955,
          0.2492867451709104,
          0.2492867451709104,
          0.5014265096581791,
          0.3103524510337844,
          0.3103524510337844,
          0.0531450498448169,
          0.0531450498448169,
          0.6365024991213987,
          0.6365024991213987 };

  double w_save[] = {
          0.0508449063702068,
          0.0508449063702068,
          0.0508449063702068,
          0.1167862757263794,
          0.1167862757263794,
          0.1167862757263794,
          0.0828510756183736,
          0.0828510756183736,
          0.0828510756183736,
          0.0828510756183736,
          0.0828510756183736,
          0.0828510756183736 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule07 ( int n, double a[], double b[], double c[], double w[] )

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

   rule07() returns the rule of precision 7.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.0337306485545878,
          0.9325387028908243,
          0.0337306485545878,
          0.2415773825954036,
          0.5168452348091929,
          0.2415773825954036,
          0.4743096925047182,
          0.0513806149905635,
          0.4743096925047182,
          0.0470366446525952,
          0.7542800405500532,
          0.1986833147973516,
          0.7542800405500532,
          0.1986833147973516,
          0.0470366446525952 };

  static double b_save[] = {
          0.9325387028908243,
          0.0337306485545878,
          0.0337306485545878,
          0.5168452348091929,
          0.2415773825954036,
          0.2415773825954036,
          0.0513806149905635,
          0.4743096925047182,
          0.4743096925047182,
          0.7542800405500532,
          0.0470366446525952,
          0.7542800405500532,
          0.1986833147973516,
          0.0470366446525952,
          0.1986833147973516 };

  static double c_save[] = {
          0.0337306485545878,
          0.0337306485545878,
          0.9325387028908243,
          0.2415773825954035,
          0.2415773825954035,
          0.5168452348091929,
          0.4743096925047183,
          0.4743096925047182,
          0.0513806149905636,
          0.1986833147973517,
          0.1986833147973517,
          0.0470366446525953,
          0.0470366446525953,
          0.7542800405500532,
          0.7542800405500532 };

  double w_save[] = {
          0.0165450501107921,
          0.0165450501107921,
          0.0165450501107921,
          0.1279441712301556,
          0.1279441712301556,
          0.1279441712301556,
          0.0770866461859861,
          0.0770866461859861,
          0.0770866461859861,
          0.0558787329031998,
          0.0558787329031998,
          0.0558787329031998,
          0.0558787329031998,
          0.0558787329031998,
          0.0558787329031998 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule08 ( int n, double a[], double b[], double c[], double w[] )

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

   rule08() returns the rule of precision 8.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.3333333333333334,
          0.4592925882927231,
          0.0814148234145537,
          0.4592925882927231,
          0.1705693077517602,
          0.6588613844964796,
          0.1705693077517602,
          0.0505472283170310,
          0.8989055433659381,
          0.0505472283170310,
          0.0083947774099576,
          0.7284923929554042,
          0.2631128296346381,
          0.7284923929554042,
          0.2631128296346381,
          0.0083947774099576 };

  static double b_save[] = {
          0.3333333333333334,
          0.0814148234145537,
          0.4592925882927231,
          0.4592925882927231,
          0.6588613844964796,
          0.1705693077517602,
          0.1705693077517602,
          0.8989055433659381,
          0.0505472283170310,
          0.0505472283170310,
          0.7284923929554042,
          0.0083947774099576,
          0.7284923929554042,
          0.2631128296346381,
          0.0083947774099576,
          0.2631128296346381 };

  static double c_save[] = {
          0.3333333333333333,
          0.4592925882927232,
          0.4592925882927231,
          0.0814148234145537,
          0.1705693077517603,
          0.1705693077517603,
          0.6588613844964796,
          0.0505472283170310,
          0.0505472283170310,
          0.8989055433659381,
          0.2631128296346381,
          0.2631128296346381,
          0.0083947774099576,
          0.0083947774099576,
          0.7284923929554044,
          0.7284923929554044 };

  double w_save[] = {
          0.1443156076777872,
          0.0950916342672846,
          0.0950916342672846,
          0.0950916342672846,
          0.1032173705347182,
          0.1032173705347182,
          0.1032173705347182,
          0.0324584976231981,
          0.0324584976231981,
          0.0324584976231981,
          0.0272303141744350,
          0.0272303141744350,
          0.0272303141744350,
          0.0272303141744350,
          0.0272303141744350,
          0.0272303141744350 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule09 ( int n, double a[], double b[], double c[], double w[] )

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

   rule09() returns the rule of precision 9.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.3333333333333334,
          0.4370895914929366,
          0.1258208170141267,
          0.4370895914929366,
          0.1882035356190327,
          0.6235929287619345,
          0.1882035356190327,
          0.4896825191987376,
          0.0206349616025248,
          0.4896825191987376,
          0.0447295133944527,
          0.9105409732110945,
          0.0447295133944527,
          0.0368384120547363,
          0.7411985987844980,
          0.2219629891607657,
          0.7411985987844980,
          0.2219629891607657,
          0.0368384120547363 };

  static double b_save[] = {
          0.3333333333333334,
          0.1258208170141267,
          0.4370895914929366,
          0.4370895914929366,
          0.6235929287619345,
          0.1882035356190327,
          0.1882035356190327,
          0.0206349616025248,
          0.4896825191987376,
          0.4896825191987376,
          0.9105409732110945,
          0.0447295133944527,
          0.0447295133944527,
          0.7411985987844980,
          0.0368384120547363,
          0.7411985987844980,
          0.2219629891607657,
          0.0368384120547363,
          0.2219629891607657 };

  static double c_save[] = {
          0.3333333333333333,
          0.4370895914929367,
          0.4370895914929367,
          0.1258208170141267,
          0.1882035356190327,
          0.1882035356190328,
          0.6235929287619345,
          0.4896825191987376,
          0.4896825191987376,
          0.0206349616025248,
          0.0447295133944527,
          0.0447295133944527,
          0.9105409732110946,
          0.2219629891607657,
          0.2219629891607657,
          0.0368384120547363,
          0.0368384120547363,
          0.7411985987844980,
          0.7411985987844980 };

  double w_save[] = {
          0.0971357962827988,
          0.0778275410047743,
          0.0778275410047743,
          0.0778275410047743,
          0.0796477389272102,
          0.0796477389272102,
          0.0796477389272102,
          0.0313347002271391,
          0.0313347002271391,
          0.0313347002271391,
          0.0255776756586980,
          0.0255776756586980,
          0.0255776756586980,
          0.0432835393772894,
          0.0432835393772894,
          0.0432835393772894,
          0.0432835393772894,
          0.0432835393772894,
          0.0432835393772894 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule10 ( int n, double a[], double b[], double c[], double w[] )

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

   rule10() returns the rule of precision 10.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.3333333333333334,
          0.0320553732169435,
          0.9358892535661130,
          0.0320553732169435,
          0.1421611010565644,
          0.7156777978868712,
          0.1421611010565644,
          0.3218129952888354,
          0.5300541189273440,
          0.1481328857838206,
          0.5300541189273440,
          0.1481328857838206,
          0.3218129952888354,
          0.0296198894887298,
          0.6012333286834592,
          0.3691467818278110,
          0.6012333286834592,
          0.3691467818278110,
          0.0296198894887298,
          0.0283676653399385,
          0.8079306009228790,
          0.1637017337371825,
          0.8079306009228790,
          0.1637017337371825,
          0.0283676653399385 };

  static double b_save[] = {
          0.3333333333333334,
          0.9358892535661130,
          0.0320553732169435,
          0.0320553732169435,
          0.7156777978868712,
          0.1421611010565644,
          0.1421611010565644,
          0.5300541189273440,
          0.3218129952888354,
          0.5300541189273440,
          0.1481328857838206,
          0.3218129952888354,
          0.1481328857838206,
          0.6012333286834592,
          0.0296198894887298,
          0.6012333286834592,
          0.3691467818278110,
          0.0296198894887298,
          0.3691467818278110,
          0.8079306009228790,
          0.0283676653399385,
          0.8079306009228790,
          0.1637017337371825,
          0.0283676653399385,
          0.1637017337371825 };

  static double c_save[] = {
          0.3333333333333333,
          0.0320553732169435,
          0.0320553732169435,
          0.9358892535661130,
          0.1421611010565644,
          0.1421611010565644,
          0.7156777978868712,
          0.1481328857838204,
          0.1481328857838206,
          0.3218129952888354,
          0.3218129952888354,
          0.5300541189273440,
          0.5300541189273440,
          0.3691467818278109,
          0.3691467818278109,
          0.0296198894887297,
          0.0296198894887297,
          0.6012333286834592,
          0.6012333286834592,
          0.1637017337371824,
          0.1637017337371824,
          0.0283676653399385,
          0.0283676653399385,
          0.8079306009228790,
          0.8079306009228790 };

  double w_save[] = {
          0.0817433291462860,
          0.0133529688131496,
          0.0133529688131496,
          0.0133529688131496,
          0.0459579636047447,
          0.0459579636047447,
          0.0459579636047447,
          0.0639049063964240,
          0.0639049063964240,
          0.0639049063964240,
          0.0639049063964240,
          0.0639049063964240,
          0.0639049063964240,
          0.0341846481629594,
          0.0341846481629594,
          0.0341846481629594,
          0.0341846481629594,
          0.0341846481629594,
          0.0341846481629594,
          0.0252977577072884,
          0.0252977577072884,
          0.0252977577072884,
          0.0252977577072884,
          0.0252977577072884,
          0.0252977577072884 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule11 ( int n, double a[], double b[], double c[], double w[] )

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

   rule11() returns the rule of precision 11.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.3333333333333334,
          0.0284854176143719,
          0.9430291647712562,
          0.0284854176143719,
          0.2102199567031783,
          0.5795600865936434,
          0.2102199567031783,
          0.1026354827122464,
          0.7947290345755071,
          0.1026354827122464,
          0.4958919009658909,
          0.0082161980682182,
          0.4958919009658909,
          0.4384659267643522,
          0.1230681464712956,
          0.4384659267643522,
          0.1493247886520824,
          0.8433497836618531,
          0.0073254276860645,
          0.8433497836618531,
          0.0073254276860645,
          0.1493247886520824,
          0.0460105001654300,
          0.6644083741968642,
          0.2895811256377059,
          0.6644083741968642,
          0.2895811256377059,
          0.0460105001654300 };

  static double b_save[] = {
          0.3333333333333334,
          0.9430291647712562,
          0.0284854176143719,
          0.0284854176143719,
          0.5795600865936434,
          0.2102199567031783,
          0.2102199567031783,
          0.7947290345755071,
          0.1026354827122464,
          0.1026354827122464,
          0.0082161980682182,
          0.4958919009658909,
          0.4958919009658909,
          0.1230681464712956,
          0.4384659267643522,
          0.4384659267643522,
          0.8433497836618531,
          0.1493247886520824,
          0.8433497836618531,
          0.0073254276860645,
          0.1493247886520824,
          0.0073254276860645,
          0.6644083741968642,
          0.0460105001654300,
          0.6644083741968642,
          0.2895811256377059,
          0.0460105001654300,
          0.2895811256377059 };

  static double c_save[] = {
          0.3333333333333333,
          0.0284854176143718,
          0.0284854176143718,
          0.9430291647712562,
          0.2102199567031782,
          0.2102199567031783,
          0.5795600865936434,
          0.1026354827122464,
          0.1026354827122464,
          0.7947290345755071,
          0.4958919009658909,
          0.4958919009658910,
          0.0082161980682182,
          0.4384659267643523,
          0.4384659267643523,
          0.1230681464712956,
          0.0073254276860645,
          0.0073254276860646,
          0.1493247886520823,
          0.1493247886520824,
          0.8433497836618531,
          0.8433497836618531,
          0.2895811256377059,
          0.2895811256377059,
          0.0460105001654300,
          0.0460105001654300,
          0.6644083741968642,
          0.6644083741968642 };

  double w_save[] = {
          0.0857611797322242,
          0.0104318705128947,
          0.0104318705128947,
          0.0104318705128947,
          0.0705156841117166,
          0.0705156841117166,
          0.0705156841117166,
          0.0386307592370193,
          0.0386307592370193,
          0.0386307592370193,
          0.0166062730545854,
          0.0166062730545854,
          0.0166062730545854,
          0.0673161540794683,
          0.0673161540794683,
          0.0673161540794683,
          0.0102902895729533,
          0.0102902895729533,
          0.0102902895729533,
          0.0102902895729533,
          0.0102902895729533,
          0.0102902895729533,
          0.0403324766405006,
          0.0403324766405006,
          0.0403324766405006,
          0.0403324766405006,
          0.0403324766405006,
          0.0403324766405006 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule12 ( int n, double a[], double b[], double c[], double w[] )

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

   rule12() returns the rule of precision 12.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.4882037509455415,
          0.0235924981089169,
          0.4882037509455415,
          0.1092578276593543,
          0.7814843446812914,
          0.1092578276593543,
          0.2714625070149261,
          0.4570749859701478,
          0.2714625070149261,
          0.0246463634363356,
          0.9507072731273288,
          0.0246463634363356,
          0.4401116486585931,
          0.1197767026828138,
          0.4401116486585931,
          0.2916556797383409,
          0.6853101639063919,
          0.0230341563552671,
          0.6853101639063919,
          0.0230341563552671,
          0.2916556797383409,
          0.1162960196779266,
          0.6282497516835561,
          0.2554542286385174,
          0.6282497516835561,
          0.2554542286385174,
          0.1162960196779266,
          0.8513377925102401,
          0.1272797172335894,
          0.0213824902561706,
          0.1272797172335894,
          0.0213824902561706,
          0.8513377925102401 };

  static double b_save[] = {
          0.0235924981089169,
          0.4882037509455415,
          0.4882037509455415,
          0.7814843446812914,
          0.1092578276593543,
          0.1092578276593543,
          0.4570749859701478,
          0.2714625070149261,
          0.2714625070149261,
          0.9507072731273288,
          0.0246463634363356,
          0.0246463634363356,
          0.1197767026828138,
          0.4401116486585931,
          0.4401116486585931,
          0.6853101639063919,
          0.2916556797383409,
          0.6853101639063919,
          0.0230341563552671,
          0.2916556797383409,
          0.0230341563552671,
          0.6282497516835561,
          0.1162960196779266,
          0.6282497516835561,
          0.2554542286385174,
          0.1162960196779266,
          0.2554542286385174,
          0.1272797172335894,
          0.8513377925102401,
          0.1272797172335894,
          0.0213824902561706,
          0.8513377925102401,
          0.0213824902561706 };

  static double c_save[] = {
          0.4882037509455416,
          0.4882037509455416,
          0.0235924981089169,
          0.1092578276593543,
          0.1092578276593543,
          0.7814843446812915,
          0.2714625070149261,
          0.2714625070149261,
          0.4570749859701478,
          0.0246463634363356,
          0.0246463634363356,
          0.9507072731273288,
          0.4401116486585932,
          0.4401116486585931,
          0.1197767026828138,
          0.0230341563552672,
          0.0230341563552672,
          0.2916556797383409,
          0.2916556797383409,
          0.6853101639063919,
          0.6853101639063919,
          0.2554542286385173,
          0.2554542286385173,
          0.1162960196779266,
          0.1162960196779266,
          0.6282497516835561,
          0.6282497516835561,
          0.0213824902561706,
          0.0213824902561706,
          0.8513377925102400,
          0.8513377925102400,
          0.1272797172335893,
          0.1272797172335893 };

  double w_save[] = {
          0.0242668380814520,
          0.0242668380814520,
          0.0242668380814520,
          0.0284860520688775,
          0.0284860520688775,
          0.0284860520688775,
          0.0625412131959028,
          0.0625412131959028,
          0.0625412131959028,
          0.0079316425099736,
          0.0079316425099736,
          0.0079316425099736,
          0.0499183349280609,
          0.0499183349280609,
          0.0499183349280609,
          0.0217835850386076,
          0.0217835850386076,
          0.0217835850386076,
          0.0217835850386076,
          0.0217835850386076,
          0.0217835850386076,
          0.0432273636594142,
          0.0432273636594142,
          0.0432273636594142,
          0.0432273636594142,
          0.0432273636594142,
          0.0432273636594142,
          0.0150836775765114,
          0.0150836775765114,
          0.0150836775765114,
          0.0150836775765114,
          0.0150836775765114,
          0.0150836775765114 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule13 ( int n, double a[], double b[], double c[], double w[] )

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

   rule13() returns the rule of precision 13.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.3333333333333334,
          0.4890769464525394,
          0.0218461070949213,
          0.4890769464525394,
          0.2213722862918329,
          0.5572554274163342,
          0.2213722862918329,
          0.4269414142598004,
          0.1461171714803992,
          0.4269414142598004,
          0.0215096811088432,
          0.9569806377823136,
          0.0215096811088432,
          0.0878954830321973,
          0.7485071158999522,
          0.1635974010678505,
          0.7485071158999522,
          0.1635974010678505,
          0.0878954830321973,
          0.1109220428034634,
          0.8647077702954428,
          0.0243701869010938,
          0.8647077702954428,
          0.0243701869010938,
          0.1109220428034634,
          0.3084417608921178,
          0.6235459955536755,
          0.0680122435542067,
          0.6235459955536755,
          0.0680122435542067,
          0.3084417608921178,
          0.0051263891023824,
          0.7223577931241880,
          0.2725158177734297,
          0.7223577931241880,
          0.2725158177734297,
          0.0051263891023824 };

  static double b_save[] = {
          0.3333333333333334,
          0.0218461070949213,
          0.4890769464525394,
          0.4890769464525394,
          0.5572554274163342,
          0.2213722862918329,
          0.2213722862918329,
          0.1461171714803992,
          0.4269414142598004,
          0.4269414142598004,
          0.9569806377823136,
          0.0215096811088432,
          0.0215096811088432,
          0.7485071158999522,
          0.0878954830321973,
          0.7485071158999522,
          0.1635974010678505,
          0.0878954830321973,
          0.1635974010678505,
          0.8647077702954428,
          0.1109220428034634,
          0.8647077702954428,
          0.0243701869010938,
          0.1109220428034634,
          0.0243701869010938,
          0.6235459955536755,
          0.3084417608921178,
          0.6235459955536755,
          0.0680122435542067,
          0.3084417608921178,
          0.0680122435542067,
          0.7223577931241880,
          0.0051263891023824,
          0.7223577931241880,
          0.2725158177734297,
          0.0051263891023824,
          0.2725158177734297 };

  static double c_save[] = {
          0.3333333333333333,
          0.4890769464525393,
          0.4890769464525394,
          0.0218461070949213,
          0.2213722862918328,
          0.2213722862918329,
          0.5572554274163342,
          0.4269414142598004,
          0.4269414142598004,
          0.1461171714803992,
          0.0215096811088432,
          0.0215096811088433,
          0.9569806377823137,
          0.1635974010678505,
          0.1635974010678505,
          0.0878954830321973,
          0.0878954830321973,
          0.7485071158999522,
          0.7485071158999522,
          0.0243701869010938,
          0.0243701869010938,
          0.1109220428034634,
          0.1109220428034634,
          0.8647077702954428,
          0.8647077702954428,
          0.0680122435542067,
          0.0680122435542068,
          0.3084417608921177,
          0.3084417608921178,
          0.6235459955536755,
          0.6235459955536755,
          0.2725158177734297,
          0.2725158177734296,
          0.0051263891023823,
          0.0051263891023823,
          0.7223577931241879,
          0.7223577931241879 };

  double w_save[] = {
          0.0679600365868316,
          0.0239944019288947,
          0.0239944019288947,
          0.0239944019288947,
          0.0582784851192000,
          0.0582784851192000,
          0.0582784851192000,
          0.0556019675304533,
          0.0556019675304533,
          0.0556019675304533,
          0.0060523371035392,
          0.0060523371035392,
          0.0060523371035392,
          0.0241790398115938,
          0.0241790398115938,
          0.0241790398115938,
          0.0241790398115938,
          0.0241790398115938,
          0.0241790398115938,
          0.0149654011051657,
          0.0149654011051657,
          0.0149654011051657,
          0.0149654011051657,
          0.0149654011051657,
          0.0149654011051657,
          0.0346412761408484,
          0.0346412761408484,
          0.0346412761408484,
          0.0346412761408484,
          0.0346412761408484,
          0.0346412761408484,
          0.0095906810035433,
          0.0095906810035433,
          0.0095906810035433,
          0.0095906810035433,
          0.0095906810035433,
          0.0095906810035433 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule14 ( int n, double a[], double b[], double c[], double w[] )

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

   rule14() returns the rule of precision 14.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.1772055324125434,
          0.6455889351749131,
          0.1772055324125434,
          0.4176447193404539,
          0.1647105613190922,
          0.4176447193404539,
          0.0617998830908726,
          0.8764002338182548,
          0.0617998830908726,
          0.4889639103621786,
          0.0220721792756427,
          0.4889639103621786,
          0.2734775283088386,
          0.4530449433823227,
          0.2734775283088386,
          0.0193909612487010,
          0.9612180775025979,
          0.0193909612487010,
          0.2983728821362578,
          0.6869801678080878,
          0.0146469500556544,
          0.6869801678080878,
          0.0146469500556544,
          0.2983728821362578,
          0.0571247574036479,
          0.7706085547749965,
          0.1722666878213556,
          0.7706085547749965,
          0.1722666878213556,
          0.0571247574036479,
          0.3368614597963450,
          0.5702222908466832,
          0.0929162493569718,
          0.5702222908466832,
          0.0929162493569718,
          0.3368614597963450,
          0.0012683309328720,
          0.8797571713701711,
          0.1189744976969568,
          0.8797571713701711,
          0.1189744976969568,
          0.0012683309328720 };

  static double b_save[] = {
          0.6455889351749131,
          0.1772055324125434,
          0.1772055324125434,
          0.1647105613190922,
          0.4176447193404539,
          0.4176447193404539,
          0.8764002338182548,
          0.0617998830908726,
          0.0617998830908726,
          0.0220721792756427,
          0.4889639103621786,
          0.4889639103621786,
          0.4530449433823227,
          0.2734775283088386,
          0.2734775283088386,
          0.9612180775025979,
          0.0193909612487010,
          0.0193909612487010,
          0.6869801678080878,
          0.2983728821362578,
          0.6869801678080878,
          0.0146469500556544,
          0.2983728821362578,
          0.0146469500556544,
          0.7706085547749965,
          0.0571247574036479,
          0.7706085547749965,
          0.1722666878213556,
          0.0571247574036479,
          0.1722666878213556,
          0.5702222908466832,
          0.3368614597963450,
          0.5702222908466832,
          0.0929162493569718,
          0.3368614597963450,
          0.0929162493569718,
          0.8797571713701711,
          0.0012683309328720,
          0.8797571713701711,
          0.1189744976969568,
          0.0012683309328720,
          0.1189744976969568 };

  static double c_save[] = {
          0.1772055324125434,
          0.1772055324125434,
          0.6455889351749131,
          0.4176447193404539,
          0.4176447193404539,
          0.1647105613190921,
          0.0617998830908726,
          0.0617998830908726,
          0.8764002338182548,
          0.4889639103621787,
          0.4889639103621787,
          0.0220721792756428,
          0.2734775283088386,
          0.2734775283088386,
          0.4530449433823227,
          0.0193909612487011,
          0.0193909612487011,
          0.9612180775025979,
          0.0146469500556543,
          0.0146469500556544,
          0.2983728821362577,
          0.2983728821362578,
          0.6869801678080878,
          0.6869801678080878,
          0.1722666878213557,
          0.1722666878213557,
          0.0571247574036480,
          0.0571247574036480,
          0.7706085547749966,
          0.7706085547749966,
          0.0929162493569718,
          0.0929162493569718,
          0.3368614597963450,
          0.3368614597963450,
          0.5702222908466832,
          0.5702222908466832,
          0.1189744976969569,
          0.1189744976969569,
          0.0012683309328720,
          0.0012683309328720,
          0.8797571713701711,
          0.8797571713701711 };

  double w_save[] = {
          0.0421625887369930,
          0.0421625887369930,
          0.0421625887369930,
          0.0327883535441253,
          0.0327883535441253,
          0.0327883535441253,
          0.0144336996697767,
          0.0144336996697767,
          0.0144336996697767,
          0.0218835813694289,
          0.0218835813694289,
          0.0218835813694289,
          0.0517741045072916,
          0.0517741045072916,
          0.0517741045072916,
          0.0049234036024001,
          0.0049234036024001,
          0.0049234036024001,
          0.0144363081135338,
          0.0144363081135338,
          0.0144363081135338,
          0.0144363081135338,
          0.0144363081135338,
          0.0144363081135338,
          0.0246657532125637,
          0.0246657532125637,
          0.0246657532125637,
          0.0246657532125637,
          0.0246657532125637,
          0.0246657532125637,
          0.0385715107870607,
          0.0385715107870607,
          0.0385715107870607,
          0.0385715107870607,
          0.0385715107870607,
          0.0385715107870607,
          0.0050102288385007,
          0.0050102288385007,
          0.0050102288385007,
          0.0050102288385007,
          0.0050102288385007,
          0.0050102288385007 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule15 ( int n, double a[], double b[], double c[], double w[] )

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

   rule15() returns the rule of precision 15.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.3333333333333334,
          0.4053622141339755,
          0.1892755717320491,
          0.4053622141339755,
          0.0701735528999860,
          0.8596528942000279,
          0.0701735528999860,
          0.4741706814380198,
          0.0516586371239604,
          0.4741706814380198,
          0.2263787134203496,
          0.5472425731593008,
          0.2263787134203496,
          0.4949969567691262,
          0.0100060864617476,
          0.4949969567691262,
          0.0158117262509886,
          0.9683765474980227,
          0.0158117262509886,
          0.0183761123856811,
          0.6669756448018681,
          0.3146482428124509,
          0.6669756448018681,
          0.3146482428124509,
          0.0183761123856811,
          0.0091392370373084,
          0.9199121577262361,
          0.0709486052364555,
          0.9199121577262361,
          0.0709486052364555,
          0.0091392370373084,
          0.1905355894763939,
          0.7152223569314506,
          0.0942420535921554,
          0.7152223569314506,
          0.0942420535921554,
          0.1905355894763939,
          0.1680686452224144,
          0.8132926410494192,
          0.0186387137281664,
          0.8132926410494192,
          0.0186387137281664,
          0.1680686452224144,
          0.3389506114752772,
          0.5652526648771142,
          0.0957967236476086,
          0.5652526648771142,
          0.0957967236476086,
          0.3389506114752772 };

  static double b_save[] = {
          0.3333333333333334,
          0.1892755717320491,
          0.4053622141339755,
          0.4053622141339755,
          0.8596528942000279,
          0.0701735528999860,
          0.0701735528999860,
          0.0516586371239604,
          0.4741706814380198,
          0.4741706814380198,
          0.5472425731593008,
          0.2263787134203496,
          0.2263787134203496,
          0.0100060864617476,
          0.4949969567691262,
          0.4949969567691262,
          0.9683765474980227,
          0.0158117262509886,
          0.0158117262509886,
          0.6669756448018681,
          0.0183761123856811,
          0.6669756448018681,
          0.3146482428124509,
          0.0183761123856811,
          0.3146482428124509,
          0.9199121577262361,
          0.0091392370373084,
          0.9199121577262361,
          0.0709486052364555,
          0.0091392370373084,
          0.0709486052364555,
          0.7152223569314506,
          0.1905355894763939,
          0.7152223569314506,
          0.0942420535921554,
          0.1905355894763939,
          0.0942420535921554,
          0.8132926410494192,
          0.1680686452224144,
          0.8132926410494192,
          0.0186387137281664,
          0.1680686452224144,
          0.0186387137281664,
          0.5652526648771142,
          0.3389506114752772,
          0.5652526648771142,
          0.0957967236476086,
          0.3389506114752772,
          0.0957967236476086 };

  static double c_save[] = {
          0.3333333333333333,
          0.4053622141339754,
          0.4053622141339754,
          0.1892755717320490,
          0.0701735528999861,
          0.0701735528999861,
          0.8596528942000280,
          0.4741706814380198,
          0.4741706814380198,
          0.0516586371239605,
          0.2263787134203497,
          0.2263787134203497,
          0.5472425731593008,
          0.4949969567691261,
          0.4949969567691263,
          0.0100060864617476,
          0.0158117262509887,
          0.0158117262509887,
          0.9683765474980227,
          0.3146482428124509,
          0.3146482428124509,
          0.0183761123856810,
          0.0183761123856810,
          0.6669756448018680,
          0.6669756448018681,
          0.0709486052364555,
          0.0709486052364554,
          0.0091392370373085,
          0.0091392370373083,
          0.9199121577262361,
          0.9199121577262361,
          0.0942420535921553,
          0.0942420535921554,
          0.1905355894763939,
          0.1905355894763940,
          0.7152223569314506,
          0.7152223569314506,
          0.0186387137281663,
          0.0186387137281664,
          0.1680686452224143,
          0.1680686452224144,
          0.8132926410494192,
          0.8132926410494192,
          0.0957967236476086,
          0.0957967236476086,
          0.3389506114752772,
          0.3389506114752772,
          0.5652526648771142,
          0.5652526648771142 };

  double w_save[] = {
          0.0443353873821841,
          0.0427137815714606,
          0.0427137815714606,
          0.0427137815714606,
          0.0164447375626252,
          0.0164447375626252,
          0.0164447375626252,
          0.0173961480007634,
          0.0173961480007634,
          0.0173961480007634,
          0.0467833617287096,
          0.0467833617287096,
          0.0467833617287096,
          0.0095738461824601,
          0.0095738461824601,
          0.0095738461824601,
          0.0029607746379054,
          0.0029607746379054,
          0.0029607746379054,
          0.0156025728305760,
          0.0156025728305760,
          0.0156025728305760,
          0.0156025728305760,
          0.0156025728305760,
          0.0156025728305760,
          0.0040298533720181,
          0.0040298533720181,
          0.0040298533720181,
          0.0040298533720181,
          0.0040298533720181,
          0.0040298533720181,
          0.0287205869252013,
          0.0287205869252013,
          0.0287205869252013,
          0.0287205869252013,
          0.0287205869252013,
          0.0287205869252013,
          0.0116726211815758,
          0.0116726211815758,
          0.0116726211815758,
          0.0116726211815758,
          0.0116726211815758,
          0.0116726211815758,
          0.0313154762849693,
          0.0313154762849693,
          0.0313154762849693,
          0.0313154762849693,
          0.0313154762849693,
          0.0313154762849693 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule16 ( int n, double a[], double b[], double c[], double w[] )

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

   rule16() returns the rule of precision 16.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.3333333333333334,
          0.2459900704671417,
          0.5080198590657166,
          0.2459900704671417,
          0.4155848968854205,
          0.1688302062291590,
          0.4155848968854205,
          0.0853555665867003,
          0.8292888668265994,
          0.0853555665867003,
          0.1619186441912712,
          0.6761627116174576,
          0.1619186441912712,
          0.5000000000000000,
          0.0000000000000000,
          0.5000000000000000,
          0.4752807275459421,
          0.0494385449081158,
          0.4752807275459421,
          0.0547551749147031,
          0.7541700614447677,
          0.1910747636405291,
          0.7541700614447677,
          0.1910747636405291,
          0.0547551749147031,
          0.0232034277688137,
          0.9682443680309587,
          0.0085522042002276,
          0.9682443680309587,
          0.0085522042002276,
          0.0232034277688137,
          0.0189317782804059,
          0.6493036982454464,
          0.3317645234741476,
          0.6493036982454464,
          0.3317645234741476,
          0.0189317782804059,
          0.0190301297436974,
          0.9002737032704295,
          0.0806961669858730,
          0.9002737032704295,
          0.0806961669858730,
          0.0190301297436974,
          0.1026061902393981,
          0.5891488405642479,
          0.3082449691963540,
          0.5891488405642479,
          0.3082449691963540,
          0.1026061902393981,
          0.0059363500168222,
          0.8066218674993957,
          0.1874417824837821,
          0.8066218674993957,
          0.1874417824837821,
          0.0059363500168222 };

  static double b_save[] = {
          0.3333333333333334,
          0.5080198590657166,
          0.2459900704671417,
          0.2459900704671417,
          0.1688302062291590,
          0.4155848968854205,
          0.4155848968854205,
          0.8292888668265994,
          0.0853555665867003,
          0.0853555665867003,
          0.6761627116174576,
          0.1619186441912712,
          0.1619186441912712,
          0.0000000000000000,
          0.5000000000000000,
          0.5000000000000000,
          0.0494385449081158,
          0.4752807275459421,
          0.4752807275459421,
          0.7541700614447677,
          0.0547551749147031,
          0.7541700614447677,
          0.1910747636405291,
          0.0547551749147031,
          0.1910747636405291,
          0.9682443680309587,
          0.0232034277688137,
          0.9682443680309587,
          0.0085522042002276,
          0.0232034277688137,
          0.0085522042002276,
          0.6493036982454464,
          0.0189317782804059,
          0.6493036982454464,
          0.3317645234741476,
          0.0189317782804059,
          0.3317645234741476,
          0.9002737032704295,
          0.0190301297436974,
          0.9002737032704295,
          0.0806961669858730,
          0.0190301297436974,
          0.0806961669858730,
          0.5891488405642479,
          0.1026061902393981,
          0.5891488405642479,
          0.3082449691963540,
          0.1026061902393981,
          0.3082449691963540,
          0.8066218674993957,
          0.0059363500168222,
          0.8066218674993957,
          0.1874417824837821,
          0.0059363500168222,
          0.1874417824837821 };

  static double c_save[] = {
          0.3333333333333333,
          0.2459900704671418,
          0.2459900704671417,
          0.5080198590657166,
          0.4155848968854206,
          0.4155848968854206,
          0.1688302062291590,
          0.0853555665867003,
          0.0853555665867003,
          0.8292888668265993,
          0.1619186441912712,
          0.1619186441912712,
          0.6761627116174576,
          0.5000000000000000,
          0.5000000000000000,
          0.0000000000000000,
          0.4752807275459421,
          0.4752807275459421,
          0.0494385449081158,
          0.1910747636405291,
          0.1910747636405292,
          0.0547551749147032,
          0.0547551749147033,
          0.7541700614447678,
          0.7541700614447678,
          0.0085522042002275,
          0.0085522042002275,
          0.0232034277688137,
          0.0232034277688137,
          0.9682443680309587,
          0.9682443680309587,
          0.3317645234741476,
          0.3317645234741476,
          0.0189317782804059,
          0.0189317782804059,
          0.6493036982454464,
          0.6493036982454464,
          0.0806961669858730,
          0.0806961669858730,
          0.0190301297436974,
          0.0190301297436974,
          0.9002737032704295,
          0.9002737032704295,
          0.3082449691963540,
          0.3082449691963540,
          0.1026061902393980,
          0.1026061902393981,
          0.5891488405642479,
          0.5891488405642479,
          0.1874417824837822,
          0.1874417824837821,
          0.0059363500168224,
          0.0059363500168222,
          0.8066218674993957,
          0.8066218674993957 };

  double w_save[] = {
          0.0452645660738188,
          0.0410929231436990,
          0.0410929231436990,
          0.0410929231436990,
          0.0407118333124254,
          0.0407118333124254,
          0.0407118333124254,
          0.0147816346902244,
          0.0147816346902244,
          0.0147816346902244,
          0.0294184096989881,
          0.0294184096989881,
          0.0294184096989881,
          0.0044185463121506,
          0.0044185463121506,
          0.0044185463121506,
          0.0259743332982772,
          0.0259743332982772,
          0.0259743332982772,
          0.0189382724644157,
          0.0189382724644157,
          0.0189382724644157,
          0.0189382724644157,
          0.0189382724644157,
          0.0189382724644157,
          0.0016544667148350,
          0.0016544667148350,
          0.0016544667148350,
          0.0016544667148350,
          0.0016544667148350,
          0.0016544667148350,
          0.0150086017842858,
          0.0150086017842858,
          0.0150086017842858,
          0.0150086017842858,
          0.0150086017842858,
          0.0150086017842858,
          0.0079475939333925,
          0.0079475939333925,
          0.0079475939333925,
          0.0079475939333925,
          0.0079475939333925,
          0.0079475939333925,
          0.0319836100793701,
          0.0319836100793701,
          0.0319836100793701,
          0.0319836100793701,
          0.0319836100793701,
          0.0319836100793701,
          0.0053911871168488,
          0.0053911871168488,
          0.0053911871168488,
          0.0053911871168488,
          0.0053911871168488,
          0.0053911871168488 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule17 ( int n, double a[], double b[], double c[], double w[] )

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

   rule17() returns the rule of precision 17.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.4171034443615992,
          0.1657931112768016,
          0.4171034443615992,
          0.0147554916607540,
          0.9704890166784921,
          0.0147554916607540,
          0.4655978716188903,
          0.0688042567622194,
          0.4655978716188903,
          0.1803581162663706,
          0.6392837674672588,
          0.1803581162663706,
          0.0666540634795969,
          0.8666918730408062,
          0.0666540634795969,
          0.2857065024365866,
          0.4285869951268267,
          0.2857065024365866,
          0.0160176423621193,
          0.8247900701650881,
          0.1591922874727927,
          0.8247900701650881,
          0.1591922874727927,
          0.0160176423621193,
          0.3062815917461865,
          0.6263690303864522,
          0.0673493778673612,
          0.6263690303864522,
          0.0673493778673612,
          0.3062815917461865,
          0.0132296727600869,
          0.5712948679446841,
          0.4154754592952291,
          0.5712948679446841,
          0.4154754592952291,
          0.0132296727600869,
          0.0780423405682824,
          0.7532351459364581,
          0.1687225134952595,
          0.7532351459364581,
          0.1687225134952595,
          0.0780423405682824,
          0.0131358708340027,
          0.7150722591106424,
          0.2717918700553549,
          0.7150722591106424,
          0.2717918700553549,
          0.0131358708340027,
          0.0115751759031806,
          0.9159193532978169,
          0.0725054707990024,
          0.9159193532978169,
          0.0725054707990024,
          0.0115751759031806,
          0.1575054779268699,
          0.5432755795961598,
          0.2992189424769703,
          0.5432755795961598,
          0.2992189424769703,
          0.1575054779268699 };

  static double b_save[] = {
          0.1657931112768016,
          0.4171034443615992,
          0.4171034443615992,
          0.9704890166784921,
          0.0147554916607540,
          0.0147554916607540,
          0.0688042567622194,
          0.4655978716188903,
          0.4655978716188903,
          0.6392837674672588,
          0.1803581162663706,
          0.1803581162663706,
          0.8666918730408062,
          0.0666540634795969,
          0.0666540634795969,
          0.4285869951268267,
          0.2857065024365866,
          0.2857065024365866,
          0.8247900701650881,
          0.0160176423621193,
          0.8247900701650881,
          0.1591922874727927,
          0.0160176423621193,
          0.1591922874727927,
          0.6263690303864522,
          0.3062815917461865,
          0.6263690303864522,
          0.0673493778673612,
          0.3062815917461865,
          0.0673493778673612,
          0.5712948679446841,
          0.0132296727600869,
          0.5712948679446841,
          0.4154754592952291,
          0.0132296727600869,
          0.4154754592952291,
          0.7532351459364581,
          0.0780423405682824,
          0.7532351459364581,
          0.1687225134952595,
          0.0780423405682824,
          0.1687225134952595,
          0.7150722591106424,
          0.0131358708340027,
          0.7150722591106424,
          0.2717918700553549,
          0.0131358708340027,
          0.2717918700553549,
          0.9159193532978169,
          0.0115751759031806,
          0.9159193532978169,
          0.0725054707990024,
          0.0115751759031806,
          0.0725054707990024,
          0.5432755795961598,
          0.1575054779268699,
          0.5432755795961598,
          0.2992189424769703,
          0.1575054779268699,
          0.2992189424769703 };

  static double c_save[] = {
          0.4171034443615993,
          0.4171034443615992,
          0.1657931112768016,
          0.0147554916607540,
          0.0147554916607540,
          0.9704890166784921,
          0.4655978716188902,
          0.4655978716188903,
          0.0688042567622194,
          0.1803581162663707,
          0.1803581162663705,
          0.6392837674672588,
          0.0666540634795969,
          0.0666540634795969,
          0.8666918730408062,
          0.2857065024365866,
          0.2857065024365866,
          0.4285869951268267,
          0.1591922874727927,
          0.1591922874727927,
          0.0160176423621192,
          0.0160176423621192,
          0.8247900701650881,
          0.8247900701650881,
          0.0673493778673613,
          0.0673493778673613,
          0.3062815917461865,
          0.3062815917461866,
          0.6263690303864523,
          0.6263690303864523,
          0.4154754592952290,
          0.4154754592952290,
          0.0132296727600869,
          0.0132296727600869,
          0.5712948679446841,
          0.5712948679446841,
          0.1687225134952595,
          0.1687225134952595,
          0.0780423405682824,
          0.0780423405682824,
          0.7532351459364581,
          0.7532351459364581,
          0.2717918700553549,
          0.2717918700553549,
          0.0131358708340027,
          0.0131358708340027,
          0.7150722591106424,
          0.7150722591106424,
          0.0725054707990024,
          0.0725054707990025,
          0.0115751759031806,
          0.0115751759031806,
          0.9159193532978169,
          0.9159193532978169,
          0.2992189424769703,
          0.2992189424769702,
          0.1575054779268699,
          0.1575054779268699,
          0.5432755795961597,
          0.5432755795961598 };

  double w_save[] = {
          0.0273109265281021,
          0.0273109265281021,
          0.0273109265281021,
          0.0027738875776376,
          0.0027738875776376,
          0.0027738875776376,
          0.0250194509504974,
          0.0250194509504974,
          0.0250194509504974,
          0.0263126305880180,
          0.0263126305880180,
          0.0263126305880180,
          0.0124590008023054,
          0.0124590008023054,
          0.0124590008023054,
          0.0377162371527953,
          0.0377162371527953,
          0.0377162371527953,
          0.0079783002059296,
          0.0079783002059296,
          0.0079783002059296,
          0.0079783002059296,
          0.0079783002059296,
          0.0079783002059296,
          0.0224877725466911,
          0.0224877725466911,
          0.0224877725466911,
          0.0224877725466911,
          0.0224877725466911,
          0.0224877725466911,
          0.0103984399558395,
          0.0103984399558395,
          0.0103984399558395,
          0.0103984399558395,
          0.0103984399558395,
          0.0103984399558395,
          0.0205578983204545,
          0.0205578983204545,
          0.0205578983204545,
          0.0205578983204545,
          0.0205578983204545,
          0.0205578983204545,
          0.0086922145010012,
          0.0086922145010012,
          0.0086922145010012,
          0.0086922145010012,
          0.0086922145010012,
          0.0086922145010012,
          0.0045843484017359,
          0.0045843484017359,
          0.0045843484017359,
          0.0045843484017359,
          0.0045843484017359,
          0.0045843484017359,
          0.0261716259353370,
          0.0261716259353370,
          0.0261716259353370,
          0.0261716259353370,
          0.0261716259353370,
          0.0261716259353370 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule18 ( int n, double a[], double b[], double c[], double w[] )

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

   rule18() returns the rule of precision 18.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.3333333333333334,
          0.3999556280675762,
          0.2000887438648475,
          0.3999556280675762,
          0.4875803015748695,
          0.0248393968502609,
          0.4875803015748695,
          0.4618095064064492,
          0.0763809871871015,
          0.4618095064064492,
          0.2422647025142720,
          0.5154705949714561,
          0.2422647025142720,
          0.0388302560886856,
          0.9223394878226288,
          0.0388302560886856,
          0.0919477421216432,
          0.8161045157567136,
          0.0919477421216432,
          0.0458049158598608,
          0.7703723762146752,
          0.1838227079254640,
          0.7703723762146752,
          0.1838227079254640,
          0.0458049158598608,
          0.2063492574338379,
          0.6709539851942345,
          0.1226967573719275,
          0.6709539851942345,
          0.1226967573719275,
          0.2063492574338379,
          0.0038976110334734,
          0.6004189546342569,
          0.3956834343322697,
          0.6004189546342569,
          0.3956834343322697,
          0.0038976110334734,
          0.0134620167414450,
          0.8783421894675217,
          0.1081957937910333,
          0.8783421894675217,
          0.1081957937910333,
          0.0134620167414450,
          0.0402602834699081,
          0.6399880920047146,
          0.3197516245253773,
          0.6399880920047146,
          0.3197516245253773,
          0.0402602834699081,
          0.0052983351866098,
          0.7589294798551984,
          0.2357721849581917,
          0.7589294798551984,
          0.2357721849581917,
          0.0052983351866098,
          0.0005483600420423,
          0.9723607289627957,
          0.0270909109951620,
          0.9723607289627957,
          0.0270909109951620,
          0.0005483600420423,
          0.1205876951639246,
          0.5459187753861946,
          0.3334935294498808,
          0.5459187753861946,
          0.3334935294498808,
          0.1205876951639246 };

  static double b_save[] = {
          0.3333333333333334,
          0.2000887438648475,
          0.3999556280675762,
          0.3999556280675762,
          0.0248393968502609,
          0.4875803015748695,
          0.4875803015748695,
          0.0763809871871015,
          0.4618095064064492,
          0.4618095064064492,
          0.5154705949714561,
          0.2422647025142720,
          0.2422647025142720,
          0.9223394878226288,
          0.0388302560886856,
          0.0388302560886856,
          0.8161045157567136,
          0.0919477421216432,
          0.0919477421216432,
          0.7703723762146752,
          0.0458049158598608,
          0.7703723762146752,
          0.1838227079254640,
          0.0458049158598608,
          0.1838227079254640,
          0.6709539851942345,
          0.2063492574338379,
          0.6709539851942345,
          0.1226967573719275,
          0.2063492574338379,
          0.1226967573719275,
          0.6004189546342569,
          0.0038976110334734,
          0.6004189546342569,
          0.3956834343322697,
          0.0038976110334734,
          0.3956834343322697,
          0.8783421894675217,
          0.0134620167414450,
          0.8783421894675217,
          0.1081957937910333,
          0.0134620167414450,
          0.1081957937910333,
          0.6399880920047146,
          0.0402602834699081,
          0.6399880920047146,
          0.3197516245253773,
          0.0402602834699081,
          0.3197516245253773,
          0.7589294798551984,
          0.0052983351866098,
          0.7589294798551984,
          0.2357721849581917,
          0.0052983351866098,
          0.2357721849581917,
          0.9723607289627957,
          0.0005483600420423,
          0.9723607289627957,
          0.0270909109951620,
          0.0005483600420423,
          0.0270909109951620,
          0.5459187753861946,
          0.1205876951639246,
          0.5459187753861946,
          0.3334935294498808,
          0.1205876951639246,
          0.3334935294498808 };

  static double c_save[] = {
          0.3333333333333333,
          0.3999556280675762,
          0.3999556280675762,
          0.2000887438648475,
          0.4875803015748696,
          0.4875803015748695,
          0.0248393968502609,
          0.4618095064064492,
          0.4618095064064492,
          0.0763809871871015,
          0.2422647025142719,
          0.2422647025142719,
          0.5154705949714561,
          0.0388302560886855,
          0.0388302560886856,
          0.9223394878226288,
          0.0919477421216433,
          0.0919477421216433,
          0.8161045157567136,
          0.1838227079254640,
          0.1838227079254640,
          0.0458049158598608,
          0.0458049158598608,
          0.7703723762146752,
          0.7703723762146752,
          0.1226967573719275,
          0.1226967573719275,
          0.2063492574338379,
          0.2063492574338379,
          0.6709539851942345,
          0.6709539851942345,
          0.3956834343322697,
          0.3956834343322697,
          0.0038976110334734,
          0.0038976110334734,
          0.6004189546342569,
          0.6004189546342569,
          0.1081957937910333,
          0.1081957937910333,
          0.0134620167414450,
          0.0134620167414450,
          0.8783421894675217,
          0.8783421894675217,
          0.3197516245253773,
          0.3197516245253773,
          0.0402602834699081,
          0.0402602834699081,
          0.6399880920047146,
          0.6399880920047146,
          0.2357721849581917,
          0.2357721849581917,
          0.0052983351866098,
          0.0052983351866098,
          0.7589294798551984,
          0.7589294798551984,
          0.0270909109951620,
          0.0270909109951620,
          0.0005483600420423,
          0.0005483600420423,
          0.9723607289627956,
          0.9723607289627956,
          0.3334935294498808,
          0.3334935294498808,
          0.1205876951639246,
          0.1205876951639246,
          0.5459187753861946,
          0.5459187753861946 };

  double w_save[] = {
          0.0363557353014267,
          0.0333044700333901,
          0.0333044700333901,
          0.0333044700333901,
          0.0120466476339997,
          0.0120466476339997,
          0.0120466476339997,
          0.0189491715067789,
          0.0189491715067789,
          0.0189491715067789,
          0.0364750894089436,
          0.0364750894089436,
          0.0364750894089436,
          0.0071293260197190,
          0.0071293260197190,
          0.0071293260197190,
          0.0165591599520032,
          0.0165591599520032,
          0.0165591599520032,
          0.0137596162349422,
          0.0137596162349422,
          0.0137596162349422,
          0.0137596162349422,
          0.0137596162349422,
          0.0137596162349422,
          0.0237819109001528,
          0.0237819109001528,
          0.0237819109001528,
          0.0237819109001528,
          0.0237819109001528,
          0.0237819109001528,
          0.0045305345022571,
          0.0045305345022571,
          0.0045305345022571,
          0.0045305345022571,
          0.0045305345022571,
          0.0045305345022571,
          0.0068401101196072,
          0.0068401101196072,
          0.0068401101196072,
          0.0068401101196072,
          0.0068401101196072,
          0.0068401101196072,
          0.0177474891020204,
          0.0177474891020204,
          0.0177474891020204,
          0.0177474891020204,
          0.0177474891020204,
          0.0177474891020204,
          0.0050106608745797,
          0.0050106608745797,
          0.0050106608745797,
          0.0050106608745797,
          0.0050106608745797,
          0.0050106608745797,
          0.0012229481269611,
          0.0012229481269611,
          0.0012229481269611,
          0.0012229481269611,
          0.0012229481269611,
          0.0012229481269611,
          0.0254821753118244,
          0.0254821753118244,
          0.0254821753118244,
          0.0254821753118244,
          0.0254821753118244,
          0.0254821753118244 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule19 ( int n, double a[], double b[], double c[], double w[] )

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

   rule19() returns the rule of precision 19.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.3333333333333334,
          0.0525238903512090,
          0.8949522192975821,
          0.0525238903512090,
          0.4925126750413369,
          0.0149746499173263,
          0.4925126750413369,
          0.1114488733230214,
          0.7771022533539573,
          0.1114488733230214,
          0.4591942010395437,
          0.0816115979209127,
          0.4591942010395437,
          0.4039697225519012,
          0.1920605548961976,
          0.4039697225519012,
          0.1781701047817643,
          0.6436597904364714,
          0.1781701047817643,
          0.0116394611837894,
          0.9767210776324211,
          0.0116394611837894,
          0.2551616329136077,
          0.4896767341727846,
          0.2551616329136077,
          0.1306976762680324,
          0.8301564644002754,
          0.0391458593316922,
          0.8301564644002754,
          0.0391458593316922,
          0.1306976762680324,
          0.3113176298095413,
          0.5593698057203009,
          0.1293125644701578,
          0.5593698057203009,
          0.1293125644701578,
          0.3113176298095413,
          0.0020689258966048,
          0.6333132931287841,
          0.3646177809746111,
          0.6333132931287841,
          0.3646177809746111,
          0.0020689258966048,
          0.0745602946016267,
          0.7040048199660421,
          0.2214348854323312,
          0.7040048199660421,
          0.2214348854323312,
          0.0745602946016267,
          0.0050072882573545,
          0.8525669543768892,
          0.1424257573657563,
          0.8525669543768892,
          0.1424257573657563,
          0.0050072882573545,
          0.0408880111960169,
          0.6050839790687079,
          0.3540280097352752,
          0.6050839790687079,
          0.3540280097352752,
          0.0408880111960169,
          0.2418945789605796,
          0.7431813689574364,
          0.0149240520819841,
          0.7431813689574364,
          0.0149240520819841,
          0.2418945789605796,
          0.0600862753223067,
          0.9301376988768051,
          0.0097760258008882,
          0.9301376988768051,
          0.0097760258008882,
          0.0600862753223067 };

  static double b_save[] = {
          0.3333333333333334,
          0.8949522192975821,
          0.0525238903512090,
          0.0525238903512090,
          0.0149746499173263,
          0.4925126750413369,
          0.4925126750413369,
          0.7771022533539573,
          0.1114488733230214,
          0.1114488733230214,
          0.0816115979209127,
          0.4591942010395437,
          0.4591942010395437,
          0.1920605548961976,
          0.4039697225519012,
          0.4039697225519012,
          0.6436597904364714,
          0.1781701047817643,
          0.1781701047817643,
          0.9767210776324211,
          0.0116394611837894,
          0.0116394611837894,
          0.4896767341727846,
          0.2551616329136077,
          0.2551616329136077,
          0.8301564644002754,
          0.1306976762680324,
          0.8301564644002754,
          0.0391458593316922,
          0.1306976762680324,
          0.0391458593316922,
          0.5593698057203009,
          0.3113176298095413,
          0.5593698057203009,
          0.1293125644701578,
          0.3113176298095413,
          0.1293125644701578,
          0.6333132931287841,
          0.0020689258966048,
          0.6333132931287841,
          0.3646177809746111,
          0.0020689258966048,
          0.3646177809746111,
          0.7040048199660421,
          0.0745602946016267,
          0.7040048199660421,
          0.2214348854323312,
          0.0745602946016267,
          0.2214348854323312,
          0.8525669543768892,
          0.0050072882573545,
          0.8525669543768892,
          0.1424257573657563,
          0.0050072882573545,
          0.1424257573657563,
          0.6050839790687079,
          0.0408880111960169,
          0.6050839790687079,
          0.3540280097352752,
          0.0408880111960169,
          0.3540280097352752,
          0.7431813689574364,
          0.2418945789605796,
          0.7431813689574364,
          0.0149240520819841,
          0.2418945789605796,
          0.0149240520819841,
          0.9301376988768051,
          0.0600862753223067,
          0.9301376988768051,
          0.0097760258008882,
          0.0600862753223067,
          0.0097760258008882 };

  static double c_save[] = {
          0.3333333333333333,
          0.0525238903512090,
          0.0525238903512090,
          0.8949522192975821,
          0.4925126750413369,
          0.4925126750413369,
          0.0149746499173262,
          0.1114488733230213,
          0.1114488733230212,
          0.7771022533539572,
          0.4591942010395437,
          0.4591942010395437,
          0.0816115979209127,
          0.4039697225519012,
          0.4039697225519012,
          0.1920605548961976,
          0.1781701047817643,
          0.1781701047817643,
          0.6436597904364714,
          0.0116394611837894,
          0.0116394611837894,
          0.9767210776324211,
          0.2551616329136077,
          0.2551616329136077,
          0.4896767341727846,
          0.0391458593316922,
          0.0391458593316922,
          0.1306976762680324,
          0.1306976762680324,
          0.8301564644002754,
          0.8301564644002754,
          0.1293125644701578,
          0.1293125644701578,
          0.3113176298095413,
          0.3113176298095413,
          0.5593698057203009,
          0.5593698057203009,
          0.3646177809746111,
          0.3646177809746111,
          0.0020689258966048,
          0.0020689258966048,
          0.6333132931287841,
          0.6333132931287841,
          0.2214348854323313,
          0.2214348854323311,
          0.0745602946016267,
          0.0745602946016266,
          0.7040048199660421,
          0.7040048199660421,
          0.1424257573657564,
          0.1424257573657564,
          0.0050072882573545,
          0.0050072882573544,
          0.8525669543768892,
          0.8525669543768892,
          0.3540280097352752,
          0.3540280097352753,
          0.0408880111960168,
          0.0408880111960169,
          0.6050839790687079,
          0.6050839790687079,
          0.0149240520819840,
          0.0149240520819840,
          0.2418945789605795,
          0.2418945789605795,
          0.7431813689574364,
          0.7431813689574364,
          0.0097760258008881,
          0.0097760258008882,
          0.0600862753223067,
          0.0600862753223068,
          0.9301376988768051,
          0.9301376988768051 };

  double w_save[] = {
          0.0344693977040123,
          0.0071092565977981,
          0.0071092565977981,
          0.0071092565977981,
          0.0103217551429443,
          0.0103217551429443,
          0.0103217551429443,
          0.0152343510930183,
          0.0152343510930183,
          0.0152343510930183,
          0.0229835900267416,
          0.0229835900267416,
          0.0229835900267416,
          0.0315375348931550,
          0.0315375348931550,
          0.0315375348931550,
          0.0246519148481909,
          0.0246519148481909,
          0.0246519148481909,
          0.0017653227764428,
          0.0017653227764428,
          0.0017653227764428,
          0.0317530193660031,
          0.0317530193660031,
          0.0317530193660031,
          0.0096954844868550,
          0.0096954844868550,
          0.0096954844868550,
          0.0096954844868550,
          0.0096954844868550,
          0.0096954844868550,
          0.0263463219773907,
          0.0263463219773907,
          0.0263463219773907,
          0.0263463219773907,
          0.0263463219773907,
          0.0263463219773907,
          0.0032820765518358,
          0.0032820765518358,
          0.0032820765518358,
          0.0032820765518358,
          0.0032820765518358,
          0.0032820765518358,
          0.0181079449312125,
          0.0181079449312125,
          0.0181079449312125,
          0.0181079449312125,
          0.0181079449312125,
          0.0181079449312125,
          0.0029263151034702,
          0.0029263151034702,
          0.0029263151034702,
          0.0029263151034702,
          0.0029263151034702,
          0.0029263151034702,
          0.0161021627640241,
          0.0161021627640241,
          0.0161021627640241,
          0.0161021627640241,
          0.0161021627640241,
          0.0161021627640241,
          0.0084558874995365,
          0.0084558874995365,
          0.0084558874995365,
          0.0084558874995365,
          0.0084558874995365,
          0.0084558874995365,
          0.0033272013628594,
          0.0033272013628594,
          0.0033272013628594,
          0.0033272013628594,
          0.0033272013628594,
          0.0033272013628594 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule20 ( int n, double a[], double b[], double c[], double w[] )

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

   rule20() returns the rule of precision 20.

   The data is given for the following triangle:

    (0,1) 
      |   \ 
      |    \ 
      |     \ 
      |      \ 
    (0,0)----(1,0)

   We suppose we are given a triangle T with vertices A, B, C.
   We call a rule with n points, returning barycentric coordinates
   a, b, c, and weights w.  Then the integral I of f(x,y) over T is 
   approximated by Q as follows:

   (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C
   Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 April 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    int n: the number of quadrature points for this rule.

  Output:

    double a[n], b[n], c[n]: the barycentric coordinates of quadrature points.

    double w[n]: the quadrature weights, which sum to 1.
*/
{
  static double a_save[] = {
          0.3333333333333334,
          0.2545792676733391,
          0.4908414646533218,
          0.2545792676733391,
          0.0109761410283978,
          0.9780477179432046,
          0.0109761410283978,
          0.1093835967117146,
          0.7812328065765708,
          0.1093835967117146,
          0.1862949977445409,
          0.6274100045109181,
          0.1862949977445409,
          0.4455510569559248,
          0.1088978860881504,
          0.4455510569559248,
          0.0373108805988847,
          0.9253782388022307,
          0.0373108805988847,
          0.3934253478170999,
          0.2131493043658003,
          0.3934253478170999,
          0.4762456115404990,
          0.0475087769190020,
          0.4762456115404990,
          0.0075707805046965,
          0.8332955118382362,
          0.1591337076570672,
          0.8332955118382362,
          0.1591337076570672,
          0.0075707805046965,
          0.0465603649076643,
          0.7549215028635474,
          0.1985181322287882,
          0.7549215028635474,
          0.1985181322287882,
          0.0465603649076643,
          0.0640905856084341,
          0.9310544767839422,
          0.0048549376076237,
          0.9310544767839422,
          0.0048549376076237,
          0.0640905856084341,
          0.0549874791429868,
          0.6118777035474257,
          0.3331348173095875,
          0.6118777035474257,
          0.3331348173095875,
          0.0549874791429868,
          0.0999522962881387,
          0.8616840189364867,
          0.0383636847753746,
          0.8616840189364867,
          0.0383636847753746,
          0.0999522962881387,
          0.1062272047202700,
          0.6781657378896355,
          0.2156070573900944,
          0.6781657378896355,
          0.2156070573900944,
          0.1062272047202700,
          0.4200237588162241,
          0.5701446928909734,
          0.0098315482928026,
          0.5701446928909734,
          0.0098315482928026,
          0.4200237588162241,
          0.3178601238357720,
          0.5423318041724281,
          0.1398080719917999,
          0.5423318041724281,
          0.1398080719917999,
          0.3178601238357720,
          0.0107372128560111,
          0.7086813757203236,
          0.2805814114236652,
          0.7086813757203236,
          0.2805814114236652,
          0.0107372128560111 };

  static double b_save[] = {
          0.3333333333333334,
          0.4908414646533218,
          0.2545792676733391,
          0.2545792676733391,
          0.9780477179432046,
          0.0109761410283978,
          0.0109761410283978,
          0.7812328065765708,
          0.1093835967117146,
          0.1093835967117146,
          0.6274100045109181,
          0.1862949977445409,
          0.1862949977445409,
          0.1088978860881504,
          0.4455510569559248,
          0.4455510569559248,
          0.9253782388022307,
          0.0373108805988847,
          0.0373108805988847,
          0.2131493043658003,
          0.3934253478170999,
          0.3934253478170999,
          0.0475087769190020,
          0.4762456115404990,
          0.4762456115404990,
          0.8332955118382362,
          0.0075707805046965,
          0.8332955118382362,
          0.1591337076570672,
          0.0075707805046965,
          0.1591337076570672,
          0.7549215028635474,
          0.0465603649076643,
          0.7549215028635474,
          0.1985181322287882,
          0.0465603649076643,
          0.1985181322287882,
          0.9310544767839422,
          0.0640905856084341,
          0.9310544767839422,
          0.0048549376076237,
          0.0640905856084341,
          0.0048549376076237,
          0.6118777035474257,
          0.0549874791429868,
          0.6118777035474257,
          0.3331348173095875,
          0.0549874791429868,
          0.3331348173095875,
          0.8616840189364867,
          0.0999522962881387,
          0.8616840189364867,
          0.0383636847753746,
          0.0999522962881387,
          0.0383636847753746,
          0.6781657378896355,
          0.1062272047202700,
          0.6781657378896355,
          0.2156070573900944,
          0.1062272047202700,
          0.2156070573900944,
          0.5701446928909734,
          0.4200237588162241,
          0.5701446928909734,
          0.0098315482928026,
          0.4200237588162241,
          0.0098315482928026,
          0.5423318041724281,
          0.3178601238357720,
          0.5423318041724281,
          0.1398080719917999,
          0.3178601238357720,
          0.1398080719917999,
          0.7086813757203236,
          0.0107372128560111,
          0.7086813757203236,
          0.2805814114236652,
          0.0107372128560111,
          0.2805814114236652 };

  static double c_save[] = {
          0.3333333333333333,
          0.2545792676733392,
          0.2545792676733392,
          0.4908414646533218,
          0.0109761410283977,
          0.0109761410283977,
          0.9780477179432044,
          0.1093835967117146,
          0.1093835967117146,
          0.7812328065765708,
          0.1862949977445409,
          0.1862949977445409,
          0.6274100045109181,
          0.4455510569559249,
          0.4455510569559248,
          0.1088978860881504,
          0.0373108805988847,
          0.0373108805988847,
          0.9253782388022306,
          0.3934253478170999,
          0.3934253478170999,
          0.2131493043658003,
          0.4762456115404990,
          0.4762456115404990,
          0.0475087769190020,
          0.1591337076570672,
          0.1591337076570672,
          0.0075707805046965,
          0.0075707805046965,
          0.8332955118382362,
          0.8332955118382362,
          0.1985181322287882,
          0.1985181322287883,
          0.0465603649076644,
          0.0465603649076645,
          0.7549215028635475,
          0.7549215028635475,
          0.0048549376076237,
          0.0048549376076238,
          0.0640905856084341,
          0.0640905856084342,
          0.9310544767839422,
          0.9310544767839422,
          0.3331348173095875,
          0.3331348173095876,
          0.0549874791429867,
          0.0549874791429869,
          0.6118777035474257,
          0.6118777035474257,
          0.0383636847753746,
          0.0383636847753746,
          0.0999522962881386,
          0.0999522962881387,
          0.8616840189364867,
          0.8616840189364867,
          0.2156070573900944,
          0.2156070573900944,
          0.1062272047202700,
          0.1062272047202701,
          0.6781657378896355,
          0.6781657378896355,
          0.0098315482928025,
          0.0098315482928025,
          0.4200237588162241,
          0.4200237588162241,
          0.5701446928909734,
          0.5701446928909732,
          0.1398080719917999,
          0.1398080719917999,
          0.3178601238357720,
          0.3178601238357720,
          0.5423318041724281,
          0.5423318041724281,
          0.2805814114236652,
          0.2805814114236653,
          0.0107372128560111,
          0.0107372128560111,
          0.7086813757203236,
          0.7086813757203236 };

  double w_save[] = {
          0.0278202214029062,
          0.0281664026150405,
          0.0281664026150405,
          0.0281664026150405,
          0.0015976815821332,
          0.0015976815821332,
          0.0015976815821332,
          0.0156604615521491,
          0.0156604615521491,
          0.0156604615521491,
          0.0183469259485058,
          0.0183469259485058,
          0.0183469259485058,
          0.0189047998664649,
          0.0189047998664649,
          0.0189047998664649,
          0.0043225508213312,
          0.0043225508213312,
          0.0043225508213312,
          0.0275761012581409,
          0.0275761012581409,
          0.0275761012581409,
          0.0142036506068169,
          0.0142036506068169,
          0.0142036506068169,
          0.0044057948371170,
          0.0044057948371170,
          0.0044057948371170,
          0.0044057948371170,
          0.0044057948371170,
          0.0044057948371170,
          0.0119727971579094,
          0.0119727971579094,
          0.0119727971579094,
          0.0119727971579094,
          0.0119727971579094,
          0.0119727971579094,
          0.0022597392042517,
          0.0022597392042517,
          0.0022597392042517,
          0.0022597392042517,
          0.0022597392042517,
          0.0022597392042517,
          0.0173344511344387,
          0.0173344511344387,
          0.0173344511344387,
          0.0173344511344387,
          0.0173344511344387,
          0.0173344511344387,
          0.0082914230552277,
          0.0082914230552277,
          0.0082914230552277,
          0.0082914230552277,
          0.0082914230552277,
          0.0082914230552277,
          0.0154452156441985,
          0.0154452156441985,
          0.0154452156441985,
          0.0154452156441985,
          0.0154452156441985,
          0.0154452156441985,
          0.0073913630005106,
          0.0073913630005106,
          0.0073913630005106,
          0.0073913630005106,
          0.0073913630005106,
          0.0073913630005106,
          0.0233834914636555,
          0.0233834914636555,
          0.0233834914636555,
          0.0233834914636555,
          0.0233834914636555,
          0.0233834914636555,
          0.0071564004769154,
          0.0071564004769154,
          0.0071564004769154,
          0.0071564004769154,
          0.0071564004769154,
          0.0071564004769154 };

  r8vec_copy ( n, a_save, a );
  r8vec_copy ( n, b_save, b );
  r8vec_copy ( n, c_save, c );
  r8vec_copy ( n, w_save, w );

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

double triangle_area ( double vert1[], double vert2[], double vert3[] )

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

    triangle_area() returns the area of a triangle.

  Licensing:

    This code is distributed under the GNU GPL license.

  Modified:

    26 June 2014

  Author:

    John Burkardt.

  Input:

    double vert1[2], vert2[2], vert3[2]: the vertices of the triangle.

  Output:

    double triangle_area: the area of the triangle.
*/
{
  double value;

  value = 0.5 * 
    ( 
        ( vert2[0] - vert1[0] ) * ( vert3[1] - vert1[1] ) 
      - ( vert3[0] - vert1[0] ) * ( vert2[1] - vert1[1] ) 
    );

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

double triangle_unit_area ( )

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

    triangle_unit_area() returns the area of a unit triangle.

  Licensing:

    This code is distributed under the GNU GPL license.

  Modified:

    24 May 2023

  Author:

    John Burkardt.

  Output:

    double triangle_unit_area: the area of the triangle.
*/
{
  double value;

  value = 0.5;

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

double triangle_unit_monomial_integral ( int expon[] )

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

    triangle_unit_monomial_integral() integrates a monomial over the unit triangle.

  Discussion:

    This routine evaluates a monomial of the form

      product ( 1 <= dim <= dim_num ) x(dim)^expon(dim)

    where the exponents are nonnegative integers.  Note that
    if the combination 0^0 is encountered, it should be treated
    as 1.

    Integral ( over unit triangle ) x^m y^n dx dy = m! * n! / ( m + n + 2 )!

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    04 July 2007

  Author:

    John Burkardt

  Input:

    int EXPON[3], the exponents.

  Output:

    double triangle_unit_monomial_integral, the value of the integral.
*/
{
  int i;
  int k;
  double value;
/*
  The first computation ends with VALUE = 1.0;
*/
  value = 1.0;

  k = 0;

  for ( i = 1; i <= expon[0]; i++ )
  {
    k = k + 1;
/*  value = value * ( double ) ( i ) / ( double ) ( k );*/
  }

  for ( i = 1; i <= expon[1]; i++ )
  {
    k = k + 1;
    value = value * ( double ) ( i ) / ( double ) ( k );
  }

  k = k + 1;
  value = value / ( double ) ( k );

  k = k + 1;
  value = value / ( double ) ( k );

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

void triangle_witherden_rule ( int p, int n, double a[], double b[], 
  double c[], double w[] )

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

    triangle_witherden_rule() returns a triangle quadrature rule of given precision.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    24 May 2023

  Author:

    John Burkardt

  Reference:

    Freddie Witherden, Peter Vincent,
    On the identification of symmetric quadrature rules for finite element methods,
    Computers and Mathematics with Applications,
    Volume 69, pages 1232-1241, 2015.

  Input:

    integer p: the precision, 0 <= p <= 20.

    integer n: the order of the rule.

  Output:

    real a(n), b(n), c(n): the barycentric coordinates of quadrature points.

    real w(n): the weights of quadrature points.
*/
{
  if ( p < 0 )
  {
    printf ( "\n" );
    printf ( "triangle_witherden_rule(): Fatal error!\n" );
    printf ( "  Input p < 0.\n" );
    exit ( 1 );
  }

  if ( 20 < p )
  {
    printf ( "\n" );
    printf ( "triangle_witherden_rule(): Fatal error!\n" );
    printf ( "  Input 20 < p.\n" );
    exit ( 1 );
  }

  if ( p == 0 )
  {
    rule00 ( n, a, b, c, w );
  }
  else if ( p == 1 )
  {
    rule01 ( n, a, b, c, w );
  }
  else if ( p == 2 )
  {
    rule02 ( n, a, b, c, w );
  }
  else if ( p == 3 )
  {
    rule03 ( n, a, b, c, w );
  }
  else if ( p == 4 )
  {
    rule04 ( n, a, b, c, w );
  }
  else if ( p == 5 )
  {
    rule05 ( n, a, b, c, w );
  }
  else if ( p == 6 )
  {
    rule06 ( n, a, b, c, w );
  }
  else if ( p == 7 )
  {
    rule07 ( n, a, b, c, w );
  }
  else if ( p == 8 )
  {
    rule08 ( n, a, b, c, w );
  }
  else if ( p == 9 )
  {
    rule09 ( n, a, b, c, w );
  }
  else if ( p == 10 )
  {
    rule10 ( n, a, b, c, w );
  }
  else if ( p == 11 )
  {
    rule11 ( n, a, b, c, w );
  }
  else if ( p == 12 )
  {
    rule12 ( n, a, b, c, w );
  }
  else if ( p == 13 )
  {
    rule13 ( n, a, b, c, w );
  }
  else if ( p == 14 )
  {
    rule14 ( n, a, b, c, w );
  }
  else if ( p == 15 )
  {
    rule15 ( n, a, b, c, w );
  }
  else if ( p == 16 )
  {
    rule16 ( n, a, b, c, w );
  }
  else if ( p == 17 )
  {
    rule17 ( n, a, b, c, w );
  }
  else if ( p == 18 )
  {
    rule18 ( n, a, b, c, w );
  }
  else if ( p == 19 )
  {
    rule19 ( n, a, b, c, w );
  }
  else if ( p == 20 )
  {
    rule20 ( n, a, b, c, w );
  }

  return;
}