libMesh
Public Member Functions | Static Public Member Functions | Public Attributes | Protected Types | Protected Member Functions | Protected Attributes | Static Protected Attributes | Private Member Functions | List of all members
libMesh::QGauss Class Reference

This class implements specific orders of Gauss quadrature. More...

#include <quadrature_gauss.h>

Inheritance diagram for libMesh::QGauss:
[legend]

Public Member Functions

 QGauss (unsigned int dim, Order order=INVALID_ORDER)
 Constructor. More...
 
 QGauss (const QGauss &)=default
 Copy/move ctor, copy/move assignment operator, and destructor are all explicitly defaulted for this simple class. More...
 
 QGauss (QGauss &&)=default
 
QGauss & operator= (const QGauss &)=default
 
QGauss & operator= (QGauss &&)=default
 
virtual ~QGauss ()=default
 
virtual QuadratureType type () const override
 
ElemType get_elem_type () const
 
unsigned int get_p_level () const
 
unsigned int n_points () const
 
unsigned int get_dim () const
 
const std::vector< Point > & get_points () const
 
std::vector< Point > & get_points ()
 
const std::vector< Real > & get_weights () const
 
std::vector< Real > & get_weights ()
 
Point qp (const unsigned int i) const
 
Real w (const unsigned int i) const
 
virtual void init (const ElemType type=INVALID_ELEM, unsigned int p_level=0)
 Initializes the data structures for a quadrature rule for an element of type type. More...
 
virtual void init (const Elem &elem, const std::vector< Real > &vertex_distance_func, unsigned int p_level=0)
 Initializes the data structures for an element potentially "cut" by a signed distance function. More...
 
Order get_order () const
 
void print_info (std::ostream &os=libMesh::out) const
 Prints information relevant to the quadrature rule, by default to libMesh::out. More...
 
void scale (std::pair< Real, Real > old_range, std::pair< Real, Real > new_range)
 Maps the points of a 1D quadrature rule defined by "old_range" to another 1D interval defined by "new_range" and scales the weights accordingly. More...
 
virtual bool shapes_need_reinit ()
 

Static Public Member Functions

static std::unique_ptr< 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. More...
 
static std::unique_ptr< QBase > build (const QuadratureType qt, const unsigned int dim, const Order order=INVALID_ORDER)
 Builds a specific quadrature rule based on the QuadratureType. More...
 
static void print_info (std::ostream &out=libMesh::out)
 Prints the reference information, by default to libMesh::out. More...
 
static std::string get_info ()
 Gets a string containing the reference information. More...
 
static unsigned int n_objects ()
 Prints the number of outstanding (created, but not yet destroyed) objects. More...
 
static void enable_print_counter_info ()
 Methods to enable/disable the reference counter output from print_info() More...
 
static void disable_print_counter_info ()
 

Public Attributes

bool allow_rules_with_negative_weights
 Flag (default true) controlling the use of quadrature rules with negative weights. More...
 

Protected Types

typedef std::map< std::string, std::pair< unsigned int, unsigned int > > Counts
 Data structure to log the information. More...
 

Protected Member Functions

virtual void init_0D (const ElemType type=INVALID_ELEM, unsigned int p_level=0)
 Initializes the 0D quadrature rule by filling the points and weights vectors with the appropriate values. More...
 
void tensor_product_quad (const QBase &q1D)
 Constructs a 2D rule from the tensor product of q1D with itself. More...
 
void tensor_product_hex (const QBase &q1D)
 Computes the tensor product quadrature rule [q1D x q1D x q1D] from the 1D rule q1D. More...
 
void tensor_product_prism (const QBase &q1D, const QBase &q2D)
 Computes the tensor product of a 1D quadrature rule and a 2D quadrature rule. More...
 
void increment_constructor_count (const std::string &name)
 Increments the construction counter. More...
 
void increment_destructor_count (const std::string &name)
 Increments the destruction counter. More...
 

Protected Attributes

unsigned int _dim
 The spatial dimension of the quadrature rule. More...
 
Order _order
 The polynomial order which the quadrature rule is capable of integrating exactly. More...
 
ElemType _type
 The type of element for which the current values have been computed. More...
 
unsigned int _p_level
 The p-level of the element for which the current values have been computed. More...
 
std::vector< Point > _points
 The locations of the quadrature points in reference element space. More...
 
std::vector< Real > _weights
 The quadrature weights. More...
 

Static Protected Attributes

static Counts _counts
 Actually holds the data. More...
 
static Threads::atomic< unsigned int > _n_objects
 The number of objects. More...
 
static Threads::spin_mutex _mutex
 Mutual exclusion object to enable thread-safe reference counting. More...
 
static bool _enable_print_counter = true
 Flag to control whether reference count information is printed when print_info is called. More...
 

Private Member Functions

virtual void init_1D (const ElemType, unsigned int) override
 Initializes the 1D quadrature rule by filling the points and weights vectors with the appropriate values. More...
 
virtual void init_2D (const ElemType, unsigned int) override
 Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate values. More...
 
virtual void init_3D (const ElemType, unsigned int) override
 Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate values. More...
 
void dunavant_rule (const Real rule_data[][4], const unsigned int n_pts)
 The Dunavant rules are for triangles. More...
 
void dunavant_rule2 (const Real *wts, const Real *a, const Real *b, const unsigned int *permutation_ids, const unsigned int n_wts)
 
void keast_rule (const Real rule_data[][4], const unsigned int n_pts)
 The Keast rules are for tets. More...
 

Detailed Description

This class implements specific orders of Gauss quadrature.

The rules of Order p are capable of integrating polynomials of degree p exactly.

Author
Benjamin Kirk
John W. Peterson
Date
2003

Implements 1, 2, and 3D "Gaussian" quadrature rules.

Definition at line 39 of file quadrature_gauss.h.

Member Typedef Documentation

◆ Counts

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > libMesh::ReferenceCounter::Counts
protectedinherited

Data structure to log the information.

The log is identified by the class name.

Definition at line 117 of file reference_counter.h.

Constructor & Destructor Documentation

◆ QGauss() [1/3]

libMesh::QGauss::QGauss ( unsigned int  dim,
Order  order = INVALID_ORDER 
)
inline

Constructor.

Declares the order of the quadrature rule. We explicitly call the init function in 1D since the other tensor-product rules require this one.

Note
The element type, EDGE2, will not be used internally, however if we called the function with INVALID_ELEM it would try to be smart and return, thinking it had already done the work.

Definition at line 52 of file quadrature_gauss.h.

53  :
54  QBase(dim, order)
55  {
56  if (dim == 1)
57  init(EDGE2);
58  }

References dim, libMesh::EDGE2, and libMesh::QBase::init().

◆ QGauss() [2/3]

libMesh::QGauss::QGauss ( const QGauss &  )
default

Copy/move ctor, copy/move assignment operator, and destructor are all explicitly defaulted for this simple class.

◆ QGauss() [3/3]

libMesh::QGauss::QGauss ( QGauss &&  )
default

◆ ~QGauss()

virtual libMesh::QGauss::~QGauss ( )
virtualdefault

Member Function Documentation

◆ build() [1/2]

std::unique_ptr< QBase > libMesh::QBase::build ( const QuadratureType  qt,
const unsigned int  dim,
const Order  order = INVALID_ORDER 
)
staticinherited

Builds a specific quadrature rule based on the QuadratureType.

This enables selection of the quadrature rule at run-time.

This function allocates memory, therefore a std::unique_ptr<QBase> is returned so that the user does not accidentally leak it.

Definition at line 52 of file quadrature_build.C.

55 {
56  switch (_qt)
57  {
58 
59  case QCLOUGH:
60  {
61 #ifdef DEBUG
62  if (_order > TWENTYTHIRD)
63  {
64  libMesh::out << "WARNING: Clough quadrature implemented" << std::endl
65  << " up to TWENTYTHIRD order." << std::endl;
66  }
67 #endif
68 
69  return libmesh_make_unique<QClough>(_dim, _order);
70  }
71 
72  case QGAUSS:
73  {
74 
75 #ifdef DEBUG
76  if (_order > FORTYTHIRD)
77  {
78  libMesh::out << "WARNING: Gauss quadrature implemented" << std::endl
79  << " up to FORTYTHIRD order." << std::endl;
80  }
81 #endif
82 
83  return libmesh_make_unique<QGauss>(_dim, _order);
84  }
85 
86  case QJACOBI_1_0:
87  {
88 
89 #ifdef DEBUG
90  if (_order > FORTYTHIRD)
91  {
92  libMesh::out << "WARNING: Jacobi(1,0) quadrature implemented" << std::endl
93  << " up to FORTYTHIRD order." << std::endl;
94  }
95 
96  if (_dim > 1)
97  {
98  libMesh::out << "WARNING: Jacobi(1,0) quadrature implemented" << std::endl
99  << " in 1D only." << std::endl;
100  }
101 #endif
102 
103  return libmesh_make_unique<QJacobi>(_dim, _order, 1, 0);
104  }
105 
106  case QJACOBI_2_0:
107  {
108 
109 #ifdef DEBUG
110  if (_order > FORTYTHIRD)
111  {
112  libMesh::out << "WARNING: Jacobi(2,0) quadrature implemented" << std::endl
113  << " up to FORTYTHIRD order." << std::endl;
114  }
115 
116  if (_dim > 1)
117  {
118  libMesh::out << "WARNING: Jacobi(2,0) quadrature implemented" << std::endl
119  << " in 1D only." << std::endl;
120  }
121 #endif
122 
123  return libmesh_make_unique<QJacobi>(_dim, _order, 2, 0);
124  }
125 
126  case QSIMPSON:
127  {
128 
129 #ifdef DEBUG
130  if (_order > THIRD)
131  {
132  libMesh::out << "WARNING: Simpson rule provides only" << std::endl
133  << " THIRD order!" << std::endl;
134  }
135 #endif
136 
137  return libmesh_make_unique<QSimpson>(_dim);
138  }
139 
140  case QTRAP:
141  {
142 
143 #ifdef DEBUG
144  if (_order > FIRST)
145  {
146  libMesh::out << "WARNING: Trapezoidal rule provides only" << std::endl
147  << " FIRST order!" << std::endl;
148  }
149 #endif
150 
151  return libmesh_make_unique<QTrap>(_dim);
152  }
153 
154  case QGRID:
155  return libmesh_make_unique<QGrid>(_dim, _order);
156 
157  case QGRUNDMANN_MOLLER:
158  return libmesh_make_unique<QGrundmann_Moller>(_dim, _order);
159 
160  case QMONOMIAL:
161  return libmesh_make_unique<QMonomial>(_dim, _order);
162 
163  case QGAUSS_LOBATTO:
164  return libmesh_make_unique<QGaussLobatto>(_dim, _order);
165 
166  case QCONICAL:
167  return libmesh_make_unique<QConical>(_dim, _order);
168 
169  case QNODAL:
170  return libmesh_make_unique<QNodal>(_dim, _order);
171 
172  default:
173  libmesh_error_msg("ERROR: Bad qt=" << _qt);
174  }
175 }

References libMesh::QBase::_dim, libMesh::QBase::_order, libMesh::FIRST, libMesh::FORTYTHIRD, libMesh::out, libMesh::QCLOUGH, libMesh::QCONICAL, libMesh::QGAUSS, libMesh::QGAUSS_LOBATTO, libMesh::QGRID, libMesh::QGRUNDMANN_MOLLER, libMesh::QJACOBI_1_0, libMesh::QJACOBI_2_0, libMesh::QMONOMIAL, libMesh::QNODAL, libMesh::QSIMPSON, libMesh::QTRAP, libMesh::THIRD, and libMesh::TWENTYTHIRD.

◆ build() [2/2]

std::unique_ptr< QBase > libMesh::QBase::build ( const std::string &  name,
const unsigned int  dim,
const Order  order = INVALID_ORDER 
)
staticinherited

Builds a specific quadrature rule based on the name string.

This enables selection of the quadrature rule at run-time. The input parameter name must be mappable through the Utility::string_to_enum<>() function.

This function allocates memory, therefore a std::unique_ptr<QBase> is returned so that the user does not accidentally leak it.

Definition at line 41 of file quadrature_build.C.

44 {
45  return QBase::build (Utility::string_to_enum<QuadratureType> (type),
46  _dim,
47  _order);
48 }

References libMesh::QBase::_dim, libMesh::QBase::_order, and libMesh::QBase::type().

Referenced by assemble_poisson(), libMesh::InfFE< Dim, T_radial, T_map >::attach_quadrature_rule(), libMesh::InfFE< Dim, T_radial, T_map >::reinit(), QuadratureTest::test1DWeights(), QuadratureTest::test2DWeights(), QuadratureTest::test3DWeights(), QuadratureTest::testBuild(), QuadratureTest::testJacobi(), QuadratureTest::testMonomialQuadrature(), QuadratureTest::testTetQuadrature(), and QuadratureTest::testTriQuadrature().

◆ disable_print_counter_info()

void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited

Definition at line 106 of file reference_counter.C.

107 {
108  _enable_print_counter = false;
109  return;
110 }

References libMesh::ReferenceCounter::_enable_print_counter.

Referenced by libMesh::LibMeshInit::LibMeshInit().

◆ dunavant_rule()

void libMesh::QGauss::dunavant_rule ( const Real  rule_data[][4],
const unsigned int  n_pts 
)
private

The Dunavant rules are for triangles.

This function takes permutation points and weights in a specific format as input and fills the _points and _weights vectors.

Definition at line 265 of file quadrature_gauss.C.

267 {
268  // The input data array has 4 columns. The first 3 are the permutation points.
269  // The last column is the weights for a given set of permutation points. A zero
270  // in two of the first 3 columns implies the point is a 1-permutation (centroid).
271  // A zero in one of the first 3 columns implies the point is a 3-permutation.
272  // Otherwise each point is assumed to be a 6-permutation.
273 
274  // Always insert into the points & weights vector relative to the offset
275  unsigned int offset=0;
276 
277 
278  for (unsigned int p=0; p<n_pts; ++p)
279  {
280 
281  // There must always be a non-zero entry to start the row
282  libmesh_assert_not_equal_to ( rule_data[p][0], static_cast<Real>(0.0) );
283 
284  // A zero weight may imply you did not set up the raw data correctly
285  libmesh_assert_not_equal_to ( rule_data[p][3], static_cast<Real>(0.0) );
286 
287  // What kind of point is this?
288  // One non-zero entry in first 3 cols ? 1-perm (centroid) point = 1
289  // Two non-zero entries in first 3 cols ? 3-perm point = 3
290  // Three non-zero entries ? 6-perm point = 6
291  unsigned int pointtype=1;
292 
293  if (rule_data[p][1] != static_cast<Real>(0.0))
294  {
295  if (rule_data[p][2] != static_cast<Real>(0.0))
296  pointtype = 6;
297  else
298  pointtype = 3;
299  }
300 
301  switch (pointtype)
302  {
303  case 1:
304  {
305  // Be sure we have enough space to insert this point
306  libmesh_assert_less (offset + 0, _points.size());
307 
308  // The point has only a single permutation (the centroid!)
309  _points[offset + 0] = Point(rule_data[p][0], rule_data[p][0]);
310 
311  // The weight is always the last entry in the row.
312  _weights[offset + 0] = rule_data[p][3];
313 
314  offset += 1;
315  break;
316  }
317 
318  case 3:
319  {
320  // Be sure we have enough space to insert these points
321  libmesh_assert_less (offset + 2, _points.size());
322 
323  // Here it's understood the second entry is to be used twice, and
324  // thus there are three possible permutations.
325  _points[offset + 0] = Point(rule_data[p][0], rule_data[p][1]);
326  _points[offset + 1] = Point(rule_data[p][1], rule_data[p][0]);
327  _points[offset + 2] = Point(rule_data[p][1], rule_data[p][1]);
328 
329  for (unsigned int j=0; j<3; ++j)
330  _weights[offset + j] = rule_data[p][3];
331 
332  offset += 3;
333  break;
334  }
335 
336  case 6:
337  {
338  // Be sure we have enough space to insert these points
339  libmesh_assert_less (offset + 5, _points.size());
340 
341  // Three individual entries with six permutations.
342  _points[offset + 0] = Point(rule_data[p][0], rule_data[p][1]);
343  _points[offset + 1] = Point(rule_data[p][0], rule_data[p][2]);
344  _points[offset + 2] = Point(rule_data[p][1], rule_data[p][0]);
345  _points[offset + 3] = Point(rule_data[p][1], rule_data[p][2]);
346  _points[offset + 4] = Point(rule_data[p][2], rule_data[p][0]);
347  _points[offset + 5] = Point(rule_data[p][2], rule_data[p][1]);
348 
349  for (unsigned int j=0; j<6; ++j)
350  _weights[offset + j] = rule_data[p][3];
351 
352  offset += 6;
353  break;
354  }
355 
356  default:
357  libmesh_error_msg("Don't know what to do with this many permutation points!");
358  }
359  }
360 }

References libMesh::QBase::_points, and libMesh::QBase::_weights.

Referenced by init_2D().

◆ dunavant_rule2()

void libMesh::QGauss::dunavant_rule2 ( const Real *  wts,
const Real *  a,
const Real *  b,
const unsigned int *  permutation_ids,
const unsigned int  n_wts 
)
private

Definition at line 190 of file quadrature_gauss.C.

195 {
196  // Figure out how many total points by summing up the entries
197  // in the permutation_ids array, and resize the _points and _weights
198  // vectors appropriately.
199  unsigned int total_pts = 0;
200  for (unsigned int p=0; p<n_wts; ++p)
201  total_pts += permutation_ids[p];
202 
203  // Resize point and weight vectors appropriately.
204  _points.resize(total_pts);
205  _weights.resize(total_pts);
206 
207  // Always insert into the points & weights vector relative to the offset
208  unsigned int offset=0;
209 
210  for (unsigned int p=0; p<n_wts; ++p)
211  {
212  switch (permutation_ids[p])
213  {
214  case 1:
215  {
216  // The point has only a single permutation (the centroid!)
217  // So we don't even need to look in the a or b arrays.
218  _points [offset + 0] = Point(Real(1)/3, Real(1)/3);
219  _weights[offset + 0] = wts[p];
220 
221  offset += 1;
222  break;
223  }
224 
225 
226  case 3:
227  {
228  // For this type of rule, don't need to look in the b array.
229  _points[offset + 0] = Point(a[p], a[p]); // (a,a)
230  _points[offset + 1] = Point(a[p], 1-2*a[p]); // (a,1-2a)
231  _points[offset + 2] = Point(1-2*a[p], a[p]); // (1-2a,a)
232 
233  for (unsigned int j=0; j<3; ++j)
234  _weights[offset + j] = wts[p];
235 
236  offset += 3;
237  break;
238  }
239 
240  case 6:
241  {
242  // This type of point uses all 3 arrays...
243  _points[offset + 0] = Point(a[p], b[p]);
244  _points[offset + 1] = Point(b[p], a[p]);
245  _points[offset + 2] = Point(a[p], 1-a[p]-b[p]);
246  _points[offset + 3] = Point(1-a[p]-b[p], a[p]);
247  _points[offset + 4] = Point(b[p], 1-a[p]-b[p]);
248  _points[offset + 5] = Point(1-a[p]-b[p], b[p]);
249 
250  for (unsigned int j=0; j<6; ++j)
251  _weights[offset + j] = wts[p];
252 
253  offset += 6;
254  break;
255  }
256 
257  default:
258  libmesh_error_msg("Unknown permutation id: " << permutation_ids[p] << "!");
259  }
260  }
261 
262 }

References libMesh::QBase::_points, libMesh::QBase::_weights, and libMesh::Real.

Referenced by init_2D().

◆ enable_print_counter_info()

void libMesh::ReferenceCounter::enable_print_counter_info ( )
staticinherited

Methods to enable/disable the reference counter output from print_info()

Definition at line 100 of file reference_counter.C.

101 {
102  _enable_print_counter = true;
103  return;
104 }

References libMesh::ReferenceCounter::_enable_print_counter.

◆ get_dim()

unsigned int libMesh::QBase::get_dim ( ) const
inlineinherited
Returns
The spatial dimension of the quadrature rule.

Definition at line 135 of file quadrature.h.

135 { return _dim; }

References libMesh::QBase::_dim.

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::attach_quadrature_rule(), libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), and libMesh::QConical::conical_product_tri().

◆ get_elem_type()

ElemType libMesh::QBase::get_elem_type ( ) const
inlineinherited
Returns
The element type we're currently using.

Definition at line 116 of file quadrature.h.

116 { return _type; }

References libMesh::QBase::_type.

◆ get_info()

std::string libMesh::ReferenceCounter::get_info ( )
staticinherited

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

48 {
49 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
50 
51  std::ostringstream oss;
52 
53  oss << '\n'
54  << " ---------------------------------------------------------------------------- \n"
55  << "| Reference count information |\n"
56  << " ---------------------------------------------------------------------------- \n";
57 
58  for (const auto & pr : _counts)
59  {
60  const std::string name(pr.first);
61  const unsigned int creations = pr.second.first;
62  const unsigned int destructions = pr.second.second;
63 
64  oss << "| " << name << " reference count information:\n"
65  << "| Creations: " << creations << '\n'
66  << "| Destructions: " << destructions << '\n';
67  }
68 
69  oss << " ---------------------------------------------------------------------------- \n";
70 
71  return oss.str();
72 
73 #else
74 
75  return "";
76 
77 #endif
78 }

References libMesh::ReferenceCounter::_counts, and libMesh::Quality::name().

Referenced by libMesh::ReferenceCounter::print_info().

◆ get_order()

Order libMesh::QBase::get_order ( ) const
inlineinherited
Returns
The current "total" order of the quadrature rule which can vary element by element, depending on the Elem::p_level(), which gets passed to us during init().

Each additional power of p increases the quadrature order required to integrate the mass matrix by 2, hence the formula below.

Definition at line 211 of file quadrature.h.

211 { return static_cast<Order>(_order + 2 * _p_level); }

References libMesh::QBase::_order, and libMesh::QBase::_p_level.

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::attach_quadrature_rule(), libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), libMesh::QGaussLobatto::init_1D(), libMesh::QConical::init_1D(), init_1D(), libMesh::QJacobi::init_1D(), libMesh::QGrundmann_Moller::init_1D(), init_2D(), libMesh::QMonomial::init_2D(), libMesh::QGrundmann_Moller::init_2D(), init_3D(), libMesh::QMonomial::init_3D(), and libMesh::QGrundmann_Moller::init_3D().

◆ get_p_level()

unsigned int libMesh::QBase::get_p_level ( ) const
inlineinherited
Returns
The p-refinement level we're currently using.

Definition at line 121 of file quadrature.h.

121 { return _p_level; }

References libMesh::QBase::_p_level.

◆ get_points() [1/2]

