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

using namespace std;

# include "hexahedron_witherden_rule.hpp"

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

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

//****************************************************************************80
//
//  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 SAVE statement to maintain the
//    variables H and T.  I have decided 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, bool &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;
}
//****************************************************************************80

double hexahedron_unit_monomial_integral ( int expon[3] )

//****************************************************************************80
//
//  Purpose:
//
//    hexahedron_unit_monomial_integral(): monomial integral in a unit hexahedron.
//
//  Discussion:
//
//    This function returns the value of the integral of X^ALPHA Y^BETA Z^GAMMA
//    over the unit hexahedron.
//
//    The unit hexahedron has vertices 
//    (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 May 2023
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int EXPON[3], the exponents.
//
//  Output:
//
//    double hexahedron_unit_monomial_integral, the integral of the monomial
//    over the hexahedron.
//
{
  double value;

  value = 1.0 / ( expon[0] + 1 ) / ( expon[1] + 1 ) / ( expon[2] + 1 );

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

double hexahedron_unit_volume ( )

//****************************************************************************80
//
//  Purpose:
//
//    hexahedron_unit_volume() returns the volume of a unit hexahedron.
//
//  Discussion:
//
//    The unit hexahedron has vertices 
//    (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1).
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    02 May 2023
//
//  Author:
//
//    John Burkardt
//
//  Output:
//
//    real value: the volume.
//
{
  double value;

  value = 1.0;

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

void hexahedron_witherden_rule ( int p, int n, double x[], double y[], 
  double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    hexahedron_witherden_rule() returns a hexahedron quadrature rule of given precision.
//
//  Discussion:
//
//    The unit hexahedron has vertices 
//    (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    22 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:
//
//    int p: the precision, 0 <= p <= 11.
//
//    int n: the order of the rule.
//
//  Output:
//
//    double x(n), y(n), z(n): the coordinates of quadrature points.
//
//    double w(n): the quadrature weights.
//
{
  if ( p < 0 ) 
  {
    cout << "\n";
    cout << "hexahedron_witherden_rule(): Fatal error!\n";
    cout << "  Input p < 0.\n";
    exit ( 1 );
  }

  if ( 11 < p )
  {
    cout << "\n";
    cout << "hexahedron_witherden_rule(): Fatal error!\n";
    cout << "  Input 11 < p.\n";
    exit ( 1 );
  }

  if ( p <= 1 )
  {
    rule01 ( n, x, y, z, w );
  }
  else if ( p <= 3 )
  {
    rule03 ( n, x, y, z, w );
  }
  else if ( p <= 5 )
  {
    rule05 ( n, x, y, z, w );
  }
  else if ( p <= 7 )
  {
    rule07 ( n, x, y, z, w );
  }
  else if ( p <= 9 )
  {
    rule09 ( n, x, y, z, w );
  }
  else if ( p <= 11 )
  {
    rule11 ( n, x, y, z, w );
  }

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

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

//****************************************************************************80
//
//  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 = new double[n];
  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;
}
//****************************************************************************80

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

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

int rule_order ( int p )

//****************************************************************************80
//
//  Purpose:
//
//    rule_order() returns the order of a hexahedron 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:
//
//    02 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:
//
//    int p: the precision, 0 <= p <= 11.
//
//  Output:
//
//    int rule_order: the order of the rule.
//
{
  int order;
  static int order_save[] = {
       1,
       1,  6,  6, 14, 14, 34, 34, 58, 58, 90,
      90 };

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

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

  order = order_save[p];

  return order;
}

void rule01 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule01() returns the hexahedron rule of precision 1.
//
//  Discussion:
//
//    We suppose we are given a hexahedron H with vertices
//      (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1),
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    H is approximated by Q as follows:
//
//    Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 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:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), x(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.5000000000000000 };

  static double y_save[] = {
          0.5000000000000000 };

  static double z_save[] = {
          0.5000000000000000 };

  double w_save[] = {
          1.0000000000000000 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule03 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule03() returns the hexahedron rule of precision 3.
//
//  Discussion:
//
//    We suppose we are given a hexahedron H with vertices
//      (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1),
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    H is approximated by Q as follows:
//
//    Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 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:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), x(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.0000000000000000,
          0.5000000000000000,
          0.5000000000000000,
          0.5000000000000000,
          1.0000000000000000,
          0.5000000000000000 };

  static double y_save[] = {
          0.5000000000000000,
          0.5000000000000000,
          1.0000000000000000,
          0.5000000000000000,
          0.5000000000000000,
          0.0000000000000000 };

  static double z_save[] = {
          0.5000000000000000,
          1.0000000000000000,
          0.5000000000000000,
          0.0000000000000000,
          0.5000000000000000,
          0.5000000000000000 };

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

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule05 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule05() returns the hexahedron rule of precision 5.
//
//  Discussion:
//
//    We suppose we are given a hexahedron H with vertices
//      (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1),
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    H is approximated by Q as follows:
//
//    Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 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:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), x(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.1020887871228893,
          0.5000000000000000,
          0.5000000000000000,
          0.5000000000000000,
          0.8979112128771107,
          0.5000000000000000,
          0.8793934553196641,
          0.1206065446803359,
          0.1206065446803359,
          0.1206065446803359,
          0.1206065446803359,
          0.8793934553196641,
          0.8793934553196641,
          0.8793934553196641 };

  static double y_save[] = {
          0.5000000000000000,
          0.5000000000000000,
          0.8979112128771107,
          0.5000000000000000,
          0.5000000000000000,
          0.1020887871228893,
          0.1206065446803359,
          0.8793934553196641,
          0.8793934553196641,
          0.1206065446803359,
          0.1206065446803359,
          0.8793934553196641,
          0.8793934553196641,
          0.1206065446803359 };

  static double z_save[] = {
          0.5000000000000000,
          0.8979112128771107,
          0.5000000000000000,
          0.1020887871228893,
          0.5000000000000000,
          0.5000000000000000,
          0.1206065446803359,
          0.8793934553196641,
          0.1206065446803359,
          0.1206065446803359,
          0.8793934553196641,
          0.1206065446803359,
          0.8793934553196641,
          0.8793934553196641 };

  double w_save[] = {
          0.1108033240997230,
          0.1108033240997230,
          0.1108033240997230,
          0.1108033240997230,
          0.1108033240997230,
          0.1108033240997230,
          0.0418975069252078,
          0.0418975069252078,
          0.0418975069252078,
          0.0418975069252078,
          0.0418975069252078,
          0.0418975069252078,
          0.0418975069252078,
          0.0418975069252078 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule07 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule07() returns the hexahedron rule of precision 7.
//
//  Discussion:
//
//    We suppose we are given a hexahedron H with vertices
//      (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1),
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    H is approximated by Q as follows:
//
//    Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 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:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), x(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.0059569194056617,
          0.5000000000000000,
          0.5000000000000000,
          0.5000000000000000,
          0.9940430805943383,
          0.5000000000000000,
          0.7039775836791597,
          0.2960224163208403,
          0.2960224163208403,
          0.2960224163208403,
          0.2960224163208403,
          0.7039775836791597,
          0.7039775836791597,
          0.7039775836791597,
          0.8905514105020593,
          0.1094485894979407,
          0.1094485894979407,
          0.1094485894979407,
          0.1094485894979407,
          0.8905514105020593,
          0.8905514105020593,
          0.8905514105020593,
          0.0759738621579806,
          0.9240261378420194,
          0.5000000000000000,
          0.9240261378420194,
          0.9240261378420194,
          0.5000000000000000,
          0.5000000000000000,
          0.0759738621579806,
          0.0759738621579806,
          0.9240261378420194,
          0.5000000000000000,
          0.0759738621579806 };

  static double y_save[] = {
          0.5000000000000000,
          0.5000000000000000,
          0.9940430805943383,
          0.5000000000000000,
          0.5000000000000000,
          0.0059569194056617,
          0.2960224163208403,
          0.7039775836791597,
          0.7039775836791597,
          0.2960224163208403,
          0.2960224163208403,
          0.7039775836791597,
          0.7039775836791597,
          0.2960224163208403,
          0.1094485894979407,
          0.8905514105020593,
          0.8905514105020593,
          0.1094485894979407,
          0.1094485894979407,
          0.8905514105020593,
          0.8905514105020593,
          0.1094485894979407,
          0.0759738621579806,
          0.5000000000000000,
          0.9240261378420194,
          0.9240261378420194,
          0.5000000000000000,
          0.0759738621579806,
          0.0759738621579806,
          0.5000000000000000,
          0.9240261378420194,
          0.0759738621579806,
          0.9240261378420194,
          0.5000000000000000 };

  static double z_save[] = {
          0.5000000000000000,
          0.9940430805943383,
          0.5000000000000000,
          0.0059569194056617,
          0.5000000000000000,
          0.5000000000000000,
          0.2960224163208403,
          0.7039775836791597,
          0.2960224163208403,
          0.2960224163208403,
          0.7039775836791597,
          0.2960224163208403,
          0.7039775836791597,
          0.7039775836791597,
          0.1094485894979407,
          0.8905514105020593,
          0.1094485894979407,
          0.1094485894979407,
          0.8905514105020593,
          0.1094485894979407,
          0.8905514105020593,
          0.8905514105020593,
          0.5000000000000000,
          0.0759738621579806,
          0.0759738621579806,
          0.5000000000000000,
          0.9240261378420194,
          0.9240261378420194,
          0.0759738621579806,
          0.9240261378420194,
          0.5000000000000000,
          0.5000000000000000,
          0.9240261378420194,
          0.0759738621579806 };

  double w_save[] = {
          0.0250162363729459,
          0.0250162363729459,
          0.0250162363729459,
          0.0250162363729459,
          0.0250162363729459,
          0.0250162363729459,
          0.0571442320087316,
          0.0571442320087316,
          0.0571442320087316,
          0.0571442320087316,
          0.0571442320087316,
          0.0571442320087316,
          0.0571442320087316,
          0.0571442320087316,
          0.0192245175082448,
          0.0192245175082448,
          0.0192245175082448,
          0.0192245175082448,
          0.0192245175082448,
          0.0192245175082448,
          0.0192245175082448,
          0.0192245175082448,
          0.0199127154688761,
          0.0199127154688761,
          0.0199127154688761,
          0.0199127154688761,
          0.0199127154688761,
          0.0199127154688761,
          0.0199127154688761,
          0.0199127154688761,
          0.0199127154688761,
          0.0199127154688761,
          0.0199127154688761,
          0.0199127154688761 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule09 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule09() returns the hexahedron rule of precision 9.
//
//  Discussion:
//
//    We suppose we are given a hexahedron H with vertices
//      (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1),
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    H is approximated by Q as follows:
//
//    Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 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:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), x(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.1931592652041455,
          0.5000000000000000,
          0.5000000000000000,
          0.5000000000000000,
          0.8068407347958545,
          0.5000000000000000,
          0.9350498923309880,
          0.0649501076690120,
          0.0649501076690120,
          0.0649501076690120,
          0.0649501076690120,
          0.9350498923309880,
          0.9350498923309880,
          0.9350498923309880,
          0.7820554035100150,
          0.2179445964899850,
          0.2179445964899850,
          0.2179445964899850,
          0.2179445964899850,
          0.7820554035100150,
          0.7820554035100150,
          0.7820554035100150,
          0.0611564383711609,
          0.9388435616288391,
          0.5000000000000000,
          0.9388435616288391,
          0.9388435616288391,
          0.5000000000000000,
          0.5000000000000000,
          0.0611564383711609,
          0.0611564383711609,
          0.9388435616288391,
          0.5000000000000000,
          0.0611564383711609,
          0.7161339513154311,
          0.2838660486845689,
          0.2838660486845689,
          0.0307347890676641,
          0.2838660486845689,
          0.9692652109323359,
          0.7161339513154311,
          0.7161339513154311,
          0.2838660486845689,
          0.7161339513154311,
          0.0307347890676641,
          0.2838660486845689,
          0.2838660486845689,
          0.2838660486845689,
          0.7161339513154311,
          0.9692652109323359,
          0.7161339513154311,
          0.9692652109323359,
          0.9692652109323359,
          0.0307347890676641,
          0.7161339513154311,
          0.2838660486845689,
          0.7161339513154311,
          0.0307347890676641 };

  static double y_save[] = {
          0.5000000000000000,
          0.5000000000000000,
          0.8068407347958545,
          0.5000000000000000,
          0.5000000000000000,
          0.1931592652041455,
          0.0649501076690120,
          0.9350498923309880,
          0.9350498923309880,
          0.0649501076690120,
          0.0649501076690120,
          0.9350498923309880,
          0.9350498923309880,
          0.0649501076690120,
          0.2179445964899850,
          0.7820554035100150,
          0.7820554035100150,
          0.2179445964899850,
          0.2179445964899850,
          0.7820554035100150,
          0.7820554035100150,
          0.2179445964899850,
          0.0611564383711609,
          0.5000000000000000,
          0.9388435616288391,
          0.9388435616288391,
          0.5000000000000000,
          0.0611564383711609,
          0.0611564383711609,
          0.5000000000000000,
          0.9388435616288391,
          0.0611564383711609,
          0.9388435616288391,
          0.5000000000000000,
          0.9692652109323359,
          0.0307347890676641,
          0.9692652109323359,
          0.2838660486845689,
          0.2838660486845689,
          0.7161339513154311,
          0.0307347890676641,
          0.0307347890676641,
          0.0307347890676641,
          0.7161339513154311,
          0.7161339513154311,
          0.7161339513154311,
          0.7161339513154311,
          0.2838660486845689,
          0.2838660486845689,
          0.2838660486845689,
          0.9692652109323359,
          0.7161339513154311,
          0.2838660486845689,
          0.7161339513154311,
          0.7161339513154311,
          0.9692652109323359,
          0.2838660486845689,
          0.2838660486845689 };

  static double z_save[] = {
          0.5000000000000000,
          0.8068407347958545,
          0.5000000000000000,
          0.1931592652041455,
          0.5000000000000000,
          0.5000000000000000,
          0.0649501076690120,
          0.9350498923309880,
          0.0649501076690120,
          0.0649501076690120,
          0.9350498923309880,
          0.0649501076690120,
          0.9350498923309880,
          0.9350498923309880,
          0.2179445964899850,
          0.7820554035100150,
          0.2179445964899850,
          0.2179445964899850,
          0.7820554035100150,
          0.2179445964899850,
          0.7820554035100150,
          0.7820554035100150,
          0.5000000000000000,
          0.0611564383711609,
          0.0611564383711609,
          0.5000000000000000,
          0.9388435616288391,
          0.9388435616288391,
          0.0611564383711609,
          0.9388435616288391,
          0.5000000000000000,
          0.5000000000000000,
          0.9388435616288391,
          0.0611564383711609,
          0.2838660486845689,
          0.7161339513154311,
          0.2838660486845689,
          0.7161339513154311,
          0.9692652109323359,
          0.2838660486845689,
          0.2838660486845689,
          0.7161339513154311,
          0.2838660486845689,
          0.0307347890676641,
          0.2838660486845689,
          0.0307347890676641,
          0.9692652109323359,
          0.0307347890676641,
          0.9692652109323359,
          0.2838660486845689,
          0.7161339513154311,
          0.7161339513154311,
          0.7161339513154311,
          0.7161339513154311,
          0.9692652109323359,
          0.7161339513154311,
          0.0307347890676641,
          0.2838660486845689 };

  double w_save[] = {
          0.0541593744687068,
          0.0541593744687068,
          0.0541593744687068,
          0.0541593744687068,
          0.0541593744687068,
          0.0541593744687068,
          0.0062685994124186,
          0.0062685994124186,
          0.0062685994124186,
          0.0062685994124186,
          0.0062685994124186,
          0.0062685994124186,
          0.0062685994124186,
          0.0062685994124186,
          0.0248574797680029,
          0.0248574797680029,
          0.0248574797680029,
          0.0248574797680029,
          0.0248574797680029,
          0.0248574797680029,
          0.0248574797680029,
          0.0248574797680029,
          0.0114737257670222,
          0.0114737257670222,
          0.0114737257670222,
          0.0114737257670222,
          0.0114737257670222,
          0.0114737257670222,
          0.0114737257670222,
          0.0114737257670222,
          0.0114737257670222,
          0.0114737257670222,
          0.0114737257670222,
          0.0114737257670222,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717,
          0.0120146004391717 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}
void rule11 ( int n, double x[], double y[], double z[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    rule11() returns the hexahedron rule of precision 11.
//
//  Discussion:
//
//    We suppose we are given a hexahedron H with vertices
//      (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1),
//    We call a rule with n points, returning coordinates
//    x, y, z, and weights w.  Then the integral I of f(x,y,z) over
//    H is approximated by Q as follows:
//
//    Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 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:
//
//    int n: the number of quadrature points for this rule.
//
//  Output:
//
//    double x(n), y(n), x(n): the coordinates of quadrature points.
//
//    double w[n]: the quadrature weights, which sum to 1.
//
{
  static double x_save[] = {
          0.0936928329501868,
          0.5000000000000000,
          0.5000000000000000,
          0.5000000000000000,
          0.9063071670498133,
          0.5000000000000000,
          0.8008376320991313,
          0.1991623679008687,
          0.1991623679008687,
          0.1991623679008687,
          0.1991623679008687,
          0.8008376320991313,
          0.8008376320991313,
          0.8008376320991313,
          0.9277278805088800,
          0.0722721194911200,
          0.0722721194911200,
          0.0722721194911200,
          0.0722721194911200,
          0.9277278805088800,
          0.9277278805088800,
          0.9277278805088800,
          0.6566967022580273,
          0.3433032977419727,
          0.3433032977419727,
          0.3433032977419727,
          0.3433032977419727,
          0.6566967022580273,
          0.6566967022580273,
          0.6566967022580273,
          0.1326658565014960,
          0.8673341434985040,
          0.5000000000000000,
          0.8673341434985040,
          0.8673341434985040,
          0.5000000000000000,
          0.5000000000000000,
          0.1326658565014960,
          0.1326658565014960,
          0.8673341434985040,
          0.5000000000000000,
          0.1326658565014960,
          0.7253999675572547,
          0.2746000324427453,
          0.2746000324427453,
          0.0174501672436448,
          0.2746000324427453,
          0.9825498327563551,
          0.7253999675572547,
          0.7253999675572547,
          0.2746000324427453,
          0.7253999675572547,
          0.0174501672436448,
          0.2746000324427453,
          0.2746000324427453,
          0.2746000324427453,
          0.7253999675572547,
          0.9825498327563551,
          0.7253999675572547,
          0.9825498327563551,
          0.9825498327563551,
          0.0174501672436448,
          0.7253999675572547,
          0.2746000324427453,
          0.7253999675572547,
          0.0174501672436448,
          0.9706224286053016,
          0.0293775713946984,
          0.0293775713946984,
          0.3230485927016850,
          0.0293775713946984,
          0.6769514072983150,
          0.9706224286053016,
          0.9706224286053016,
          0.0293775713946984,
          0.9706224286053016,
          0.3230485927016850,
          0.0293775713946984,
          0.0293775713946984,
          0.0293775713946984,
          0.9706224286053016,
          0.6769514072983150,
          0.9706224286053016,
          0.6769514072983150,
          0.6769514072983150,
          0.3230485927016850,
          0.9706224286053016,
          0.0293775713946984,
          0.9706224286053016,
          0.3230485927016850 };

  static double y_save[] = {
          0.5000000000000000,
          0.5000000000000000,
          0.9063071670498133,
          0.5000000000000000,
          0.5000000000000000,
          0.0936928329501868,
          0.1991623679008687,
          0.8008376320991313,
          0.8008376320991313,
          0.1991623679008687,
          0.1991623679008687,
          0.8008376320991313,
          0.8008376320991313,
          0.1991623679008687,
          0.0722721194911200,
          0.9277278805088800,
          0.9277278805088800,
          0.0722721194911200,
          0.0722721194911200,
          0.9277278805088800,
          0.9277278805088800,
          0.0722721194911200,
          0.3433032977419727,
          0.6566967022580273,
          0.6566967022580273,
          0.3433032977419727,
          0.3433032977419727,
          0.6566967022580273,
          0.6566967022580273,
          0.3433032977419727,
          0.1326658565014960,
          0.5000000000000000,
          0.8673341434985040,
          0.8673341434985040,
          0.5000000000000000,
          0.1326658565014960,
          0.1326658565014960,
          0.5000000000000000,
          0.8673341434985040,
          0.1326658565014960,
          0.8673341434985040,
          0.5000000000000000,
          0.9825498327563551,
          0.0174501672436448,
          0.9825498327563551,
          0.2746000324427453,
          0.2746000324427453,
          0.7253999675572547,
          0.0174501672436448,
          0.0174501672436448,
          0.0174501672436448,
          0.7253999675572547,
          0.7253999675572547,
          0.7253999675572547,
          0.7253999675572547,
          0.2746000324427453,
          0.2746000324427453,
          0.2746000324427453,
          0.9825498327563551,
          0.7253999675572547,
          0.2746000324427453,
          0.7253999675572547,
          0.7253999675572547,
          0.9825498327563551,
          0.2746000324427453,
          0.2746000324427453,
          0.6769514072983150,
          0.3230485927016850,
          0.6769514072983150,
          0.0293775713946984,
          0.0293775713946984,
          0.9706224286053016,
          0.3230485927016850,
          0.3230485927016850,
          0.3230485927016850,
          0.9706224286053016,
          0.9706224286053016,
          0.9706224286053016,
          0.9706224286053016,
          0.0293775713946984,
          0.0293775713946984,
          0.0293775713946984,
          0.6769514072983150,
          0.9706224286053016,
          0.0293775713946984,
          0.9706224286053016,
          0.9706224286053016,
          0.6769514072983150,
          0.0293775713946984,
          0.0293775713946984 };

  static double z_save[] = {
          0.5000000000000000,
          0.9063071670498133,
          0.5000000000000000,
          0.0936928329501868,
          0.5000000000000000,
          0.5000000000000000,
          0.1991623679008687,
          0.8008376320991313,
          0.1991623679008687,
          0.1991623679008687,
          0.8008376320991313,
          0.1991623679008687,
          0.8008376320991313,
          0.8008376320991313,
          0.0722721194911200,
          0.9277278805088800,
          0.0722721194911200,
          0.0722721194911200,
          0.9277278805088800,
          0.0722721194911200,
          0.9277278805088800,
          0.9277278805088800,
          0.3433032977419727,
          0.6566967022580273,
          0.3433032977419727,
          0.3433032977419727,
          0.6566967022580273,
          0.3433032977419727,
          0.6566967022580273,
          0.6566967022580273,
          0.5000000000000000,
          0.1326658565014960,
          0.1326658565014960,
          0.5000000000000000,
          0.8673341434985040,
          0.8673341434985040,
          0.1326658565014960,
          0.8673341434985040,
          0.5000000000000000,
          0.5000000000000000,
          0.8673341434985040,
          0.1326658565014960,
          0.2746000324427453,
          0.7253999675572547,
          0.2746000324427453,
          0.7253999675572547,
          0.9825498327563551,
          0.2746000324427453,
          0.2746000324427453,
          0.7253999675572547,
          0.2746000324427453,
          0.0174501672436448,
          0.2746000324427453,
          0.0174501672436448,
          0.9825498327563551,
          0.0174501672436448,
          0.9825498327563551,
          0.2746000324427453,
          0.7253999675572547,
          0.7253999675572547,
          0.7253999675572547,
          0.7253999675572547,
          0.9825498327563551,
          0.7253999675572547,
          0.0174501672436448,
          0.2746000324427453,
          0.0293775713946984,
          0.9706224286053016,
          0.0293775713946984,
          0.9706224286053016,
          0.6769514072983150,
          0.0293775713946984,
          0.0293775713946984,
          0.9706224286053016,
          0.0293775713946984,
          0.3230485927016850,
          0.0293775713946984,
          0.3230485927016850,
          0.6769514072983150,
          0.3230485927016850,
          0.6769514072983150,
          0.0293775713946984,
          0.9706224286053016,
          0.9706224286053016,
          0.9706224286053016,
          0.9706224286053016,
          0.6769514072983150,
          0.9706224286053016,
          0.3230485927016850,
          0.0293775713946984 };

  double w_save[] = {
          0.0253096342016000,
          0.0253096342016000,
          0.0253096342016000,
          0.0253096342016000,
          0.0253096342016000,
          0.0253096342016000,
          0.0146922934945570,
          0.0146922934945570,
          0.0146922934945570,
          0.0146922934945570,
          0.0146922934945570,
          0.0146922934945570,
          0.0146922934945570,
          0.0146922934945570,
          0.0055804890098537,
          0.0055804890098537,
          0.0055804890098537,
          0.0055804890098537,
          0.0055804890098537,
          0.0055804890098537,
          0.0055804890098537,
          0.0055804890098537,
          0.0269990056568711,
          0.0269990056568711,
          0.0269990056568711,
          0.0269990056568711,
          0.0269990056568711,
          0.0269990056568711,
          0.0269990056568711,
          0.0269990056568711,
          0.0181499182325145,
          0.0181499182325145,
          0.0181499182325145,
          0.0181499182325145,
          0.0181499182325145,
          0.0181499182325145,
          0.0181499182325145,
          0.0181499182325145,
          0.0181499182325145,
          0.0181499182325145,
          0.0181499182325145,
          0.0181499182325145,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0076802492622294,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527,
          0.0028267870173527 };

  r8vec_copy ( n, x_save, x );
  r8vec_copy ( n, y_save, y );
  r8vec_copy ( n, z_save, z );
  r8vec_copy ( n, w_save, w );

  return;
}