QuadratureTest Class Reference

Inheritance diagram for QuadratureTest:

Public Member Functions
	LIBMESH_CPPUNIT_TEST_SUITE (QuadratureTest)
	TEST_TWENTIETH_ORDER (QGAUSS, 9999)
	TEST_ONE_ORDER (QSIMPSON, FIRST, 1)
	TEST_ONE_ORDER (QSIMPSON, SECOND, 2)
	TEST_ONE_ORDER (QSIMPSON, THIRD, 3)
	TEST_ONE_ORDER (QTRAP, FIRST, 1)
	TEST_NINTH_ORDER (QGRID, 1)
	TEST_ONE_ORDER (QNODAL, FIRST, 1)
	TEST_NINTH_ORDER (QGAUSS_LOBATTO, 9)
	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)
	CPPUNIT_TEST (test1DWeights< QGAUSS MACROCOMMA TWENTYFIRST MACROCOMMA 21 >)
	CPPUNIT_TEST (test1DWeights< QGAUSS MACROCOMMA TWENTYFIFTH MACROCOMMA 25 >)
	CPPUNIT_TEST (test1DWeights< QGAUSS MACROCOMMA TWENTYSEVENTH MACROCOMMA 27 >)
	CPPUNIT_TEST (test1DWeights< QGAUSS MACROCOMMA TWENTYNINTH MACROCOMMA 29 >)
	CPPUNIT_TEST (test1DWeights< QGAUSS MACROCOMMA THIRTYFIRST MACROCOMMA 31 >)
	CPPUNIT_TEST (test1DWeights< QGAUSS MACROCOMMA THIRTYTHIRD MACROCOMMA 33 >)
	CPPUNIT_TEST (test1DWeights< QGAUSS MACROCOMMA THIRTYFIFTH MACROCOMMA 35 >)
	CPPUNIT_TEST (test1DWeights< QGAUSS MACROCOMMA THIRTYSEVENTH MACROCOMMA 37 >)
	CPPUNIT_TEST (test1DWeights< QGAUSS MACROCOMMA THIRTYNINTH MACROCOMMA 39 >)
	CPPUNIT_TEST (test1DWeights< QGAUSS MACROCOMMA FORTYFIRST MACROCOMMA 41 >)
	CPPUNIT_TEST (test1DWeights< QGAUSS MACROCOMMA FORTYTHIRD MACROCOMMA 43 >)
	CPPUNIT_TEST (testMonomialQuadrature)
	CPPUNIT_TEST (testNodalQuadrature)
	CPPUNIT_TEST (testTriQuadrature)
	CPPUNIT_TEST (testTetQuadrature)
	CPPUNIT_TEST (testJacobi)
	CPPUNIT_TEST_SUITE_END ()
void	setUp ()
void	tearDown ()
void	testNodalQuadrature ()
void	testMonomialQuadrature ()
void	testTetQuadrature ()
void	testTriQuadrature ()
void	testJacobi ()
template<QuadratureType qtype, Order order>
void	testBuild ()
template<QuadratureType qtype, Order order, unsigned int exactorder>
void	test1DWeights ()
template<QuadratureType qtype, Order order, unsigned int exactorder>
void	test2DWeights ()
template<QuadratureType qtype, Order order, unsigned int exactorder>
void	test3DWeights ()

Private Member Functions
void	testPolynomial (QBase &qrule, int xp, int yp, int zp, Real true_value)
void	testPolynomials (QuadratureType qtype, int order, ElemType elem_type, const std::function< Real(int, int, int)> &true_value, int exactorder)

Private Attributes
Real	quadrature_tolerance
const std::function< Real(int, int, int)>	edge_integrals
const std::function< Real(int, int, int)>	quad_integrals
const std::function< Real(int, int, int)>	tri_integrals
const std::function< Real(int, int, int)>	hex_integrals
const std::function< Real(int, int, int)>	tet_integrals
const std::function< Real(int, int, int)>	prism_integrals
const std::function< Real(int, int, int)>	pyramid_integrals

Detailed Description

Definition at line 72 of file quadrature_test.C.

Member Function Documentation

CPPUNIT_TEST() [1/16]

QuadratureTest::CPPUNIT_TEST

(

test1DWeights< QGAUSS MACROCOMMA FORTYFIRST MACROCOMMA 41 >

)

CPPUNIT_TEST() [2/16]

QuadratureTest::CPPUNIT_TEST

(

test1DWeights< QGAUSS MACROCOMMA FORTYTHIRD MACROCOMMA 43 >

)

CPPUNIT_TEST() [3/16]

QuadratureTest::CPPUNIT_TEST