std::vector<Point>& libMesh::QBase::get_points ( )
inlineinherited
Returns
A std::vector containing the quadrature point locations in reference element space as a writable reference.

Definition at line 147 of file quadrature.h.

147 { return _points; }

References libMesh::QBase::_points.

◆ get_points() [2/2]

const std::vector<Point>& libMesh::QBase::get_points ( ) const
inlineinherited
Returns
A std::vector containing the quadrature point locations in reference element space.

Definition at line 141 of file quadrature.h.

141 { return _points; }

References libMesh::QBase::_points.

Referenced by assemble_ellipticdg(), libMesh::QClough::init_1D(), libMesh::QConical::init_1D(), libMesh::QNodal::init_1D(), libMesh::QMonomial::init_1D(), libMesh::QGrundmann_Moller::init_1D(), libMesh::QClough::init_2D(), init_2D(), libMesh::QNodal::init_2D(), libMesh::QMonomial::init_2D(), init_3D(), libMesh::QNodal::init_3D(), libMesh::QMonomial::init_3D(), EIM_F::interior_assembly(), and AssemblyEIM::interior_assembly().

◆ get_weights() [1/2]

std::vector<Real>& libMesh::QBase::get_weights ( )
inlineinherited
Returns
A writable references to a std::vector containing the quadrature weights.

Definition at line 159 of file quadrature.h.

159 { return _weights; }

References libMesh::QBase::_weights.

◆ get_weights() [2/2]

const std::vector<Real>& libMesh::QBase::get_weights ( ) const
inlineinherited
Returns
A constant reference to a std::vector containing the quadrature weights.

Definition at line 153 of file quadrature.h.

153 { return _weights; }

References libMesh::QBase::_weights.

Referenced by libMesh::QClough::init_1D(), libMesh::QConical::init_1D(), libMesh::QNodal::init_1D(), libMesh::QMonomial::init_1D(), libMesh::QGrundmann_Moller::init_1D(), libMesh::QClough::init_2D(), init_2D(), libMesh::QNodal::init_2D(), libMesh::QMonomial::init_2D(), init_3D(), libMesh::QNodal::init_3D(), and libMesh::QMonomial::init_3D().

◆ increment_constructor_count()

void libMesh::ReferenceCounter::increment_constructor_count ( const std::string &  name)
inlineprotectedinherited

Increments the construction counter.

Should be called in the constructor of any derived class that will be reference counted.

Definition at line 181 of file reference_counter.h.

182 {
183  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
184  std::pair<unsigned int, unsigned int> & p = _counts[name];
185 
186  p.first++;
187 }

