<commit_before>#include <libmesh/quadrature.h>
#include <libmesh/string_to_enum.h>
#include <libmesh/utility.h>
#include <libmesh/enum_quadrature_type.h>

#include <iomanip>
#include <numeric> // std::iota

#include "libmesh_cppunit.h"


using namespace libMesh;

#define MACROCOMMA ,

#if LIBMESH_DIM > 2
#define TEST_ONE_ORDER(qtype, order, maxorder)                          \
  CPPUNIT_TEST( testBuild<qtype MACROCOMMA order> );                    \
  CPPUNIT_TEST( test1DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> ); \
  CPPUNIT_TEST( test2DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> ); \
  CPPUNIT_TEST( test3DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> );
#elif LIBMESH_DIM > 1
#define TEST_ONE_ORDER(qtype, order, maxorder)                          \
  CPPUNIT_TEST( testBuild<qtype MACROCOMMA order> );                    \
  CPPUNIT_TEST( test1DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> ); \
  CPPUNIT_TEST( test2DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> );
#else
#define TEST_ONE_ORDER(qtype, order, maxorder)                          \
  CPPUNIT_TEST( testBuild<qtype MACROCOMMA order> );                    \
  CPPUNIT_TEST( test1DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> );
#endif

// std::min isn't constexpr, and C++03 lacks constexpr anyway
#define mymin(a, b) (a < b ? a : b)

#define TEST_ALL_ORDERS(qtype, maxorder)                \
  TEST_ONE_ORDER(qtype, FIRST, mymin(1,maxorder));      \
  TEST_ONE_ORDER(qtype, SECOND, mymin(2,maxorder));     \
  TEST_ONE_ORDER(qtype, THIRD, mymin(3,maxorder));      \
  TEST_ONE_ORDER(qtype, FOURTH, mymin(4,maxorder));     \
  TEST_ONE_ORDER(qtype, FIFTH, mymin(5,maxorder));      \
  TEST_ONE_ORDER(qtype, SIXTH, mymin(6,maxorder));      \
  TEST_ONE_ORDER(qtype, SEVENTH, mymin(7,maxorder));    \
  TEST_ONE_ORDER(qtype, EIGHTH, mymin(8,maxorder));     \
  TEST_ONE_ORDER(qtype, NINTH, mymin(9,maxorder));

#define LIBMESH_ASSERT_REALS_EQUAL(first, second, tolerance)            \
  if (std::abs(first-second) >= tolerance)                              \
    {                                                                   \
      std::cerr << "first = " << first << std::endl;                    \
      std::cerr << "second = " << second << std::endl;                  \
      std::cerr << "error = " << std::abs(first-second) << std::endl;   \
      std::cerr << "tolerance = " << tolerance << std::endl;            \
    }                                                                   \
  CPPUNIT_ASSERT (std::abs(first-second) < tolerance)

class QuadratureTest : public CppUnit::TestCase {
public:
  CPPUNIT_TEST_SUITE( QuadratureTest );

  TEST_ALL_ORDERS(QGAUSS, 9999);
  TEST_ONE_ORDER(QSIMPSON, FIRST,  3);
  TEST_ONE_ORDER(QSIMPSON, SECOND, 3);
  TEST_ONE_ORDER(QSIMPSON, THIRD,  3);
  TEST_ONE_ORDER(QTRAP, FIRST, 1);
  TEST_ALL_ORDERS(QGRID, 1);

  // In general, QNodal rules (e.g. QTRAP) are only exact for linears.
  // QSIMPSON is a special case of a nodal quadrature which obtains
  // higher accuracy.
  TEST_ONE_ORDER(QNODAL, /*ignored*/FIRST, /*max order=*/1);

  // The TEST_ALL_ORDERS macro only goes up to 9th-order
  TEST_ALL_ORDERS(QGAUSS_LOBATTO, 9);

  // The Gauss-Lobatto quadrature rules passed all these tests during
  // development, but we don't run them with every 'make check'
  // because it takes too long.
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, ELEVENTH,    11);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTEENTH,  13);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, FIFTEENTH,   15);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, SEVENTEENTH, 17);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, NINETEENTH,  19);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYFIRST, 21);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYTHIRD, 23);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYFIFTH, 25);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYSEVENTH, 27);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYNINTH, 29);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYFIRST, 31);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYTHIRD, 33);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYFIFTH, 35);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYSEVENTH, 37);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYNINTH, 39);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, FORTYFIRST, 41);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, FORTYTHIRD, 43);

  // Test monomial quadrature rules on quads and hexes
  CPPUNIT_TEST( testMonomialQuadrature );

  // Test quadrature rules on Triangles