(

test1DWeights< QGAUSS MACROCOMMA THIRTYFIFTH MACROCOMMA 35 >

)

CPPUNIT_TEST() [4/16]

QuadratureTest::CPPUNIT_TEST

(

test1DWeights< QGAUSS MACROCOMMA THIRTYFIRST MACROCOMMA 31 >

)

CPPUNIT_TEST() [5/16]

QuadratureTest::CPPUNIT_TEST

(

test1DWeights< QGAUSS MACROCOMMA THIRTYNINTH MACROCOMMA 39 >

)

CPPUNIT_TEST() [6/16]

QuadratureTest::CPPUNIT_TEST

(

test1DWeights< QGAUSS MACROCOMMA THIRTYSEVENTH MACROCOMMA 37 >

)

CPPUNIT_TEST() [7/16]

QuadratureTest::CPPUNIT_TEST

(

test1DWeights< QGAUSS MACROCOMMA THIRTYTHIRD MACROCOMMA 33 >

)

CPPUNIT_TEST() [8/16]

QuadratureTest::CPPUNIT_TEST

(

test1DWeights< QGAUSS MACROCOMMA TWENTYFIFTH MACROCOMMA 25 >

)

CPPUNIT_TEST() [9/16]

QuadratureTest::CPPUNIT_TEST

(

test1DWeights< QGAUSS MACROCOMMA TWENTYFIRST MACROCOMMA 21 >

)

CPPUNIT_TEST() [10/16]

QuadratureTest::CPPUNIT_TEST

(

test1DWeights< QGAUSS MACROCOMMA TWENTYNINTH MACROCOMMA 29 >

)

CPPUNIT_TEST() [11/16]

QuadratureTest::CPPUNIT_TEST

(

test1DWeights< QGAUSS MACROCOMMA TWENTYSEVENTH MACROCOMMA 27 >

)

CPPUNIT_TEST() [12/16]

QuadratureTest::CPPUNIT_TEST

(

testJacobi

)

CPPUNIT_TEST() [13/16]

QuadratureTest::CPPUNIT_TEST

(

testMonomialQuadrature

)

CPPUNIT_TEST() [14/16]

QuadratureTest::CPPUNIT_TEST

(

testNodalQuadrature

)

CPPUNIT_TEST() [15/16]

QuadratureTest::CPPUNIT_TEST

(

testTetQuadrature

)

CPPUNIT_TEST() [16/16]

QuadratureTest::CPPUNIT_TEST

(

testTriQuadrature

)

CPPUNIT_TEST_SUITE_END()

QuadratureTest::CPPUNIT_TEST_SUITE_END

(

)

LIBMESH_CPPUNIT_TEST_SUITE()

QuadratureTest::LIBMESH_CPPUNIT_TEST_SUITE

(

QuadratureTest

)

setUp()

void QuadratureTest::setUp ( )

inline

Definition at line 321 of file quadrature_test.C.

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

References quadrature_tolerance, and libMesh::TOLERANCE.

tearDown()

void QuadratureTest::tearDown ( )

inline

Definition at line 324 of file quadrature_test.C.

  {}

test1DWeights()

template<QuadratureType qtype, Order order, unsigned int exactorder>

void QuadratureTest::test1DWeights ( )

inline

Definition at line 546 of file quadrature_test.C.

  {
    LOG_UNIT_TEST;
 
    testPolynomials(qtype, order, EDGE3, edge_integrals, exactorder);
  }

References libMesh::EDGE3, edge_integrals, and testPolynomials().

test2DWeights()

template<QuadratureType qtype, Order order, unsigned int exactorder>

void QuadratureTest::test2DWeights ( )

inline

Definition at line 558 of file quadrature_test.C.

  {
    LOG_UNIT_TEST;
 
    testPolynomials(qtype, order, QUAD8, quad_integrals, exactorder);
 
    // We may eventually support Gauss-Lobatto type quadrature on triangles...
    if (qtype == QGAUSS_LOBATTO)
      return;
 
    // QGrid needs to be changed to use symmetric offsets on triangles
    // so it can at *least* get linears right...
    if (qtype == QGRID)
      testPolynomials(qtype, order, TRI6, tri_integrals, 0);
    // QSimpson doesn't even get all quadratics on a triangle
    else if (qtype == QSIMPSON)
      testPolynomials(qtype, order, TRI6, tri_integrals, std::min(1u,exactorder));
    else
      {
        // Our 18th order and higher triangle rules were only computed
        // in double precision, so we use a lower tolerance here.
        if ((sizeof(Real) > sizeof(double)) && (order > 17))
          quadrature_tolerance = TOLERANCE;
 
        testPolynomials(qtype, order, TRI6, tri_integrals, exactorder);
      }
  }