References libMesh::ReferenceCounter::_counts, libMesh::Quality::name(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::ReferenceCountedObject().

◆ increment_destructor_count()

void libMesh::ReferenceCounter::increment_destructor_count ( const std::string &  name)
inlineprotectedinherited

Increments the destruction counter.

Should be called in the destructor of any derived class that will be reference counted.

Definition at line 194 of file reference_counter.h.

195 {
196  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
197  std::pair<unsigned int, unsigned int> & p = _counts[name];
198 
199  p.second++;
200 }

References libMesh::ReferenceCounter::_counts, libMesh::Quality::name(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::~ReferenceCountedObject().

◆ init() [1/2]

void libMesh::QBase::init ( const Elem &  elem,
const std::vector< Real > &  vertex_distance_func,
unsigned int  p_level = 0 
)
virtualinherited

Initializes the data structures for an element potentially "cut" by a signed distance function.

The array vertex_distance_func contains vertex values of the signed distance function. If the signed distance function changes sign on the vertices, then the element is considered to be cut.) This interface can be extended by derived classes in order to subdivide the element and construct a composite quadrature rule.

Definition at line 103 of file quadrature.C.

106 {
107  // dispatch generic implementation
108  this->init(elem.type(), p_level);
109 }

References libMesh::QBase::init(), and libMesh::Elem::type().

◆ init() [2/2]

void libMesh::QBase::init ( const ElemType  type = INVALID_ELEM,
unsigned int  p_level = 0 
)
virtualinherited

Initializes the data structures for a quadrature rule for an element of type type.

Definition at line 59 of file quadrature.C.

61 {
62  // check to see if we have already
63  // done the work for this quadrature rule
64  if (t == _type && p == _p_level)
65  return;
66  else
67  {
68  _type = t;
69  _p_level = p;
70  }
71 
72 
73 
74  switch(_dim)
75  {
76  case 0:
77  this->init_0D();
78 
79  return;
80 
81  case 1:
82  this->init_1D();
83 
84  return;
85 
86  case 2:
87  this->init_2D();
88 
89  return;
90 
91  case 3:
92  this->init_3D();
93 
94  return;
95 
96  default:
97  libmesh_error_msg("Invalid dimension _dim = " << _dim);
98  }
99 }

References libMesh::QBase::_dim, libMesh::QBase::_p_level, libMesh::QBase::_type, libMesh::QBase::init_0D(), libMesh::QBase::init_1D(), libMesh::QBase::init_2D(), and libMesh::QBase::init_3D().

Referenced by libMesh::QBase::init(), libMesh::QClough::init_1D(), libMesh::QNodal::init_1D(), libMesh::QMonomial::init_1D(), libMesh::QClough::init_2D(), libMesh::QGaussLobatto::init_2D(), libMesh::QGrid::init_2D(), libMesh::QTrap::init_2D(), init_2D(), libMesh::QSimpson::init_2D(), libMesh::QNodal::init_2D(), libMesh::QMonomial::init_2D(), libMesh::QGaussLobatto::init_3D(), libMesh::QGrid::init_3D(), libMesh::QTrap::init_3D(), libMesh::QSimpson::init_3D(), init_3D(), libMesh::QNodal::init_3D(), libMesh::QMonomial::init_3D(), QGauss(), libMesh::QGaussLobatto::QGaussLobatto(), libMesh::QJacobi::QJacobi(), libMesh::QNodal::QNodal(), libMesh::QSimpson::QSimpson(), and libMesh::QTrap::QTrap().

◆ init_0D()

void libMesh::QBase::init_0D ( const ElemType  type = INVALID_ELEM,
unsigned int  p_level = 0 
)
protectedvirtualinherited

Initializes the 0D quadrature rule by filling the points and weights vectors with the appropriate values.

Generally this is just one point with weight 1.

Note
The arguments should no longer be used for anything in derived classes, they are only maintained for backwards compatibility and will eventually be removed.

Definition at line 113 of file quadrature.C.

114 {
115  _points.resize(1);
116  _weights.resize(1);
117  _points[0] = Point(0.);
118  _weights[0] = 1.0;
119 }

References libMesh::QBase::_points, and libMesh::QBase::_weights.

Referenced by libMesh::QBase::init().

◆ init_1D()

void libMesh::QGauss::init_1D ( const  type,
unsigned int  p_level 
)
overrideprivatevirtual

Initializes the 1D quadrature rule by filling the points and weights vectors with the appropriate values.

The order of the rule will be defined by the implementing class. It is assumed that derived quadrature rules will at least define the init_1D function, therefore it is pure virtual.

Note
The arguments should no longer be used for anything in derived classes, they are only maintained for backwards compatibility and will eventually be removed.

Implements libMesh::QBase.

Definition at line 30 of file quadrature_gauss_1D.C.

31 {
32  //----------------------------------------------------------------------
33  // 1D quadrature rules
34  switch(get_order())
35  {
36  case CONSTANT:
37  case FIRST:
38  {
39  _points.resize (1);
40  _weights.resize(1);
41 
42  _points[0](0) = 0.;
43 
44  _weights[0] = 2.;
45 
46  return;
47  }
48  case SECOND:
49  case THIRD:
50  {
51  _points.resize (2);
52  _weights.resize(2);
53 
54  _points[0](0) = Real(-5.7735026918962576450914878050196e-01L); // -sqrt(3)/3
55  _points[1] = -_points[0];
56 
57  _weights[0] = 1.;
58  _weights[1] = _weights[0];
59 
60  return;
61  }
62  case FOURTH:
63  case FIFTH:
64  {
65  _points.resize (3);
66  _weights.resize(3);
67 
68  _points[ 0](0) = Real(-7.7459666924148337703585307995648e-01L);
69  _points[ 1](0) = 0.;
70  _points[ 2] = -_points[0];
71 
72  _weights[ 0] = Real(5.5555555555555555555555555555556e-01L);
73  _weights[ 1] = Real(8.8888888888888888888888888888889e-01L);
74  _weights[ 2] = _weights[0];
75 
76  return;
77  }
78  case SIXTH:
79  case SEVENTH:
80  {
81  _points.resize (4);
82  _weights.resize(4);
83 
84  _points[ 0](0) = Real(-8.6113631159405257522394648889281e-01L);
85  _points[ 1](0) = Real(-3.3998104358485626480266575910324e-01L);
86  _points[ 2] = -_points[1];
87  _points[ 3] = -_points[0];
88 
89  _weights[ 0] = Real(3.4785484513745385737306394922200e-01L);
90  _weights[ 1] = Real(6.5214515486254614262693605077800e-01L);
91  _weights[ 2] = _weights[1];
92  _weights[ 3] = _weights[0];
93 
94  return;
95  }
96  case EIGHTH:
97  case NINTH:
98  {
99  _points.resize (5);
100  _weights.resize(5);
101 
102  _points[ 0](0) = Real(-9.0617984593866399279762687829939e-01L);
103  _points[ 1](0) = Real(-5.3846931010568309103631442070021e-01L);
104  _points[ 2](0) = 0.;
105  _points[ 3] = -_points[1];
106  _points[ 4] = -_points[0];
107 
108  _weights[ 0] = Real(2.3692688505618908751426404071992e-01L);
109  _weights[ 1] = Real(4.7862867049936646804129151483564e-01L);
110  _weights[ 2] = Real(5.6888888888888888888888888888889e-01L);
111  _weights[ 3] = _weights[1];
112  _weights[ 4] = _weights[0];
113 
114  return;
115  }
116  case TENTH:
117  case ELEVENTH:
118  {
119  _points.resize (6);
120  _weights.resize(6);
121 
122  _points[ 0](0) = Real(-9.3246951420315202781230155449399e-01L);
123  _points[ 1](0) = Real(-6.6120938646626451366139959501991e-01L);
124  _points[ 2](0) = Real(-2.3861918608319690863050172168071e-01L);
125  _points[ 3] = -_points[2];
126  _points[ 4] = -_points[1];
127  _points[ 5] = -_points[0];
128 
129  _weights[ 0] = Real(1.7132449237917034504029614217273e-01L);
130  _weights[ 1] = Real(3.6076157304813860756983351383772e-01L);
131  _weights[ 2] = Real(4.6791393457269104738987034398955e-01L);
132  _weights[ 3] = _weights[2];
133  _weights[ 4] = _weights[1];
134  _weights[ 5] = _weights[0];
135 
136  return;
137  }
138  case TWELFTH:
139  case THIRTEENTH:
140  {
141  _points.resize (7);
142  _weights.resize(7);
143 
144  _points[ 0](0) = Real(-9.4910791234275852452618968404785e-01L);
145  _points[ 1](0) = Real(-7.4153118559939443986386477328079e-01L);
146  _points[ 2](0) = Real(-4.0584515137739716690660641207696e-01L);
147  _points[ 3](0) = 0.;
148  _points[ 4] = -_points[2];
149  _points[ 5] = -_points[1];
150  _points[ 6] = -_points[0];
151 
152  _weights[ 0] = Real(1.2948496616886969327061143267908e-01L);
153  _weights[ 1] = Real(2.7970539148927666790146777142378e-01L);
154  _weights[ 2] = Real(3.8183005050511894495036977548898e-01L);
155  _weights[ 3] = Real(4.1795918367346938775510204081633e-01L);
156  _weights[ 4] = _weights[2];
157  _weights[ 5] = _weights[1];
158  _weights[ 6] = _weights[0];
159 
160  return;
161  }
162  case FOURTEENTH:
163  case FIFTEENTH:
164  {
165  _points.resize (8);
166  _weights.resize(8);
167 
168  _points[ 0](0) = Real(-9.6028985649753623168356086856947e-01L);
169  _points[ 1](0) = Real(-7.9666647741362673959155393647583e-01L);
170  _points[ 2](0) = Real(-5.2553240991632898581773904918925e-01L);
171  _points[ 3](0) = Real(-1.8343464249564980493947614236018e-01L);
172  _points[ 4] = -_points[3];
173  _points[ 5] = -_points[2];
174  _points[ 6] = -_points[1];
175  _points[ 7] = -_points[0];
176 
177  _weights[ 0] = Real(1.0122853629037625915253135430996e-01L);
178  _weights[ 1] = Real(2.2238103445337447054435599442624e-01L);
179  _weights[ 2] = Real(3.1370664587788728733796220198660e-01L);
180  _weights[ 3] = Real(3.6268378337836198296515044927720e-01L);
181  _weights[ 4] = _weights[3];
182  _weights[ 5] = _weights[2];
183  _weights[ 6] = _weights[1];
184  _weights[ 7] = _weights[0];
185 
186  return;
187  }
188  case SIXTEENTH:
189  case SEVENTEENTH:
190  {
191  _points.resize (9);
192  _weights.resize(9);
193 
194  _points[ 0](0) = Real(-9.6816023950762608983557620290367e-01L);
195  _points[ 1](0) = Real(-8.3603110732663579429942978806973e-01L);
196  _points[ 2](0) = Real(-6.1337143270059039730870203934147e-01L);
197  _points[ 3](0) = Real(-3.2425342340380892903853801464334e-01L);
198  _points[ 4](0) = 0.;
199  _points[ 5] = -_points[3];
200  _points[ 6] = -_points[2];
201  _points[ 7] = -_points[1];
202  _points[ 8] = -_points[0];
203 
204  _weights[ 0] = Real(8.1274388361574411971892158110524e-02L);
205  _weights[ 1] = Real(1.8064816069485740405847203124291e-01L);
206  _weights[ 2] = Real(2.6061069640293546231874286941863e-01L);
207  _weights[ 3] = Real(3.1234707704000284006863040658444e-01L);
208  _weights[ 4] = Real(3.3023935500125976316452506928697e-01L);
209  _weights[ 5] = _weights[3];
210  _weights[ 6] = _weights[2];
211  _weights[ 7] = _weights[1];
212  _weights[ 8] = _weights[0];
213 
214  return;
215  }
216  case EIGHTTEENTH:
217  case NINETEENTH:
218  {
219  _points.resize (10);
220  _weights.resize(10);
221 
222  _points[ 0](0) = Real(-9.7390652851717172007796401208445e-01L);
223  _points[ 1](0) = Real(-8.6506336668898451073209668842349e-01L);
224  _points[ 2](0) = Real(-6.7940956829902440623432736511487e-01L);
225  _points[ 3](0) = Real(-4.3339539412924719079926594316578e-01L);
226  _points[ 4](0) = Real(-1.4887433898163121088482600112972e-01L);
227  _points[ 5] = -_points[4];
228  _points[ 6] = -_points[3];
229  _points[ 7] = -_points[2];
230  _points[ 8] = -_points[1];
231  _points[ 9] = -_points[0];
232 
233  _weights[ 0] = Real(6.6671344308688137593568809893332e-02L);
234  _weights[ 1] = Real(1.4945134915058059314577633965770e-01L);
235  _weights[ 2] = Real(2.1908636251598204399553493422816e-01L);
236  _weights[ 3] = Real(2.6926671930999635509122692156947e-01L);
237  _weights[ 4] = Real(2.9552422471475287017389299465134e-01L);
238  _weights[ 5] = _weights[4];
239  _weights[ 6] = _weights[3];
240  _weights[ 7] = _weights[2];
241  _weights[ 8] = _weights[1];
242  _weights[ 9] = _weights[0];
243 
244  return;
245  }
246 
247  case TWENTIETH:
248  case TWENTYFIRST:
249  {
250  _points.resize (11);
251  _weights.resize(11);
252 
253  _points[ 0](0) = Real(-9.7822865814605699280393800112286e-01L);
254  _points[ 1](0) = Real(-8.8706259976809529907515776930393e-01L);
255  _points[ 2](0) = Real(-7.3015200557404932409341625203115e-01L);
256  _points[ 3](0) = Real(-5.1909612920681181592572566945861e-01L);
257  _points[ 4](0) = Real(-2.6954315595234497233153198540086e-01L);
258  _points[ 5](0) = 0.;
259  _points[ 6] = -_points[4];
260  _points[ 7] = -_points[3];
261  _points[ 8] = -_points[2];
262  _points[ 9] = -_points[1];
263  _points[10] = -_points[0];
264 
265  _weights[ 0] = Real(5.5668567116173666482753720442549e-02L);
266  _weights[ 1] = Real(1.2558036946490462463469429922394e-01L);
267  _weights[ 2] = Real(1.8629021092773425142609764143166e-01L);
268  _weights[ 3] = Real(2.3319376459199047991852370484318e-01L);
269  _weights[ 4] = Real(2.6280454451024666218068886989051e-01L);
270  _weights[ 5] = Real(2.7292508677790063071448352833634e-01L);
271  _weights[ 6] = _weights[4];
272  _weights[ 7] = _weights[3];
273  _weights[ 8] = _weights[2];
274  _weights[ 9] = _weights[1];
275  _weights[10] = _weights[0];
276 
277  return;
278  }
279 
280  case TWENTYSECOND:
281  case TWENTYTHIRD:
282  {
283  _points.resize (12);
284  _weights.resize(12);
285 
286  _points[ 0](0) = Real(-9.8156063424671925069054909014928e-01L);
287  _points[ 1](0) = Real(-9.0411725637047485667846586611910e-01L);
288  _points[ 2](0) = Real(-7.6990267419430468703689383321282e-01L);
289  _points[ 3](0) = Real(-5.8731795428661744729670241894053e-01L);
290  _points[ 4](0) = Real(-3.6783149899818019375269153664372e-01L);
291  _points[ 5](0) = Real(-1.2523340851146891547244136946385e-01L);
292  _points[ 6] = -_points[5];
293  _points[ 7] = -_points[4];
294  _points[ 8] = -_points[3];
295  _points[ 9] = -_points[2];
296  _points[10] = -_points[1];
297  _points[11] = -_points[0];
298 
299  _weights[ 0] = Real(4.7175336386511827194615961485017e-02L);
300  _weights[ 1] = Real(1.0693932599531843096025471819400e-01L);
301  _weights[ 2] = Real(1.6007832854334622633465252954336e-01L);
302  _weights[ 3] = Real(2.0316742672306592174906445580980e-01L);
303  _weights[ 4] = Real(2.3349253653835480876084989892483e-01L);
304  _weights[ 5] = Real(2.4914704581340278500056243604295e-01L);
305  _weights[ 6] = _weights[5];
306  _weights[ 7] = _weights[4];
307  _weights[ 8] = _weights[3];
308  _weights[ 9] = _weights[2];
309  _weights[10] = _weights[1];
310  _weights[11] = _weights[0];
311 
312  return;
313  }
314 
315  case TWENTYFOURTH:
316  case TWENTYFIFTH:
317  {
318  _points.resize (13);
319  _weights.resize(13);
320 
321  _points[ 0](0) = Real(-9.8418305471858814947282944880711e-01L);
322  _points[ 1](0) = Real(-9.1759839922297796520654783650072e-01L);
323  _points[ 2](0) = Real(-8.0157809073330991279420648958286e-01L);
324  _points[ 3](0) = Real(-6.4234933944034022064398460699552e-01L);
325  _points[ 4](0) = Real(-4.4849275103644685287791285212764e-01L);
326  _points[ 5](0) = Real(-2.3045831595513479406552812109799e-01L);
327  _points[ 6](0) = 0.;
328  _points[ 7] = -_points[5];
329  _points[ 8] = -_points[4];
330  _points[ 9] = -_points[3];
331  _points[10] = -_points[2];
332  _points[11] = -_points[1];
333  _points[12] = -_points[0];
334 
335  _weights[ 0] = Real(4.0484004765315879520021592200986e-02L);
336  _weights[ 1] = Real(9.2121499837728447914421775953797e-02L);
337  _weights[ 2] = Real(1.3887351021978723846360177686887e-01L);
338  _weights[ 3] = Real(1.7814598076194573828004669199610e-01L);
339  _weights[ 4] = Real(2.0781604753688850231252321930605e-01L);
340  _weights[ 5] = Real(2.2628318026289723841209018603978e-01L);
341  _weights[ 6] = Real(2.3255155323087391019458951526884e-01L);
342  _weights[ 7] = _weights[5];
343  _weights[ 8] = _weights[4];
344  _weights[ 9] = _weights[3];
345  _weights[10] = _weights[2];
346  _weights[11] = _weights[1];
347  _weights[12] = _weights[0];
348 
349  return;
350  }
351 
352  case TWENTYSIXTH:
353  case TWENTYSEVENTH:
354  {
355  _points.resize (14);
356  _weights.resize(14);
357 
358  _points[ 0](0) = Real(-9.8628380869681233884159726670405e-01L);
359  _points[ 1](0) = Real(-9.2843488366357351733639113937787e-01L);
360  _points[ 2](0) = Real(-8.2720131506976499318979474265039e-01L);
361  _points[ 3](0) = Real(-6.8729290481168547014801980301933e-01L);
362  _points[ 4](0) = Real(-5.1524863635815409196529071855119e-01L);
363  _points[ 5](0) = Real(-3.1911236892788976043567182416848e-01L);
364  _points[ 6](0) = Real(-1.0805494870734366206624465021983e-01L);
365  _points[ 7] = -_points[6];
366  _points[ 8] = -_points[5];
367  _points[ 9] = -_points[4];
368  _points[10] = -_points[3];
369  _points[11] = -_points[2];
370  _points[12] = -_points[1];
371  _points[13] = -_points[0];
372 
373  _weights[ 0] = Real(3.5119460331751863031832876138192e-02L);
374  _weights[ 1] = Real(8.0158087159760209805633277062854e-02L);
375  _weights[ 2] = Real(1.2151857068790318468941480907248e-01L);
376  _weights[ 3] = Real(1.5720316715819353456960193862384e-01L);
377  _weights[ 4] = Real(1.8553839747793781374171659012516e-01L);
378  _weights[ 5] = Real(2.0519846372129560396592406566122e-01L);
379  _weights[ 6] = Real(2.1526385346315779019587644331626e-01L);
380  _weights[ 7] = _weights[6];
381  _weights[ 8] = _weights[5];
382  _weights[ 9] = _weights[4];
383  _weights[10] = _weights[3];
384  _weights[11] = _weights[2];
385  _weights[12] = _weights[1];
386  _weights[13] = _weights[0];
387 
388  return;
389  }
390 
391  case TWENTYEIGHTH:
392  case TWENTYNINTH:
393  {
394  _points.resize (15);
395  _weights.resize(15);
396 
397  _points[ 0](0) = Real(-9.8799251802048542848956571858661e-01L);
398  _points[ 1](0) = Real(-9.3727339240070590430775894771021e-01L);
399  _points[ 2](0) = Real(-8.4820658341042721620064832077422e-01L);
400  _points[ 3](0) = Real(-7.2441773136017004741618605461394e-01L);
401  _points[ 4](0) = Real(-5.7097217260853884753722673725391e-01L);
402  _points[ 5](0) = Real(-3.9415134707756336989720737098105e-01L);
403  _points[ 6](0) = Real(-2.0119409399743452230062830339460e-01L);
404  _points[ 7](0) = 0.;
405  _points[ 8] = -_points[6];
406  _points[ 9] = -_points[5];
407  _points[10] = -_points[4];
408  _points[11] = -_points[3];
409  _points[12] = -_points[2];
410  _points[13] = -_points[1];
411  _points[14] = -_points[0];
412 
413  _weights[ 0] = Real(3.0753241996117268354628393577204e-02L);
414  _weights[ 1] = Real(7.0366047488108124709267416450667e-02L);
415  _weights[ 2] = Real(1.0715922046717193501186954668587e-01L);
416  _weights[ 3] = Real(1.3957067792615431444780479451103e-01L);
417  _weights[ 4] = Real(1.6626920581699393355320086048121e-01L);
418  _weights[ 5] = Real(1.8616100001556221102680056186642e-01L);
419  _weights[ 6] = Real(1.9843148532711157645611832644384e-01L);
420  _weights[ 7] = Real(2.0257824192556127288062019996752e-01L);
421  _weights[ 8] = _weights[6];
422  _weights[ 9] = _weights[5];
423  _weights[10] = _weights[4];
424  _weights[11] = _weights[3];
425  _weights[12] = _weights[2];
426  _weights[13] = _weights[1];
427  _weights[14] = _weights[0];
428 
429  return;
430  }
431 
432  case THIRTIETH:
433  case THIRTYFIRST:
434  {
435  _points.resize (16);
436  _weights.resize(16);
437 
438  _points[ 0](0) = Real(-9.8940093499164993259615417345033e-01L);
439  _points[ 1](0) = Real(-9.4457502307323257607798841553461e-01L);
440  _points[ 2](0) = Real(-8.6563120238783174388046789771239e-01L);
441  _points[ 3](0) = Real(-7.5540440835500303389510119484744e-01L);
442  _points[ 4](0) = Real(-6.1787624440264374844667176404879e-01L);
443  _points[ 5](0) = Real(-4.5801677765722738634241944298358e-01L);
444  _points[ 6](0) = Real(-2.8160355077925891323046050146050e-01L);
445  _points[ 7](0) = Real(-9.5012509837637440185319335424958e-02L);
446  _points[ 8] = -_points[7];
447  _points[ 9] = -_points[6];
448  _points[10] = -_points[5];
449  _points[11] = -_points[4];
450  _points[12] = -_points[3];
451  _points[13] = -_points[2];
452  _points[14] = -_points[1];
453  _points[15] = -_points[0];
454 
455  _weights[ 0] = Real(2.7152459411754094851780572456018e-02L);
456  _weights[ 1] = Real(6.2253523938647892862843836994378e-02L);
457  _weights[ 2] = Real(9.5158511682492784809925107602246e-02L);
458  _weights[ 3] = Real(1.2462897125553387205247628219202e-01L);
459  _weights[ 4] = Real(1.4959598881657673208150173054748e-01L);
460  _weights[ 5] = Real(1.6915651939500253818931207903033e-01L);
461  _weights[ 6] = Real(1.8260341504492358886676366796922e-01L);
462  _weights[ 7] = Real(1.8945061045506849628539672320828e-01L);
463  _weights[ 8] = _weights[7];
464  _weights[ 9] = _weights[6];
465  _weights[10] = _weights[5];
466  _weights[11] = _weights[4];
467  _weights[12] = _weights[3];
468  _weights[13] = _weights[2];
469  _weights[14] = _weights[1];
470  _weights[15] = _weights[0];
471 
472  return;
473  }
474 
475  case THIRTYSECOND:
476  case THIRTYTHIRD:
477  {
478  _points.resize (17);
479  _weights.resize(17);
480 
481  _points[ 0](0) = Real(-9.9057547531441733567543401994067e-01L);
482  _points[ 1](0) = Real(-9.5067552176876776122271695789580e-01L);
483  _points[ 2](0) = Real(-8.8023915372698590212295569448816e-01L);
484  _points[ 3](0) = Real(-7.8151400389680140692523005552048e-01L);
485  _points[ 4](0) = Real(-6.5767115921669076585030221664300e-01L);
486  _points[ 5](0) = Real(-5.1269053708647696788624656862955e-01L);
487  _points[ 6](0) = Real(-3.5123176345387631529718551709535e-01L);
488  _points[ 7](0) = Real(-1.7848418149584785585067749365407e-01L);
489  _points[ 8](0) = 0.;
490  _points[ 9] = -_points[7];
491  _points[10] = -_points[6];
492  _points[11] = -_points[5];
493  _points[12] = -_points[4];
494  _points[13] = -_points[3];
495  _points[14] = -_points[2];
496  _points[15] = -_points[1];
497  _points[16] = -_points[0];
498 
499  _weights[ 0] = Real(2.4148302868547931960110026287565e-02L);
500  _weights[ 1] = Real(5.5459529373987201129440165358245e-02L);
501  _weights[ 2] = Real(8.5036148317179180883535370191062e-02L);
502  _weights[ 3] = Real(1.1188384719340397109478838562636e-01L);
503  _weights[ 4] = Real(1.3513636846852547328631998170235e-01L);
504  _weights[ 5] = Real(1.5404576107681028808143159480196e-01L);
505  _weights[ 6] = Real(1.6800410215645004450997066378832e-01L);
506  _weights[ 7] = Real(1.7656270536699264632527099011320e-01L);
507  _weights[ 8] = Real(1.7944647035620652545826564426189e-01L);
508  _weights[ 9] = _weights[7];
509  _weights[10] = _weights[6];
510  _weights[11] = _weights[5];
511  _weights[12] = _weights[4];
512  _weights[13] = _weights[3];
513  _weights[14] = _weights[2];
514  _weights[15] = _weights[1];
515  _weights[16] = _weights[0];
516 
517  return;
518  }
519 
520  case THIRTYFOURTH:
521  case THIRTYFIFTH:
522  {
523  _points.resize (18);
524  _weights.resize(18);
525 
526  _points[ 0](0) = Real(-9.9156516842093094673001600470615e-01L);
527  _points[ 1](0) = Real(-9.5582394957139775518119589292978e-01L);
528  _points[ 2](0) = Real(-8.9260246649755573920606059112715e-01L);
529  _points[ 3](0) = Real(-8.0370495897252311568241745501459e-01L);
530  _points[ 4](0) = Real(-6.9168704306035320787489108128885e-01L);
531  _points[ 5](0) = Real(-5.5977083107394753460787154852533e-01L);
532  _points[ 6](0) = Real(-4.1175116146284264603593179383305e-01L);
533  _points[ 7](0) = Real(-2.5188622569150550958897285487791e-01L);
534  _points[ 8](0) = Real(-8.4775013041735301242261852935784e-02L);
535  _points[ 9] = -_points[8];
536  _points[10] = -_points[7];
537  _points[11] = -_points[6];
538  _points[12] = -_points[5];
539  _points[13] = -_points[4];
540  _points[14] = -_points[3];
541  _points[15] = -_points[2];
542  _points[16] = -_points[1];
543  _points[17] = -_points[0];
544 
545  _weights[ 0] = Real(2.1616013526483310313342710266452e-02L);
546  _weights[ 1] = Real(4.9714548894969796453334946202639e-02L);
547  _weights[ 2] = Real(7.6425730254889056529129677616637e-02L);
548  _weights[ 3] = Real(1.0094204410628716556281398492483e-01L);
549  _weights[ 4] = Real(1.2255520671147846018451912680020e-01L);
550  _weights[ 5] = Real(1.4064291467065065120473130375195e-01L);
551  _weights[ 6] = Real(1.5468467512626524492541800383637e-01L);
552  _weights[ 7] = Real(1.6427648374583272298605377646593e-01L);
553  _weights[ 8] = Real(1.6914238296314359184065647013499e-01L);
554  _weights[ 9] = _weights[8];
555  _weights[10] = _weights[7];
556  _weights[11] = _weights[6];
557  _weights[12] = _weights[5];
558  _weights[13] = _weights[4];
559  _weights[14] = _weights[3];
560  _weights[15] = _weights[2];
561  _weights[16] = _weights[1];
562  _weights[17] = _weights[0];
563 
564  return;
565  }
566 
567  case THIRTYSIXTH:
568  case THIRTYSEVENTH:
569  {
570  _points.resize (19);
571  _weights.resize(19);
572 
573  _points[ 0](0) = Real(-9.9240684384358440318901767025326e-01L);
574  _points[ 1](0) = Real(-9.6020815213483003085277884068765e-01L);
575  _points[ 2](0) = Real(-9.0315590361481790164266092853231e-01L);
576  _points[ 3](0) = Real(-8.2271465653714282497892248671271e-01L);
577  _points[ 4](0) = Real(-7.2096617733522937861709586082378e-01L);
578  _points[ 5](0) = Real(-6.0054530466168102346963816494624e-01L);
579  _points[ 6](0) = Real(-4.6457074137596094571726714810410e-01L);
580  _points[ 7](0) = Real(-3.1656409996362983199011732884984e-01L);
581  _points[ 8](0) = Real(-1.6035864564022537586809611574074e-01L);
582  _points[ 9](0) = 0.;
583  _points[10] = -_points[8];
584  _points[11] = -_points[7];
585  _points[12] = -_points[6];
586  _points[13] = -_points[5];
587  _points[14] = -_points[4];
588  _points[15] = -_points[3];
589  _points[16] = -_points[2];
590  _points[17] = -_points[1];
591  _points[18] = -_points[0];
592 
593  _weights[ 0] = Real(1.9461788229726477036312041464438e-02L);
594  _weights[ 1] = Real(4.4814226765699600332838157401994e-02L);
595  _weights[ 2] = Real(6.9044542737641226580708258006013e-02L);
596  _weights[ 3] = Real(9.1490021622449999464462094123840e-02L);
597  _weights[ 4] = Real(1.1156664554733399471602390168177e-01L);
598  _weights[ 5] = Real(1.2875396253933622767551578485688e-01L);
599  _weights[ 6] = Real(1.4260670217360661177574610944190e-01L);
600  _weights[ 7] = Real(1.5276604206585966677885540089766e-01L);
601  _weights[ 8] = Real(1.5896884339395434764995643946505e-01L);
602  _weights[ 9] = Real(1.6105444984878369597916362532092e-01L);
603  _weights[10] = _weights[8];
604  _weights[11] = _weights[7];
605  _weights[12] = _weights[6];
606  _weights[13] = _weights[5];
607  _weights[14] = _weights[4];
608  _weights[15] = _weights[3];
609  _weights[16] = _weights[2];
610  _weights[17] = _weights[1];
611  _weights[18] = _weights[0];
612 
613  return;
614  }
615 
616  case THIRTYEIGHTH:
617  case THIRTYNINTH:
618  {
619  _points.resize (20);
620  _weights.resize(20);
621 
622  _points[ 0](0) = Real(-9.9312859918509492478612238847132e-01L);
623  _points[ 1](0) = Real(-9.6397192727791379126766613119728e-01L);
624  _points[ 2](0) = Real(-9.1223442825132590586775244120330e-01L);
625  _points[ 3](0) = Real(-8.3911697182221882339452906170152e-01L);
626  _points[ 4](0) = Real(-7.4633190646015079261430507035564e-01L);
627  _points[ 5](0) = Real(-6.3605368072651502545283669622629e-01L);
628  _points[ 6](0) = Real(-5.1086700195082709800436405095525e-01L);
629  _points[ 7](0) = Real(-3.7370608871541956067254817702493e-01L);
630  _points[ 8](0) = Real(-2.2778585114164507808049619536857e-01L);
631  _points[ 9](0) = Real(-7.6526521133497333754640409398838e-02L);
632  _points[10] = -_points[9];
633  _points[11] = -_points[8];
634  _points[12] = -_points[7];
635  _points[13] = -_points[6];
636  _points[14] = -_points[5];
637  _points[15] = -_points[4];
638  _points[16] = -_points[3];
639  _points[17] = -_points[2];
640  _points[18] = -_points[1];
641  _points[19] = -_points[0];
642 
643  _weights[ 0] = Real(1.7614007139152118311861962351853e-02L);
644  _weights[ 1] = Real(4.0601429800386941331039952274932e-02L);
645  _weights[ 2] = Real(6.2672048334109063569506535187042e-02L);
646  _weights[ 3] = Real(8.3276741576704748724758143222046e-02L);
647  _weights[ 4] = Real(1.0193011981724043503675013548035e-01L);
648  _weights[ 5] = Real(1.1819453196151841731237737771138e-01L);
649  _weights[ 6] = Real(1.3168863844917662689849449974816e-01L);
650  _weights[ 7] = Real(1.4209610931838205132929832506716e-01L);
651  _weights[ 8] = Real(1.4917298647260374678782873700197e-01L);
652  _weights[ 9] = Real(1.5275338713072585069808433195510e-01L);
653  _weights[10] = _weights[9];
654  _weights[11] = _weights[8];
655  _weights[12] = _weights[7];
656  _weights[13] = _weights[6];
657  _weights[14] = _weights[5];
658  _weights[15] = _weights[4];
659  _weights[16] = _weights[3];
660  _weights[17] = _weights[2];
661  _weights[18] = _weights[1];
662  _weights[19] = _weights[0];
663 
664  return;
665  }
666 
667  case FORTIETH:
668  case FORTYFIRST:
669  {
670  _points.resize (21);
671  _weights.resize(21);
672 
673  _points[ 0](0) = Real(-9.9375217062038950026024203593794e-01L);
674  _points[ 1](0) = Real(-9.6722683856630629431662221490770e-01L);
675  _points[ 2](0) = Real(-9.2009933415040082879018713371497e-01L);
676  _points[ 3](0) = Real(-8.5336336458331728364725063858757e-01L);
677  _points[ 4](0) = Real(-7.6843996347567790861587785130623e-01L);
678  _points[ 5](0) = Real(-6.6713880419741231930596666999034e-01L);
679  _points[ 6](0) = Real(-5.5161883588721980705901879672431e-01L);
680  _points[ 7](0) = Real(-4.2434212020743878357366888854379e-01L);
681  _points[ 8](0) = Real(-2.8802131680240109660079251606460e-01L);
682  _points[ 9](0) = Real(-1.4556185416089509093703098233869e-01L);
683  _points[10](0) = 0.;
684  _points[11] = -_points[9];
685  _points[12] = -_points[8];
686  _points[13] = -_points[7];
687  _points[14] = -_points[6];
688  _points[15] = -_points[5];
689  _points[16] = -_points[4];
690  _points[17] = -_points[3];
691  _points[18] = -_points[2];
692  _points[19] = -_points[1];
693  _points[20] = -_points[0];
694 
695  _weights[ 0] = Real(1.6017228257774333324224616858471e-02L);
696  _weights[ 1] = Real(3.6953789770852493799950668299330e-02L);
697  _weights[ 2] = Real(5.7134425426857208283635826472448e-02L);
698  _weights[ 3] = Real(7.6100113628379302017051653300183e-02L);
699  _weights[ 4] = Real(9.3444423456033861553289741113932e-02L);
700  _weights[ 5] = Real(1.0879729916714837766347457807011e-01L);
701  _weights[ 6] = Real(1.2183141605372853419536717712572e-01L);
702  _weights[ 7] = Real(1.3226893863333746178105257449678e-01L);
703  _weights[ 8] = Real(1.3988739479107315472213342386758e-01L);
704  _weights[ 9] = Real(1.4452440398997005906382716655375e-01L);
705  _weights[10] = Real(1.4608113364969042719198514768337e-01L);
706  _weights[11] = _weights[9];
707  _weights[12] = _weights[8];
708  _weights[13] = _weights[7];
709  _weights[14] = _weights[6];
710  _weights[15] = _weights[5];
711  _weights[16] = _weights[4];
712  _weights[17] = _weights[3];
713  _weights[18] = _weights[2];
714  _weights[19] = _weights[1];
715  _weights[20] = _weights[0];
716 
717  return;
718  }
719 
720  case FORTYSECOND:
721  case FORTYTHIRD:
722  {
723  _points.resize (22);
724  _weights.resize(22);
725 
726  _points[ 0](0) = Real(-9.9429458548239929207303142116130e-01L);
727  _points[ 1](0) = Real(-9.7006049783542872712395098676527e-01L);
728  _points[ 2](0) = Real(-9.2695677218717400052069293925905e-01L);
729  _points[ 3](0) = Real(-8.6581257772030013653642563701938e-01L);
730  _points[ 4](0) = Real(-7.8781680597920816200427795540835e-01L);
731  _points[ 5](0) = Real(-6.9448726318668278005068983576226e-01L);
732  _points[ 6](0) = Real(-5.8764040350691159295887692763865e-01L);
733  _points[ 7](0) = Real(-4.6935583798675702640633071096641e-01L);
734  _points[ 8](0) = Real(-3.4193582089208422515814742042738e-01L);
735  _points[ 9](0) = Real(-2.0786042668822128547884653391955e-01L);
736  _points[10](0) = Real(-6.9739273319722221213841796118628e-02L);
737  _points[11] = -_points[10];
738  _points[12] = -_points[9];
739  _points[13] = -_points[8];
740  _points[14] = -_points[7];
741  _points[15] = -_points[6];
742  _points[16] = -_points[5];
743  _points[17] = -_points[4];
744  _points[18] = -_points[3];
745  _points[19] = -_points[2];
746  _points[20] = -_points[1];
747  _points[21] = -_points[0];
748 
749  _weights[ 0] = Real(1.4627995298272200684991098047185e-02L);
750  _weights[ 1] = Real(3.3774901584814154793302246865913e-02L);
751  _weights[ 2] = Real(5.2293335152683285940312051273211e-02L);
752  _weights[ 3] = Real(6.9796468424520488094961418930218e-02L);
753  _weights[ 4] = Real(8.5941606217067727414443681372703e-02L);
754  _weights[ 5] = Real(1.0041414444288096493207883783054e-01L);
755  _weights[ 6] = Real(1.1293229608053921839340060742175e-01L);
756  _weights[ 7] = Real(1.2325237681051242428556098615481e-01L);
757  _weights[ 8] = Real(1.3117350478706237073296499253031e-01L);
758  _weights[ 9] = Real(1.3654149834601517135257383123152e-01L);
759  _weights[10] = Real(1.3925187285563199337541024834181e-01L);
760  _weights[11] = _weights[10];
761  _weights[12] = _weights[9];
762  _weights[13] = _weights[8];
763  _weights[14] = _weights[7];
764  _weights[15] = _weights[6];
765  _weights[16] = _weights[5];
766  _weights[17] = _weights[4];
767  _weights[18] = _weights[3];
768  _weights[19] = _weights[2];
769  _weights[20] = _weights[1];
770  _weights[21] = _weights[0];
771 
772  return;
773  }
774 
775 
776  default:
777  libmesh_error_msg("Quadrature rule " << _order << " not supported!");
778  }
779 }

References libMesh::QBase::_order, libMesh::QBase::_points, libMesh::QBase::_weights, libMesh::CONSTANT, libMesh::EIGHTH, libMesh::EIGHTTEENTH, libMesh::ELEVENTH, libMesh::FIFTEENTH, libMesh::FIFTH, libMesh::FIRST, libMesh::FORTIETH, libMesh::FORTYFIRST, libMesh::FORTYSECOND, libMesh::FORTYTHIRD, libMesh::FOURTEENTH, libMesh::FOURTH, libMesh::QBase::get_order(), libMesh::NINETEENTH, libMesh::NINTH, libMesh::Real, libMesh::SECOND, libMesh::SEVENTEENTH, libMesh::SEVENTH, libMesh::SIXTEENTH, libMesh::SIXTH, libMesh::TENTH, libMesh::THIRD, libMesh::THIRTEENTH, libMesh::THIRTIETH, libMesh::THIRTYEIGHTH, libMesh::THIRTYFIFTH, libMesh::THIRTYFIRST, libMesh::THIRTYFOURTH, libMesh::THIRTYNINTH, libMesh::THIRTYSECOND, libMesh::THIRTYSEVENTH, libMesh::THIRTYSIXTH, libMesh::THIRTYTHIRD, libMesh::TWELFTH, libMesh::TWENTIETH, libMesh::TWENTYEIGHTH, libMesh::TWENTYFIFTH, libMesh::TWENTYFIRST, libMesh::TWENTYFOURTH, libMesh::TWENTYNINTH, libMesh::TWENTYSECOND, libMesh::TWENTYSEVENTH, libMesh::TWENTYSIXTH, and libMesh::TWENTYTHIRD.

◆ init_2D()

void libMesh::QGauss::init_2D ( const  type,
unsigned int  p_level 
)
overrideprivatevirtual

Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate values.

The order of the rule will be defined by the implementing class. Should not be pure virtual since a derived quadrature rule may only be defined in 1D. If not overridden, throws an error.

Note
The arguments should no longer be used for anything in derived classes, they are only maintained for backwards compatibility and will eventually be removed.

Reimplemented from libMesh::QBase.

Definition at line 29 of file quadrature_gauss_2D.C.

30 {
31 #if LIBMESH_DIM > 1
32 
33  //-----------------------------------------------------------------------
34  // 2D quadrature rules
35  switch (_type)
36  {
37 
38 
39  //---------------------------------------------
40  // Quadrilateral quadrature rules
41  case QUAD4:
42  case QUADSHELL4:
43  case QUAD8:
44  case QUADSHELL8:
45  case QUAD9:
46  {
47  // We compute the 2D quadrature rule as a tensor
48  // product of the 1D quadrature rule.
49  //
50  // For QUADs, a quadrature rule of order 'p' must be able to integrate
51  // bilinear (p=1), biquadratic (p=2), bicubic (p=3), etc. polynomials of the form
52  //
53  // (x^p + x^{p-1} + ... + 1) * (y^p + y^{p-1} + ... + 1)
54  //
55  // These polynomials have terms *up to* degree 2p but they are *not* complete
56  // polynomials of degree 2p. For example, when p=2 we have
57  // 1
58  // x y
59  // x^2 xy y^2
60  // yx^2 xy^2
61  // x^2y^2
62  QGauss q1D(1,_order);
63  q1D.init(EDGE2, _p_level);
64  tensor_product_quad( q1D );
65  return;
66  }
67 
68 
69  //---------------------------------------------
70  // Triangle quadrature rules
71  case TRI3:
72  case TRISHELL3:
73  case TRI3SUBDIVISION:
74  case TRI6:
75  {
76  switch(get_order())
77  {
78  case CONSTANT:
79  case FIRST:
80  {
81  // Exact for linears
82  _points.resize(1);
83  _weights.resize(1);
84 
85  _points[0](0) = Real(1)/3;
86  _points[0](1) = Real(1)/3;
87 
88  _weights[0] = 0.5;
89 
90  return;
91  }
92  case SECOND:
93  {
94  // Exact for quadratics
95  _points.resize(3);
96  _weights.resize(3);
97 
98  // Alternate rule with points on ref. elt. boundaries.
99  // Not ideal for problems with material coefficient discontinuities
100  // aligned along element boundaries.
101  // _points[0](0) = .5;
102  // _points[0](1) = .5;
103  // _points[1](0) = 0.;
104  // _points[1](1) = .5;
105  // _points[2](0) = .5;
106  // _points[2](1) = .0;
107 
108  _points[0](0) = Real(2)/3;
109  _points[0](1) = Real(1)/6;
110 
111  _points[1](0) = Real(1)/6;
112  _points[1](1) = Real(2)/3;
113 
114  _points[2](0) = Real(1)/6;
115  _points[2](1) = Real(1)/6;
116 
117 
118  _weights[0] = Real(1)/6;
119  _weights[1] = Real(1)/6;
120  _weights[2] = Real(1)/6;
121 
122  return;
123  }
124  case THIRD:
125  {
126  // Exact for cubics
127  _points.resize(4);
128  _weights.resize(4);
129 
130  // This rule is formed from a tensor product of
131  // appropriately-scaled Gauss and Jacobi rules. (See
132  // also the QConical quadrature class, this is a
133  // hard-coded version of one of those rules.) For high
134  // orders these rules generally have too many points,
135  // but at extremely low order they are competitive and
136  // have the additional benefit of having all positive
137  // weights.
138  _points[0](0) = Real(1.5505102572168219018027159252941e-01L);
139  _points[0](1) = Real(1.7855872826361642311703513337422e-01L);
140  _points[1](0) = Real(6.4494897427831780981972840747059e-01L);
141  _points[1](1) = Real(7.5031110222608118177475598324603e-02L);
142  _points[2](0) = Real(1.5505102572168219018027159252941e-01L);
143  _points[2](1) = Real(6.6639024601470138670269327409637e-01L);
144  _points[3](0) = Real(6.4494897427831780981972840747059e-01L);
145  _points[3](1) = Real(2.8001991549907407200279599420481e-01L);
146 
147  _weights[0] = Real(1.5902069087198858469718450103758e-01L);
148  _weights[1] = Real(9.0979309128011415302815498962418e-02L);
149  _weights[2] = Real(1.5902069087198858469718450103758e-01L);
150  _weights[3] = Real(9.0979309128011415302815498962418e-02L);
151 
152  return;
153 
154 
155  // The following third-order rule is quite commonly cited
156  // in the literature and most likely works fine. However,
157  // we generally prefer a rule with all positive weights
158  // and an equal number of points, when available.
159  //
160  // (allow_rules_with_negative_weights)
161  // {
162  // // Exact for cubics
163  // _points.resize(4);
164  // _weights.resize(4);
165  //
166  // _points[0](0) = .33333333333333333333333333333333;
167  // _points[0](1) = .33333333333333333333333333333333;
168  //
169  // _points[1](0) = .2;
170  // _points[1](1) = .6;
171  //
172  // _points[2](0) = .2;
173  // _points[2](1) = .2;
174  //
175  // _points[3](0) = .6;
176  // _points[3](1) = .2;
177  //
178  //
179  // _weights[0] = -27./96.;
180  // _weights[1] = 25./96.;
181  // _weights[2] = 25./96.;
182  // _weights[3] = 25./96.;
183  //
184  // return;
185  // } // end if (allow_rules_with_negative_weights)
186  // Note: if !allow_rules_with_negative_weights, fall through to next case.
187  }
188 
189 
190 
191  // A degree 4 rule with six points. This rule can be found in many places
192  // including:
193  //
194  // J.N. Lyness and D. Jespersen, Moderate degree symmetric
195  // quadrature rules for the triangle, J. Inst. Math. Appl. 15 (1975),
196  // 19--32.
197  //
198  // We used the code in:
199  // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
200  // on triangles and tetrahedra" Journal of Computational Mathematics,
201  // v. 27, no. 1, 2009, pp. 89-96.
202  // to generate additional precision.
203  case FOURTH:
204  {
205  const unsigned int n_wts = 2;
206  const Real wts[n_wts] =
207  {
208  Real(1.1169079483900573284750350421656140e-01L),
209  Real(5.4975871827660933819163162450105264e-02L)
210  };
211 
212  const Real a[n_wts] =
213  {
214  Real(4.4594849091596488631832925388305199e-01L),
215  Real(9.1576213509770743459571463402201508e-02L)
216  };
217 
218  const Real b[n_wts] = {0., 0.}; // not used
219  const unsigned int permutation_ids[n_wts] = {3, 3};
220 
221  dunavant_rule2(wts, a, b, permutation_ids, n_wts); // 6 total points
222 
223  return;
224  }
225 
226 
227 
228  // Exact for quintics
229  // Can be found in "Quadrature on Simplices of Arbitrary
230  // Dimension" by Walkington.
231  case FIFTH:
232  {
233  const unsigned int n_wts = 3;
234  const Real wts[n_wts] =
235  {
236  Real(9)/80,
237  Real(31)/480 + std::sqrt(Real(15.0L))/2400,
238  Real(31)/480 - std::sqrt(Real(15.0L))/2400
239  };
240 
241  const Real a[n_wts] =
242  {
243  0., // 'a' parameter not used for origin permutation
244  Real(2)/7 + std::sqrt(Real(15.0L))/21,
245  Real(2)/7 - std::sqrt(Real(15.0L))/21
246  };
247 
248  const Real b[n_wts] = {0., 0., 0.}; // not used
249  const unsigned int permutation_ids[n_wts] = {1, 3, 3};
250 
251  dunavant_rule2(wts, a, b, permutation_ids, n_wts); // 7 total points
252 
253  return;
254  }
255 
256 
257 
258  // A degree 6 rule with 12 points. This rule can be found in many places
259  // including:
260  //
261  // J.N. Lyness and D. Jespersen, Moderate degree symmetric
262  // quadrature rules for the triangle, J. Inst. Math. Appl. 15 (1975),
263  // 19--32.
264  //
265  // We used the code in:
266  // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
267  // on triangles and tetrahedra" Journal of Computational Mathematics,
268  // v. 27, no. 1, 2009, pp. 89-96.
269  // to generate additional precision.
270  //
271  // Note that the following 7th-order Ro3-invariant rule also has only 12 points,
272  // which technically makes it the superior rule. This one is here for completeness.
273  case SIXTH:
274  {
275  const unsigned int n_wts = 3;
276  const Real wts[n_wts] =
277  {
278  Real(5.8393137863189683012644805692789721e-02L),
279  Real(2.5422453185103408460468404553434492e-02L),
280  Real(4.1425537809186787596776728210221227e-02L)
281  };
282 
283  const Real a[n_wts] =
284  {
285  Real(2.4928674517091042129163855310701908e-01L),
286  Real(6.3089014491502228340331602870819157e-02L),
287  Real(3.1035245103378440541660773395655215e-01L)
288  };
289 
290  const Real b[n_wts] =
291  {
292  0.,
293  0.,
294  Real(6.3650249912139864723014259441204970e-01L)
295  };
296 
297  const unsigned int permutation_ids[n_wts] = {3, 3, 6}; // 12 total points
298 
299  dunavant_rule2(wts, a, b, permutation_ids, n_wts);
300 
301  return;
302  }
303 
304 
305  // A degree 7 rule with 12 points. This rule can be found in:
306  //
307  // K. Gatermann, The construction of symmetric cubature
308  // formulas for the square and the triangle, Computing 40
309  // (1988), 229--240.
310  //
311  // This rule, which is provably minimal in the number of
312  // integration points, is said to be 'Ro3 invariant' which
313  // means that a given set of barycentric coordinates
314  // (z1,z2,z3) implies the quadrature points (z1,z2),
315  // (z3,z1), (z2,z3) which are formed by taking the first
316  // two entries in cyclic permutations of the barycentric
317  // point. Barycentric coordinates are related in the
318  // sense that: z3 = 1 - z1 - z2.
319  //
320  // The 12-point sixth-order rule for triangles given in
321  // Flaherty's (http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps)
322  // lecture notes has been removed in favor of this rule
323  // which is higher-order (for the same number of
324  // quadrature points) and has a few more digits of
325  // precision in the points and weights. Some 10-point
326  // degree 6 rules exist for the triangle but they have
327  // quadrature points outside the region of integration.
328  case SEVENTH:
329  {
330  _points.resize (12);
331  _weights.resize(12);
332 
333  const unsigned int nrows=4;
334 
335  // In each of the rows below, the first two entries are (z1, z2) which imply
336  // z3. The third entry is the weight for each of the points in the cyclic permutation.
337  // The original publication tabulated about 16 decimal digits for each point and weight
338  // parameter. The additional digits shown here were obtained using a code in the
339  // mp-quadrature library, https://github.com/jwpeterson/mp-quadrature
340  const Real rule_data[nrows][3] = {
341  {Real(6.2382265094402118173683000996350e-02L), Real(6.7517867073916085442557131050869e-02L), Real(2.6517028157436251428754180460739e-02L)}, // group A
342  {Real(5.5225456656926611737479190275645e-02L), Real(3.2150249385198182266630784919920e-01L), Real(4.3881408714446055036769903139288e-02L)}, // group B
343  {Real(3.4324302945097146469630642483938e-02L), Real(6.6094919618673565761198031019780e-01L), Real(2.8775042784981585738445496900219e-02L)}, // group C
344  {Real(5.1584233435359177925746338682643e-01L), Real(2.7771616697639178256958187139372e-01L), Real(6.7493187009802774462697086166421e-02L)} // group D
345  };
346 
347  for (unsigned int i=0, offset=0; i<nrows; ++i)
348  {
349  _points[offset + 0] = Point(rule_data[i][0], rule_data[i][1]); // (z1,z2)
350  _points[offset + 1] = Point(1.-rule_data[i][0]-rule_data[i][1], rule_data[i][0]); // (z3,z1)
351  _points[offset + 2] = Point(rule_data[i][1], 1.-rule_data[i][0]-rule_data[i][1]); // (z2,z3)
352 
353  // All these points get the same weight
354  _weights[offset + 0] = rule_data[i][2];
355  _weights[offset + 1] = rule_data[i][2];
356  _weights[offset + 2] = rule_data[i][2];
357 
358  // Increment offset
359  offset += 3;
360  }
361 
362  return;
363 
364 
365  // // The following is an inferior 7th-order Lyness-style rule with 15 points.
366  // // It's here only for completeness and the Ro3-invariant rule above should
367  // // be used instead!
368  // const unsigned int n_wts = 3;
369  // const Real wts[n_wts] =
370  // {
371  // Real(2.6538900895116205835977487499847719e-02L),
372  // Real(3.5426541846066783659206291623201826e-02L),
373  // Real(3.4637341039708446756138297960207647e-02L)
374  // };
375  //
376  // const Real a[n_wts] =
377  // {
378  // Real(6.4930513159164863078379776030396538e-02L),
379  // Real(2.8457558424917033519741605734978046e-01L),
380  // Real(3.1355918438493150795585190219862865e-01L)
381  // };
382  //
383  // const Real b[n_wts] =
384  // {
385  // 0.,
386  // Real(1.9838447668150671917987659863332941e-01L),
387  // Real(4.3863471792372471511798695971295936e-02L)
388  // };
389  //
390  // const unsigned int permutation_ids[n_wts] = {3, 6, 6}; // 15 total points
391  //
392  // dunavant_rule2(wts, a, b, permutation_ids, n_wts);
393  //
394  // return;
395  }
396 
397 
398 
399 
400  // Another Dunavant rule. This one has all positive weights. This rule has
401  // 16 points while a comparable conical product rule would have 5*5=25.
402  //
403  // It was copied 23rd June 2008 from:
404  // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
405  //
406  // Additional precision obtained from the code in:
407  // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
408  // on triangles and tetrahedra" Journal of Computational Mathematics,
409  // v. 27, no. 1, 2009, pp. 89-96.
410  case EIGHTH:
411  {
412  const unsigned int n_wts = 5;
413  const Real wts[n_wts] =
414  {
415  Real(7.2157803838893584125545555244532310e-02L),
416  Real(4.7545817133642312396948052194292159e-02L),
417  Real(5.1608685267359125140895775146064515e-02L),
418  Real(1.6229248811599040155462964170890299e-02L),
419  Real(1.3615157087217497132422345036954462e-02L)
420  };
421 
422  const Real a[n_wts] =
423  {
424  0.0, // 'a' parameter not used for origin permutation
425  Real(4.5929258829272315602881551449416932e-01L),
426  Real(1.7056930775176020662229350149146450e-01L),
427  Real(5.0547228317030975458423550596598947e-02L),
428  Real(2.6311282963463811342178578628464359e-01L),
429  };
430 
431  const Real b[n_wts] =
432  {
433  0.,
434  0.,
435  0.,
436  0.,
437  Real(7.2849239295540428124100037917606196e-01L)
438  };
439 
440  const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 6}; // 16 total points
441 
442  dunavant_rule2(wts, a, b, permutation_ids, n_wts);
443 
444  return;
445  }
446 
447 
448 
449  // Another Dunavant rule. This one has all positive weights. This rule has 19
450  // points. The comparable conical product rule would have 25.
451  // It was copied 23rd June 2008 from:
452  // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
453  //
454  // Additional precision obtained from the code in:
455  // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
456  // on triangles and tetrahedra" Journal of Computational Mathematics,
457  // v. 27, no. 1, 2009, pp. 89-96.
458  case NINTH:
459  {
460  const unsigned int n_wts = 6;
461  const Real wts[n_wts] =
462  {
463  Real(4.8567898141399416909620991253644315e-02L),
464  Real(1.5667350113569535268427415643604658e-02L),
465  Real(1.2788837829349015630839399279499912e-02L),
466  Real(3.8913770502387139658369678149701978e-02L),
467  Real(3.9823869463605126516445887132022637e-02L),
468  Real(2.1641769688644688644688644688644689e-02L)
469  };
470 
471  const Real a[n_wts] =
472  {
473  0.0, // 'a' parameter not used for origin permutation
474  Real(4.8968251919873762778370692483619280e-01L),
475  Real(4.4729513394452709865106589966276365e-02L),
476  Real(4.3708959149293663726993036443535497e-01L),
477  Real(1.8820353561903273024096128046733557e-01L),
478  Real(2.2196298916076569567510252769319107e-01L)
479  };
480 
481  const Real b[n_wts] =
482  {
483  0.,
484  0.,
485  0.,
486  0.,
487  0.,
488  Real(7.4119859878449802069007987352342383e-01L)
489  };
490 
491  const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 3, 6}; // 19 total points
492 
493  dunavant_rule2(wts, a, b, permutation_ids, n_wts);
494 
495  return;
496  }
497 
498 
499  // Another Dunavant rule with all positive weights. This rule has 25
500  // points. The comparable conical product rule would have 36.
501  // It was copied 23rd June 2008 from:
502  // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
503  //
504  // Additional precision obtained from the code in:
505  // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
506  // on triangles and tetrahedra" Journal of Computational Mathematics,
507  // v. 27, no. 1, 2009, pp. 89-96.
508  case TENTH:
509  {
510  const unsigned int n_wts = 6;
511  const Real wts[n_wts] =
512  {
513  Real(4.5408995191376790047643297550014267e-02L),
514  Real(1.8362978878233352358503035945683300e-02L),
515  Real(2.2660529717763967391302822369298659e-02L),
516  Real(3.6378958422710054302157588309680344e-02L),
517  Real(1.4163621265528742418368530791049552e-02L),
518  Real(4.7108334818664117299637354834434138e-03L)
519  };
520 
521  const Real a[n_wts] =
522  {
523  0.0, // 'a' parameter not used for origin permutation
524  Real(4.8557763338365737736750753220812615e-01L),
525  Real(1.0948157548503705479545863134052284e-01L),
526  Real(3.0793983876412095016515502293063162e-01L),
527  Real(2.4667256063990269391727646541117681e-01L),
528  Real(6.6803251012200265773540212762024737e-02L)
529  };
530 
531  const Real b[n_wts] =
532  {
533  0.,
534  0.,
535  0.,
536  Real(5.5035294182099909507816172659300821e-01L),
537  Real(7.2832390459741092000873505358107866e-01L),
538  Real(9.2365593358750027664630697761508843e-01L)
539  };
540 
541  const unsigned int permutation_ids[n_wts] = {1, 3, 3, 6, 6, 6}; // 25 total points
542 
543  dunavant_rule2(wts, a, b, permutation_ids, n_wts);
544 
545  return;
546  }
547 
548 
549  // Dunavant's 11th-order rule contains points outside the region of
550  // integration, and is thus unacceptable for our FEM calculations.
551  //
552  // This 30-point, 11th-order rule was obtained by me [JWP] using the code in
553  //
554  // Additional precision obtained from the code in:
555  // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
556  // on triangles and tetrahedra" Journal of Computational Mathematics,
557  // v. 27, no. 1, 2009, pp. 89-96.
558  //
559  // Note: the 28-point 11th-order rule obtained by Zhang in the paper above
560  // does not appear to be unique. It is a solution in the sense that it
561  // minimizes the error in the least-squares minimization problem, but
562  // it involves too many unknowns and the Jacobian is therefore singular
563  // when attempting to improve the solution via Newton's method.
564  case ELEVENTH:
565  {
566  const unsigned int n_wts = 6;
567  const Real wts[n_wts] =
568  {
569  Real(3.6089021198604635216985338480426484e-02L),
570  Real(2.1607717807680420303346736867931050e-02L),
571  Real(3.1144524293927978774861144478241807e-03L),
572  Real(2.9086855161081509446654185084988077e-02L),
573  Real(8.4879241614917017182977532679947624e-03L),
574  Real(1.3795732078224796530729242858347546e-02L)
575  };
576 
577  const Real a[n_wts] =
578  {
579  Real(3.9355079629947969884346551941969960e-01L),
580  Real(4.7979065808897448654107733982929214e-01L),
581  Real(5.1003445645828061436081405648347852e-03L),
582  Real(2.6597620190330158952732822450744488e-01L),
583  Real(2.8536418538696461608233522814483715e-01L),
584  Real(1.3723536747817085036455583801851025e-01L)
585  };
586 
587  const Real b[n_wts] =
588  {
589  0.,
590  0.,
591  Real(5.6817155788572446538150614865768991e-02L),
592  Real(1.2539956353662088473247489775203396e-01L),
593  Real(1.2409970153698532116262152247041742e-02L),
594  Real(5.2792057988217708934207928630851643e-02L)
595  };
596 
597  const unsigned int permutation_ids[n_wts] = {3, 3, 6, 6, 6, 6}; // 30 total points
598 
599  dunavant_rule2(wts, a, b, permutation_ids, n_wts);
600 
601  return;
602  }
603 
604 
605 
606 
607  // Another Dunavant rule with all positive weights. This rule has 33
608  // points. The comparable conical product rule would have 36 (ELEVENTH) or 49 (TWELFTH).
609  //
610  // It was copied 23rd June 2008 from:
611  // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
612  //
613  // Additional precision obtained from the code in:
614  // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
615  // on triangles and tetrahedra" Journal of Computational Mathematics,
616  // v. 27, no. 1, 2009, pp. 89-96.
617  case TWELFTH:
618  {
619  const unsigned int n_wts = 8;
620  const Real wts[n_wts] =
621  {
622  Real(3.0831305257795086169332418926151771e-03L),
623  Real(3.1429112108942550177135256546441273e-02L),
624  Real(1.7398056465354471494664198647499687e-02L),
625  Real(2.1846272269019201067728631278737487e-02L),
626  Real(1.2865533220227667708895461535782215e-02L),
627  Real(1.1178386601151722855919538351159995e-02L),
628  Real(8.6581155543294461858210504055170332e-03L),
629  Real(2.0185778883190464758914349626118386e-02L)
630  };
631 
632  const Real a[n_wts] =
633  {
634  Real(2.1317350453210370246856975515728246e-02L),
635  Real(2.7121038501211592234595134039689474e-01L),
636  Real(1.2757614554158592467389632515428357e-01L),
637  Real(4.3972439229446027297973662348436108e-01L),
638  Real(4.8821738977380488256466206525881104e-01L),
639  Real(2.8132558098993954824813069297455275e-01L),
640  Real(1.1625191590759714124135414784260182e-01L),
641  Real(2.7571326968551419397479634607976398e-01L)
642  };
643 
644  const Real b[n_wts] =
645  {
646  0.,
647  0.,
648  0.,
649  0.,
650  0.,
651  Real(6.9583608678780342214163552323607254e-01L),
652  Real(8.5801403354407263059053661662617818e-01L),
653  Real(6.0894323577978780685619243776371007e-01L)
654  };
655 
656  const unsigned int permutation_ids[n_wts] = {3, 3, 3, 3, 3, 6, 6, 6}; // 33 total points
657 
658  dunavant_rule2(wts, a, b, permutation_ids, n_wts);
659 
660  return;
661  }
662 
663 
664  // Another Dunavant rule with all positive weights. This rule has 37
665  // points. The comparable conical product rule would have 49 points.
666  //
667  // It was copied 23rd June 2008 from:
668  // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
669  //
670  // A second rule with additional precision obtained from the code in:
671  // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
672  // on triangles and tetrahedra" Journal of Computational Mathematics,
673  // v. 27, no. 1, 2009, pp. 89-96.
674  case THIRTEENTH:
675  {
676  const unsigned int n_wts = 9;
677  const Real wts[n_wts] =
678  {
679  Real(3.3980018293415822140887212340442440e-02L),
680  Real(2.7800983765226664353628733005230734e-02L),
681  Real(2.9139242559599990702383541756669905e-02L),
682  Real(3.0261685517695859208964000161454122e-03L),
683  Real(1.1997200964447365386855399725479827e-02L),
684  Real(1.7320638070424185232993414255459110e-02L),
685  Real(7.4827005525828336316229285664517190e-03L),
686  Real(1.2089519905796909568722872786530380e-02L),
687  Real(4.7953405017716313612975450830554457e-03L)
688  };
689 
690  const Real a[n_wts] =
691  {
692  0., // 'a' parameter not used for origin permutation
693  Real(4.2694141425980040602081253503137421e-01L),
694  Real(2.2137228629183290065481255470507908e-01L),
695  Real(2.1509681108843183869291313534052083e-02L),
696  Real(4.8907694645253934990068971909020439e-01L),
697  Real(3.0844176089211777465847185254124531e-01L),
698  Real(1.1092204280346339541286954522167452e-01L),
699  Real(1.6359740106785048023388790171095725e-01L),
700  Real(2.7251581777342966618005046435408685e-01L)
701  };
702 
703  const Real b[n_wts] =
704  {
705  0.,
706  0.,
707  0.,
708  0.,
709  0.,
710  Real(6.2354599555367557081585435318623659e-01L),
711  Real(8.6470777029544277530254595089569318e-01L),
712  Real(7.4850711589995219517301859578870965e-01L),
713  Real(7.2235779312418796526062013230478405e-01L)
714  };
715 
716  const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 3, 6, 6, 6, 6}; // 37 total points
717 
718  dunavant_rule2(wts, a, b, permutation_ids, n_wts);
719 
720  return;
721  }
722 
723 
724  // Another Dunavant rule. This rule has 42 points, while
725  // a comparable conical product rule would have 64.
726  //
727  // It was copied 23rd June 2008 from:
728  // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
729  //
730  // Additional precision obtained from the code in:
731  // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
732  // on triangles and tetrahedra" Journal of Computational Mathematics,
733  // v. 27, no. 1, 2009, pp. 89-96.
734  case FOURTEENTH:
735  {
736  const unsigned int n_wts = 10;
737  const Real wts[n_wts] =
738  {
739  Real(1.0941790684714445320422472981662986e-02L),
740  Real(1.6394176772062675320655489369312672e-02L),
741  Real(2.5887052253645793157392455083198201e-02L),
742  Real(2.1081294368496508769115218662093065e-02L),
743  Real(7.2168498348883338008549607403266583e-03L),
744  Real(2.4617018012000408409130117545210774e-03L),
745  Real(1.2332876606281836981437622591818114e-02L),
746  Real(1.9285755393530341614244513905205430e-02L),
747  Real(7.2181540567669202480443459995079017e-03L),
748  Real(2.5051144192503358849300465412445582e-03L)
749  };
750 
751  const Real a[n_wts] =
752  {
753  Real(4.8896391036217863867737602045239024e-01L),
754  Real(4.1764471934045392250944082218564344e-01L),
755  Real(2.7347752830883865975494428326269856e-01L),
756  Real(1.7720553241254343695661069046505908e-01L),
757  Real(6.1799883090872601267478828436935788e-02L),
758  Real(1.9390961248701048178250095054529511e-02L),
759  Real(1.7226668782135557837528960161365733e-01L),
760  Real(3.3686145979634500174405519708892539e-01L),
761  Real(2.9837288213625775297083151805961273e-01L),
762  Real(1.1897449769695684539818196192990548e-01L)
763  };
764 
765  const Real b[n_wts] =
766  {
767  0.,
768  0.,
769  0.,
770  0.,
771  0.,
772  0.,
773  Real(7.7060855477499648258903327416742796e-01L),
774  Real(5.7022229084668317349769621336235426e-01L),
775  Real(6.8698016780808783735862715402031306e-01L),
776  Real(8.7975717137017112951457163697460183e-01L)
777  };
778 
779  const unsigned int permutation_ids[n_wts]
780  = {3, 3, 3, 3, 3, 3, 6, 6, 6, 6}; // 42 total points
781 
782  dunavant_rule2(wts, a, b, permutation_ids, n_wts);
783 
784  return;
785  }
786 
787 
788  // This 49-point rule was found by me [JWP] using the code in:
789  //
790  // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
791  // on triangles and tetrahedra" Journal of Computational Mathematics,
792  // v. 27, no. 1, 2009, pp. 89-96.
793  //
794  // A 54-point, 15th-order rule is reported by
795  //
796  // Stephen Wandzura, Hong Xiao,
797  // Symmetric Quadrature Rules on a Triangle,
798  // Computers and Mathematics with Applications,
799  // Volume 45, Number 12, June 2003, pages 1829-1840.
800  //
801  // can be found here:
802  // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
803  //
804  // but this 49-point rule is superior.
805  case FIFTEENTH:
806  {
807  const unsigned int n_wts = 11;
808  const Real wts[n_wts] =
809  {
810  Real(2.4777380743035579804788826970198951e-02L),
811  Real(9.2433943023307730591540642828347660e-03L),
812  Real(2.2485768962175402793245929133296627e-03L),
813  Real(6.7052581900064143760518398833360903e-03L),
814  Real(1.9011381726930579256700190357527956e-02L),
815  Real(1.4605445387471889398286155981802858e-02L),
816  Real(1.5087322572773133722829435011138258e-02L),
817  Real(1.5630213780078803020711746273129099e-02L),
818  Real(6.1808086085778203192616856133701233e-03L),
819  Real(3.2209366452594664857296985751120513e-03L),
820  Real(5.8747373242569702667677969985668817e-03L)
821  };
822 
823  const Real a[n_wts] =
824  {
825  0.0, // 'a' parameter not used for origin
826  Real(7.9031013655541635005816956762252155e-02L),
827  Real(1.8789501810770077611247984432284226e-02L),
828  Real(4.9250168823249670532514526605352905e-01L),
829  Real(4.0886316907744105975059040108092775e-01L),
830  Real(5.3877851064220142445952549348423733e-01L),
831  Real(2.0250549804829997692885033941362673e-01L),
832  Real(5.5349674918711643207148086558288110e-01L),
833  Real(7.8345022567320812359258882143250181e-01L),
834  Real(8.9514624528794883409864566727625002e-01L),
835  Real(3.2515745241110782862789881780746490e-01L)
836  };
837 
838  const Real b[n_wts] =
839  {
840  0.,
841  0.,
842  0.,
843  0.,
844  0.,
845  Real(1.9412620368774630292701241080996842e-01L),
846  Real(9.8765911355712115933807754318089099e-02L),
847  Real(7.7663767064308164090246588765178087e-02L),
848  Real(2.1594628433980258573654682690950798e-02L),
849  Real(1.2563596287784997705599005477153617e-02L),
850  Real(1.5082654870922784345283124845552190e-02L)
851  };
852 
853  const unsigned int permutation_ids[n_wts]
854  = {1, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6}; // 49 total points
855 
856  dunavant_rule2(wts, a, b, permutation_ids, n_wts);
857 
858  return;
859  }
860 
861 
862 
863 
864  // Dunavant's 16th-order rule contains points outside the region of
865  // integration, and is thus unacceptable for our FEM calculations.
866  //
867  // This 55-point, 16th-order rule was obtained by me [JWP] using the code in
868  //
869  // Additional precision obtained from the code in:
870  // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
871  // on triangles and tetrahedra" Journal of Computational Mathematics,
872  // v. 27, no. 1, 2009, pp. 89-96.
873  //
874  // Note: the 55-point 16th-order rule obtained by Zhang in the paper above
875  // does not appear to be unique. It is a solution in the sense that it
876  // minimizes the error in the least-squares minimization problem, but
877  // it involves too many unknowns and the Jacobian is therefore singular
878  // when attempting to improve the solution via Newton's method.
879  case SIXTEENTH:
880  {
881  const unsigned int n_wts = 12;
882  const Real wts[n_wts] =
883  {
884  Real(2.2668082505910087151996321171534230e-02L),
885  Real(8.4043060714818596159798961899306135e-03L),
886  Real(1.0850949634049747713966288634484161e-03L),
887  Real(7.2252773375423638869298219383808751e-03L),
888  Real(1.2997715227338366024036316182572871e-02L),
889  Real(2.0054466616677715883228810959112227e-02L),
890  Real(9.7299841600417010281624372720122710e-03L),
891  Real(1.1651974438298104227427176444311766e-02L),
892  Real(9.1291185550484450744725847363097389e-03L),
893  Real(3.5568614040947150231712567900113671e-03L),
894  Real(5.8355861686234326181790822005304303e-03L),
895  Real(4.7411314396804228041879331486234396e-03L)
896  };
897 
898  const Real a[n_wts] =
899  {
900  0.0, // 'a' parameter not used for centroid weight
901  Real(8.5402539407933203673769900926355911e-02L),
902  Real(1.2425572001444092841183633409631260e-02L),
903  Real(4.9174838341891594024701017768490960e-01L),
904  Real(4.5669426695387464162068900231444462e-01L),
905  Real(4.8506759880447437974189793537259677e-01L),
906  Real(2.0622099278664205707909858461264083e-01L),
907  Real(3.2374950270039093446805340265853956e-01L),
908  Real(7.3834330556606586255186213302750029e-01L),
909  Real(9.1210673061680792565673823935174611e-01L),
910  Real(6.6129919222598721544966837350891531e-01L),
911  Real(1.7807138906021476039088828811346122e-01L)
912  };
913 
914  const Real b[n_wts] =
915  {
916  0.0,
917  0.0,
918  0.0,
919  0.0,
920  0.0,
921  Real(3.2315912848634384647700266402091638e-01L),
922  Real(1.5341553679414688425981898952416987e-01L),
923  Real(7.4295478991330687632977899141707872e-02L),
924  Real(7.1278762832147862035977841733532020e-02L),
925  Real(1.6623223223705792825395256602140459e-02L),
926  Real(1.4160772533794791868984026749196156e-02L),
927  Real(1.4539694958941854654807449467759690e-02L)
928  };
929 
930  const unsigned int permutation_ids[n_wts]
931  = {1, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6, 6}; // 55 total points
932 
933  dunavant_rule2(wts, a, b, permutation_ids, n_wts);
934 
935  return;
936  }
937 
938 
939  // Dunavant's 17th-order rule has 61 points, while a
940  // comparable conical product rule would have 81 (16th and 17th orders).
941  //
942  // It can be found here:
943  // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
944  //
945  // Zhang reports an identical rule in:
946  // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
947  // on triangles and tetrahedra" Journal of Computational Mathematics,
948  // v. 27, no. 1, 2009, pp. 89-96.
949  //
950  // Note: the 61-point 17th-order rule obtained by Dunavant and Zhang
951  // does not appear to be unique. It is a solution in the sense that it
952  // minimizes the error in the least-squares minimization problem, but
953  // it involves too many unknowns and the Jacobian is therefore singular
954  // when attempting to improve the solution via Newton's method.
955  //
956  // Therefore, we prefer the following 63-point rule which
957  // I [JWP] found. It appears to be more accurate than the
958  // rule reported by Dunavant and Zhang, even though it has
959  // a few more points.
960  case SEVENTEENTH:
961  {
962  const unsigned int n_wts = 12;
963  const Real wts[n_wts] =
964  {
965  Real(1.7464603792572004485690588092246146e-02L),
966  Real(5.9429003555801725246549713984660076e-03L),
967  Real(1.2490753345169579649319736639588729e-02L),
968  Real(1.5386987188875607593083456905596468e-02L),
969  Real(1.1185807311917706362674684312990270e-02L),
970  Real(1.0301845740670206831327304917180007e-02L),
971  Real(1.1767783072977049696840016810370464e-02L),
972  Real(3.8045312849431209558329128678945240e-03L),
973  Real(4.5139302178876351271037137230354382e-03L),
974  Real(2.2178812517580586419412547665472893e-03L),
975  Real(5.2216271537483672304731416553063103e-03L),
976  Real(9.8381136389470256422419930926212114e-04L)
977  };
978 
979  const Real a[n_wts] =
980  {
981  Real(2.8796825754667362165337965123570514e-01L),
982  Real(4.9216175986208465345536805750663939e-01L),
983  Real(4.6252866763171173685916780827044612e-01L),
984  Real(1.6730292951631792248498303276090273e-01L),
985  Real(1.5816335500814652972296428532213019e-01L),
986  Real(1.6352252138387564873002458959679529e-01L),
987  Real(6.2447680488959768233910286168417367e-01L),
988  Real(8.7317249935244454285263604347964179e-01L),
989  Real(3.4428164322282694677972239461699271e-01L),
990  Real(9.1584484467813674010523309855340209e-02L),
991  Real(2.0172088013378989086826623852040632e-01L),
992  Real(9.6538762758254643474731509845084691e-01L)
993  };
994 
995  const Real b[n_wts] =
996  {
997  0.0,
998  0.0,
999  0.0,
1000  Real(3.4429160695501713926320695771253348e-01L),
1001  Real(2.2541623431550639817203145525444726e-01L),
1002  Real(8.0670083153531811694942222940484991e-02L),
1003  Real(6.5967451375050925655738829747288190e-02L),
1004  Real(4.5677879890996762665044366994439565e-02L),
1005  Real(1.1528411723154215812386518751976084e-02L),
1006  Real(9.3057714323900610398389176844165892e-03L),
1007  Real(1.5916814107619812717966560404970160e-02L),
1008  Real(1.0734733163764032541125434215228937e-02L)
1009  };
1010 
1011  const unsigned int permutation_ids[n_wts]
1012  = {3, 3, 3, 6, 6, 6, 6, 6, 6, 6, 6, 6}; // 63 total points
1013 
1014  dunavant_rule2(wts, a, b, permutation_ids, n_wts);
1015 
1016  return;
1017 
1018  // _points.resize (61);
1019  // _weights.resize(61);
1020 
1021  // // The raw data for the quadrature rule.
1022  // const Real p[15][4] = {
1023  // { 1./3., 0., 0., 0.033437199290803e+00 / 2.0}, // 1-perm
1024  // {0.005658918886452e+00, 0.497170540556774e+00, 0., 0.005093415440507e+00 / 2.0}, // 3-perm
1025  // {0.035647354750751e+00, 0.482176322624625e+00, 0., 0.014670864527638e+00 / 2.0}, // 3-perm
1026  // {0.099520061958437e+00, 0.450239969020782e+00, 0., 0.024350878353672e+00 / 2.0}, // 3-perm
1027  // {0.199467521245206e+00, 0.400266239377397e+00, 0., 0.031107550868969e+00 / 2.0}, // 3-perm
1028  // {0.495717464058095e+00, 0.252141267970953e+00, 0., 0.031257111218620e+00 / 2.0}, // 3-perm
1029  // {0.675905990683077e+00, 0.162047004658461e+00, 0., 0.024815654339665e+00 / 2.0}, // 3-perm
1030  // {0.848248235478508e+00, 0.075875882260746e+00, 0., 0.014056073070557e+00 / 2.0}, // 3-perm
1031  // {0.968690546064356e+00, 0.015654726967822e+00, 0., 0.003194676173779e+00 / 2.0}, // 3-perm
1032  // {0.010186928826919e+00, 0.334319867363658e+00, 0.655493203809423e+00, 0.008119655318993e+00 / 2.0}, // 6-perm
1033  // {0.135440871671036e+00, 0.292221537796944e+00, 0.572337590532020e+00, 0.026805742283163e+00 / 2.0}, // 6-perm
1034  // {0.054423924290583e+00, 0.319574885423190e+00, 0.626001190286228e+00, 0.018459993210822e+00 / 2.0}, // 6-perm
1035  // {0.012868560833637e+00, 0.190704224192292e+00, 0.796427214974071e+00, 0.008476868534328e+00 / 2.0}, // 6-perm
1036  // {0.067165782413524e+00, 0.180483211648746e+00, 0.752351005937729e+00, 0.018292796770025e+00 / 2.0}, // 6-perm
1037  // {0.014663182224828e+00, 0.080711313679564e+00, 0.904625504095608e+00, 0.006665632004165e+00 / 2.0} // 6-perm
1038  // };
1039 
1040 
1041  // // Now call the dunavant routine to generate _points and _weights
1042  // dunavant_rule(p, 15);
1043 
1044  // return;
1045  }
1046 
1047 
1048 
1049  // Dunavant's 18th-order rule contains points outside the region and is therefore unsuitable
1050  // for our FEM calculations. His 19th-order rule has 73 points, compared with 100 points for
1051  // a comparable-order conical product rule.
1052  //
1053  // It was copied 23rd June 2008 from:
1054  // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
1055  case EIGHTTEENTH:
1056  case NINETEENTH:
1057  {
1058  _points.resize (73);
1059  _weights.resize(73);
1060 
1061  // The raw data for the quadrature rule.
1062  const Real rule_data[17][4] = {
1063  { 1./3., 0., 0., 0.032906331388919e+00 / 2.0}, // 1-perm
1064  {0.020780025853987e+00, 0.489609987073006e+00, 0., 0.010330731891272e+00 / 2.0}, // 3-perm
1065  {0.090926214604215e+00, 0.454536892697893e+00, 0., 0.022387247263016e+00 / 2.0}, // 3-perm
1066  {0.197166638701138e+00, 0.401416680649431e+00, 0., 0.030266125869468e+00 / 2.0}, // 3-perm
1067  {0.488896691193805e+00, 0.255551654403098e+00, 0., 0.030490967802198e+00 / 2.0}, // 3-perm
1068  {0.645844115695741e+00, 0.177077942152130e+00, 0., 0.024159212741641e+00 / 2.0}, // 3-perm
1069  {0.779877893544096e+00, 0.110061053227952e+00, 0., 0.016050803586801e+00 / 2.0}, // 3-perm
1070  {0.888942751496321e+00, 0.055528624251840e+00, 0., 0.008084580261784e+00 / 2.0}, // 3-perm
1071  {0.974756272445543e+00, 0.012621863777229e+00, 0., 0.002079362027485e+00 / 2.0}, // 3-perm
1072  {0.003611417848412e+00, 0.395754787356943e+00, 0.600633794794645e+00, 0.003884876904981e+00 / 2.0}, // 6-perm
1073  {0.134466754530780e+00, 0.307929983880436e+00, 0.557603261588784e+00, 0.025574160612022e+00 / 2.0}, // 6-perm
1074  {0.014446025776115e+00, 0.264566948406520e+00, 0.720987025817365e+00, 0.008880903573338e+00 / 2.0}, // 6-perm
1075  {0.046933578838178e+00, 0.358539352205951e+00, 0.594527068955871e+00, 0.016124546761731e+00 / 2.0}, // 6-perm
1076  {0.002861120350567e+00, 0.157807405968595e+00, 0.839331473680839e+00, 0.002491941817491e+00 / 2.0}, // 6-perm
1077  {0.223861424097916e+00, 0.075050596975911e+00, 0.701087978926173e+00, 0.018242840118951e+00 / 2.0}, // 6-perm
1078  {0.034647074816760e+00, 0.142421601113383e+00, 0.822931324069857e+00, 0.010258563736199e+00 / 2.0}, // 6-perm
1079  {0.010161119296278e+00, 0.065494628082938e+00, 0.924344252620784e+00, 0.003799928855302e+00 / 2.0} // 6-perm
1080  };
1081 
1082 
1083  // Now call the dunavant routine to generate _points and _weights
1084  dunavant_rule(rule_data, 17);
1085 
1086  return;
1087  }
1088 
1089 
1090  // 20th-order rule by Wandzura.
1091  //
1092  // Stephen Wandzura, Hong Xiao,
1093  // Symmetric Quadrature Rules on a Triangle,
1094  // Computers and Mathematics with Applications,
1095  // Volume 45, Number 12, June 2003, pages 1829-1840.
1096  //
1097  // Wandzura's work extends the work of Dunavant by providing degree
1098  // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
1099  //
1100  // Copied on 3rd July 2008 from:
1101  // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
1102  case TWENTIETH:
1103  {
1104  // The equivalent conical product rule would have 121 points
1105  _points.resize (85);
1106  _weights.resize(85);
1107 
1108  // The raw data for the quadrature rule.
1109  const Real rule_data[19][4] = {
1110  {0.33333333333333e+00, 0.0, 0.0, 0.2761042699769952e-01 / 2.0}, // 1-perm
1111  {0.00150064932443e+00, 0.49924967533779e+00, 0.0, 0.1779029547326740e-02 / 2.0}, // 3-perm
1112  {0.09413975193895e+00, 0.45293012403052e+00, 0.0, 0.2011239811396117e-01 / 2.0}, // 3-perm
1113  {0.20447212408953e+00, 0.39776393795524e+00, 0.0, 0.2681784725933157e-01 / 2.0}, // 3-perm
1114  {0.47099959493443e+00, 0.26450020253279e+00, 0.0, 0.2452313380150201e-01 / 2.0}, // 3-perm
1115  {0.57796207181585e+00, 0.21101896409208e+00, 0.0, 0.1639457841069539e-01 / 2.0}, // 3-perm
1116  {0.78452878565746e+00, 0.10773560717127e+00, 0.0, 0.1479590739864960e-01 / 2.0}, // 3-perm
1117  {0.92186182432439e+00, 0.03906908783780e+00, 0.0, 0.4579282277704251e-02 / 2.0}, // 3-perm
1118  {0.97765124054134e+00, 0.01117437972933e+00, 0.0, 0.1651826515576217e-02 / 2.0}, // 3-perm
1119  {0.00534961818734e+00, 0.06354966590835e+00, 0.93110071590431e+00, 0.2349170908575584e-02 / 2.0}, // 6-perm
1120  {0.00795481706620e+00, 0.15710691894071e+00, 0.83493826399309e+00, 0.4465925754181793e-02 / 2.0}, // 6-perm
1121  {0.01042239828126e+00, 0.39564211436437e+00, 0.59393548735436e+00, 0.6099566807907972e-02 / 2.0}, // 6-perm
1122  {0.01096441479612e+00, 0.27316757071291e+00, 0.71586801449097e+00, 0.6891081327188203e-02 / 2.0}, // 6-perm
1123  {0.03856671208546e+00, 0.10178538248502e+00, 0.85964790542952e+00, 0.7997475072478163e-02 / 2.0}, // 6-perm
1124  {0.03558050781722e+00, 0.44665854917641e+00, 0.51776094300637e+00, 0.7386134285336024e-02 / 2.0}, // 6-perm
1125  {0.04967081636276e+00, 0.19901079414950e+00, 0.75131838948773e+00, 0.1279933187864826e-01 / 2.0}, // 6-perm
1126  {0.05851972508433e+00, 0.32426118369228e+00, 0.61721909122339e+00, 0.1725807117569655e-01 / 2.0}, // 6-perm
1127  {0.12149778700439e+00, 0.20853136321013e+00, 0.66997084978547e+00, 0.1867294590293547e-01 / 2.0}, // 6-perm
1128  {0.14071084494394e+00, 0.32317056653626e+00, 0.53611858851980e+00, 0.2281822405839526e-01 / 2.0} // 6-perm
1129  };
1130 
1131 
1132  // Now call the dunavant routine to generate _points and _weights
1133  dunavant_rule(rule_data, 19);
1134 
1135  return;
1136  }
1137 
1138 
1139 
1140  // 25th-order rule by Wandzura.
1141  //
1142  // Stephen Wandzura, Hong Xiao,
1143  // Symmetric Quadrature Rules on a Triangle,
1144  // Computers and Mathematics with Applications,
1145  // Volume 45, Number 12, June 2003, pages 1829-1840.
1146  //
1147  // Wandzura's work extends the work of Dunavant by providing degree
1148  // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
1149  //
1150  // Copied on 3rd July 2008 from:
1151  // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
1152  // case TWENTYFIRST: // fall through to 121 point conical product rule below
1153  case TWENTYSECOND:
1154  case TWENTYTHIRD:
1155  case TWENTYFOURTH:
1156  case TWENTYFIFTH:
1157  {
1158  // The equivalent conical product rule would have 169 points
1159  _points.resize (126);
1160  _weights.resize(126);
1161 
1162  // The raw data for the quadrature rule.
1163  const Real rule_data[26][4] = {
1164  {0.02794648307317e+00, 0.48602675846341e+00, 0.0, 0.8005581880020417e-02 / 2.0}, // 3-perm
1165  {0.13117860132765e+00, 0.43441069933617e+00, 0.0, 0.1594707683239050e-01 / 2.0}, // 3-perm
1166  {0.22022172951207e+00, 0.38988913524396e+00, 0.0, 0.1310914123079553e-01 / 2.0}, // 3-perm
1167  {0.40311353196039e+00, 0.29844323401980e+00, 0.0, 0.1958300096563562e-01 / 2.0}, // 3-perm
1168  {0.53191165532526e+00, 0.23404417233737e+00, 0.0, 0.1647088544153727e-01 / 2.0}, // 3-perm
1169  {0.69706333078196e+00, 0.15146833460902e+00, 0.0, 0.8547279074092100e-02 / 2.0}, // 3-perm
1170  {0.77453221290801e+00, 0.11273389354599e+00, 0.0, 0.8161885857226492e-02 / 2.0}, // 3-perm
1171  {0.84456861581695e+00, 0.07771569209153e+00, 0.0, 0.6121146539983779e-02 / 2.0}, // 3-perm
1172  {0.93021381277141e+00, 0.03489309361430e+00, 0.0, 0.2908498264936665e-02 / 2.0}, // 3-perm
1173  {0.98548363075813e+00, 0.00725818462093e+00, 0.0, 0.6922752456619963e-03 / 2.0}, // 3-perm
1174  {0.00129235270444e+00, 0.22721445215336e+00, 0.77149319514219e+00, 0.1248289199277397e-02 / 2.0}, // 6-perm
1175  {0.00539970127212e+00, 0.43501055485357e+00, 0.55958974387431e+00, 0.3404752908803022e-02 / 2.0}, // 6-perm
1176  {0.00638400303398e+00, 0.32030959927220e+00, 0.67330639769382e+00, 0.3359654326064051e-02 / 2.0}, // 6-perm
1177  {0.00502821150199e+00, 0.09175032228001e+00, 0.90322146621800e+00, 0.1716156539496754e-02 / 2.0}, // 6-perm
1178  {0.00682675862178e+00, 0.03801083585872e+00, 0.95516240551949e+00, 0.1480856316715606e-02 / 2.0}, // 6-perm
1179  {0.01001619963993e+00, 0.15742521848531e+00, 0.83255858187476e+00, 0.3511312610728685e-02 / 2.0}, // 6-perm
1180  {0.02575781317339e+00, 0.23988965977853e+00, 0.73435252704808e+00, 0.7393550149706484e-02 / 2.0}, // 6-perm
1181  {0.03022789811992e+00, 0.36194311812606e+00, 0.60782898375402e+00, 0.7983087477376558e-02 / 2.0}, // 6-perm
1182  {0.03050499010716e+00, 0.08355196095483e+00, 0.88594304893801e+00, 0.4355962613158041e-02 / 2.0}, // 6-perm
1183  {0.04595654736257e+00, 0.14844322073242e+00, 0.80560023190501e+00, 0.7365056701417832e-02 / 2.0}, // 6-perm
1184  {0.06744280054028e+00, 0.28373970872753e+00, 0.64881749073219e+00, 0.1096357284641955e-01 / 2.0}, // 6-perm
1185  {0.07004509141591e+00, 0.40689937511879e+00, 0.52305553346530e+00, 0.1174996174354112e-01 / 2.0}, // 6-perm
1186  {0.08391152464012e+00, 0.19411398702489e+00, 0.72197448833499e+00, 0.1001560071379857e-01 / 2.0}, // 6-perm
1187  {0.12037553567715e+00, 0.32413434700070e+00, 0.55549011732214e+00, 0.1330964078762868e-01 / 2.0}, // 6-perm
1188  {0.14806689915737e+00, 0.22927748355598e+00, 0.62265561728665e+00, 0.1415444650522614e-01 / 2.0}, // 6-perm
1189  {0.19177186586733e+00, 0.32561812259598e+00, 0.48261001153669e+00, 0.1488137956116801e-01 / 2.0} // 6-perm
1190  };
1191 
1192 
1193  // Now call the dunavant routine to generate _points and _weights
1194  dunavant_rule(rule_data, 26);
1195 
1196  return;
1197  }
1198 
1199 
1200 
1201  // 30th-order rule by Wandzura.
1202  //
1203  // Stephen Wandzura, Hong Xiao,
1204  // Symmetric Quadrature Rules on a Triangle,
1205  // Computers and Mathematics with Applications,
1206  // Volume 45, Number 12, June 2003, pages 1829-1840.
1207  //
1208  // Wandzura's work extends the work of Dunavant by providing degree
1209  // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
1210  //
1211  // Copied on 3rd July 2008 from:
1212  // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
1213  case TWENTYSIXTH:
1214  case TWENTYSEVENTH:
1215  case TWENTYEIGHTH:
1216  case TWENTYNINTH:
1217  case THIRTIETH:
1218  {
1219  // The equivalent conical product rule would have 256 points
1220  _points.resize (175);
1221  _weights.resize(175);
1222 
1223  // The raw data for the quadrature rule.
1224  const Real rule_data[36][4] = {
1225  {0.33333333333333e+00, 0.0, 0.0, 0.1557996020289920e-01 / 2.0}, // 1-perm
1226  {0.00733011643277e+00, 0.49633494178362e+00, 0.0, 0.3177233700534134e-02 / 2.0}, // 3-perm
1227  {0.08299567580296e+00, 0.45850216209852e+00, 0.0, 0.1048342663573077e-01 / 2.0}, // 3-perm
1228  {0.15098095612541e+00, 0.42450952193729e+00, 0.0, 0.1320945957774363e-01 / 2.0}, // 3-perm
1229  {0.23590585989217e+00, 0.38204707005392e+00, 0.0, 0.1497500696627150e-01 / 2.0}, // 3-perm
1230  {0.43802430840785e+00, 0.28098784579608e+00, 0.0, 0.1498790444338419e-01 / 2.0}, // 3-perm
1231  {0.54530204829193e+00, 0.22734897585403e+00, 0.0, 0.1333886474102166e-01 / 2.0}, // 3-perm
1232  {0.65088177698254e+00, 0.17455911150873e+00, 0.0, 0.1088917111390201e-01 / 2.0}, // 3-perm
1233  {0.75348314559713e+00, 0.12325842720144e+00, 0.0, 0.8189440660893461e-02 / 2.0}, // 3-perm
1234  {0.83983154221561e+00, 0.08008422889220e+00, 0.0, 0.5575387588607785e-02 / 2.0}, // 3-perm
1235  {0.90445106518420e+00, 0.04777446740790e+00, 0.0, 0.3191216473411976e-02 / 2.0}, // 3-perm
1236  {0.95655897063972e+00, 0.02172051468014e+00, 0.0, 0.1296715144327045e-02 / 2.0}, // 3-perm
1237  {0.99047064476913e+00, 0.00476467761544e+00, 0.0, 0.2982628261349172e-03 / 2.0}, // 3-perm
1238  {0.00092537119335e+00, 0.41529527091331e+00, 0.58377935789334e+00, 0.9989056850788964e-03 / 2.0}, // 6-perm
1239  {0.00138592585556e+00, 0.06118990978535e+00, 0.93742416435909e+00, 0.4628508491732533e-03 / 2.0}, // 6-perm
1240  {0.00368241545591e+00, 0.16490869013691e+00, 0.83140889440718e+00, 0.1234451336382413e-02 / 2.0}, // 6-perm
1241  {0.00390322342416e+00, 0.02503506223200e+00, 0.97106171434384e+00, 0.5707198522432062e-03 / 2.0}, // 6-perm
1242  {0.00323324815501e+00, 0.30606446515110e+00, 0.69070228669389e+00, 0.1126946125877624e-02 / 2.0}, // 6-perm
1243  {0.00646743211224e+00, 0.10707328373022e+00, 0.88645928415754e+00, 0.1747866949407337e-02 / 2.0}, // 6-perm
1244  {0.00324747549133e+00, 0.22995754934558e+00, 0.76679497516308e+00, 0.1182818815031657e-02 / 2.0}, // 6-perm
1245  {0.00867509080675e+00, 0.33703663330578e+00, 0.65428827588746e+00, 0.1990839294675034e-02 / 2.0}, // 6-perm
1246  {0.01559702646731e+00, 0.05625657618206e+00, 0.92814639735063e+00, 0.1900412795035980e-02 / 2.0}, // 6-perm
1247  {0.01797672125369e+00, 0.40245137521240e+00, 0.57957190353391e+00, 0.4498365808817451e-02 / 2.0}, // 6-perm
1248  {0.01712424535389e+00, 0.24365470201083e+00, 0.73922105263528e+00, 0.3478719460274719e-02 / 2.0}, // 6-perm
1249  {0.02288340534658e+00, 0.16538958561453e+00, 0.81172700903888e+00, 0.4102399036723953e-02 / 2.0}, // 6-perm
1250  {0.03273759728777e+00, 0.09930187449585e+00, 0.86796052821639e+00, 0.4021761549744162e-02 / 2.0}, // 6-perm
1251  {0.03382101234234e+00, 0.30847833306905e+00, 0.65770065458860e+00, 0.6033164660795066e-02 / 2.0}, // 6-perm
1252  {0.03554761446002e+00, 0.46066831859211e+00, 0.50378406694787e+00, 0.3946290302129598e-02 / 2.0}, // 6-perm
1253  {0.05053979030687e+00, 0.21881529945393e+00, 0.73064491023920e+00, 0.6644044537680268e-02 / 2.0}, // 6-perm
1254  {0.05701471491573e+00, 0.37920955156027e+00, 0.56377573352399e+00, 0.8254305856078458e-02 / 2.0}, // 6-perm
1255  {0.06415280642120e+00, 0.14296081941819e+00, 0.79288637416061e+00, 0.6496056633406411e-02 / 2.0}, // 6-perm
1256  {0.08050114828763e+00, 0.28373128210592e+00, 0.63576756960645e+00, 0.9252778144146602e-02 / 2.0}, // 6-perm
1257  {0.10436706813453e+00, 0.19673744100444e+00, 0.69889549086103e+00, 0.9164920726294280e-02 / 2.0}, // 6-perm
1258  {0.11384489442875e+00, 0.35588914121166e+00, 0.53026596435959e+00, 0.1156952462809767e-01 / 2.0}, // 6-perm
1259  {0.14536348771552e+00, 0.25981868535191e+00, 0.59481782693256e+00, 0.1176111646760917e-01 / 2.0}, // 6-perm
1260  {0.18994565282198e+00, 0.32192318123130e+00, 0.48813116594672e+00, 0.1382470218216540e-01 / 2.0} // 6-perm
1261  };
1262 
1263 
1264  // Now call the dunavant routine to generate _points and _weights
1265  dunavant_rule(rule_data, 36);
1266 
1267  return;
1268  }
1269 
1270 
1271  // By default, we fall back on the conical product rules. If the user
1272  // requests an order higher than what is currently available in the 1D
1273  // rules, an error will be thrown from the respective 1D code.
1274  default:
1275  {
1276  // The following quadrature rules are generated as
1277  // conical products. These tend to be non-optimal
1278  // (use too many points, cluster points in certain
1279  // regions of the domain) but they are quite easy to
1280  // automatically generate using a 1D Gauss rule on
1281  // [0,1] and two 1D Jacobi-Gauss rules on [0,1].
1282  QConical conical_rule(2, _order);
1283  conical_rule.init(_type, _p_level);
1284 
1285  // Swap points and weights with the about-to-be destroyed rule.
1286  _points.swap (conical_rule.get_points() );
1287  _weights.swap(conical_rule.get_weights());
1288 
1289  return;
1290  }
1291  }
1292  }
1293 
1294 
1295  //---------------------------------------------
1296  // Unsupported type
1297  default:
1298  libmesh_error_msg("Element type not supported:" << Utility::enum_to_string(_type));
1299  }
1300 #endif
1301 }