#if LIBMESH_DIM > 1
  CPPUNIT_TEST( testTriQuadrature );
#endif

  // Test quadrature rules on Tetrahedra
#if LIBMESH_DIM > 2
  CPPUNIT_TEST( testTetQuadrature );
#endif

  // Test Jacobi quadrature rules with special weighting function
  CPPUNIT_TEST( testJacobi );

  CPPUNIT_TEST_SUITE_END();

private:

  Real quadrature_tolerance;


public:
  void setUp ()
  { quadrature_tolerance = TOLERANCE * std::sqrt(TOLERANCE); }

  void tearDown ()
  {}

  void testMonomialQuadrature ()
  {
    ElemType elem_type[2] = {QUAD4, HEX8};
    int dims[2]           = {2, 3};

    for (int i=0; i<(LIBMESH_DIM-1); ++i)
      {
        // std::cout << "\nChecking monomial quadrature on element type " << elem_type[i] << std::endl;

        for (int order=0; order<7; ++order)
          {
            std::unique_ptr<QBase> qrule = QBase::build(QMONOMIAL,
                                                        dims[i],
                                                        static_cast<Order>(order));
            qrule->init(elem_type[i]);

            // In 3D, max(z_power)==order, in 2D max(z_power)==0
            int max_z_power = dims[i]==2 ? 0 : order;

            int xyz_power[3];
            for (xyz_power[0]=0; xyz_power[0]<=order; ++xyz_power[0])
              for (xyz_power[1]=0; xyz_power[1]<=order; ++xyz_power[1])
                for (xyz_power[2]=0; xyz_power[2]<=max_z_power; ++xyz_power[2])
                  {
                    // Only try to integrate polynomials we can integrate exactly
                    if (xyz_power[0] + xyz_power[1] + xyz_power[2] > order)
                      continue;

                    // Compute the integral via quadrature.  Note that
                    // std::pow(0,0) returns 1 in the 2D case.
                    Real sumq = 0.;
                    for (unsigned int qp=0; qp<qrule->n_points(); qp++)
                      {
                        Real term = qrule->w(qp);
                        for (unsigned int d=0; d != LIBMESH_DIM; ++d)
                          term *= std::pow(qrule->qp(qp)(d), xyz_power[d]);
                        sumq += term;
                      }

                    // std::cout << "Quadrature of x^" << xyz_power[0]
                    //           << " y^" << xyz_power[1]
                    //           << " z^" << xyz_power[2]
                    //           << " = " << sumq << std::endl;

                    // Copy-pasted code from test3DWeights()
                    Real exact_x = (xyz_power[0] % 2) ? 0 : (Real(2.0) / (xyz_power[0]+1));
                    Real exact_y = (xyz_power[1] % 2) ? 0 : (Real(2.0) / (xyz_power[1]+1));
                    Real exact_z = (xyz_power[2] % 2) ? 0 : (Real(2.0) / (xyz_power[2]+1));

                    // Handle 2D
                    if (dims[i]==2)
                      exact_z = 1.0;

                    Real exact = exact_x*exact_y*exact_z;

                    // Make sure that the quadrature solution matches the exact solution
                    LIBMESH_ASSERT_REALS_EQUAL(exact, sumq, quadrature_tolerance);
                  }
          } // end for (order)
      } // end for (i)
  }

  void testTetQuadrature ()
  {
    // There are 3 different families of quadrature rules for tetrahedra
    QuadratureType qtype[3] = {QCONICAL, QGRUNDMANN_MOLLER, QGAUSS};

    for (int qt=0; qt<3; ++qt)
      for (int order=0; order<7; ++order)
        {
          std::unique_ptr<QBase> qrule = QBase::build(qtype[qt],
                                                      /*dim=*/3,
                                                      static_cast<Order>(order));

          // Initialize on a TET element
          qrule->init (TET4);

          // Test the sum of the weights for this order
          Real sumw = 0.;
          for (unsigned int qp=0; qp<qrule->n_points(); qp++)
            sumw += qrule->w(qp);
<commit_msg>More 64-bit-precision-only quadrature unit tests<commit_after>#include <libmesh/quadrature.h>
#include <libmesh/string_to_enum.h>
#include <libmesh/utility.h>
#include <libmesh/enum_quadrature_type.h>

#include <iomanip>
#include <numeric> // std::iota

#include "libmesh_cppunit.h"


using namespace libMesh;

#define MACROCOMMA ,

#if LIBMESH_DIM > 2
#define TEST_ONE_ORDER(qtype, order, maxorder)                          \
  CPPUNIT_TEST( testBuild<qtype MACROCOMMA order> );                    \
  CPPUNIT_TEST( test1DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> ); \
  CPPUNIT_TEST( test2DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> ); \
  CPPUNIT_TEST( test3DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> );
#elif LIBMESH_DIM > 1
#define TEST_ONE_ORDER(qtype, order, maxorder)                          \
  CPPUNIT_TEST( testBuild<qtype MACROCOMMA order> );                    \
  CPPUNIT_TEST( test1DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> ); \
  CPPUNIT_TEST( test2DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> );
#else
#define TEST_ONE_ORDER(qtype, order, maxorder)                          \
  CPPUNIT_TEST( testBuild<qtype MACROCOMMA order> );                    \
  CPPUNIT_TEST( test1DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> );
#endif

// std::min isn't constexpr, and C++03 lacks constexpr anyway
#define mymin(a, b) (a < b ? a : b)

#define TEST_ALL_ORDERS(qtype, maxorder)                \
  TEST_ONE_ORDER(qtype, FIRST, mymin(1,maxorder));      \
  TEST_ONE_ORDER(qtype, SECOND, mymin(2,maxorder));     \
  TEST_ONE_ORDER(qtype, THIRD, mymin(3,maxorder));      \
  TEST_ONE_ORDER(qtype, FOURTH, mymin(4,maxorder));     \
  TEST_ONE_ORDER(qtype, FIFTH, mymin(5,maxorder));      \
  TEST_ONE_ORDER(qtype, SIXTH, mymin(6,maxorder));      \
  TEST_ONE_ORDER(qtype, SEVENTH, mymin(7,maxorder));    \
  TEST_ONE_ORDER(qtype, EIGHTH, mymin(8,maxorder));     \
  TEST_ONE_ORDER(qtype, NINTH, mymin(9,maxorder));

#define LIBMESH_ASSERT_REALS_EQUAL(first, second, tolerance)            \
  if (std::abs(first-second) >= tolerance)                              \
    {                                                                   \
      std::cerr << "first = " << first << std::endl;                    \
      std::cerr << "second = " << second << std::endl;                  \
      std::cerr << "error = " << std::abs(first-second) << std::endl;   \
      std::cerr << "tolerance = " << tolerance << std::endl;            \
    }                                                                   \
  CPPUNIT_ASSERT (std::abs(first-second) < tolerance)

class QuadratureTest : public CppUnit::TestCase {
public:
  CPPUNIT_TEST_SUITE( QuadratureTest );

  TEST_ALL_ORDERS(QGAUSS, 9999);
  TEST_ONE_ORDER(QSIMPSON, FIRST,  3);
  TEST_ONE_ORDER(QSIMPSON, SECOND, 3);
  TEST_ONE_ORDER(QSIMPSON, THIRD,  3);
  TEST_ONE_ORDER(QTRAP, FIRST, 1);
  TEST_ALL_ORDERS(QGRID, 1);

  // In general, QNodal rules (e.g. QTRAP) are only exact for linears.
  // QSIMPSON is a special case of a nodal quadrature which obtains
  // higher accuracy.
  TEST_ONE_ORDER(QNODAL, /*ignored*/FIRST, /*max order=*/1);

  // The TEST_ALL_ORDERS macro only goes up to 9th-order
  TEST_ALL_ORDERS(QGAUSS_LOBATTO, 9);

  // The Gauss-Lobatto quadrature rules passed all these tests during
  // development, but we don't run them with every 'make check'
  // because it takes too long.
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, ELEVENTH,    11);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTEENTH,  13);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, FIFTEENTH,   15);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, SEVENTEENTH, 17);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, NINETEENTH,  19);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYFIRST, 21);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYTHIRD, 23);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYFIFTH, 25);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYSEVENTH, 27);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYNINTH, 29);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYFIRST, 31);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYTHIRD, 33);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYFIFTH, 35);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYSEVENTH, 37);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYNINTH, 39);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, FORTYFIRST, 41);
  // TEST_ONE_ORDER(QGAUSS_LOBATTO, FORTYTHIRD, 43);

  // Test monomial quadrature rules on quads and hexes
  CPPUNIT_TEST( testMonomialQuadrature );

  // Test quadrature rules on Triangles