References libMesh::QGAUSS_LOBATTO, libMesh::QGRID, libMesh::QSIMPSON, libMesh::QUAD8, quad_integrals, quadrature_tolerance, libMesh::Real, testPolynomials(), libMesh::TOLERANCE, libMesh::TRI6, and tri_integrals.

test3DWeights()

template<QuadratureType qtype, Order order, unsigned int exactorder>

void QuadratureTest::test3DWeights ( )

inline

Definition at line 591 of file quadrature_test.C.

  {
    LOG_UNIT_TEST;
 
    testPolynomials(qtype, order, HEX20, hex_integrals, exactorder);
 
    // We may eventually support Gauss-Lobatto type quadrature on tets and prisms...
    if (qtype == QGAUSS_LOBATTO)
      return;
 
    // QGrid needs to be changed to use symmetric offsets on
    // non-tensor product elements so it can at *least* get linears
    // right...
    if (qtype == QGRID)
      {
        testPolynomials(qtype, order, TET10, tet_integrals, 0);
        testPolynomials(qtype, order, PYRAMID14, pyramid_integrals, 0);
      }
    // QSimpson doesn't get all quadratics on a simplex or its extrusion
    else if (qtype == QSIMPSON)
      {
        testPolynomials(qtype, order, TET10, tet_integrals, std::min(1u,exactorder));
        testPolynomials(qtype, order, PRISM15, prism_integrals, std::min(1u,exactorder));
 
        // And on pyramids we gave up and redid QTrap
        testPolynomials(qtype, order, PYRAMID14, pyramid_integrals, std::min(1u,exactorder));
      }
    else
      {
        testPolynomials(qtype, order, TET10, tet_integrals, exactorder);
        testPolynomials(qtype, order, PYRAMID14, pyramid_integrals, exactorder);
 
        // Our 18th order and higher triangle rules (used to construct
        // prism rules) were only computed in double precision, so we
        // use a lower tolerance here.
        if ((sizeof(Real) > sizeof(double)) && (order > 17))
          quadrature_tolerance = TOLERANCE;
 
        testPolynomials(qtype, order, PRISM15, prism_integrals, exactorder);
      }
  }

References libMesh::HEX20, hex_integrals, libMesh::PRISM15, prism_integrals, libMesh::PYRAMID14, pyramid_integrals, libMesh::QGAUSS_LOBATTO, libMesh::QGRID, libMesh::QSIMPSON, quadrature_tolerance, libMesh::Real, testPolynomials(), libMesh::TET10, tet_integrals, and libMesh::TOLERANCE.

TEST_NINTH_ORDER() [1/2]

QuadratureTest::TEST_NINTH_ORDER	(	QGAUSS_LOBATTO	,
		9
	)

TEST_NINTH_ORDER() [2/2]

QuadratureTest::TEST_NINTH_ORDER	(	QGRID	,
		1
	)

TEST_ONE_ORDER() [1/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		ELEVENTH	,
		11
	)

TEST_ONE_ORDER() [2/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		FIFTEENTH	,
		15
	)

TEST_ONE_ORDER() [3/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		FORTYFIRST	,
		41
	)

TEST_ONE_ORDER() [4/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		FORTYTHIRD	,
		43
	)

TEST_ONE_ORDER() [5/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		NINETEENTH	,
		19
	)

TEST_ONE_ORDER() [6/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		SEVENTEENTH	,
		17
	)

TEST_ONE_ORDER() [7/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		THIRTEENTH	,
		13
	)

TEST_ONE_ORDER() [8/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		THIRTYFIFTH	,
		35
	)

TEST_ONE_ORDER() [9/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		THIRTYFIRST	,
		31
	)

TEST_ONE_ORDER() [10/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		THIRTYNINTH	,
		39
	)

TEST_ONE_ORDER() [11/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		THIRTYSEVENTH	,
		37
	)

TEST_ONE_ORDER() [12/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		THIRTYTHIRD	,
		33
	)

TEST_ONE_ORDER() [13/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		TWENTYFIFTH	,
		25
	)

TEST_ONE_ORDER() [14/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		TWENTYFIRST	,
		21
	)

TEST_ONE_ORDER() [15/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		TWENTYNINTH	,
		29
	)

TEST_ONE_ORDER() [16/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		TWENTYSEVENTH	,
		27
	)

TEST_ONE_ORDER() [17/22]

QuadratureTest::TEST_ONE_ORDER	(	QGAUSS_LOBATTO	,
		TWENTYTHIRD	,
		23
	)