References libMesh::QBase::_order, libMesh::QBase::_p_level, libMesh::QBase::_points, libMesh::QBase::_type, libMesh::QBase::_weights, libMesh::CONSTANT, dunavant_rule(), dunavant_rule2(), libMesh::EDGE2, libMesh::EIGHTH, libMesh::EIGHTTEENTH, libMesh::ELEVENTH, libMesh::Utility::enum_to_string(), libMesh::FIFTEENTH, libMesh::FIFTH, libMesh::FIRST, libMesh::FOURTEENTH, libMesh::FOURTH, libMesh::QBase::get_order(), libMesh::QBase::get_points(), libMesh::QBase::get_weights(), libMesh::QBase::init(), libMesh::NINETEENTH, libMesh::NINTH, libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::QUADSHELL4, libMesh::QUADSHELL8, libMesh::Real, libMesh::SECOND, libMesh::SEVENTEENTH, libMesh::SEVENTH, libMesh::SIXTEENTH, libMesh::SIXTH, std::sqrt(), libMesh::QBase::tensor_product_quad(), libMesh::TENTH, libMesh::THIRD, libMesh::THIRTEENTH, libMesh::THIRTIETH, libMesh::TRI3, libMesh::TRI3SUBDIVISION, libMesh::TRI6, libMesh::TRISHELL3, libMesh::TWELFTH, libMesh::TWENTIETH, libMesh::TWENTYEIGHTH, libMesh::TWENTYFIFTH, libMesh::TWENTYFOURTH, libMesh::TWENTYNINTH, libMesh::TWENTYSECOND, libMesh::TWENTYSEVENTH, libMesh::TWENTYSIXTH, and libMesh::TWENTYTHIRD.