#if LIBMESH_DIM > 1
  CPPUNIT_TEST( testTriQuadrature );
#endif

  // Test quadrature rules on Tetrahedra
#if LIBMESH_DIM > 2
  CPPUNIT_TEST( testTetQuadrature );
#endif

  // Test Jacobi quadrature rules with special weighting function
  CPPUNIT_TEST( testJacobi );

  CPPUNIT_TEST_SUITE_END();

private:

  Real quadrature_tolerance;


public:
  void setUp ()
  { quadrature_tolerance = TOLERANCE * std::sqrt(TOLERANCE); }

  void tearDown ()
  {}

  void testMonomialQuadrature ()
  {
    ElemType elem_type[2] = {QUAD4, HEX8};
    int dims[2]           = {2, 3};

    for (int i=0; i<(LIBMESH_DIM-1); ++i)
      {
        // std::cout << "\nChecking monomial quadrature on element type " << elem_type[i] << std::endl;

        for (int order=0; order<7; ++order)
          {
            std::unique_ptr<QBase> qrule = QBase::build(QMONOMIAL,
                                                        dims[i],
                                                        static_cast<Order>(order));
            qrule->init(elem_type[i]);

            // In 3D, max(z_power)==order, in 2D max(z_power)==0
            int max_z_power = dims[i]==2 ? 0 : order;

            int xyz_power[3];
            for (xyz_power[0]=0; xyz_power[0]<=order; ++xyz_power[0])
              for (xyz_power[1]=0; xyz_power[1]<=order; ++xyz_power[1])
                for (xyz_power[2]=0; xyz_power[2]<=max_z_power; ++xyz_power[2])
                  {
                    // Only try to integrate polynomials we can integrate exactly
                    if (xyz_power[0] + xyz_power[1] + xyz_power[2] > order)
                      continue;

                    // Compute the integral via quadrature.  Note that
                    // std::pow(0,0) returns 1 in the 2D case.
                    Real sumq = 0.;
                    for (unsigned int qp=0; qp<qrule->n_points(); qp++)
                      {
                        Real term = qrule->w(qp);
                        for (unsigned int d=0; d != LIBMESH_DIM; ++d)
                          term *= std::pow(qrule->qp(qp)(d), xyz_power[d]);
                        sumq += term;
                      }

                    // std::cout << "Quadrature of x^" << xyz_power[0]
                    //           << " y^" << xyz_power[1]
                    //           << " z^" << xyz_power[2]
                    //           << " = " << sumq << std::endl;

                    // Copy-pasted code from test3DWeights()
                    Real exact_x = (xyz_power[0] % 2) ? 0 : (Real(2.0) / (xyz_power[0]+1));
                    Real exact_y = (xyz_power[1] % 2) ? 0 : (Real(2.0) / (xyz_power[1]+1));
                    Real exact_z = (xyz_power[2] % 2) ? 0 : (Real(2.0) / (xyz_power[2]+1));

                    // Handle 2D
                    if (dims[i]==2)
                      exact_z = 1.0;

                    Real exact = exact_x*exact_y*exact_z;

                    // Make sure that the quadrature solution matches the exact solution
                    LIBMESH_ASSERT_REALS_EQUAL(exact, sumq, quadrature_tolerance);
                  }
          } // end for (order)
      } // end for (i)
  }

  void testTetQuadrature ()
  {
    // There are 3 different families of quadrature rules for tetrahedra
    QuadratureType qtype[3] = {QCONICAL, QGRUNDMANN_MOLLER, QGAUSS};

    int end_order = 7;
    // Our higher order tet rules were only computed to double precision
    if (quadrature_tolerance < 1e-16)
      end_order = 2;

    for (int qt=0; qt<3; ++qt)
      for (int order=0; order<end_order; ++order)
        {
          std::unique_ptr<QBase> qrule = QBase::build(qtype[qt],
                                                      /*dim=*/3,
                                                      static_cast<Order>(order));

          // Initialize on a TET element
          qrule->init (TET4);

          // Test the sum of the weights for this order
          Real sumw = 0.;
          for (unsigned int qp=0; qp<qrule->n_points(); qp++)
            sumw += qrule->w(qp);

{'fraction_non_alphanumeric': 0.12151642997207202, 'fraction_numerical': 0.017469843722146296, 'mean_word_length': 2.090909090909091, 'pattern_counts': {'":': 0, '<': 160, '<?xml version=': 0, '>': 67, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 34, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1865258', 'n_tokens_mistral': 6185, 'n_tokens_neox': 5744, 'n_words': 1193}