TEST_ONE_ORDER() [18/22]

QuadratureTest::TEST_ONE_ORDER	(	QNODAL	,
		FIRST	,
		1
	)

TEST_ONE_ORDER() [19/22]

QuadratureTest::TEST_ONE_ORDER	(	QSIMPSON	,
		FIRST	,
		1
	)

TEST_ONE_ORDER() [20/22]

QuadratureTest::TEST_ONE_ORDER	(	QSIMPSON	,
		SECOND	,
		2
	)

TEST_ONE_ORDER() [21/22]

QuadratureTest::TEST_ONE_ORDER	(	QSIMPSON	,
		THIRD	,
		3
	)

TEST_ONE_ORDER() [22/22]

QuadratureTest::TEST_ONE_ORDER	(	QTRAP	,
		FIRST	,
		1
	)

TEST_TWENTIETH_ORDER()

QuadratureTest::TEST_TWENTIETH_ORDER	(	QGAUSS	,
		9999
	)

testBuild()

template<QuadratureType qtype, Order order>

void QuadratureTest::testBuild ( )

inline

Definition at line 512 of file quadrature_test.C.

  {
    LOG_UNIT_TEST;
 
    std::unique_ptr<QBase> qrule1D = QBase::build (qtype, 1, order);
    std::unique_ptr<QBase> qrule2D = QBase::build (qtype, 2, order);
    std::unique_ptr<QBase> qrule3D = QBase::build (qtype, 3, order);
 
    CPPUNIT_ASSERT_EQUAL ( qtype, qrule1D->type() );
    CPPUNIT_ASSERT_EQUAL ( qtype, qrule2D->type() );
    CPPUNIT_ASSERT_EQUAL ( qtype, qrule3D->type() );
 
    // Simpson's Rule is inherently third-order
    if (qtype != QSIMPSON)
      {
        CPPUNIT_ASSERT_EQUAL ( order, qrule1D->get_base_order() );
        CPPUNIT_ASSERT_EQUAL ( order, qrule2D->get_base_order() );
        CPPUNIT_ASSERT_EQUAL ( order, qrule3D->get_base_order() );
      }
 
    CPPUNIT_ASSERT_EQUAL ( static_cast<unsigned int>(1) , qrule1D->get_dim() );
    CPPUNIT_ASSERT_EQUAL ( static_cast<unsigned int>(2) , qrule2D->get_dim() );
    CPPUNIT_ASSERT_EQUAL ( static_cast<unsigned int>(3) , qrule3D->get_dim() );
 
    // Test the enum-to-string conversion for this qtype is
    // implemented, but not what the actual value is.
    Utility::enum_to_string(qtype);
  }

References libMesh::QBase::build(), libMesh::Utility::enum_to_string(), and libMesh::QSIMPSON.

testJacobi()

void QuadratureTest::testJacobi ( )

inline

Definition at line 434 of file quadrature_test.C.

  {
    LOG_UNIT_TEST;
 
    constexpr int max_order = 43;  // We seriously implemented that?!
 
    // LibMesh supports two different types of Jacobi quadrature
    QuadratureType qtype[2] = {QJACOBI_1_0, QJACOBI_2_0};
 
    // The weights of the Jacobi quadrature rules in libmesh have been
    // scaled based on their intended use:
    // (alpha=1, beta=0) rule weights sum to 1/2.
    // (alpha=2, beta=0) rule weights sum to 1/3.
    Real initial_sum_weights[2] = {.5, 1./3.};
 
    // The points of the Jacobi rules in LibMesh are also defined on
    // [0,1]... this has to be taken into account when computing the
    // exact integral values in Maple!  Also note: if you scale the
    // points to a different interval, you need to also compute what
    // the sum of the weights should be for that interval, it will not
    // simply be the element length for weighted quadrature rules like
    // Jacobi.  For general alpha and beta=0, the sum of the weights
    // on the interval [-1,1] is 2^(alpha+1) / (alpha+1).
    std::vector<std::vector<Real>> true_integrals(2);
 
    // alpha=1 integral values
    // int((1-x)*x^p, x=0..1) = 1 / (p^2 + 3p + 2)
    true_integrals[0].resize(max_order);
    for (std::size_t p=0; p<true_integrals[0].size(); ++p)
      true_integrals[0][p] = 1. / (p*p + 3.*p + 2.);
 
    // alpha=2 integral values
    // int((1-x)^2*x^p, x=0..1) = 2 / (p^3 + 6*p^2 + 11*p + 6)
    true_integrals[1].resize(max_order);
    for (std::size_t p=0; p<true_integrals[1].size(); ++p)
      true_integrals[1][p] = 2. / (p*p*p + 6.*p*p + 11.*p + 6.);
 
    // Test both types of Jacobi quadrature rules
    for (int qt=0; qt<2; ++qt)
      {
        for (int order=0; order<max_order; ++order)
          {
            std::unique_ptr<QBase> qrule = QBase::build(qtype[qt],
                                                        /*dim=*/1,
                                                        static_cast<Order>(order));
 
            // Initialize on a 1D element, EDGE2/3/4 should not matter...
            qrule->init (EDGE2, 0, true);
 
            // 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);
 
            // Make sure that the weights add up to the value we expect
            LIBMESH_ASSERT_REALS_EQUAL(initial_sum_weights[qt], sumw, quadrature_tolerance);
 
            // Test integrating different polynomial powers
            for (int testpower=0; testpower<=order; ++testpower)
              {
                // Note that the weighting function, (1-x)^alpha *
                // (1+x)^beta, is built into these quadrature rules;
                // the polynomials we actually integrate are just the
                // usual monomials.
                Real sumq = 0.;
                for (unsigned int qp=0; qp<qrule->n_points(); qp++)
                  sumq += qrule->w(qp) * std::pow(qrule->qp(qp)(0), testpower);
 
                // Make sure that the computed integral agrees with the "true" value
                LIBMESH_ASSERT_REALS_EQUAL(true_integrals[qt][testpower], sumq, quadrature_tolerance);
              } // end for(testpower)
          } // end for(order)
      } // end for(qt)
  } // testJacobi