◆ init_3D()

void libMesh::QGauss::init_3D ( const  type,
unsigned int  p_level 
)
overrideprivatevirtual

Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate values.

The order of the rule will be defined by the implementing class. Should not be pure virtual since a derived quadrature rule may only be defined in 1D. If not overridden, throws an error.

Note
The arguments should no longer be used for anything in derived classes, they are only maintained for backwards compatibility and will eventually be removed.

Reimplemented from libMesh::QBase.

Definition at line 29 of file quadrature_gauss_3D.C.

30 {
31 #if LIBMESH_DIM == 3
32 
33  //-----------------------------------------------------------------------
34  // 3D quadrature rules
35  switch (_type)
36  {
37  //---------------------------------------------
38  // Hex quadrature rules
39  case HEX8:
40  case HEX20:
41  case HEX27:
42  {
43  // We compute the 3D quadrature rule as a tensor
44  // product of the 1D quadrature rule.
45  QGauss q1D(1, _order);
46  q1D.init(EDGE2, _p_level);
47  tensor_product_hex( q1D );
48  return;
49  }
50 
51 
52 
53  //---------------------------------------------
54  // Tetrahedral quadrature rules
55  case TET4:
56  case TET10:
57  {
58  switch(get_order())
59  {
60  // Taken from pg. 222 of "The finite element method," vol. 1
61  // ed. 5 by Zienkiewicz & Taylor
62  case CONSTANT:
63  case FIRST:
64  {
65  // Exact for linears
66  _points.resize(1);
67  _weights.resize(1);
68 
69 
70  _points[0](0) = .25;
71  _points[0](1) = .25;
72  _points[0](2) = .25;
73 
74  _weights[0] = Real(1)/6;
75 
76  return;
77  }
78  case SECOND:
79  {
80  // Exact for quadratics
81  _points.resize(4);
82  _weights.resize(4);
83 
84 
85  const Real a = .585410196624969;
86  const Real b = .138196601125011;
87 
88  _points[0](0) = a;
89  _points[0](1) = b;
90  _points[0](2) = b;
91 
92  _points[1](0) = b;
93  _points[1](1) = a;
94  _points[1](2) = b;
95 
96  _points[2](0) = b;
97  _points[2](1) = b;
98  _points[2](2) = a;
99 
100  _points[3](0) = b;
101  _points[3](1) = b;
102  _points[3](2) = b;
103 
104 
105 
106  _weights[0] = Real(1)/24;
107  _weights[1] = _weights[0];
108  _weights[2] = _weights[0];
109  _weights[3] = _weights[0];
110 
111  return;
112  }
113 
114 
115 
116  // Can be found in the class notes
117  // http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps
118  // by Flaherty.
119  //
120  // Caution: this rule has a negative weight and may be
121  // unsuitable for some problems.
122  // Exact for cubics.
123  //
124  // Note: Keast (see ref. elsewhere in this file) also gives
125  // a third-order rule with positive weights, but it contains points
126  // on the ref. elt. boundary, making it less suitable for FEM calculations.
127  case THIRD:
128  {
129  if (allow_rules_with_negative_weights)
130  {
131  _points.resize(5);
132  _weights.resize(5);
133 
134 
135  _points[0](0) = .25;
136  _points[0](1) = .25;
137  _points[0](2) = .25;
138 
139  _points[1](0) = .5;
140  _points[1](1) = Real(1)/6;
141  _points[1](2) = Real(1)/6;
142 
143  _points[2](0) = Real(1)/6;
144  _points[2](1) = .5;
145  _points[2](2) = Real(1)/6;
146 
147  _points[3](0) = Real(1)/6;
148  _points[3](1) = Real(1)/6;
149  _points[3](2) = .5;
150 
151  _points[4](0) = Real(1)/6;
152  _points[4](1) = Real(1)/6;
153  _points[4](2) = Real(1)/6;
154 
155 
156  _weights[0] = Real(-2)/15;
157  _weights[1] = .075;
158  _weights[2] = _weights[1];
159  _weights[3] = _weights[1];
160  _weights[4] = _weights[1];
161 
162  return;
163  } // end if (allow_rules_with_negative_weights)
164  else
165  {
166  // If a rule with positive weights is required, the 2x2x2 conical
167  // product rule is third-order accurate and has less points than
168  // the next-available positive-weight rule at FIFTH order.
169  QConical conical_rule(3, _order);
170  conical_rule.init(_type, _p_level);
171 
172  // Swap points and weights with the about-to-be destroyed rule.
173  _points.swap (conical_rule.get_points() );
174  _weights.swap(conical_rule.get_weights());
175 
176  return;
177  }
178  // Note: if !allow_rules_with_negative_weights, fall through to next case.
179  }
180 
181 
182 
183  // Originally a Keast rule,
184  // Patrick Keast,
185  // Moderate Degree Tetrahedral Quadrature Formulas,
186  // Computer Methods in Applied Mechanics and Engineering,
187  // Volume 55, Number 3, May 1986, pages 339-348.
188  //
189  // Can also be found the class notes
190  // http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps
191  // by Flaherty.
192  //
193  // Caution: this rule has a negative weight and may be
194  // unsuitable for some problems.
195  case FOURTH:
196  {
197  if (allow_rules_with_negative_weights)
198  {
199  _points.resize(11);
200  _weights.resize(11);
201 
202  // The raw data for the quadrature rule.
203  const Real rule_data[3][4] = {
204  {0.250000000000000000e+00, 0., 0., -0.131555555555555556e-01}, // 1
205  {0.785714285714285714e+00, 0.714285714285714285e-01, 0., 0.762222222222222222e-02}, // 4
206  {0.399403576166799219e+00, 0., 0.100596423833200785e+00, 0.248888888888888889e-01} // 6
207  };
208 
209 
210  // Now call the keast routine to generate _points and _weights
211  keast_rule(rule_data, 3);
212 
213  return;
214  } // end if (allow_rules_with_negative_weights)
215  // Note: if !allow_rules_with_negative_weights, fall through to next case.
216  }
217 
218  libmesh_fallthrough();
219 
220 
221  // Walkington's fifth-order 14-point rule from
222  // "Quadrature on Simplices of Arbitrary Dimension"
223  //
224  // We originally had a Keast rule here, but this rule had
225  // more points than an equivalent rule by Walkington and
226  // also contained points on the boundary of the ref. elt,
227  // making it less suitable for FEM calculations.
228  case FIFTH:
229  {
230  _points.resize(14);
231  _weights.resize(14);
232 
233  // permutations of these points and suitably-modified versions of
234  // these points are the quadrature point locations
235  const Real a[3] = {0.31088591926330060980, // a1 from the paper
236  0.092735250310891226402, // a2 from the paper
237  0.045503704125649649492}; // a3 from the paper
238 
239  // weights. a[] and wt[] are the only floating-point inputs required
240  // for this rule.
241  const Real wt[3] = {0.018781320953002641800, // w1 from the paper
242  0.012248840519393658257, // w2 from the paper
243  0.0070910034628469110730}; // w3 from the paper
244 
245  // The first two sets of 4 points are formed in a similar manner
246  for (unsigned int i=0; i<2; ++i)
247  {
248  // Where we will insert values into _points and _weights
249  const unsigned int offset=4*i;
250 
251  // Stuff points and weights values into their arrays
252  const Real b = 1. - 3.*a[i];
253 
254  // Here are the permutations. Order of these is not important,
255  // all have the same weight
256  _points[offset + 0] = Point(a[i], a[i], a[i]);
257  _points[offset + 1] = Point(a[i], b, a[i]);
258  _points[offset + 2] = Point( b, a[i], a[i]);
259  _points[offset + 3] = Point(a[i], a[i], b);
260 
261  // These 4 points all have the same weights
262  for (unsigned int j=0; j<4; ++j)
263  _weights[offset + j] = wt[i];
264  } // end for
265 
266 
267  {
268  // The third set contains 6 points and is formed a little differently
269  const unsigned int offset = 8;
270  const Real b = 0.5*(1. - 2.*a[2]);
271 
272  // Here are the permutations. Order of these is not important,
273  // all have the same weight
274  _points[offset + 0] = Point(b , b, a[2]);
275  _points[offset + 1] = Point(b , a[2], a[2]);
276  _points[offset + 2] = Point(a[2], a[2], b);
277  _points[offset + 3] = Point(a[2], b, a[2]);
278  _points[offset + 4] = Point( b, a[2], b);
279  _points[offset + 5] = Point(a[2], b, b);
280 
281  // These 6 points all have the same weights
282  for (unsigned int j=0; j<6; ++j)
283  _weights[offset + j] = wt[2];
284  }
285 
286 
287  // Original rule by Keast, unsuitable because it has points on the
288  // reference element boundary!
289  // _points.resize(15);
290  // _weights.resize(15);
291 
292  // _points[0](0) = 0.25;
293  // _points[0](1) = 0.25;
294  // _points[0](2) = 0.25;
295 
296  // {
297  // const Real a = 0.;
298  // const Real b = Real(1)/3;
299 
300  // _points[1](0) = a;
301  // _points[1](1) = b;
302  // _points[1](2) = b;
303 
304  // _points[2](0) = b;
305  // _points[2](1) = a;
306  // _points[2](2) = b;
307 
308  // _points[3](0) = b;
309  // _points[3](1) = b;
310  // _points[3](2) = a;
311 
312  // _points[4](0) = b;
313  // _points[4](1) = b;
314  // _points[4](2) = b;
315  // }
316  // {
317  // const Real a = Real(8)/11;
318  // const Real b = Real(1)/11;
319 
320  // _points[5](0) = a;
321  // _points[5](1) = b;
322  // _points[5](2) = b;
323 
324  // _points[6](0) = b;
325  // _points[6](1) = a;
326  // _points[6](2) = b;
327 
328  // _points[7](0) = b;
329  // _points[7](1) = b;
330  // _points[7](2) = a;
331 
332  // _points[8](0) = b;
333  // _points[8](1) = b;
334  // _points[8](2) = b;
335  // }
336  // {
337  // const Real a = 0.066550153573664;
338  // const Real b = 0.433449846426336;
339 
340  // _points[9](0) = b;
341  // _points[9](1) = a;
342  // _points[9](2) = a;
343 
344  // _points[10](0) = a;
345  // _points[10](1) = a;
346  // _points[10](2) = b;
347 
348  // _points[11](0) = a;
349  // _points[11](1) = b;
350  // _points[11](2) = b;
351 
352  // _points[12](0) = b;
353  // _points[12](1) = a;
354  // _points[12](2) = b;
355 
356  // _points[13](0) = b;
357  // _points[13](1) = b;
358  // _points[13](2) = a;
359 
360  // _points[14](0) = a;
361  // _points[14](1) = b;
362  // _points[14](2) = a;
363  // }
364 
365  // _weights[0] = 0.030283678097089;
366  // _weights[1] = 0.006026785714286;
367  // _weights[2] = _weights[1];
368  // _weights[3] = _weights[1];
369  // _weights[4] = _weights[1];
370  // _weights[5] = 0.011645249086029;
371  // _weights[6] = _weights[5];
372  // _weights[7] = _weights[5];
373  // _weights[8] = _weights[5];
374  // _weights[9] = 0.010949141561386;
375  // _weights[10] = _weights[9];
376  // _weights[11] = _weights[9];
377  // _weights[12] = _weights[9];
378  // _weights[13] = _weights[9];
379  // _weights[14] = _weights[9];
380 
381  return;
382  }
383 
384 
385 
386 
387  // This rule is originally from Keast:
388  // Patrick Keast,
389  // Moderate Degree Tetrahedral Quadrature Formulas,
390  // Computer Methods in Applied Mechanics and Engineering,
391  // Volume 55, Number 3, May 1986, pages 339-348.
392  //
393  // It is accurate on 6th-degree polynomials and has 24 points
394  // vs. 64 for the comparable conical product rule.
395  //
396  // Values copied 24th June 2008 from:
397  // http://people.scs.fsu.edu/~burkardt/f_src/keast/keast.f90
398  case SIXTH:
399  {
400  _points.resize (24);
401  _weights.resize(24);
402 
403  // The raw data for the quadrature rule.
404  const Real rule_data[4][4] = {
405  {0.356191386222544953e+00 , 0.214602871259151684e+00 , 0., 0.00665379170969464506e+00}, // 4
406  {0.877978124396165982e+00 , 0.0406739585346113397e+00, 0., 0.00167953517588677620e+00}, // 4
407  {0.0329863295731730594e+00, 0.322337890142275646e+00 , 0., 0.00922619692394239843e+00}, // 4
408  {0.0636610018750175299e+00, 0.269672331458315867e+00 , 0.603005664791649076e+00, 0.00803571428571428248e+00} // 12
409  };
410 
411 
412  // Now call the keast routine to generate _points and _weights
413  keast_rule(rule_data, 4);
414 
415  return;
416  }
417 
418 
419  // Keast's 31 point, 7th-order rule contains points on the reference
420  // element boundary, so we've decided not to include it here.
421  //
422  // Keast's 8th-order rule has 45 points. and a negative
423  // weight, so if you've explicitly disallowed such rules
424  // you will fall through to the conical product rules
425  // below.
426  case SEVENTH:
427  case EIGHTH:
428  {
429  if (allow_rules_with_negative_weights)
430  {
431  _points.resize (45);
432  _weights.resize(45);
433 
434  // The raw data for the quadrature rule.
435  const Real rule_data[7][4] = {
436  {0.250000000000000000e+00, 0., 0., -0.393270066412926145e-01}, // 1
437  {0.617587190300082967e+00, 0.127470936566639015e+00, 0., 0.408131605934270525e-02}, // 4
438  {0.903763508822103123e+00, 0.320788303926322960e-01, 0., 0.658086773304341943e-03}, // 4
439  {0.497770956432810185e-01, 0., 0.450222904356718978e+00, 0.438425882512284693e-02}, // 6
440  {0.183730447398549945e+00, 0., 0.316269552601450060e+00, 0.138300638425098166e-01}, // 6
441  {0.231901089397150906e+00, 0.229177878448171174e-01, 0.513280033360881072e+00, 0.424043742468372453e-02}, // 12
442  {0.379700484718286102e-01, 0.730313427807538396e+00, 0.193746475248804382e+00, 0.223873973961420164e-02} // 12
443  };
444 
445 
446  // Now call the keast routine to generate _points and _weights
447  keast_rule(rule_data, 7);
448 
449  return;
450  } // end if (allow_rules_with_negative_weights)
451  // Note: if !allow_rules_with_negative_weights, fall through to next case.
452  }
453 
454  libmesh_fallthrough();
455 
456 
457  // Fall back on Grundmann-Moller or Conical Product rules at high orders.
458  default:
459  {
460  if ((allow_rules_with_negative_weights) && (get_order() < 34))
461  {
462  // The Grundmann-Moller rules are defined to arbitrary order and
463  // can have significantly fewer evaluation points than conical product
464  // rules. If you allow rules with negative weights, the GM rules
465  // will be more efficient for degree up to 33 (for degree 34 and
466  // higher, CP is more efficient!) but may be more susceptible
467  // to round-off error. Safest is to disallow rules with negative
468  // weights, but this decision should be made on a case-by-case basis.
469  QGrundmann_Moller gm_rule(3, _order);
470  gm_rule.init(_type, _p_level);
471 
472  // Swap points and weights with the about-to-be destroyed rule.
473  _points.swap (gm_rule.get_points() );
474  _weights.swap(gm_rule.get_weights());
475 
476  return;
477  }
478 
479  else
480  {
481  // The following quadrature rules are generated as
482  // conical products. These tend to be non-optimal
483  // (use too many points, cluster points in certain
484  // regions of the domain) but they are quite easy to
485  // automatically generate using a 1D Gauss rule on
486  // [0,1] and two 1D Jacobi-Gauss rules on [0,1].
487  QConical conical_rule(3, _order);
488  conical_rule.init(_type, _p_level);
489 
490  // Swap points and weights with the about-to-be destroyed rule.
491  _points.swap (conical_rule.get_points() );
492  _weights.swap(conical_rule.get_weights());
493 
494  return;
495  }
496  }
497  }
498  } // end case TET4,TET10
499 
500 
501 
502  //---------------------------------------------
503  // Prism quadrature rules
504  case PRISM6:
505  case PRISM15:
506  case PRISM18:
507  {
508  // We compute the 3D quadrature rule as a tensor
509  // product of the 1D quadrature rule and a 2D
510  // triangle quadrature rule
511 
512  QGauss q1D(1,_order);
513  QGauss q2D(2,_order);
514 
515  // Initialize
516  q1D.init(EDGE2, _p_level);
517  q2D.init(TRI3, _p_level);
518 
519  tensor_product_prism(q1D, q2D);
520 
521  return;
522  }
523 
524 
525 
526  //---------------------------------------------
527  // Pyramid
528  case PYRAMID5:
529  case PYRAMID13:
530  case PYRAMID14:
531  {
532  // We compute the Pyramid rule as a conical product of a
533  // Jacobi rule with alpha==2 on the interval [0,1] two 1D
534  // Gauss rules with suitably modified points. The idea comes
535  // from: Stroud, A.H. "Approximate Calculation of Multiple
536  // Integrals."
537  QConical conical_rule(3, _order);
538  conical_rule.init(_type, _p_level);
539 
540  // Swap points and weights with the about-to-be destroyed rule.
541  _points.swap (conical_rule.get_points() );
542  _weights.swap(conical_rule.get_weights());
543 
544  return;
545 
546  }
547 
548 
549 
550  //---------------------------------------------
551  // Unsupported type
552  default:
553  libmesh_error_msg("ERROR: Unsupported type: " << Utility::enum_to_string(_type));
554  }
555 #endif
556 }

