libMesh
quadrature_test.C
Go to the documentation of this file.
1 // Ignore unused parameter warnings coming from cppunit headers
2 #include <libmesh/ignore_warnings.h>
3 #include <cppunit/extensions/HelperMacros.h>
4 #include <cppunit/TestCase.h>
5 #include <libmesh/restore_warnings.h>
6 
7 #include <libmesh/quadrature.h>
8 #include <libmesh/string_to_enum.h>
9 #include <libmesh/utility.h>
10 
11 #include <iomanip>
12 
13 // THE CPPUNIT_TEST_SUITE_END macro expands to code that involves
14 // std::auto_ptr, which in turn produces -Wdeprecated-declarations
15 // warnings. These can be ignored in GCC as long as we wrap the
16 // offending code in appropriate pragmas. We can't get away with a
17 // single ignore_warnings.h inclusion at the beginning of this file,
18 // since the libmesh headers pull in a restore_warnings.h at some
19 // point. We also don't bother restoring warnings at the end of this
20 // file since it's not a header.
21 #include <libmesh/ignore_warnings.h>
22 
23 using namespace libMesh;
24 
25 #define MACROCOMMA ,
26 
27 #define TEST_ONE_ORDER(qtype, order, maxorder) \
28  CPPUNIT_TEST( testBuild<qtype MACROCOMMA order> ); \
29  CPPUNIT_TEST( test1DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> ); \
30  CPPUNIT_TEST( test2DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> ); \
31  CPPUNIT_TEST( test3DWeights<qtype MACROCOMMA order MACROCOMMA maxorder> );
32 
33 // std::min isn't constexpr, and C++03 lacks constexpr anyway
34 #define mymin(a, b) (a < b ? a : b)
35 
36 #define TEST_ALL_ORDERS(qtype, maxorder) \
37  TEST_ONE_ORDER(qtype, FIRST, mymin(1,maxorder)); \
38  TEST_ONE_ORDER(qtype, SECOND, mymin(2,maxorder)); \
39  TEST_ONE_ORDER(qtype, THIRD, mymin(3,maxorder)); \
40  TEST_ONE_ORDER(qtype, FOURTH, mymin(4,maxorder)); \
41  TEST_ONE_ORDER(qtype, FIFTH, mymin(5,maxorder)); \
42  TEST_ONE_ORDER(qtype, SIXTH, mymin(6,maxorder)); \
43  TEST_ONE_ORDER(qtype, SEVENTH, mymin(7,maxorder)); \
44  TEST_ONE_ORDER(qtype, EIGHTH, mymin(8,maxorder)); \
45  TEST_ONE_ORDER(qtype, NINTH, mymin(9,maxorder));
46 
47 #define LIBMESH_ASSERT_REALS_EQUAL(first, second, tolerance) \
48  if (std::abs(first-second) >= tolerance) \
49  { \
50  std::cerr << "first = " << first << std::endl; \
51  std::cerr << "second = " << second << std::endl; \
52  std::cerr << "error = " << std::abs(first-second) << std::endl; \
53  std::cerr << "tolerance = " << tolerance << std::endl; \
54  } \
55  CPPUNIT_ASSERT (std::abs(first-second) < tolerance)
56 
57 class QuadratureTest : public CppUnit::TestCase {
58 public:
59  CPPUNIT_TEST_SUITE( QuadratureTest );
60 
61  TEST_ALL_ORDERS(QGAUSS, 9999);
62  TEST_ONE_ORDER(QSIMPSON, FIRST, 3);
63  TEST_ONE_ORDER(QSIMPSON, SECOND, 3);
64  TEST_ONE_ORDER(QSIMPSON, THIRD, 3);
65  TEST_ONE_ORDER(QTRAP, FIRST, 1);
66  TEST_ALL_ORDERS(QGRID, 1);
67 
68  // The TEST_ALL_ORDERS macro only goes up to 9th-order
69  TEST_ALL_ORDERS(QGAUSS_LOBATTO, 9);
70 
71  // The Gauss-Lobatto quadrature rules passed all these tests during
72  // development, but we don't run them with every 'make check'
73  // because it takes too long.
74  // TEST_ONE_ORDER(QGAUSS_LOBATTO, ELEVENTH, 11);
75  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTEENTH, 13);
76  // TEST_ONE_ORDER(QGAUSS_LOBATTO, FIFTEENTH, 15);
77  // TEST_ONE_ORDER(QGAUSS_LOBATTO, SEVENTEENTH, 17);
78  // TEST_ONE_ORDER(QGAUSS_LOBATTO, NINETEENTH, 19);
79  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYFIRST, 21);
80  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYTHIRD, 23);
81  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYFIFTH, 25);
82  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYSEVENTH, 27);
83  // TEST_ONE_ORDER(QGAUSS_LOBATTO, TWENTYNINTH, 29);
84  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYFIRST, 31);
85  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYTHIRD, 33);
86  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYFIFTH, 35);
87  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYSEVENTH, 37);
88  // TEST_ONE_ORDER(QGAUSS_LOBATTO, THIRTYNINTH, 39);
89  // TEST_ONE_ORDER(QGAUSS_LOBATTO, FORTYFIRST, 41);
90  // TEST_ONE_ORDER(QGAUSS_LOBATTO, FORTYTHIRD, 43);
91 
92  // Test monomial quadrature rules on quads and hexes
93  CPPUNIT_TEST( testMonomialQuadrature );
94 
95  // Test quadrature rules on Triangles
96  CPPUNIT_TEST( testTriQuadrature );
97 
98  // Test quadrature rules on Tetrahedra
99  CPPUNIT_TEST( testTetQuadrature );
100 
101  // Test Jacobi quadrature rules with special weighting function
102  CPPUNIT_TEST( testJacobi );
103 
104  CPPUNIT_TEST_SUITE_END();
105 
106 private:
107 
109 
110 
111 public:
112  void setUp ()
113  { quadrature_tolerance = TOLERANCE * std::sqrt(TOLERANCE); }
114 
115  void tearDown ()
116  {}
117 
119  {
120  ElemType elem_type[2] = {QUAD4, HEX8};
121  int dims[2] = {2, 3};
122 
123  for (int i=0; i<2; ++i)
124  {
125  // std::cout << "\nChecking monomial quadrature on element type " << elem_type[i] << std::endl;
126 
127  for (int order=0; order<7; ++order)
128  {
130  dims[i],
131  static_cast<Order>(order));
132  qrule->init(elem_type[i]);
133 
134  // In 3D, max(z_power)==order, in 2D max(z_power)==0
135  int max_z_power = dims[i]==2 ? 0 : order;
136 
137  for (int x_power=0; x_power<=order; ++x_power)
138  for (int y_power=0; y_power<=order; ++y_power)
139  for (int z_power=0; z_power<=max_z_power; ++z_power)
140  {
141  // Only try to integrate polynomials we can integrate exactly
142  if (x_power + y_power + z_power > order)
143  continue;
144 
145  // Compute the integral via quadrature. Note that
146  // std::pow(0,0) returns 1 in the 2D case.
147  Real sumq = 0.;
148  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
149  sumq += qrule->w(qp)
150  * std::pow(qrule->qp(qp)(0), x_power)
151  * std::pow(qrule->qp(qp)(1), y_power)
152  * std::pow(qrule->qp(qp)(2), z_power);
153 
154  // std::cout << "Quadrature of x^" << x_power
155  // << " y^" << y_power
156  // << " z^" << z_power
157  // << " = " << sumq << std::endl;
158 
159  // Copy-pasted code from test3DWeights()
160  Real exact_x = (x_power % 2) ? 0 : (Real(2.0) / (x_power+1));
161  Real exact_y = (y_power % 2) ? 0 : (Real(2.0) / (y_power+1));
162  Real exact_z = (z_power % 2) ? 0 : (Real(2.0) / (z_power+1));
163 
164  // Handle 2D
165  if (dims[i]==2)
166  exact_z = 1.0;
167 
168  Real exact = exact_x*exact_y*exact_z;
169 
170  // Make sure that the quadrature solution matches the exact solution
171  LIBMESH_ASSERT_REALS_EQUAL(exact, sumq, quadrature_tolerance);
172  }
173  } // end for (order)
174  } // end for (i)
175  }
176 
178  {
179  // There are 3 different families of quadrature rules for tetrahedra
181 
182  for (int qt=0; qt<3; ++qt)
183  for (int order=0; order<7; ++order)
184  {
185  UniquePtr<QBase> qrule = QBase::build(qtype[qt],
186  /*dim=*/3,
187  static_cast<Order>(order));
188 
189  // Initialize on a TET element
190  qrule->init (TET4);
191 
192  // Test the sum of the weights for this order
193  Real sumw = 0.;
194  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
195  sumw += qrule->w(qp);
196 
197  // Make sure that the weights add up to the value we expect
198  LIBMESH_ASSERT_REALS_EQUAL(1./6., sumw, quadrature_tolerance);
199 
200  // Test integrating different polynomial powers
201  for (int x_power=0; x_power<=order; ++x_power)
202  for (int y_power=0; y_power<=order; ++y_power)
203  for (int z_power=0; z_power<=order; ++z_power)
204  {
205  // Only try to integrate polynomials we can integrate exactly
206  if (x_power + y_power + z_power > order)
207  continue;
208 
209  // Compute the integral via quadrature
210  Real sumq = 0.;
211  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
212  sumq += qrule->w(qp)
213  * std::pow(qrule->qp(qp)(0), x_power)
214  * std::pow(qrule->qp(qp)(1), y_power)
215  * std::pow(qrule->qp(qp)(2), z_power);
216 
217  // std::cout << "sumq = " << sumq << std::endl;
218 
219  // Compute the true integral, a! b! c! / (a + b + c + 3)!
220  Real analytical = 1.0;
221  {
222  // Sort the a, b, c values
223  int sorted_powers[3] = {x_power, y_power, z_power};
224  std::sort(sorted_powers, sorted_powers+3);
225 
226  // Cancel the largest power with the denominator, fill in the
227  // entries for the remaining numerator terms and the denominator.
228  std::vector<int>
229  numerator_1(sorted_powers[0] > 1 ? sorted_powers[0]-1 : 0),
230  numerator_2(sorted_powers[1] > 1 ? sorted_powers[1]-1 : 0),
231  denominator(3 + sorted_powers[0] + sorted_powers[1]);
232 
233  // Fill up the vectors with sequences starting at the right values.
234  Utility::iota(numerator_1.begin(), numerator_1.end(), 2);
235  Utility::iota(numerator_2.begin(), numerator_2.end(), 2);
236  Utility::iota(denominator.begin(), denominator.end(), sorted_powers[2]+1);
237 
238  // The denominator is guaranteed to have the most terms...
239  for (std::size_t i=0; i<denominator.size(); ++i)
240  {
241  if (i < numerator_1.size())
242  analytical *= numerator_1[i];
243 
244  if (i < numerator_2.size())
245  analytical *= numerator_2[i];
246 
247  analytical /= denominator[i];
248  }
249  }
250 
251  // std::cout << "analytical = " << analytical << std::endl;
252 
253  // Make sure that the computed integral agrees with the "true" value
254  LIBMESH_ASSERT_REALS_EQUAL(analytical, sumq, quadrature_tolerance);
255  } // end for(testpower)
256  } // end for(order)
257  }
258 
260  {
262 
263  for (int qt=0; qt<4; ++qt)
264  for (int order=0; order<10; ++order)
265  {
266  UniquePtr<QBase> qrule = QBase::build(qtype[qt],
267  /*dim=*/2,
268  static_cast<Order>(order));
269 
270  // Initialize on a TRI element
271  qrule->init (TRI3);
272 
273  // Test the sum of the weights for this order
274  Real sumw = 0.;
275  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
276  sumw += qrule->w(qp);
277 
278  // Make sure that the weights add up to the value we expect
279  LIBMESH_ASSERT_REALS_EQUAL(0.5, sumw, quadrature_tolerance);
280 
281  // Test integrating different polynomial powers
282  for (int x_power=0; x_power<=order; ++x_power)
283  for (int y_power=0; y_power<=order; ++y_power)
284  {
285  // Only try to integrate polynomials we can integrate exactly
286  if (x_power + y_power > order)
287  continue;
288 
289  // Compute the integral via quadrature
290  Real sumq = 0.;
291  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
292  sumq += qrule->w(qp)
293  * std::pow(qrule->qp(qp)(0), x_power)
294  * std::pow(qrule->qp(qp)(1), y_power);
295 
296  // std::cout << "sumq = " << sumq << std::endl;
297 
298  // Compute the true integral, a! b! / (a + b + 2)!
299  Real analytical = 1.0;
300  {
301  unsigned
302  larger_power = std::max(x_power, y_power),
303  smaller_power = std::min(x_power, y_power);
304 
305  // Cancel the larger of the two numerator terms with the
306  // denominator, and fill in the remaining entries.
307  std::vector<unsigned>
308  numerator(smaller_power > 1 ? smaller_power-1 : 0),
309  denominator(2+smaller_power);
310 
311  // Fill up the vectors with sequences starting at the right values.
312  Utility::iota(numerator.begin(), numerator.end(), 2);
313  Utility::iota(denominator.begin(), denominator.end(), larger_power+1);
314 
315  // The denominator is guaranteed to have more terms...
316  for (std::size_t i=0; i<denominator.size(); ++i)
317  {
318  if (i < numerator.size())
319  analytical *= numerator[i];
320  analytical /= denominator[i];
321  }
322  }
323 
324  // std::cout << "analytical = " << analytical << std::endl;
325 
326  // Make sure that the computed integral agrees with the "true" value
327  LIBMESH_ASSERT_REALS_EQUAL(analytical, sumq, quadrature_tolerance);
328  } // end for(testpower)
329  } // end for(order)
330  }
331 
332  void testJacobi ()
333  {
334  // LibMesh supports two different types of Jacobi quadrature
335  QuadratureType qtype[2] = {QJACOBI_1_0, QJACOBI_2_0};
336 
337  // The weights of the Jacobi quadrature rules in libmesh have been
338  // scaled based on their intended use:
339  // (alpha=1, beta=0) rule weights sum to 1/2.
340  // (alpha=2, beta=0) rule weights sum to 1/3.
341  Real initial_sum_weights[2] = {.5, 1./3.};
342 
343  // The points of the Jacobi rules in LibMesh are also defined on
344  // [0,1]... this has to be taken into account when computing the
345  // exact integral values in Maple! Also note: if you scale the
346  // points to a different interval, you need to also compute what
347  // the sum of the weights should be for that interval, it will not
348  // simply be the element length for weighted quadrature rules like
349  // Jacobi. For general alpha and beta=0, the sum of the weights
350  // on the interval [-1,1] is 2^(alpha+1) / (alpha+1).
351  std::vector<std::vector<Real>> true_integrals(2);
352 
353  // alpha=1 integral values
354  // int((1-x)*x^p, x=0..1) = 1 / (p^2 + 3p + 2)
355  true_integrals[0].resize(10);
356  for (std::size_t p=0; p<true_integrals[0].size(); ++p)
357  true_integrals[0][p] = 1. / (p*p + 3.*p + 2.);
358 
359  // alpha=2 integral values
360  // int((1-x)^2*x^p, x=0..1) = 2 / (p^3 + 6*p^2 + 11*p + 6)
361  true_integrals[1].resize(10);
362  for (std::size_t p=0; p<true_integrals[1].size(); ++p)
363  true_integrals[1][p] = 2. / (p*p*p + 6.*p*p + 11.*p + 6.);
364 
365  // Test both types of Jacobi quadrature rules
366  for (int qt=0; qt<2; ++qt)
367  {
368  for (int order=0; order<10; ++order)
369  {
370  UniquePtr<QBase> qrule = QBase::build(qtype[qt],
371  /*dim=*/1,
372  static_cast<Order>(order));
373 
374  // Initialize on a 1D element, EDGE2/3/4 should not matter...
375  qrule->init (EDGE2);
376 
377  // Test the sum of the weights for this order
378  Real sumw = 0.;
379  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
380  sumw += qrule->w(qp);
381 
382  // Make sure that the weights add up to the value we expect
383  LIBMESH_ASSERT_REALS_EQUAL(initial_sum_weights[qt], sumw, quadrature_tolerance);
384 
385  // Test integrating different polynomial powers
386  for (int testpower=0; testpower<=order; ++testpower)
387  {
388  // Note that the weighting function, (1-x)^alpha *
389  // (1+x)^beta, is built into these quadrature rules;
390  // the polynomials we actually integrate are just the
391  // usual monomials.
392  Real sumq = 0.;
393  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
394  sumq += qrule->w(qp) * std::pow(qrule->qp(qp)(0), testpower);
395 
396  // Make sure that the computed integral agrees with the "true" value
397  LIBMESH_ASSERT_REALS_EQUAL(true_integrals[qt][testpower], sumq, quadrature_tolerance);
398  } // end for(testpower)
399  } // end for(order)
400  } // end for(qt)
401  } // testJacobi
402 
403 
404 
405  template <QuadratureType qtype, Order order>
406  void testBuild ()
407  {
408  UniquePtr<QBase> qrule1D = QBase::build (qtype, 1, order);
409  UniquePtr<QBase> qrule2D = QBase::build (qtype, 2, order);
410  UniquePtr<QBase> qrule3D = QBase::build (qtype, 3, order);
411 
412  CPPUNIT_ASSERT_EQUAL ( static_cast<unsigned int>(1) , qrule1D->get_dim() );
413  CPPUNIT_ASSERT_EQUAL ( static_cast<unsigned int>(2) , qrule2D->get_dim() );
414  CPPUNIT_ASSERT_EQUAL ( static_cast<unsigned int>(3) , qrule3D->get_dim() );
415 
416  // Test the enum-to-string conversion for this qtype is
417  // implemented, but not what the actual value is.
419  }
420 
421 
422 
423  //-------------------------------------------------------
424  // 1D Quadrature Rule Test
425  template <QuadratureType qtype, Order order, unsigned int exactorder>
427  {
428  UniquePtr<QBase> qrule = QBase::build(qtype , 1, order);
429  qrule->init (EDGE3);
430 
431  for (unsigned int mode=0; mode <= exactorder; ++mode)
432  {
433  Real sum = 0;
434 
435  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
436  sum += qrule->w(qp) * std::pow(qrule->qp(qp)(0), static_cast<Real>(mode));
437 
438  const Real exact = (mode % 2) ?
439  0 : (Real(2.0) / (mode+1));
440 
441  if (std::abs(exact - sum) >= quadrature_tolerance)
442  {
443  std::cout << "qtype = " << qtype << std::endl;
444  std::cout << "order = " << order << std::endl;
445  std::cout << "exactorder = " << exactorder << std::endl;
446  std::cout << "mode = " << mode << std::endl;
447  std::cout << "exact = " << exact << std::endl;
448  std::cout << "sum = " << sum << std::endl << std::endl;
449  }
450 
451  LIBMESH_ASSERT_REALS_EQUAL( exact , sum , quadrature_tolerance );
452  }
453  }
454 
455 
456 
457  //-------------------------------------------------------
458  // 2D Quadrature Rule Test
459  template <QuadratureType qtype, Order order, unsigned int exactorder>
461  {
462  UniquePtr<QBase> qrule = QBase::build(qtype, 2, order);
463  qrule->init (QUAD8);
464 
465  for (unsigned int modex=0; modex <= exactorder; ++modex)
466  for (unsigned int modey=0; modey <= exactorder; ++modey)
467  {
468  Real sum = 0;
469 
470  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
471  sum += qrule->w(qp)
472  * std::pow(qrule->qp(qp)(0), static_cast<Real>(modex))
473  * std::pow(qrule->qp(qp)(1), static_cast<Real>(modey));
474 
475  const Real exactx = (modex % 2) ?
476  0 : (Real(2.0) / (modex+1));
477 
478  const Real exacty = (modey % 2) ?
479  0 : (Real(2.0) / (modey+1));
480 
481  const Real exact = exactx*exacty;
482 
483  LIBMESH_ASSERT_REALS_EQUAL( exact , sum , quadrature_tolerance );
484  }
485 
486  // We may eventually support Gauss-Lobatto type quadrature on triangles...
487  if (qtype != QGAUSS_LOBATTO)
488  {
489  qrule->init (TRI6);
490 
491  Real sum = 0;
492 
493  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
494  sum += qrule->w(qp);
495 
496  LIBMESH_ASSERT_REALS_EQUAL( 0.5 , sum , quadrature_tolerance );
497  }
498  }
499 
500 
501 
502  //-------------------------------------------------------
503  // 3D Gauss Rule Test
504  template <QuadratureType qtype, Order order, unsigned int exactorder>
506  {
507  UniquePtr<QBase> qrule = QBase::build(qtype, 3, order);
508  qrule->init (HEX20);
509 
510  for (unsigned int modex=0; modex <= exactorder; ++modex)
511  for (unsigned int modey=0; modey <= exactorder; ++modey)
512  for (unsigned int modez=0; modez <= exactorder; ++modez)
513  {
514  Real sum = 0;
515 
516  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
517  sum += qrule->w(qp)
518  * std::pow(qrule->qp(qp)(0), static_cast<Real>(modex))
519  * std::pow(qrule->qp(qp)(1), static_cast<Real>(modey))
520  * std::pow(qrule->qp(qp)(2), static_cast<Real>(modez));
521 
522  const Real exactx = (modex % 2) ?
523  0 : (Real(2.0) / (modex+1));
524 
525  const Real exacty = (modey % 2) ?
526  0 : (Real(2.0) / (modey+1));
527 
528  const Real exactz = (modez % 2) ?
529  0 : (Real(2.0) / (modez+1));
530 
531  const Real exact = exactx*exacty*exactz;
532 
533  LIBMESH_ASSERT_REALS_EQUAL( exact , sum , quadrature_tolerance );
534  }
535 
536  // We may eventually support Gauss-Lobatto type quadrature on tets and prisms...
537  if (qtype != QGAUSS_LOBATTO)
538  {
539  qrule->init (TET10);
540 
541  Real sum = 0;
542 
543  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
544  sum += qrule->w(qp);
545 
546  LIBMESH_ASSERT_REALS_EQUAL( 1./6., sum , quadrature_tolerance );
547 
548  qrule->init (PRISM15);
549 
550  sum = 0;
551 
552  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
553  sum += qrule->w(qp);
554 
555  LIBMESH_ASSERT_REALS_EQUAL( 1., sum , quadrature_tolerance );
556  }
557  }
558 };
559 
double abs(double a)
QuadratureType
Defines an enum for currently available quadrature rules.
void testTetQuadrature()
ElemType
Defines an enum for geometric element types.
static const Real TOLERANCE
The libMesh namespace provides an interface to certain functionality in the library.
void sum(T &r, const Communicator &comm=Communicator_World)
long double max(long double a, double b)
void iota(ForwardIter first, ForwardIter last, T value)
Utility::iota is a duplication of the SGI STL extension std::iota.
Definition: utility.h:57
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
void testMonomialQuadrature()
double pow(double a, int b)
std::string enum_to_string(const T e)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
CPPUNIT_TEST_SUITE_REGISTRATION(QuadratureTest)
static UniquePtr< QBase > build(const std::string &name, const unsigned int dim, const Order order=INVALID_ORDER)
Builds a specific quadrature rule based on the name string.
void testTriQuadrature()
long double min(long double a, double b)