References libMesh::QBase::build(), libMesh::EDGE2, libMesh::QJACOBI_1_0, libMesh::QJACOBI_2_0, quadrature_tolerance, and libMesh::Real.

testMonomialQuadrature()

void QuadratureTest::testMonomialQuadrature ( )

inline

Definition at line 385 of file quadrature_test.C.

  {
    LOG_UNIT_TEST;
 
    const std::vector<std::vector<ElemType>> all_types =
      {{EDGE2}, {QUAD4, TRI3}, {HEX8, TET4, PRISM6, PYRAMID5}};
    const std::vector<std::vector<std::function<Real(int,int,int)>>>
      true_values =
      {{edge_integrals},
       {quad_integrals, tri_integrals},
       {hex_integrals, tet_integrals, prism_integrals, pyramid_integrals}};
 
    for (auto i : index_range(all_types))
      for (auto j : index_range(all_types[i]))
        for (int order=0; order<17; ++order)
          testPolynomials(QMONOMIAL, order, all_types[i][j], true_values[i][j], order);
  }

References libMesh::EDGE2, edge_integrals, libMesh::HEX8, hex_integrals, libMesh::index_range(), libMesh::PRISM6, prism_integrals, libMesh::PYRAMID5, pyramid_integrals, libMesh::QMONOMIAL, libMesh::QUAD4, quad_integrals, libMesh::Real, testPolynomials(), libMesh::TET4, tet_integrals, libMesh::TRI3, and tri_integrals.

testNodalQuadrature()

void QuadratureTest::testNodalQuadrature ( )

inline

Definition at line 327 of file quadrature_test.C.

  {
    LOG_UNIT_TEST;
 
    const std::vector<std::vector<ElemType>> all_types =
      {{EDGE2, EDGE3, EDGE4},
       {TRI3, TRI6, TRI7},
       {QUAD4, QUAD8, QUAD9},
       {TET4, TET10, TET14},
       {HEX8, HEX20, HEX27},
       {PRISM6, PRISM15, PRISM18, PRISM20, PRISM21},
       {PYRAMID5, PYRAMID13, PYRAMID14, PYRAMID18}};
 
    const std::function<Real(int,int,int)> true_values[] =
      {edge_integrals,
       tri_integrals,
       quad_integrals,
       tet_integrals,
       hex_integrals,
       prism_integrals,
       pyramid_integrals};
 
    for (auto i : index_range(all_types))
      for (ElemType elem_type : all_types[i])
        {
          const unsigned int order = Elem::build(elem_type)->default_order();
 
          unsigned int exactorder = order;
          // There are some sad exceptions to the "positive nodal
          // quadrature can integrate polynomials up to the element
          // order" rule.
          //
          // Some "quadratic" elements can only exactly integrate
          // linears when mass lumping.
          if (elem_type == TRI6 || elem_type == QUAD8 ||
              elem_type == TET10 || elem_type == HEX20 ||
              elem_type == PRISM15 || elem_type == PRISM18 ||
              elem_type == PYRAMID13 || elem_type == PYRAMID14 ||
          // And some partially-cubic elements can only do quadratics
              elem_type == TET14 || elem_type == PRISM20)
            exactorder--;
 
          // And one partially-cubic element can only do linears.
          // Seriously.
          if (elem_type == PYRAMID18)
            exactorder=1;
 
          // Our PRISM20 rule was only computed in double precision.
          if ((sizeof(Real) > sizeof(double)) && elem_type == PRISM20)
            quadrature_tolerance = TOLERANCE;
          else // restore the default
            quadrature_tolerance = TOLERANCE * std::sqrt(TOLERANCE);
 
          testPolynomials(QNODAL, order, elem_type, true_values[i], exactorder);
        }
  }