References libMesh::QBase::_order, libMesh::QBase::_p_level, libMesh::QBase::_points, libMesh::QBase::_type, libMesh::QBase::_weights, libMesh::QBase::allow_rules_with_negative_weights, libMesh::CONSTANT, libMesh::EDGE2, libMesh::EIGHTH, libMesh::Utility::enum_to_string(), libMesh::FIFTH, libMesh::FIRST, libMesh::FOURTH, libMesh::QBase::get_order(), libMesh::QBase::get_points(), libMesh::QBase::get_weights(), libMesh::HEX20, libMesh::HEX27, libMesh::HEX8, libMesh::QBase::init(), keast_rule(), libMesh::PRISM15, libMesh::PRISM18, libMesh::PRISM6, libMesh::PYRAMID13, libMesh::PYRAMID14, libMesh::PYRAMID5, libMesh::Real, libMesh::SECOND, libMesh::SEVENTH, libMesh::SIXTH, libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), libMesh::TET10, libMesh::TET4, libMesh::THIRD, and libMesh::TRI3.

◆ keast_rule()

void libMesh::QGauss::keast_rule ( const Real  rule_data[][4],
const unsigned int  n_pts 
)
private

The Keast rules are for tets.

This function takes permutation points and weights in a specific format as input and fills the _points and _weights vectors.