References libMesh::Elem::build(), libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, edge_integrals, libMesh::HEX20, libMesh::HEX27, libMesh::HEX8, hex_integrals, libMesh::index_range(), libMesh::PRISM15, libMesh::PRISM18, libMesh::PRISM20, libMesh::PRISM21, libMesh::PRISM6, prism_integrals, libMesh::PYRAMID13, libMesh::PYRAMID14, libMesh::PYRAMID18, libMesh::PYRAMID5, pyramid_integrals, libMesh::QNODAL, libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, quad_integrals, quadrature_tolerance, libMesh::Real, testPolynomials(), libMesh::TET10, libMesh::TET14, libMesh::TET4, tet_integrals, libMesh::TOLERANCE, libMesh::TRI3, libMesh::TRI6, libMesh::TRI7, and tri_integrals.

testPolynomial()

void QuadratureTest::testPolynomial ( QBase &  qrule, int  xp, int  yp, int  zp, Real  true_value  )

inlineprivate

Definition at line 150 of file quadrature_test.C.

  {
    Real sum = 0.;
    for (unsigned int qp=0; qp<qrule.n_points(); qp++)
      {
        const Point p = qrule.qp(qp);
        Real val = qrule.w(qp) * std::pow(p(0),xp);
        if (qrule.get_dim() > 1)
          val *= std::pow(p(1),yp);
        if (qrule.get_dim() > 2)
          val *= std::pow(p(2),zp);
        sum += val;
      }
 
    // Make sure that the weighted evaluations add up to the true
    // integral value we expect
    LIBMESH_ASSERT_REALS_EQUAL(true_value, sum, quadrature_tolerance);
  }

References libMesh::QBase::get_dim(), libMesh::QBase::n_points(), libMesh::QBase::qp(), quadrature_tolerance, libMesh::Real, and libMesh::QBase::w().

Referenced by testPolynomials().

testPolynomials()

void QuadratureTest::testPolynomials ( QuadratureType  qtype, int  order, ElemType  elem_type, const std::function< Real(int, int, int)> &  true_value, int  exactorder  )

inlineprivate

Definition at line 171 of file quadrature_test.C.

  {
    auto elem = Elem::build(elem_type);
    unsigned int dim = elem->dim();
    std::unique_ptr<QBase> qrule =
      QBase::build(qtype, dim, static_cast<Order>(order));
 
    // We are testing integration on the reference elements here, so
    // the element map is not relevant, and we can safely use nodal
    // Pyramid quadrature.
    if (elem_type == PYRAMID5 || elem_type == PYRAMID13 ||
        elem_type == PYRAMID14 || elem_type == PYRAMID18)
      qrule->allow_nodal_pyramid_quadrature = true;
 
    qrule->init (*elem);
 
    const int max_y_order = dim>1 ? 0 : exactorder;
    const int max_z_order = dim>2 ? 0 : exactorder;
 
    for (int xp=0; xp <= exactorder; ++xp)
      for (int yp=0; yp <= max_y_order; ++yp)
        for (int zp=0; zp <= max_z_order; ++zp)
          {
            // Only try to integrate polynomials we can integrate exactly
            if (xp + yp + zp > exactorder)
              continue;
 
            const Real exact = true_value(xp, yp, zp);
            testPolynomial(*qrule, xp, yp, zp, exact);
          }
  }

References libMesh::Elem::build(), libMesh::QBase::build(), dim, libMesh::PYRAMID13, libMesh::PYRAMID14, libMesh::PYRAMID18, libMesh::PYRAMID5, libMesh::Real, and testPolynomial().

Referenced by test1DWeights(), test2DWeights(), test3DWeights(), testMonomialQuadrature(), testNodalQuadrature(), testTetQuadrature(), and testTriQuadrature().

testTetQuadrature()

void QuadratureTest::testTetQuadrature ( )

inline

Definition at line 403 of file quadrature_test.C.

  {
    LOG_UNIT_TEST;
 
    // There are 3 different families of quadrature rules for tetrahedra
    QuadratureType qtype[3] = {QCONICAL, QGRUNDMANN_MOLLER, QGAUSS};
 
    int end_order = 7;
 
    for (int qt=0; qt<3; ++qt)
      for (int order=0; order<end_order; ++order)
        {
          // Our higher order tet rules were only computed to double
          // precision
          if ((sizeof(Real) > sizeof(double)) && order > 2)
            quadrature_tolerance = TOLERANCE;
          testPolynomials(qtype[qt], order, TET4, tet_integrals, order);
        }
  }

References libMesh::QCONICAL, libMesh::QGAUSS, libMesh::QGRUNDMANN_MOLLER, quadrature_tolerance, libMesh::Real, testPolynomials(), libMesh::TET4, tet_integrals, and libMesh::TOLERANCE.

testTriQuadrature()

void QuadratureTest::testTriQuadrature ( )

inline

Definition at line 423 of file quadrature_test.C.

  {
    LOG_UNIT_TEST;
 
    QuadratureType qtype[4] = {QCONICAL, QCLOUGH, QGAUSS, QGRUNDMANN_MOLLER};
 
    for (int qt=0; qt<4; ++qt)
      for (int order=0; order<10; ++order)
        testPolynomials(qtype[qt], order, TRI3, tri_integrals, order);
  }

References libMesh::QCLOUGH, libMesh::QCONICAL, libMesh::QGAUSS, libMesh::QGRUNDMANN_MOLLER, testPolynomials(), libMesh::TRI3, and tri_integrals.

Member Data Documentation

edge_integrals

const std::function<Real(int,int,int)> QuadratureTest::edge_integrals

private

Initial value:

=
  [](int mode, int, int) {
    return (mode % 2) ?  0 : (Real(2.0) / (mode+1));
  }

Definition at line 206 of file quadrature_test.C.

                         {
    return (mode % 2) ?  0 : (Real(2.0) / (mode+1));
  };

Referenced by test1DWeights(), testMonomialQuadrature(), and testNodalQuadrature().

hex_integrals

const std::function<Real(int,int,int)> QuadratureTest::hex_integrals

private

Initial value:

=
  [](int modex, int modey, int modez) {
    const Real exactx = (modex % 2) ?
      0 : (Real(2.0) / (modex+1));
 
    const Real exacty = (modey % 2) ?
      0 : (Real(2.0) / (modey+1));
 
    const Real exactz = (modez % 2) ?
      0 : (Real(2.0) / (modez+1));
 
    return exactx*exacty*exactz;
  }

Definition at line 251 of file quadrature_test.C.

                                      {
    const Real exactx = (modex % 2) ?
      0 : (Real(2.0) / (modex+1));
 
    const Real exacty = (modey % 2) ?
      0 : (Real(2.0) / (modey+1));
 
    const Real exactz = (modez % 2) ?
      0 : (Real(2.0) / (modez+1));
 
    return exactx*exacty*exactz;
  };

Referenced by test3DWeights(), testMonomialQuadrature(), and testNodalQuadrature().

prism_integrals

const std::function<Real(int,int,int)> QuadratureTest::prism_integrals

private

Initial value:

=
  [this](int modex, int modey, int modez) {
    const Real exactz = (modez % 2) ?
      0 : (Real(2.0) / (modez+1));
 
    return exactz * tri_integrals(modex, modey, 0);
  }

Definition at line 300 of file quadrature_test.C.

                                          {
    const Real exactz = (modez % 2) ?
      0 : (Real(2.0) / (modez+1));
 
    return exactz * tri_integrals(modex, modey, 0);
  };

Referenced by test3DWeights(), testMonomialQuadrature(), and testNodalQuadrature().

pyramid_integrals

const std::function<Real(int,int,int)> QuadratureTest::pyramid_integrals

private

Initial value:

=
  [](int modex, int modey, int modez) {
 
    const int binom = Utility::binomial(modex+modey+modez+3, modez);
 
    if (modex%2 || modey%2)
      return Real(0);
 
    return Real(4)/((modex+1)*(modey+1)*binom*(modex+modey+modez+3));
  }

Definition at line 308 of file quadrature_test.C.

                                      {
 
    const int binom = Utility::binomial(modex+modey+modez+3, modez);
 
    if (modex%2 || modey%2)
      return Real(0);
 
    return Real(4)/((modex+1)*(modey+1)*binom*(modex+modey+modez+3));
  };

Referenced by test3DWeights(), testMonomialQuadrature(), and testNodalQuadrature().

quad_integrals