Definition at line 38 of file quadrature_gauss.C.

40 {
41  // Like the Dunavant rule, the input data should have 4 columns. These columns
42  // have the following format and implied permutations (w=weight).
43  // {a, 0, 0, w} = 1-permutation (a,a,a)
44  // {a, b, 0, w} = 4-permutation (a,b,b), (b,a,b), (b,b,a), (b,b,b)
45  // {a, 0, b, w} = 6-permutation (a,a,b), (a,b,b), (b,b,a), (b,a,b), (b,a,a), (a,b,a)
46  // {a, b, c, w} = 12-permutation (a,a,b), (a,a,c), (b,a,a), (c,a,a), (a,b,a), (a,c,a)
47  // (a,b,c), (a,c,b), (b,a,c), (b,c,a), (c,a,b), (c,b,a)
48 
49  // Always insert into the points & weights vector relative to the offset
50  unsigned int offset=0;
51 
52 
53  for (unsigned int p=0; p<n_pts; ++p)
54  {
55 
56  // There must always be a non-zero entry to start the row
57  libmesh_assert_not_equal_to (rule_data[p][0], static_cast<Real>(0.0));
58 
59  // A zero weight may imply you did not set up the raw data correctly
60  libmesh_assert_not_equal_to (rule_data[p][3], static_cast<Real>(0.0));
61 
62  // What kind of point is this?
63  // One non-zero entry in first 3 cols ? 1-perm (centroid) point = 1
64  // Two non-zero entries in first 3 cols ? 3-perm point = 3
65  // Three non-zero entries ? 6-perm point = 6
66  unsigned int pointtype=1;
67 
68  if (rule_data[p][1] != static_cast<Real>(0.0))
69  {
70  if (rule_data[p][2] != static_cast<Real>(0.0))
71  pointtype = 12;
72  else
73  pointtype = 4;
74  }
75  else
76  {
77  // The second entry is zero. What about the third?
78  if (rule_data[p][2] != static_cast<Real>(0.0))
79  pointtype = 6;
80  }
81 
82 
83  switch (pointtype)
84  {
85  case 1:
86  {
87  // Be sure we have enough space to insert this point
88  libmesh_assert_less (offset + 0, _points.size());
89 
90  const Real a = rule_data[p][0];
91 
92  // The point has only a single permutation (the centroid!)
93  _points[offset + 0] = Point(a,a,a);
94 
95  // The weight is always the last entry in the row.
96  _weights[offset + 0] = rule_data[p][3];
97 
98  offset += pointtype;
99  break;
100  }
101 
102  case 4:
103  {
104  // Be sure we have enough space to insert these points
105  libmesh_assert_less (offset + 3, _points.size());
106 
107  const Real a = rule_data[p][0];
108  const Real b = rule_data[p][1];
109  const Real wt = rule_data[p][3];
110 
111  // Here it's understood the second entry is to be used twice, and
112  // thus there are three possible permutations.
113  _points[offset + 0] = Point(a,b,b);
114  _points[offset + 1] = Point(b,a,b);
115  _points[offset + 2] = Point(b,b,a);
116  _points[offset + 3] = Point(b,b,b);
117 
118  for (unsigned int j=0; j<pointtype; ++j)
119  _weights[offset + j] = wt;
120 
121  offset += pointtype;
122  break;
123  }
124 
125  case 6:
126  {
127  // Be sure we have enough space to insert these points
128  libmesh_assert_less (offset + 5, _points.size());
129 
130  const Real a = rule_data[p][0];
131  const Real b = rule_data[p][2];
132  const Real wt = rule_data[p][3];
133 
134  // Three individual entries with six permutations.
135  _points[offset + 0] = Point(a,a,b);
136  _points[offset + 1] = Point(a,b,b);
137  _points[offset + 2] = Point(b,b,a);
138  _points[offset + 3] = Point(b,a,b);
139  _points[offset + 4] = Point(b,a,a);
140  _points[offset + 5] = Point(a,b,a);
141 
142  for (unsigned int j=0; j<pointtype; ++j)
143  _weights[offset + j] = wt;
144 
145  offset += pointtype;
146  break;
147  }
148 
149 
150  case 12:
151  {
152  // Be sure we have enough space to insert these points
153  libmesh_assert_less (offset + 11, _points.size());
154 
155  const Real a = rule_data[p][0];
156  const Real b = rule_data[p][1];
157  const Real c = rule_data[p][2];
158  const Real wt = rule_data[p][3];
159 
160  // Three individual entries with six permutations.
161  _points[offset + 0] = Point(a,a,b); _points[offset + 6] = Point(a,b,c);
162  _points[offset + 1] = Point(a,a,c); _points[offset + 7] = Point(a,c,b);
163  _points[offset + 2] = Point(b,a,a); _points[offset + 8] = Point(b,a,c);
164  _points[offset + 3] = Point(c,a,a); _points[offset + 9] = Point(b,c,a);
165  _points[offset + 4] = Point(a,b,a); _points[offset + 10] = Point(c,a,b);
166  _points[offset + 5] = Point(a,c,a); _points[offset + 11] = Point(c,b,a);
167 
168  for (unsigned int j=0; j<pointtype; ++j)
169  _weights[offset + j] = wt;
170 
171  offset += pointtype;
172  break;
173  }
174 
175  default:
176  libmesh_error_msg("Don't know what to do with this many permutation points!");
177  }
178 
179  }
180 
181 }

References libMesh::QBase::_points, libMesh::QBase::_weights, and libMesh::Real.

Referenced by init_3D().

◆ n_objects()

static unsigned int libMesh::ReferenceCounter::n_objects ( )
inlinestaticinherited

Prints the number of outstanding (created, but not yet destroyed) objects.

Definition at line 83 of file reference_counter.h.

84  { return _n_objects; }

References libMesh::ReferenceCounter::_n_objects.

◆ n_points()

unsigned int libMesh::QBase::n_points ( ) const
inlineinherited
Returns
The number of points associated with the quadrature rule.

Definition at line 126 of file quadrature.h.

127  {
128  libmesh_assert (!_points.empty());
129  return cast_int<unsigned int>(_points.size());
130  }

References libMesh::QBase::_points, and libMesh::libmesh_assert().

Referenced by libMesh::ExactSolution::_compute_error(), assemble(), assemble_1D(), assemble_ellipticdg(), assemble_helmholtz(), assemble_poisson(), assemble_shell(), assemble_wave(), assembly_with_dg_fem_context(), AssemblyA0::boundary_assembly(), AssemblyF0::boundary_assembly(), AssemblyA1::boundary_assembly(), AssemblyF1::boundary_assembly(), AssemblyF2::boundary_assembly(), A2::boundary_assembly(), AssemblyA2::boundary_assembly(), A3::boundary_assembly(), F0::boundary_assembly(), Output0::boundary_assembly(), libMesh::System::calculate_norm(), libMesh::FirstOrderUnsteadySolver::compute_second_order_eqns(), libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), SecondOrderScalarSystemSecondOrderTimeSolverBase::damping_residual(), SecondOrderScalarSystemFirstOrderTimeSolverBase::damping_residual(), NavierSystem::element_constraint(), CoupledSystem::element_constraint(), PoissonSystem::element_postprocess(), LaplaceSystem::element_postprocess(), LaplaceQoI::element_qoi(), LaplaceQoI::element_qoi_derivative(), LaplaceSystem::element_qoi_derivative(), HeatSystem::element_qoi_derivative(), NavierSystem::element_time_derivative(), SolidSystem::element_time_derivative(), PoissonSystem::element_time_derivative(), LaplaceSystem::element_time_derivative(), CurlCurlSystem::element_time_derivative(), ElasticitySystem::element_time_derivative(), L2System::element_time_derivative(), CoupledSystem::element_time_derivative(), FirstOrderScalarSystemBase::element_time_derivative(), SecondOrderScalarSystemFirstOrderTimeSolverBase::element_time_derivative(), libMesh::RBEIMConstruction::enrich_RB_space(), B::interior_assembly(), A0::interior_assembly(), M0::interior_assembly(), A1::interior_assembly(), AssemblyA0::interior_assembly(), EIM_IP_assembly::interior_assembly(), AcousticsInnerProduct::interior_assembly(), AssemblyA1::interior_assembly(), A2::interior_assembly(), AssemblyA2::interior_assembly(), EIM_F::interior_assembly(), F0::interior_assembly(), OutputAssembly::interior_assembly(), InnerProductAssembly::interior_assembly(), AssemblyEIM::interior_assembly(), AssemblyF0::interior_assembly(), AssemblyF1::interior_assembly(), Ex6InnerProduct::interior_assembly(), Ex6EIMInnerProduct::interior_assembly(), LaplaceYoung::jacobian(), NavierSystem::mass_residual(), ElasticitySystem::mass_residual(), libMesh::FEMPhysics::mass_residual(), FirstOrderScalarSystemBase::mass_residual(), SecondOrderScalarSystemSecondOrderTimeSolverBase::mass_residual(), SecondOrderScalarSystemFirstOrderTimeSolverBase::mass_residual(), libMesh::QBase::print_info(), LaplaceYoung::residual(), LaplaceSystem::side_constraint(), LaplaceSystem::side_postprocess(), CoupledSystemQoI::side_qoi(), CoupledSystemQoI::side_qoi_derivative(), LaplaceSystem::side_qoi_derivative(), SolidSystem::side_time_derivative(), CurlCurlSystem::side_time_derivative(), ElasticitySystem::side_time_derivative(), libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), libMesh::QBase::tensor_product_quad(), and libMesh::RBEIMConstruction::truth_solve().

◆ operator=() [1/2]

QGauss& libMesh::QGauss::operator= ( const QGauss &  )
default

◆ operator=() [2/2]

QGauss& libMesh::QGauss::operator= ( QGauss &&  )
default

◆ print_info() [1/2]

void libMesh::QBase::print_info ( std::ostream &  os = libMesh::out) const
inherited

Prints information relevant to the quadrature rule, by default to libMesh::out.

Definition at line 37 of file quadrature.C.

38 {
39  libmesh_assert(!_points.empty());
40  libmesh_assert(!_weights.empty());
41 
42  Real summed_weights=0;
43  os << "N_Q_Points=" << this->n_points() << std::endl << std::endl;
44  for (auto qpoint: index_range(_points))
45  {
46  os << " Point " << qpoint << ":\n"
47  << " "
48  << _points[qpoint]
49  << "\n Weight:\n "
50  << " w=" << _weights[qpoint] << "\n" << std::endl;
51 
52  summed_weights += _weights[qpoint];
53  }
54  os << "Summed Weights: " << summed_weights << std::endl;
55 }

References libMesh::QBase::_points, libMesh::QBase::_weights, libMesh::index_range(), libMesh::libmesh_assert(), libMesh::QBase::n_points(), and libMesh::Real.

Referenced by libMesh::operator<<().

◆ print_info() [2/2]

void libMesh::ReferenceCounter::print_info ( std::ostream &  out = libMesh::out)
staticinherited

Prints the reference information, by default to libMesh::out.

Definition at line 87 of file reference_counter.C.

88 {
89  if (_enable_print_counter)
90  out_stream << ReferenceCounter::get_info();
91 }

References libMesh::ReferenceCounter::_enable_print_counter, and libMesh::ReferenceCounter::get_info().

◆ qp()

Point libMesh::QBase::qp ( const unsigned int  i) const
inlineinherited
Returns
The \( i^{th} \) quadrature point in reference element space.

Definition at line 164 of file quadrature.h.

165  {
166  libmesh_assert_less (i, _points.size());
167  return _points[i];
168  }

References libMesh::QBase::_points.

Referenced by libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), and libMesh::QBase::tensor_product_quad().

◆ scale()

void libMesh::QBase::scale ( std::pair< Real, Real >  old_range,
std::pair< Real, Real >  new_range 
)
inherited

Maps the points of a 1D quadrature rule defined by "old_range" to another 1D interval defined by "new_range" and scales the weights accordingly.

Definition at line 137 of file quadrature.C.

139 {
140  // Make sure we are in 1D
141  libmesh_assert_equal_to (_dim, 1);
142 
143  Real
144  h_new = new_range.second - new_range.first,
145  h_old = old_range.second - old_range.first;
146 
147  // Make sure that we have sane ranges
148  libmesh_assert_greater (h_new, 0.);
149  libmesh_assert_greater (h_old, 0.);
150 
151  // Make sure there are some points
152  libmesh_assert_greater (_points.size(), 0);
153 
154  // Compute the scale factor
155  Real scfact = h_new/h_old;
156 
157  // We're mapping from old_range -> new_range
158  for (auto i : index_range(_points))
159  {
160  _points[i](0) = new_range.first +
161  (_points[i](0) - old_range.first) * scfact;
162 
163  // Scale the weights
164  _weights[i] *= scfact;
165  }
166 }

References libMesh::QBase::_dim, libMesh::QBase::_points, libMesh::QBase::_weights, libMesh::index_range(), and libMesh::Real.