const std::function<Real(int,int,int)> QuadratureTest::quad_integrals

private

Initial value:

=
  [](int modex, int modey, int) {
    const Real exactx = (modex % 2) ?
      0 : (Real(2.0) / (modex+1));
 
    const Real exacty = (modey % 2) ?
      0 : (Real(2.0) / (modey+1));
 
    return exactx*exacty;
  }

Definition at line 211 of file quadrature_test.C.

                                {
    const Real exactx = (modex % 2) ?
      0 : (Real(2.0) / (modex+1));
 
    const Real exacty = (modey % 2) ?
      0 : (Real(2.0) / (modey+1));
 
    return exactx*exacty;
  };

Referenced by test2DWeights(), testMonomialQuadrature(), and testNodalQuadrature().

quadrature_tolerance

Real QuadratureTest::quadrature_tolerance

private

Definition at line 148 of file quadrature_test.C.

Referenced by setUp(), test2DWeights(), test3DWeights(), testJacobi(), testNodalQuadrature(), testPolynomial(), and testTetQuadrature().

tet_integrals

const std::function<Real(int,int,int)> QuadratureTest::tet_integrals

private

Definition at line 265 of file quadrature_test.C.

                                            {
    // Compute the true integral, a! b! c! / (a + b + c + 3)!
    Real analytical = 1.0;
 
    // Sort the a, b, c values
    int sorted_powers[3] = {x_power, y_power, z_power};
    std::sort(sorted_powers, sorted_powers+3);
 
    // Cancel the largest power with the denominator, fill in the
    // entries for the remaining numerator terms and the denominator.
    std::vector<int>
      numerator_1(sorted_powers[0] > 1 ? sorted_powers[0]-1 : 0),
      numerator_2(sorted_powers[1] > 1 ? sorted_powers[1]-1 : 0),
      denominator(3 + sorted_powers[0] + sorted_powers[1]);
 
    // Fill up the vectors with sequences starting at the right values.
    std::iota(numerator_1.begin(), numerator_1.end(), 2);
    std::iota(numerator_2.begin(), numerator_2.end(), 2);
    std::iota(denominator.begin(), denominator.end(), sorted_powers[2]+1);
 
    // The denominator is guaranteed to have the most terms...
    for (std::size_t i=0; i<denominator.size(); ++i)
      {
        if (i < numerator_1.size())
          analytical *= numerator_1[i];
 
        if (i < numerator_2.size())
          analytical *= numerator_2[i];
 
        analytical /= denominator[i];
      }
    return analytical;
  };

Referenced by test3DWeights(), testMonomialQuadrature(), testNodalQuadrature(), and testTetQuadrature().

tri_integrals

const std::function<Real(int,int,int)> QuadratureTest::tri_integrals

private

Initial value:

=
  [](int x_power, int y_power, int) {
    
    Real analytical = 1.0;
 
    unsigned
      larger_power = std::max(x_power, y_power),
      smaller_power = std::min(x_power, y_power);
 
    
    
    std::vector<unsigned>
      numerator(smaller_power > 1 ? smaller_power-1 : 0),
      denominator(2+smaller_power);
 
    
    std::iota(numerator.begin(), numerator.end(), 2);
    std::iota(denominator.begin(), denominator.end(), larger_power+1);
 
    
    for (std::size_t i=0; i<denominator.size(); ++i)
      {
        if (i < numerator.size())
          analytical *= numerator[i];
        analytical /= denominator[i];
      }
    return analytical;
  }

Definition at line 222 of file quadrature_test.C.

                                    {
    // Compute the true integral, a! b! / (a + b + 2)!
    Real analytical = 1.0;
 
    unsigned
      larger_power = std::max(x_power, y_power),
      smaller_power = std::min(x_power, y_power);
 
    // Cancel the larger of the two numerator terms with the
    // denominator, and fill in the remaining entries.
    std::vector<unsigned>
      numerator(smaller_power > 1 ? smaller_power-1 : 0),
      denominator(2+smaller_power);
 
    // Fill up the vectors with sequences starting at the right values.
    std::iota(numerator.begin(), numerator.end(), 2);
    std::iota(denominator.begin(), denominator.end(), larger_power+1);
 
    // The denominator is guaranteed to have more terms...
    for (std::size_t i=0; i<denominator.size(); ++i)
      {
        if (i < numerator.size())
          analytical *= numerator[i];
        analytical /= denominator[i];
      }
    return analytical;
  };

Referenced by test2DWeights(), testMonomialQuadrature(), testNodalQuadrature(), and testTriQuadrature().

---

The documentation for this class was generated from the following file:

• tests/quadrature/quadrature_test.C