Referenced by libMesh::QConical::conical_product_tet(), and libMesh::QConical::conical_product_tri().

◆ shapes_need_reinit()

virtual bool libMesh::QBase::shapes_need_reinit ( )
inlinevirtualinherited
Returns
true if the shape functions need to be recalculated, false otherwise.

This may be required if the number of quadrature points or their position changes.

Definition at line 239 of file quadrature.h.

239 { return false; }

◆ tensor_product_hex()

void libMesh::QBase::tensor_product_hex ( const QBase &  q1D)
protectedinherited

Computes the tensor product quadrature rule [q1D x q1D x q1D] from the 1D rule q1D.

Used in the init_3D routines for hexahedral element types.

Definition at line 198 of file quadrature.C.

199 {
200  const unsigned int np = q1D.n_points();
201 
202  _points.resize(np * np * np);
203 
204  _weights.resize(np * np * np);
205 
206  unsigned int q=0;
207 
208  for (unsigned int k=0; k<np; k++)
209  for (unsigned int j=0; j<np; j++)
210  for (unsigned int i=0; i<np; i++)
211  {
212  _points[q](0) = q1D.qp(i)(0);
213  _points[q](1) = q1D.qp(j)(0);
214  _points[q](2) = q1D.qp(k)(0);
215 
216  _weights[q] = q1D.w(i) * q1D.w(j) * q1D.w(k);
217 
218  q++;
219  }
220 }

References libMesh::QBase::_points, libMesh::QBase::_weights, libMesh::QBase::n_points(), libMesh::QBase::qp(), and libMesh::QBase::w().

Referenced by libMesh::QGaussLobatto::init_3D(), libMesh::QTrap::init_3D(), libMesh::QGrid::init_3D(), libMesh::QSimpson::init_3D(), and init_3D().

◆ tensor_product_prism()

void libMesh::QBase::tensor_product_prism ( const QBase &  q1D,
const QBase &  q2D 
)
protectedinherited

Computes the tensor product of a 1D quadrature rule and a 2D quadrature rule.

Used in the init_3D routines for prismatic element types.

Definition at line 225 of file quadrature.C.

226 {
227  const unsigned int n_points1D = q1D.n_points();
228  const unsigned int n_points2D = q2D.n_points();
229 
230  _points.resize (n_points1D * n_points2D);
231  _weights.resize (n_points1D * n_points2D);
232 
233  unsigned int q=0;
234 
235  for (unsigned int j=0; j<n_points1D; j++)
236  for (unsigned int i=0; i<n_points2D; i++)
237  {
238  _points[q](0) = q2D.qp(i)(0);
239  _points[q](1) = q2D.qp(i)(1);
240  _points[q](2) = q1D.qp(j)(0);
241 
242  _weights[q] = q2D.w(i) * q1D.w(j);
243 
244  q++;
245  }
246 
247 }

References libMesh::QBase::_points, libMesh::QBase::_weights, libMesh::QBase::n_points(), libMesh::QBase::qp(), and libMesh::QBase::w().

Referenced by libMesh::QGrid::init_3D(), libMesh::QTrap::init_3D(), libMesh::QSimpson::init_3D(), and init_3D().

◆ tensor_product_quad()

void libMesh::QBase::tensor_product_quad ( const QBase &  q1D)
protectedinherited

Constructs a 2D rule from the tensor product of q1D with itself.

Used in the init_2D() routines for quadrilateral element types.

Definition at line 171 of file quadrature.C.

172 {
173 
174  const unsigned int np = q1D.n_points();
175 
176  _points.resize(np * np);
177 
178  _weights.resize(np * np);
179 
180  unsigned int q=0;
181 
182  for (unsigned int j=0; j<np; j++)
183  for (unsigned int i=0; i<np; i++)
184  {
185  _points[q](0) = q1D.qp(i)(0);
186  _points[q](1) = q1D.qp(j)(0);
187 
188  _weights[q] = q1D.w(i)*q1D.w(j);
189 
190  q++;
191  }
192 }

References libMesh::QBase::_points, libMesh::QBase::_weights, libMesh::QBase::n_points(), libMesh::QBase::qp(), and libMesh::QBase::w().

Referenced by libMesh::QGaussLobatto::init_2D(), libMesh::QTrap::init_2D(), libMesh::QGrid::init_2D(), libMesh::QSimpson::init_2D(), and init_2D().

◆ type()

QuadratureType libMesh::QGauss::type ( ) const
overridevirtual
Returns
QGAUSS.

Implements libMesh::QBase.

Definition at line 33 of file quadrature_gauss.C.

34 {
35  return QGAUSS;
36 }

References libMesh::QGAUSS.

◆ w()

Real libMesh::QBase::w ( const unsigned int  i) const
inlineinherited
Returns
The \( i^{th} \) quadrature weight.

Definition at line 173 of file quadrature.h.

174  {
175  libmesh_assert_less (i, _weights.size());
176  return _weights[i];
177  }

References libMesh::QBase::_weights.

Referenced by libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), and libMesh::QBase::tensor_product_quad().

Member Data Documentation

◆ _counts

ReferenceCounter::Counts libMesh::ReferenceCounter::_counts
staticprotectedinherited

Actually holds the data.

Definition at line 122 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::get_info(), libMesh::ReferenceCounter::increment_constructor_count(), and libMesh::ReferenceCounter::increment_destructor_count().

◆ _dim

unsigned int libMesh::QBase::_dim
protectedinherited

The spatial dimension of the quadrature rule.

Definition at line 339 of file quadrature.h.

Referenced by libMesh::QBase::build(), libMesh::QBase::get_dim(), libMesh::QBase::init(), libMesh::QGaussLobatto::QGaussLobatto(), libMesh::QJacobi::QJacobi(), libMesh::QNodal::QNodal(), libMesh::QSimpson::QSimpson(), libMesh::QTrap::QTrap(), and libMesh::QBase::scale().

◆ _enable_print_counter

bool libMesh::ReferenceCounter::_enable_print_counter = true
staticprotectedinherited

Flag to control whether reference count information is printed when print_info is called.

Definition at line 141 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::disable_print_counter_info(), libMesh::ReferenceCounter::enable_print_counter_info(), and libMesh::ReferenceCounter::print_info().

◆ _mutex

Threads::spin_mutex libMesh::ReferenceCounter::_mutex
staticprotectedinherited

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 135 of file reference_counter.h.

◆ _n_objects

Threads::atomic< unsigned int > libMesh::ReferenceCounter::_n_objects
staticprotectedinherited

The number of objects.

Print the reference count information when the number returns to 0.

Definition at line 130 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::n_objects(), libMesh::ReferenceCounter::ReferenceCounter(), and libMesh::ReferenceCounter::~ReferenceCounter().

◆ _order

Order libMesh::QBase::_order
protectedinherited

The polynomial order which the quadrature rule is capable of integrating exactly.

Definition at line 345 of file quadrature.h.

Referenced by libMesh::QBase::build(), libMesh::QBase::get_order(), libMesh::QClough::init_1D(), libMesh::QGaussLobatto::init_1D(), libMesh::QGrid::init_1D(), init_1D(), libMesh::QNodal::init_1D(), libMesh::QJacobi::init_1D(), libMesh::QMonomial::init_1D(), libMesh::QClough::init_2D(), libMesh::QGaussLobatto::init_2D(), libMesh::QGrid::init_2D(), init_2D(), libMesh::QNodal::init_2D(), libMesh::QMonomial::init_2D(), libMesh::QGaussLobatto::init_3D(), libMesh::QGrid::init_3D(), init_3D(), libMesh::QNodal::init_3D(), and libMesh::QMonomial::init_3D().

◆ _p_level

unsigned int libMesh::QBase::_p_level
protectedinherited

The p-level of the element for which the current values have been computed.

Definition at line 357 of file quadrature.h.

Referenced by libMesh::QBase::get_order(), libMesh::QBase::get_p_level(), libMesh::QBase::init(), libMesh::QClough::init_1D(), libMesh::QNodal::init_1D(), libMesh::QMonomial::init_1D(), libMesh::QClough::init_2D(), libMesh::QGaussLobatto::init_2D(), init_2D(), libMesh::QNodal::init_2D(), libMesh::QMonomial::init_2D(), libMesh::QGaussLobatto::init_3D(), init_3D(), libMesh::QNodal::init_3D(), and libMesh::QMonomial::init_3D().

◆ _points

std::vector<Point> libMesh::QBase::_points
protectedinherited

The locations of the quadrature points in reference element space.

Definition at line 363 of file quadrature.h.

Referenced by libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), dunavant_rule(), dunavant_rule2(), libMesh::QBase::get_points(), libMesh::QGrundmann_Moller::gm_rule(), libMesh::QBase::init_0D(), libMesh::QClough::init_1D(), libMesh::QGaussLobatto::init_1D(), libMesh::QTrap::init_1D(), libMesh::QConical::init_1D(), libMesh::QGrid::init_1D(), init_1D(), libMesh::QSimpson::init_1D(), libMesh::QNodal::init_1D(), libMesh::QJacobi::init_1D(), libMesh::QMonomial::init_1D(), libMesh::QGrundmann_Moller::init_1D(), libMesh::QClough::init_2D(), libMesh::QTrap::init_2D(), libMesh::QGrid::init_2D(), libMesh::QSimpson::init_2D(), init_2D(), libMesh::QNodal::init_2D(), libMesh::QMonomial::init_2D(), libMesh::QGrid::init_3D(), libMesh::QTrap::init_3D(), libMesh::QSimpson::init_3D(), init_3D(), libMesh::QNodal::init_3D(), libMesh::QMonomial::init_3D(), libMesh::QGrundmann_Moller::init_3D(), keast_rule(), libMesh::QMonomial::kim_rule(), libMesh::QBase::n_points(), libMesh::QBase::print_info(), libMesh::QBase::qp(), libMesh::QBase::scale(), libMesh::QMonomial::stroud_rule(), libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), libMesh::QBase::tensor_product_quad(), and libMesh::QMonomial::wissmann_rule().

◆ _type

ElemType libMesh::QBase::_type
protectedinherited

The type of element for which the current values have been computed.

Definition at line 351 of file quadrature.h.

Referenced by libMesh::QBase::get_elem_type(), libMesh::QBase::init(), libMesh::QClough::init_1D(), libMesh::QNodal::init_1D(), libMesh::QMonomial::init_1D(), libMesh::QClough::init_2D(), libMesh::QGaussLobatto::init_2D(), libMesh::QConical::init_2D(), libMesh::QTrap::init_2D(), libMesh::QGrid::init_2D(), init_2D(), libMesh::QSimpson::init_2D(), libMesh::QNodal::init_2D(), libMesh::QMonomial::init_2D(), libMesh::QGrundmann_Moller::init_2D(), libMesh::QClough::init_3D(), libMesh::QGaussLobatto::init_3D(), libMesh::QTrap::init_3D(), libMesh::QGrid::init_3D(), libMesh::QConical::init_3D(), libMesh::QSimpson::init_3D(), init_3D(), libMesh::QNodal::init_3D(), libMesh::QMonomial::init_3D(), and libMesh::QGrundmann_Moller::init_3D().

◆ _weights

std::vector<Real> libMesh::QBase::_weights
protectedinherited

The quadrature weights.

The order of the weights matches the ordering of the _points vector.

Definition at line 369 of file quadrature.h.

Referenced by libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), dunavant_rule(), dunavant_rule2(), libMesh::QBase::get_weights(), libMesh::QGrundmann_Moller::gm_rule(), libMesh::QBase::init_0D(), libMesh::QClough::init_1D(), libMesh::QGaussLobatto::init_1D(), libMesh::QConical::init_1D(), libMesh::QTrap::init_1D(), libMesh::QGrid::init_1D(), init_1D(), libMesh::QSimpson::init_1D(), libMesh::QNodal::init_1D(), libMesh::QJacobi::init_1D(), libMesh::QMonomial::init_1D(), libMesh::QGrundmann_Moller::init_1D(), libMesh::QClough::init_2D(), libMesh::QTrap::init_2D(), libMesh::QGrid::init_2D(), init_2D(), libMesh::QSimpson::init_2D(), libMesh::QNodal::init_2D(), libMesh::QMonomial::init_2D(), libMesh::QGrid::init_3D(), libMesh::QTrap::init_3D(), libMesh::QSimpson::init_3D(), init_3D(), libMesh::QNodal::init_3D(), libMesh::QMonomial::init_3D(), libMesh::QGrundmann_Moller::init_3D(), keast_rule(), libMesh::QMonomial::kim_rule(), libMesh::QBase::print_info(), libMesh::QBase::scale(), libMesh::QMonomial::stroud_rule(), libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), libMesh::QBase::tensor_product_quad(), libMesh::QBase::w(), and libMesh::QMonomial::wissmann_rule().

◆ allow_rules_with_negative_weights

bool libMesh::QBase::allow_rules_with_negative_weights
inherited

Flag (default true) controlling the use of quadrature rules with negative weights.

Set this to false to require rules with all positive weights.

Rules with negative weights can be unsuitable for some problems. For example, it is possible for a rule with negative weights to obtain a negative result when integrating a positive function.

A particular example: if rules with negative weights are not allowed, a request for TET,THIRD (5 points) will return the TET,FIFTH (14 points) rule instead, nearly tripling the computational effort required!

Definition at line 254 of file quadrature.h.

Referenced by libMesh::QGrundmann_Moller::init_2D(), init_3D(), libMesh::QMonomial::init_3D(), and libMesh::QGrundmann_Moller::init_3D().

The documentation for this class was generated from the following files:
libMesh::HEX20
Definition: enum_elem_type.h:48
libMesh::THIRTYFOURTH
Definition: enum_order.h:75
libMesh::PRISM6
Definition: enum_elem_type.h:50
libMesh::QBase::_p_level
unsigned int _p_level
The p-level of the element for which the current values have been computed.
Definition: quadrature.h:357
libMesh::QMONOMIAL
Definition: enum_quadrature_type.h:41
libMesh::TWENTYFIFTH
Definition: enum_order.h:66
libMesh::THIRTYSECOND
Definition: enum_order.h:73
libMesh::THIRTYFIFTH
Definition: enum_order.h:76
libMesh::QCLOUGH
Definition: enum_quadrature_type.h:44
libMesh::QBase::tensor_product_hex
void tensor_product_hex(const QBase &q1D)
Computes the tensor product quadrature rule [q1D x q1D x q1D] from the 1D rule q1D.
Definition: quadrature.C:198
libMesh::HEX8
Definition: enum_elem_type.h:47
libMesh::NINETEENTH
Definition: enum_order.h:60
libMesh::SIXTH
Definition: enum_order.h:47
libMesh::QBase::n_points
unsigned int n_points() const
Definition: quadrature.h:126
libMesh::FOURTEENTH
Definition: enum_order.h:55
libMesh::FORTYSECOND
Definition: enum_order.h:83
libMesh::index_range
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:106
libMesh::TET10
Definition: enum_elem_type.h:46
libMesh::ReferenceCounter::_counts
static Counts _counts
Actually holds the data.
Definition: reference_counter.h:122
std::sqrt
MetaPhysicL::DualNumber< T, D > sqrt(const MetaPhysicL::DualNumber< T, D > &in)
libMesh::QGauss::QGauss
QGauss(unsigned int dim, Order order=INVALID_ORDER)
Constructor.
Definition: quadrature_gauss.h:52
libMesh::QBase::init_1D
virtual void init_1D(const ElemType type=INVALID_ELEM, unsigned int p_level=0)=0
Initializes the 1D quadrature rule by filling the points and weights vectors with the appropriate val...
libMesh::ReferenceCounter::_n_objects
static Threads::atomic< unsigned int > _n_objects
The number of objects.
Definition: reference_counter.h:130
libMesh::FIFTH
Definition: enum_order.h:46
libMesh::THIRTEENTH
Definition: enum_order.h:54
libMesh::THIRTYSEVENTH
Definition: enum_order.h:78
libMesh::SECOND
Definition: enum_order.h:43
libMesh::QGAUSS
Definition: enum_quadrature_type.h:34
libMesh::THIRTIETH
Definition: enum_order.h:71
libMesh::ReferenceCounter::get_info
static std::string get_info()
Gets a string containing the reference information.
Definition: reference_counter.C:47
libMesh::QGauss::dunavant_rule
void dunavant_rule(const Real rule_data[][4], const unsigned int n_pts)
The Dunavant rules are for triangles.
Definition: quadrature_gauss.C:265
libMesh::TWENTIETH
Definition: enum_order.h:61
libMesh::QBase::init_3D
virtual void init_3D(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:130
libMesh::QBase::init
virtual void init(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Initializes the data structures for a quadrature rule for an element of type type.
Definition: quadrature.C:59
libMesh::TWELFTH
Definition: enum_order.h:53
dim
unsigned int dim
Definition: adaptivity_ex3.C:113
libMesh::TWENTYEIGHTH
Definition: enum_order.h:69
libMesh::QBase::build
static std::unique_ptr< 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.
Definition: quadrature_build.C:41
libMesh::FORTYTHIRD
Definition: enum_order.h:84
libMesh::TET4
Definition: enum_elem_type.h:45
libMesh::QSIMPSON
Definition: enum_quadrature_type.h:37
libMesh::PRISM15
Definition: enum_elem_type.h:51
libMesh::THIRTYEIGHTH
Definition: enum_order.h:79
libMesh::QBase::_order
Order _order
The polynomial order which the quadrature rule is capable of integrating exactly.
Definition: quadrature.h:345
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::ELEVENTH
Definition: enum_order.h:52
libMesh::EIGHTTEENTH
Definition: enum_order.h:59
libMesh::HEX27
Definition: enum_elem_type.h:49
libMesh::TWENTYSECOND
Definition: enum_order.h:63
libMesh::TWENTYTHIRD
Definition: enum_order.h:64
libMesh::QGauss::dunavant_rule2
void dunavant_rule2(const Real *wts, const Real *a, const Real *b, const unsigned int *permutation_ids, const unsigned int n_wts)
Definition: quadrature_gauss.C:190
libMesh::QTRAP
Definition: enum_quadrature_type.h:38
libMesh::QNODAL
Definition: enum_quadrature_type.h:46
libMesh::TWENTYFIRST
Definition: enum_order.h:62
libMesh::FORTYFIRST
Definition: enum_order.h:82
libMesh::Threads::spin_mtx
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:29
libMesh::TWENTYSIXTH
Definition: enum_order.h:67
libMesh::QUAD4
Definition: enum_elem_type.h:41
libMesh::THIRTYFIRST
Definition: enum_order.h:72
libMesh::QBase::QBase
QBase(unsigned int dim, Order order=INVALID_ORDER)
Constructor.
Definition: quadrature.C:27
libMesh::QBase::_type
ElemType _type
The type of element for which the current values have been computed.
Definition: quadrature.h:351
libMesh::FORTIETH
Definition: enum_order.h:81
libMesh::TRI3
Definition: enum_elem_type.h:39
libMesh::CONSTANT
Definition: enum_order.h:41
libMesh::TWENTYSEVENTH
Definition: enum_order.h:68
libMesh::EIGHTH
Definition: enum_order.h:49
libMesh::TENTH
Definition: enum_order.h:51
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::QBase::tensor_product_prism
void tensor_product_prism(const QBase &q1D, const QBase &q2D)
Computes the tensor product of a 1D quadrature rule and a 2D quadrature rule.
Definition: quadrature.C:225
libMesh::FOURTH
Definition: enum_order.h:45
libMesh::TWENTYNINTH
Definition: enum_order.h:70
libMesh::QBase::_weights
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:369
libMesh::TRI6
Definition: enum_elem_type.h:40
libMesh::QBase::tensor_product_quad
void tensor_product_quad(const QBase &q1D)
Constructs a 2D rule from the tensor product of q1D with itself.
Definition: quadrature.C:171
libMesh::QGAUSS_LOBATTO
Definition: enum_quadrature_type.h:43
libMesh::QBase::type
virtual QuadratureType type() const =0
libMesh::QJACOBI_1_0
Definition: enum_quadrature_type.h:35
libMesh::SEVENTH
Definition: enum_order.h:48
libMesh::QUADSHELL8
Definition: enum_elem_type.h:73
libMesh::PYRAMID5
Definition: enum_elem_type.h:53
libMesh::QUADSHELL4
Definition: enum_elem_type.h:72
libMesh::FIFTEENTH
Definition: enum_order.h:56
libMesh::QBase::allow_rules_with_negative_weights
bool allow_rules_with_negative_weights
Flag (default true) controlling the use of quadrature rules with negative weights.
Definition: quadrature.h:254
libMesh::THIRTYNINTH
Definition: enum_order.h:80
libMesh::QJACOBI_2_0
Definition: enum_quadrature_type.h:36
libMesh::THIRTYTHIRD
Definition: enum_order.h:74
libMesh::SIXTEENTH
Definition: enum_order.h:57
libMesh::THIRD
Definition: enum_order.h:44
libMesh::QUAD9
Definition: enum_elem_type.h:43
libMesh::NINTH
Definition: enum_order.h:50
libMesh::TRI3SUBDIVISION
Definition: enum_elem_type.h:69
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::QBase::get_order
Order get_order() const
Definition: quadrature.h:211
libMesh::QCONICAL
Definition: enum_quadrature_type.h:42
libMesh::ReferenceCounter::_enable_print_counter
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called.
Definition: reference_counter.h:141
libMesh::TWENTYFOURTH
Definition: enum_order.h:65
libMesh::PRISM18
Definition: enum_elem_type.h:52
libMesh::TRISHELL3
Definition: enum_elem_type.h:71
libMesh::QBase::init_0D
virtual void init_0D(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Initializes the 0D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:113
libMesh::QBase::_dim
unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:339
libMesh::out
OStreamProxy out
libMesh::QBase::_points
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:363
libMesh::FIRST
Definition: enum_order.h:42
libMesh::THIRTYSIXTH
Definition: enum_order.h:77
libMesh::PYRAMID13
Definition: enum_elem_type.h:54
libMesh::Quality::name
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
libMesh::PYRAMID14
Definition: enum_elem_type.h:55
libMesh::EDGE2
Definition: enum_elem_type.h:35
libMesh::QUAD8
Definition: enum_elem_type.h:42
libMesh::QGauss::keast_rule
void keast_rule(const Real rule_data[][4], const unsigned int n_pts)
The Keast rules are for tets.
Definition: quadrature_gauss.C:38
libMesh::QGRUNDMANN_MOLLER
Definition: enum_quadrature_type.h:40
libMesh::QGRID
Definition: enum_quadrature_type.h:39
libMesh::QBase::init_2D
virtual void init_2D(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:123
libMesh::SEVENTEENTH
Definition: enum_order.h:58
Generated on Sat Jan 25 2020 12:08:40 for libMesh by   doxygen 1.8.16