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
 
QGaussoperator= (const QGauss &)=default
 
QGaussoperator= (QGauss &&)=default
 
virtual ~QGauss ()=default
 
virtual QuadratureType type () const override
 
virtual std::unique_ptr< QBaseclone () const override
 
ElemType get_elem_type () const
 
unsigned int get_p_level () const
 
unsigned int n_points () const
 
unsigned int size () const
 Alias for n_points() to enable use in index_range. More...
 
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 Elem &e, unsigned int p_level=invalid_uint)
 Initializes the data structures for a quadrature rule for the element e. More...
 
virtual void init (const ElemType type=INVALID_ELEM, unsigned int p_level=0, bool simple_type_only=false)
 Initializes the data structures for a quadrature rule for an element of type type. More...
 
virtual void init (const QBase &other_rule)
 Initializes the data structures for a quadrature rule based on the element, element type, and p_level settings of other_rule. 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
 
Order get_base_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< QBasebuild (std::string_view 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< QBasebuild (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_stream=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...
 
bool allow_nodal_pyramid_quadrature
 The flag's value defaults to false so that one does not accidentally use a nodal quadrature rule on Pyramid elements, since evaluating the inverse element Jacobian (e.g. 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 ()
 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) noexcept
 Increments the construction counter. More...
 
void increment_destructor_count (const std::string &name) noexcept
 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...
 
const Elem_elem
 The element for which the current values were computed, or nullptr if values were computed without a specific element. 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 () override
 Initializes the 1D quadrature rule by filling the points and weights vectors with the appropriate values. More...
 
virtual void init_2D () override
 Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate values. More...
 
virtual void init_3D () 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 119 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.

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

53  :
54  QBase(dim, order)
55  {
56  if (dim == 1)
57  init(EDGE2);
58  }
unsigned int dim
QBase(unsigned int dim, Order order=INVALID_ORDER)
Constructor.
Definition: quadrature.C:27
virtual void init(const Elem &e, unsigned int p_level=invalid_uint)
Initializes the data structures for a quadrature rule for the element e.
Definition: quadrature.C:65

◆ 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 ( std::string_view  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 43 of file quadrature_build.C.

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::QBase::clone(), main(), libMesh::OverlapCoupling::OverlapCoupling(), libMesh::InfFE< Dim, T_radial, T_map >::reinit(), QuadratureTest::testBuild(), QuadratureTest::testJacobi(), and QuadratureTest::testPolynomials().

46 {
47  return QBase::build (Utility::string_to_enum<QuadratureType> (type),
48  _dim,
49  _order);
50 }
unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:379
virtual QuadratureType type() const =0
Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:385
static std::unique_ptr< QBase > build(std::string_view name, const unsigned int dim, const Order order=INVALID_ORDER)
Builds a specific quadrature rule based on the name string.

◆ build() [2/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 54 of file quadrature_build.C.

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.

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

◆ clone()

std::unique_ptr< QBase > libMesh::QGauss::clone ( ) const
overridevirtual
Returns
A copy of this quadrature rule wrapped in a smart pointer.

Reimplemented from libMesh::QBase.

Definition at line 38 of file quadrature_gauss.C.

39 {
40  return std::make_unique<QGauss>(*this);
41 }

◆ disable_print_counter_info()

void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited

Definition at line 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

101 {
102  _enable_print_counter = false;
103  return;
104 }
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called...

◆ 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 270 of file quadrature_gauss.C.

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

Referenced by init_2D().

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

◆ 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 195 of file quadrature_gauss.C.

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

Referenced by init_2D().

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

◆ 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 94 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

95 {
96  _enable_print_counter = true;
97  return;
98 }
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called...

◆ get_base_order()

Order libMesh::QBase::get_base_order ( ) const
inlineinherited
Returns
The "base" order of the quadrature rule, independent of element.

This function should be used when comparing quadrature objects independently of their last initialization.

Definition at line 258 of file quadrature.h.

References libMesh::QBase::_order.

258 { return static_cast<Order>(_order); }
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:385

◆ get_dim()

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

Definition at line 150 of file quadrature.h.

References libMesh::QBase::_dim.

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

150 { return _dim; }
unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:379

◆ get_elem_type()

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

Definition at line 121 of file quadrature.h.

References libMesh::QBase::_type.

121 { return _type; }
ElemType _type
The type of element for which the current values have been computed.
Definition: quadrature.h:391

◆ 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.

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

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

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 & [name, cd] : _counts)
59  oss << "| " << name << " reference count information:\n"
60  << "| Creations: " << cd.first << '\n'
61  << "| Destructions: " << cd.second << '\n';
62 
63  oss << " ---------------------------------------------------------------------------- \n";
64 
65  return oss.str();
66 
67 #else
68 
69  return "";
70 
71 #endif
72 }
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
static Counts _counts
Actually holds the data.

◆ 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.

Todo:
This function should also be used in all of the Order switch statements in the rules themselves.

Definition at line 249 of file quadrature.h.

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

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

249 { return static_cast<Order>(_order + 2 * _p_level); }
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
unsigned int _p_level
The p-level of the element for which the current values have been computed.
Definition: quadrature.h:403
Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:385

◆ get_p_level()

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

Definition at line 126 of file quadrature.h.

References libMesh::QBase::_p_level.

126 { return _p_level; }
unsigned int _p_level
The p-level of the element for which the current values have been computed.
Definition: quadrature.h:403

◆ get_points() [1/2]

const std::vector<Point>& libMesh::QBase::get_points ( ) const
inlineinherited

◆ get_points() [2/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 162 of file quadrature.h.

References libMesh::QBase::_points.

162 { return _points; }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409

◆ get_weights() [1/2]

const std::vector<Real>& libMesh::QBase::get_weights ( ) const
inlineinherited

◆ get_weights() [2/2]

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

Definition at line 174 of file quadrature.h.

References libMesh::QBase::_weights.

174 { return _weights; }
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415

◆ increment_constructor_count()

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

Increments the construction counter.

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

Definition at line 183 of file reference_counter.h.

References libMesh::err, libMesh::BasicOStreamProxy< charT, traits >::get(), libMesh::Quality::name(), and libMesh::Threads::spin_mtx.

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

184 {
185  libmesh_try
186  {
187  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
188  std::pair<unsigned int, unsigned int> & p = _counts[name];
189  p.first++;
190  }
191  libmesh_catch (...)
192  {
193  auto stream = libMesh::err.get();
194  stream->exceptions(stream->goodbit); // stream must not throw
195  libMesh::err << "Encountered unrecoverable error while calling "
196  << "ReferenceCounter::increment_constructor_count() "
197  << "for a(n) " << name << " object." << std::endl;
198  std::terminate();
199  }
200 }
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
OStreamProxy err
static Counts _counts
Actually holds the data.
streamT * get()
Rather than implement every ostream/ios/ios_base function, we&#39;ll be lazy and make esoteric uses go th...
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:30

◆ increment_destructor_count()

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

Increments the destruction counter.

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

Definition at line 207 of file reference_counter.h.

References libMesh::err, libMesh::BasicOStreamProxy< charT, traits >::get(), libMesh::Quality::name(), and libMesh::Threads::spin_mtx.

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

208 {
209  libmesh_try
210  {
211  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
212  std::pair<unsigned int, unsigned int> & p = _counts[name];
213  p.second++;
214  }
215  libmesh_catch (...)
216  {
217  auto stream = libMesh::err.get();
218  stream->exceptions(stream->goodbit); // stream must not throw
219  libMesh::err << "Encountered unrecoverable error while calling "
220  << "ReferenceCounter::increment_destructor_count() "
221  << "for a(n) " << name << " object." << std::endl;
222  std::terminate();
223  }
224 }
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
OStreamProxy err
static Counts _counts
Actually holds the data.
streamT * get()
Rather than implement every ostream/ios/ios_base function, we&#39;ll be lazy and make esoteric uses go th...
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:30

◆ init() [1/4]

void libMesh::QBase::init ( const Elem e,
unsigned int  p_level = invalid_uint 
)
virtualinherited

Initializes the data structures for a quadrature rule for the element e.

If p_level is specified it overrides the element p_level() elevation to use.

Definition at line 65 of file quadrature.C.

References libMesh::QBase::_dim, libMesh::QBase::_elem, libMesh::QBase::_p_level, libMesh::QBase::_type, libMesh::Elem::dim(), libMesh::QBase::init_0D(), libMesh::QBase::init_1D(), libMesh::QBase::init_2D(), libMesh::QBase::init_3D(), libMesh::invalid_uint, libMesh::Elem::p_level(), libMesh::Elem::runtime_topology(), and libMesh::Elem::type().

Referenced by libMesh::QBase::init(), libMesh::QClough::init_1D(), libMesh::QConical::init_1D(), libMesh::QNodal::init_1D(), libMesh::QMonomial::init_1D(), libMesh::QGaussLobatto::init_2D(), libMesh::QClough::init_2D(), init_2D(), libMesh::QNodal::init_2D(), libMesh::QMonomial::init_2D(), libMesh::QGaussLobatto::init_3D(), init_3D(), libMesh::QSimpson::init_3D(), libMesh::QTrap::init_3D(), libMesh::QGrid::init_3D(), libMesh::QNodal::init_3D(), libMesh::QMonomial::init_3D(), libMesh::OverlapCoupling::operator()(), libMesh::QClough::QClough(), libMesh::QConical::QConical(), QGauss(), libMesh::QGaussLobatto::QGaussLobatto(), libMesh::QGrid::QGrid(), libMesh::QGrundmann_Moller::QGrundmann_Moller(), libMesh::QJacobi::QJacobi(), libMesh::QMonomial::QMonomial(), libMesh::QNodal::QNodal(), libMesh::QSimpson::QSimpson(), libMesh::QTrap::QTrap(), and libMesh::FE< Dim, LAGRANGE_VEC >::reinit_default_dual_shape_coeffs().

67 {
68  libmesh_assert_equal_to(elem.dim(), _dim);
69 
70  // Default to the element p_level() value
71  if (p == invalid_uint)
72  p = elem.p_level();
73 
74  ElemType t = elem.type();
75 
76  // check to see if we have already
77  // done the work for this quadrature rule
78  //
79  // If we have something like a Polygon subclass then we're going to
80  // need to recompute to be safe; even if we're using the same
81  // element, it might have been distorted enough that its subtriangle
82  // triangulation has been changed.
83  if (t == _type && p == _p_level && !elem.runtime_topology())
84  return;
85  else
86  {
87  _elem = &elem;
88  _type = t;
89  _p_level = p;
90  }
91 
92  switch(_elem->dim())
93  {
94  case 0:
95  this->init_0D();
96 
97  return;
98 
99  case 1:
100  this->init_1D();
101 
102  return;
103 
104  case 2:
105  this->init_2D();
106 
107  return;
108 
109  case 3:
110  this->init_3D();
111 
112  return;
113 
114  default:
115  libmesh_error_msg("Invalid dimension _dim = " << _dim);
116  }
117 }
ElemType
Defines an enum for geometric element types.
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
virtual void init_2D()
Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:208
ElemType _type
The type of element for which the current values have been computed.
Definition: quadrature.h:391
unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:379
const Elem * _elem
The element for which the current values were computed, or nullptr if values were computed without a ...
Definition: quadrature.h:397
unsigned int _p_level
The p-level of the element for which the current values have been computed.
Definition: quadrature.h:403
virtual void init_0D()
Initializes the 0D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:198
virtual void init_1D()=0
Initializes the 1D quadrature rule by filling the points and weights vectors with the appropriate val...
virtual unsigned short dim() const =0
virtual void init_3D()
Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:215

◆ init() [2/4]

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

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

Some types, such as Polygon subclasses, might require more detailed element information and so might not be compatible with this API.

New code should use the Elem-based API for most use cases, but some code may initialize quadrature rules with simple ElemType values like triangles and edges for use in tensor product or conical product constructions; this code can set the simple_type_only flag to avoid being identified as deprecated.

Definition at line 121 of file quadrature.C.

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

124 {
125  // Some element types require data from a specific element, so can
126  // only be used with newer APIs.
127  if (t == C0POLYGON || t == C0POLYHEDRON)
128  libmesh_error_msg("Code (see stack trace) used an outdated quadrature function overload.\n"
129  "Quadrature rules on a C0Polygon are not defined by its ElemType alone.");
130 
131  // This API is dangerous to use on general meshes, which may include
132  // element types where the desired quadrature depends on the
133  // physical element, but we still want to be able to initialize
134  // based on only a type for certain simple cases
135  if (t != EDGE2 && !simple_type_only)
136  libmesh_deprecated();
137 
138  // check to see if we have already
139  // done the work for this quadrature rule
140  if (t == _type && p == _p_level)
141  return;
142  else
143  {
144  _elem = nullptr;
145  _type = t;
146  _p_level = p;
147  }
148 
149  switch(_dim)
150  {
151  case 0:
152  this->init_0D();
153 
154  return;
155 
156  case 1:
157  this->init_1D();
158 
159  return;
160 
161  case 2:
162  this->init_2D();
163 
164  return;
165 
166  case 3:
167  this->init_3D();
168 
169  return;
170 
171  default:
172  libmesh_error_msg("Invalid dimension _dim = " << _dim);
173  }
174 }
virtual void init_2D()
Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:208
ElemType _type
The type of element for which the current values have been computed.
Definition: quadrature.h:391
unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:379
const Elem * _elem
The element for which the current values were computed, or nullptr if values were computed without a ...
Definition: quadrature.h:397
unsigned int _p_level
The p-level of the element for which the current values have been computed.
Definition: quadrature.h:403
virtual void init_0D()
Initializes the 0D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:198
virtual void init_1D()=0
Initializes the 1D quadrature rule by filling the points and weights vectors with the appropriate val...
virtual void init_3D()
Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:215

◆ init() [3/4]

void libMesh::QBase::init ( const QBase other_rule)
virtualinherited

Initializes the data structures for a quadrature rule based on the element, element type, and p_level settings of other_rule.

Definition at line 178 of file quadrature.C.

References libMesh::QBase::_elem, libMesh::QBase::_p_level, libMesh::QBase::_type, and libMesh::QBase::init().

179 {
180  if (other_rule._elem)
181  this->init(*other_rule._elem, other_rule._p_level);
182  else
183  this->init(other_rule._type, other_rule._p_level, true);
184 }
virtual void init(const Elem &e, unsigned int p_level=invalid_uint)
Initializes the data structures for a quadrature rule for the element e.
Definition: quadrature.C:65

◆ init() [4/4]

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 188 of file quadrature.C.

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

191 {
192  // dispatch generic implementation
193  this->init(elem.type(), p_level);
194 }
virtual void init(const Elem &e, unsigned int p_level=invalid_uint)
Initializes the data structures for a quadrature rule for the element e.
Definition: quadrature.C:65

◆ init_0D()

void libMesh::QBase::init_0D ( )
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.

Definition at line 198 of file quadrature.C.

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

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

199 {
200  _points.resize(1);
201  _weights.resize(1);
202  _points[0] = Point(0.);
203  _weights[0] = 1.0;
204 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415

◆ init_1D()

void libMesh::QGauss::init_1D ( )
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.

Implements libMesh::QBase.

Definition at line 30 of file quadrature_gauss_1D.C.

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.

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) = -std::sqrt(Real(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) = -7.7459666924148337703585307995648e-01_R;
69  _points[ 1](0) = 0.;
70  _points[ 2] = -_points[0];
71 
72  _weights[ 0] = Real(5)/9;
73  _weights[ 1] = Real(8)/9;
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) = -8.6113631159405257522394648889281e-01_R;
85  _points[ 1](0) = -3.3998104358485626480266575910324e-01_R;
86  _points[ 2] = -_points[1];
87  _points[ 3] = -_points[0];
88 
89  _weights[ 0] = 3.4785484513745385737306394922200e-01_R;
90  _weights[ 1] = 6.5214515486254614262693605077800e-01_R;
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) = -9.0617984593866399279762687829939e-01_R;
103  _points[ 1](0) = -5.3846931010568309103631442070021e-01_R;
104  _points[ 2](0) = 0.;
105  _points[ 3] = -_points[1];
106  _points[ 4] = -_points[0];
107 
108  _weights[ 0] = 2.3692688505618908751426404071992e-01_R;
109  _weights[ 1] = 4.7862867049936646804129151483564e-01_R;
110  _weights[ 2] = 5.6888888888888888888888888888889e-01_R;
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) = -9.3246951420315202781230155449399e-01_R;
123  _points[ 1](0) = -6.6120938646626451366139959501991e-01_R;
124  _points[ 2](0) = -2.3861918608319690863050172168071e-01_R;
125  _points[ 3] = -_points[2];
126  _points[ 4] = -_points[1];
127  _points[ 5] = -_points[0];
128 
129  _weights[ 0] = 1.7132449237917034504029614217273e-01_R;
130  _weights[ 1] = 3.6076157304813860756983351383772e-01_R;
131  _weights[ 2] = 4.6791393457269104738987034398955e-01_R;
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) = -9.4910791234275852452618968404785e-01_R;
145  _points[ 1](0) = -7.4153118559939443986386477328079e-01_R;
146  _points[ 2](0) = -4.0584515137739716690660641207696e-01_R;
147  _points[ 3](0) = 0.;
148  _points[ 4] = -_points[2];
149  _points[ 5] = -_points[1];
150  _points[ 6] = -_points[0];
151 
152  _weights[ 0] = 1.2948496616886969327061143267908e-01_R;
153  _weights[ 1] = 2.7970539148927666790146777142378e-01_R;
154  _weights[ 2] = 3.8183005050511894495036977548898e-01_R;
155  _weights[ 3] = 4.1795918367346938775510204081633e-01_R;
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) = -9.6028985649753623168356086856947e-01_R;
169  _points[ 1](0) = -7.9666647741362673959155393647583e-01_R;
170  _points[ 2](0) = -5.2553240991632898581773904918925e-01_R;
171  _points[ 3](0) = -1.8343464249564980493947614236018e-01_R;
172  _points[ 4] = -_points[3];
173  _points[ 5] = -_points[2];
174  _points[ 6] = -_points[1];
175  _points[ 7] = -_points[0];
176 
177  _weights[ 0] = 1.0122853629037625915253135430996e-01_R;
178  _weights[ 1] = 2.2238103445337447054435599442624e-01_R;
179  _weights[ 2] = 3.1370664587788728733796220198660e-01_R;
180  _weights[ 3] = 3.6268378337836198296515044927720e-01_R;
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) = -9.6816023950762608983557620290367e-01_R;
195  _points[ 1](0) = -8.3603110732663579429942978806973e-01_R;
196  _points[ 2](0) = -6.1337143270059039730870203934147e-01_R;
197  _points[ 3](0) = -3.2425342340380892903853801464334e-01_R;
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] = 8.1274388361574411971892158110524e-02_R;
205  _weights[ 1] = 1.8064816069485740405847203124291e-01_R;
206  _weights[ 2] = 2.6061069640293546231874286941863e-01_R;
207  _weights[ 3] = 3.1234707704000284006863040658444e-01_R;
208  _weights[ 4] = 3.3023935500125976316452506928697e-01_R;
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) = -9.7390652851717172007796401208445e-01_R;
223  _points[ 1](0) = -8.6506336668898451073209668842349e-01_R;
224  _points[ 2](0) = -6.7940956829902440623432736511487e-01_R;
225  _points[ 3](0) = -4.3339539412924719079926594316578e-01_R;
226  _points[ 4](0) = -1.4887433898163121088482600112972e-01_R;
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] = 6.6671344308688137593568809893332e-02_R;
234  _weights[ 1] = 1.4945134915058059314577633965770e-01_R;
235  _weights[ 2] = 2.1908636251598204399553493422816e-01_R;
236  _weights[ 3] = 2.6926671930999635509122692156947e-01_R;
237  _weights[ 4] = 2.9552422471475287017389299465134e-01_R;
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) = -9.7822865814605699280393800112286e-01_R;
254  _points[ 1](0) = -8.8706259976809529907515776930393e-01_R;
255  _points[ 2](0) = -7.3015200557404932409341625203115e-01_R;
256  _points[ 3](0) = -5.1909612920681181592572566945861e-01_R;
257  _points[ 4](0) = -2.6954315595234497233153198540086e-01_R;
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] = 5.5668567116173666482753720442549e-02_R;
266  _weights[ 1] = 1.2558036946490462463469429922394e-01_R;
267  _weights[ 2] = 1.8629021092773425142609764143166e-01_R;
268  _weights[ 3] = 2.3319376459199047991852370484318e-01_R;
269  _weights[ 4] = 2.6280454451024666218068886989051e-01_R;
270  _weights[ 5] = 2.7292508677790063071448352833634e-01_R;
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) = -9.8156063424671925069054909014928e-01_R;
287  _points[ 1](0) = -9.0411725637047485667846586611910e-01_R;
288  _points[ 2](0) = -7.6990267419430468703689383321282e-01_R;
289  _points[ 3](0) = -5.8731795428661744729670241894053e-01_R;
290  _points[ 4](0) = -3.6783149899818019375269153664372e-01_R;
291  _points[ 5](0) = -1.2523340851146891547244136946385e-01_R;
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] = 4.7175336386511827194615961485017e-02_R;
300  _weights[ 1] = 1.0693932599531843096025471819400e-01_R;
301  _weights[ 2] = 1.6007832854334622633465252954336e-01_R;
302  _weights[ 3] = 2.0316742672306592174906445580980e-01_R;
303  _weights[ 4] = 2.3349253653835480876084989892483e-01_R;
304  _weights[ 5] = 2.4914704581340278500056243604295e-01_R;
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) = -9.8418305471858814947282944880711e-01_R;
322  _points[ 1](0) = -9.1759839922297796520654783650072e-01_R;
323  _points[ 2](0) = -8.0157809073330991279420648958286e-01_R;
324  _points[ 3](0) = -6.4234933944034022064398460699552e-01_R;
325  _points[ 4](0) = -4.4849275103644685287791285212764e-01_R;
326  _points[ 5](0) = -2.3045831595513479406552812109799e-01_R;
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] = 4.0484004765315879520021592200986e-02_R;
336  _weights[ 1] = 9.2121499837728447914421775953797e-02_R;
337  _weights[ 2] = 1.3887351021978723846360177686887e-01_R;
338  _weights[ 3] = 1.7814598076194573828004669199610e-01_R;
339  _weights[ 4] = 2.0781604753688850231252321930605e-01_R;
340  _weights[ 5] = 2.2628318026289723841209018603978e-01_R;
341  _weights[ 6] = 2.3255155323087391019458951526884e-01_R;
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) = -9.8628380869681233884159726670405e-01_R;
359  _points[ 1](0) = -9.2843488366357351733639113937787e-01_R;
360  _points[ 2](0) = -8.2720131506976499318979474265039e-01_R;
361  _points[ 3](0) = -6.8729290481168547014801980301933e-01_R;
362  _points[ 4](0) = -5.1524863635815409196529071855119e-01_R;
363  _points[ 5](0) = -3.1911236892788976043567182416848e-01_R;
364  _points[ 6](0) = -1.0805494870734366206624465021983e-01_R;
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] = 3.5119460331751863031832876138192e-02_R;
374  _weights[ 1] = 8.0158087159760209805633277062854e-02_R;
375  _weights[ 2] = 1.2151857068790318468941480907248e-01_R;
376  _weights[ 3] = 1.5720316715819353456960193862384e-01_R;
377  _weights[ 4] = 1.8553839747793781374171659012516e-01_R;
378  _weights[ 5] = 2.0519846372129560396592406566122e-01_R;
379  _weights[ 6] = 2.1526385346315779019587644331626e-01_R;
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) = -9.8799251802048542848956571858661e-01_R;
398  _points[ 1](0) = -9.3727339240070590430775894771021e-01_R;
399  _points[ 2](0) = -8.4820658341042721620064832077422e-01_R;
400  _points[ 3](0) = -7.2441773136017004741618605461394e-01_R;
401  _points[ 4](0) = -5.7097217260853884753722673725391e-01_R;
402  _points[ 5](0) = -3.9415134707756336989720737098105e-01_R;
403  _points[ 6](0) = -2.0119409399743452230062830339460e-01_R;
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] = 3.0753241996117268354628393577204e-02_R;
414  _weights[ 1] = 7.0366047488108124709267416450667e-02_R;
415  _weights[ 2] = 1.0715922046717193501186954668587e-01_R;
416  _weights[ 3] = 1.3957067792615431444780479451103e-01_R;
417  _weights[ 4] = 1.6626920581699393355320086048121e-01_R;
418  _weights[ 5] = 1.8616100001556221102680056186642e-01_R;
419  _weights[ 6] = 1.9843148532711157645611832644384e-01_R;
420  _weights[ 7] = 2.0257824192556127288062019996752e-01_R;
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) = -9.8940093499164993259615417345033e-01_R;
439  _points[ 1](0) = -9.4457502307323257607798841553461e-01_R;
440  _points[ 2](0) = -8.6563120238783174388046789771239e-01_R;
441  _points[ 3](0) = -7.5540440835500303389510119484744e-01_R;
442  _points[ 4](0) = -6.1787624440264374844667176404879e-01_R;
443  _points[ 5](0) = -4.5801677765722738634241944298358e-01_R;
444  _points[ 6](0) = -2.8160355077925891323046050146050e-01_R;
445  _points[ 7](0) = -9.5012509837637440185319335424958e-02_R;
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] = 2.7152459411754094851780572456018e-02_R;
456  _weights[ 1] = 6.2253523938647892862843836994378e-02_R;
457  _weights[ 2] = 9.5158511682492784809925107602246e-02_R;
458  _weights[ 3] = 1.2462897125553387205247628219202e-01_R;
459  _weights[ 4] = 1.4959598881657673208150173054748e-01_R;
460  _weights[ 5] = 1.6915651939500253818931207903033e-01_R;
461  _weights[ 6] = 1.8260341504492358886676366796922e-01_R;
462  _weights[ 7] = 1.8945061045506849628539672320828e-01_R;
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) = -9.9057547531441733567543401994067e-01_R;
482  _points[ 1](0) = -9.5067552176876776122271695789580e-01_R;
483  _points[ 2](0) = -8.8023915372698590212295569448816e-01_R;
484  _points[ 3](0) = -7.8151400389680140692523005552048e-01_R;
485  _points[ 4](0) = -6.5767115921669076585030221664300e-01_R;
486  _points[ 5](0) = -5.1269053708647696788624656862955e-01_R;
487  _points[ 6](0) = -3.5123176345387631529718551709535e-01_R;
488  _points[ 7](0) = -1.7848418149584785585067749365407e-01_R;
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] = 2.4148302868547931960110026287565e-02_R;
500  _weights[ 1] = 5.5459529373987201129440165358245e-02_R;
501  _weights[ 2] = 8.5036148317179180883535370191062e-02_R;
502  _weights[ 3] = 1.1188384719340397109478838562636e-01_R;
503  _weights[ 4] = 1.3513636846852547328631998170235e-01_R;
504  _weights[ 5] = 1.5404576107681028808143159480196e-01_R;
505  _weights[ 6] = 1.6800410215645004450997066378832e-01_R;
506  _weights[ 7] = 1.7656270536699264632527099011320e-01_R;
507  _weights[ 8] = 1.7944647035620652545826564426189e-01_R;
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) = -9.9156516842093094673001600470615e-01_R;
527  _points[ 1](0) = -9.5582394957139775518119589292978e-01_R;
528  _points[ 2](0) = -8.9260246649755573920606059112715e-01_R;
529  _points[ 3](0) = -8.0370495897252311568241745501459e-01_R;
530  _points[ 4](0) = -6.9168704306035320787489108128885e-01_R;
531  _points[ 5](0) = -5.5977083107394753460787154852533e-01_R;
532  _points[ 6](0) = -4.1175116146284264603593179383305e-01_R;
533  _points[ 7](0) = -2.5188622569150550958897285487791e-01_R;
534  _points[ 8](0) = -8.4775013041735301242261852935784e-02_R;
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] = 2.1616013526483310313342710266452e-02_R;
546  _weights[ 1] = 4.9714548894969796453334946202639e-02_R;
547  _weights[ 2] = 7.6425730254889056529129677616637e-02_R;
548  _weights[ 3] = 1.0094204410628716556281398492483e-01_R;
549  _weights[ 4] = 1.2255520671147846018451912680020e-01_R;
550  _weights[ 5] = 1.4064291467065065120473130375195e-01_R;
551  _weights[ 6] = 1.5468467512626524492541800383637e-01_R;
552  _weights[ 7] = 1.6427648374583272298605377646593e-01_R;
553  _weights[ 8] = 1.6914238296314359184065647013499e-01_R;
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) = -9.9240684384358440318901767025326e-01_R;
574  _points[ 1](0) = -9.6020815213483003085277884068765e-01_R;
575  _points[ 2](0) = -9.0315590361481790164266092853231e-01_R;
576  _points[ 3](0) = -8.2271465653714282497892248671271e-01_R;
577  _points[ 4](0) = -7.2096617733522937861709586082378e-01_R;
578  _points[ 5](0) = -6.0054530466168102346963816494624e-01_R;
579  _points[ 6](0) = -4.6457074137596094571726714810410e-01_R;
580  _points[ 7](0) = -3.1656409996362983199011732884984e-01_R;
581  _points[ 8](0) = -1.6035864564022537586809611574074e-01_R;
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] = 1.9461788229726477036312041464438e-02_R;
594  _weights[ 1] = 4.4814226765699600332838157401994e-02_R;
595  _weights[ 2] = 6.9044542737641226580708258006013e-02_R;
596  _weights[ 3] = 9.1490021622449999464462094123840e-02_R;
597  _weights[ 4] = 1.1156664554733399471602390168177e-01_R;
598  _weights[ 5] = 1.2875396253933622767551578485688e-01_R;
599  _weights[ 6] = 1.4260670217360661177574610944190e-01_R;
600  _weights[ 7] = 1.5276604206585966677885540089766e-01_R;
601  _weights[ 8] = 1.5896884339395434764995643946505e-01_R;
602  _weights[ 9] = 1.6105444984878369597916362532092e-01_R;
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) = -9.9312859918509492478612238847132e-01_R;
623  _points[ 1](0) = -9.6397192727791379126766613119728e-01_R;
624  _points[ 2](0) = -9.1223442825132590586775244120330e-01_R;
625  _points[ 3](0) = -8.3911697182221882339452906170152e-01_R;
626  _points[ 4](0) = -7.4633190646015079261430507035564e-01_R;
627  _points[ 5](0) = -6.3605368072651502545283669622629e-01_R;
628  _points[ 6](0) = -5.1086700195082709800436405095525e-01_R;
629  _points[ 7](0) = -3.7370608871541956067254817702493e-01_R;
630  _points[ 8](0) = -2.2778585114164507808049619536857e-01_R;
631  _points[ 9](0) = -7.6526521133497333754640409398838e-02_R;
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] = 1.7614007139152118311861962351853e-02_R;
644  _weights[ 1] = 4.0601429800386941331039952274932e-02_R;
645  _weights[ 2] = 6.2672048334109063569506535187042e-02_R;
646  _weights[ 3] = 8.3276741576704748724758143222046e-02_R;
647  _weights[ 4] = 1.0193011981724043503675013548035e-01_R;
648  _weights[ 5] = 1.1819453196151841731237737771138e-01_R;
649  _weights[ 6] = 1.3168863844917662689849449974816e-01_R;
650  _weights[ 7] = 1.4209610931838205132929832506716e-01_R;
651  _weights[ 8] = 1.4917298647260374678782873700197e-01_R;
652  _weights[ 9] = 1.5275338713072585069808433195510e-01_R;
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) = -9.9375217062038950026024203593794e-01_R;
674  _points[ 1](0) = -9.6722683856630629431662221490770e-01_R;
675  _points[ 2](0) = -9.2009933415040082879018713371497e-01_R;
676  _points[ 3](0) = -8.5336336458331728364725063858757e-01_R;
677  _points[ 4](0) = -7.6843996347567790861587785130623e-01_R;
678  _points[ 5](0) = -6.6713880419741231930596666999034e-01_R;
679  _points[ 6](0) = -5.5161883588721980705901879672431e-01_R;
680  _points[ 7](0) = -4.2434212020743878357366888854379e-01_R;
681  _points[ 8](0) = -2.8802131680240109660079251606460e-01_R;
682  _points[ 9](0) = -1.4556185416089509093703098233869e-01_R;
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] = 1.6017228257774333324224616858471e-02_R;
696  _weights[ 1] = 3.6953789770852493799950668299330e-02_R;
697  _weights[ 2] = 5.7134425426857208283635826472448e-02_R;
698  _weights[ 3] = 7.6100113628379302017051653300183e-02_R;
699  _weights[ 4] = 9.3444423456033861553289741113932e-02_R;
700  _weights[ 5] = 1.0879729916714837766347457807011e-01_R;
701  _weights[ 6] = 1.2183141605372853419536717712572e-01_R;
702  _weights[ 7] = 1.3226893863333746178105257449678e-01_R;
703  _weights[ 8] = 1.3988739479107315472213342386758e-01_R;
704  _weights[ 9] = 1.4452440398997005906382716655375e-01_R;
705  _weights[10] = 1.4608113364969042719198514768337e-01_R;
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) = -9.9429458548239929207303142116130e-01_R;
727  _points[ 1](0) = -9.7006049783542872712395098676527e-01_R;
728  _points[ 2](0) = -9.2695677218717400052069293925905e-01_R;
729  _points[ 3](0) = -8.6581257772030013653642563701938e-01_R;
730  _points[ 4](0) = -7.8781680597920816200427795540835e-01_R;
731  _points[ 5](0) = -6.9448726318668278005068983576226e-01_R;
732  _points[ 6](0) = -5.8764040350691159295887692763865e-01_R;
733  _points[ 7](0) = -4.6935583798675702640633071096641e-01_R;
734  _points[ 8](0) = -3.4193582089208422515814742042738e-01_R;
735  _points[ 9](0) = -2.0786042668822128547884653391955e-01_R;
736  _points[10](0) = -6.9739273319722221213841796118628e-02_R;
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] = 1.4627995298272200684991098047185e-02_R;
750  _weights[ 1] = 3.3774901584814154793302246865913e-02_R;
751  _weights[ 2] = 5.2293335152683285940312051273211e-02_R;
752  _weights[ 3] = 6.9796468424520488094961418930218e-02_R;
753  _weights[ 4] = 8.5941606217067727414443681372703e-02_R;
754  _weights[ 5] = 1.0041414444288096493207883783054e-01_R;
755  _weights[ 6] = 1.1293229608053921839340060742175e-01_R;
756  _weights[ 7] = 1.2325237681051242428556098615481e-01_R;
757  _weights[ 8] = 1.3117350478706237073296499253031e-01_R;
758  _weights[ 9] = 1.3654149834601517135257383123152e-01_R;
759  _weights[10] = 1.3925187285563199337541024834181e-01_R;
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 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415
Order get_order() const
Definition: quadrature.h:249
Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:385
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ init_2D()

void libMesh::QGauss::init_2D ( )
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.

Reimplemented from libMesh::QBase.

Definition at line 31 of file quadrature_gauss_2D.C.

References libMesh::QBase::_elem, libMesh::QBase::_order, libMesh::QBase::_p_level, libMesh::QBase::_points, libMesh::QBase::_type, libMesh::QBase::_weights, libMesh::C0POLYGON, libMesh::CONSTANT, dunavant_rule(), dunavant_rule2(), 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::libmesh_assert(), libMesh::Polygon::n_subtriangles(), libMesh::NINETEENTH, libMesh::NINTH, libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::QUADSHELL4, libMesh::QUADSHELL8, libMesh::QUADSHELL9, libMesh::Real, libMesh::SECOND, libMesh::SEVENTEENTH, libMesh::SEVENTH, libMesh::SIXTEENTH, libMesh::SIXTH, libMesh::QBase::tensor_product_quad(), libMesh::TENTH, libMesh::THIRD, libMesh::THIRTEENTH, libMesh::THIRTIETH, libMesh::TRI3, libMesh::TRI3SUBDIVISION, libMesh::TRI6, libMesh::TRI7, libMesh::TRISHELL3, libMesh::TWELFTH, libMesh::TWENTIETH, libMesh::TWENTYEIGHTH, libMesh::TWENTYFIFTH, libMesh::TWENTYFOURTH, libMesh::TWENTYNINTH, libMesh::TWENTYSECOND, libMesh::TWENTYSEVENTH, libMesh::TWENTYSIXTH, libMesh::TWENTYTHIRD, and libMesh::Elem::type().

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

◆ init_3D()

void libMesh::QGauss::init_3D ( )
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.

Reimplemented from libMesh::QBase.

Definition at line 30 of file quadrature_gauss_3D.C.

References libMesh::QBase::_elem, libMesh::QBase::_order, libMesh::QBase::_p_level, libMesh::QBase::_points, libMesh::QBase::_type, libMesh::QBase::_weights, libMesh::QBase::allow_rules_with_negative_weights, libMesh::C0POLYHEDRON, libMesh::CONSTANT, 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::libmesh_assert(), libMesh::Polyhedron::n_subelements(), libMesh::PRISM15, libMesh::PRISM18, libMesh::PRISM20, libMesh::PRISM21, libMesh::PRISM6, libMesh::PYRAMID13, libMesh::PYRAMID14, libMesh::PYRAMID18, libMesh::PYRAMID5, libMesh::Real, libMesh::SECOND, libMesh::SEVENTH, libMesh::SIXTH, libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), libMesh::TET10, libMesh::TET14, libMesh::TET4, libMesh::THIRD, libMesh::TRI3, libMesh::triple_product(), and libMesh::Elem::type().

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

◆ 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 43 of file quadrature_gauss.C.

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

Referenced by init_3D().

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

◆ n_objects()

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

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

Definition at line 85 of file reference_counter.h.

References libMesh::ReferenceCounter::_n_objects.

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

86  { return _n_objects; }
static Threads::atomic< unsigned int > _n_objects
The number of objects.

◆ n_points()

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

Definition at line 131 of file quadrature.h.

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(), AssemblyA1::boundary_assembly(), AssemblyF0::boundary_assembly(), AssemblyF1::boundary_assembly(), AssemblyA2::boundary_assembly(), AssemblyF2::boundary_assembly(), A2::boundary_assembly(), A3::boundary_assembly(), F0::boundary_assembly(), Output0::boundary_assembly(), compute_enriched_soln(), libMesh::FirstOrderUnsteadySolver::compute_second_order_eqns(), libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), libMesh::HDGProblem::create_identity_jacobian(), libMesh::HDGProblem::create_identity_residual(), SecondOrderScalarSystemSecondOrderTimeSolverBase::damping_residual(), SecondOrderScalarSystemFirstOrderTimeSolverBase::damping_residual(), NavierSystem::element_constraint(), CoupledSystem::element_constraint(), PoissonSystem::element_postprocess(), LaplaceSystem::element_postprocess(), LaplaceQoI::element_qoi(), HeatSystem::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(), CoupledSystem::element_time_derivative(), HilbertSystem::element_time_derivative(), HeatSystem::element_time_derivative(), FirstOrderScalarSystemBase::element_time_derivative(), SecondOrderScalarSystemFirstOrderTimeSolverBase::element_time_derivative(), fe_assembly(), A0::interior_assembly(), B::interior_assembly(), M0::interior_assembly(), A1::interior_assembly(), AssemblyA0::interior_assembly(), AcousticsInnerProduct::interior_assembly(), AssemblyA1::interior_assembly(), A2::interior_assembly(), AssemblyA2::interior_assembly(), F0::interior_assembly(), OutputAssembly::interior_assembly(), EIM_IP_assembly::interior_assembly(), EIM_F::interior_assembly(), AssemblyEIM::interior_assembly(), InnerProductAssembly::interior_assembly(), AssemblyF0::interior_assembly(), AssemblyF1::interior_assembly(), Ex6InnerProduct::interior_assembly(), LaplaceYoung::jacobian(), NavierSystem::mass_residual(), ElasticitySystem::mass_residual(), libMesh::FEMPhysics::mass_residual(), FirstOrderScalarSystemBase::mass_residual(), SecondOrderScalarSystemSecondOrderTimeSolverBase::mass_residual(), SecondOrderScalarSystemFirstOrderTimeSolverBase::mass_residual(), libMesh::FEAbstract::n_quadrature_points(), 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::size(), libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), libMesh::QBase::tensor_product_quad(), and QuadratureTest::testPolynomial().

132  {
133  libmesh_assert (!_points.empty());
134  return cast_int<unsigned int>(_points.size());
135  }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
libmesh_assert(ctx)

◆ 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::ReferenceCounter::print_info ( std::ostream &  out_stream = libMesh::out)
staticinherited

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

Definition at line 81 of file reference_counter.C.

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

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

82 {
84  out_stream << ReferenceCounter::get_info();
85 }
static std::string get_info()
Gets a string containing the reference information.
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called...

◆ print_info() [2/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 43 of file quadrature.C.

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

Referenced by libMesh::operator<<().

44 {
45  libmesh_assert(!_points.empty());
46  libmesh_assert(!_weights.empty());
47 
48  Real summed_weights=0;
49  os << "N_Q_Points=" << this->n_points() << std::endl << std::endl;
50  for (auto qpoint: index_range(_points))
51  {
52  os << " Point " << qpoint << ":\n"
53  << " "
54  << _points[qpoint]
55  << "\n Weight:\n "
56  << " w=" << _weights[qpoint] << "\n" << std::endl;
57 
58  summed_weights += _weights[qpoint];
59  }
60  os << "Summed Weights: " << summed_weights << std::endl;
61 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415
libmesh_assert(ctx)
unsigned int n_points() const
Definition: quadrature.h:131
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117

◆ qp()

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

Definition at line 179 of file quadrature.h.

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(), libMesh::QBase::tensor_product_quad(), and QuadratureTest::testPolynomial().

180  {
181  libmesh_assert_less (i, _points.size());
182  return _points[i];
183  }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409

◆ 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 222 of file quadrature.C.

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().

224 {
225  // Make sure we are in 1D
226  libmesh_assert_equal_to (_dim, 1);
227 
228  Real
229  h_new = new_range.second - new_range.first,
230  h_old = old_range.second - old_range.first;
231 
232  // Make sure that we have sane ranges
233  libmesh_assert_greater (h_new, 0.);
234  libmesh_assert_greater (h_old, 0.);
235 
236  // Make sure there are some points
237  libmesh_assert_greater (_points.size(), 0);
238 
239  // Compute the scale factor
240  Real scfact = h_new/h_old;
241 
242  // We're mapping from old_range -> new_range
243  for (auto i : index_range(_points))
244  {
245  _points[i](0) = new_range.first +
246  (_points[i](0) - old_range.first) * scfact;
247 
248  // Scale the weights
249  _weights[i] *= scfact;
250  }
251 }
unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:379
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117

◆ 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 286 of file quadrature.h.

286 { return false; }

◆ size()

unsigned int libMesh::QBase::size ( ) const
inlineinherited

Alias for n_points() to enable use in index_range.

Returns
The number of points associated with the quadrature rule.

Definition at line 142 of file quadrature.h.

References libMesh::QBase::n_points().

Referenced by libMesh::Elem::true_centroid().

143  {
144  return n_points();
145  }
unsigned int n_points() const
Definition: quadrature.h:131

◆ 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 283 of file quadrature.C.

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(), init_3D(), libMesh::QSimpson::init_3D(), and libMesh::QGrid::init_3D().

284 {
285  const unsigned int np = q1D.n_points();
286 
287  _points.resize(np * np * np);
288 
289  _weights.resize(np * np * np);
290 
291  unsigned int q=0;
292 
293  for (unsigned int k=0; k<np; k++)
294  for (unsigned int j=0; j<np; j++)
295  for (unsigned int i=0; i<np; i++)
296  {
297  _points[q](0) = q1D.qp(i)(0);
298  _points[q](1) = q1D.qp(j)(0);
299  _points[q](2) = q1D.qp(k)(0);
300 
301  _weights[q] = q1D.w(i) * q1D.w(j) * q1D.w(k);
302 
303  q++;
304  }
305 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415

◆ 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 310 of file quadrature.C.

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

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

311 {
312  const unsigned int n_points1D = q1D.n_points();
313  const unsigned int n_points2D = q2D.n_points();
314 
315  _points.resize (n_points1D * n_points2D);
316  _weights.resize (n_points1D * n_points2D);
317 
318  unsigned int q=0;
319 
320  for (unsigned int j=0; j<n_points1D; j++)
321  for (unsigned int i=0; i<n_points2D; i++)
322  {
323  _points[q](0) = q2D.qp(i)(0);
324  _points[q](1) = q2D.qp(i)(1);
325  _points[q](2) = q1D.qp(j)(0);
326 
327  _weights[q] = q2D.w(i) * q1D.w(j);
328 
329  q++;
330  }
331 
332 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415

◆ 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 256 of file quadrature.C.

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(), init_2D(), libMesh::QSimpson::init_2D(), and libMesh::QGrid::init_2D().

257 {
258 
259  const unsigned int np = q1D.n_points();
260 
261  _points.resize(np * np);
262 
263  _weights.resize(np * np);
264 
265  unsigned int q=0;
266 
267  for (unsigned int j=0; j<np; j++)
268  for (unsigned int i=0; i<np; i++)
269  {
270  _points[q](0) = q1D.qp(i)(0);
271  _points[q](1) = q1D.qp(j)(0);
272 
273  _weights[q] = q1D.w(i)*q1D.w(j);
274 
275  q++;
276  }
277 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415

◆ type()

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

Implements libMesh::QBase.

Definition at line 33 of file quadrature_gauss.C.

References libMesh::QGAUSS.

34 {
35  return QGAUSS;
36 }

◆ w()

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

Definition at line 188 of file quadrature.h.

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(), libMesh::QBase::tensor_product_quad(), and QuadratureTest::testPolynomial().

189  {
190  libmesh_assert_less (i, _weights.size());
191  return _weights[i];
192  }
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415

Member Data Documentation

◆ _counts

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

Actually holds the data.

Definition at line 124 of file reference_counter.h.

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

◆ _dim

unsigned int libMesh::QBase::_dim
protectedinherited

◆ _elem

const Elem* libMesh::QBase::_elem
protectedinherited

The element for which the current values were computed, or nullptr if values were computed without a specific element.

Definition at line 397 of file quadrature.h.

Referenced by libMesh::QBase::init(), init_2D(), libMesh::QNodal::init_2D(), init_3D(), and libMesh::QNodal::init_3D().

◆ _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 143 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 137 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 132 of file reference_counter.h.

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

◆ _order

Order libMesh::QBase::_order
protectedinherited

◆ _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 403 of file quadrature.h.

Referenced by libMesh::QBase::get_order(), libMesh::QBase::get_p_level(), libMesh::QBase::init(), libMesh::QClough::init_2D(), init_2D(), init_3D(), libMesh::QSimpson::init_3D(), and libMesh::QTrap::init_3D().

◆ _points

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

◆ _type

ElemType libMesh::QBase::_type
protectedinherited

◆ _weights

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

◆ allow_nodal_pyramid_quadrature

bool libMesh::QBase::allow_nodal_pyramid_quadrature
inherited

The flag's value defaults to false so that one does not accidentally use a nodal quadrature rule on Pyramid elements, since evaluating the inverse element Jacobian (e.g.

dphi) is not well-defined at the Pyramid apex because the element Jacobian is zero there.

We do not want to completely prevent someone from using a nodal quadrature rule on Pyramids, however, since there are legitimate use cases (lumped mass matrix) so the flag can be set to true to override this behavior.

Definition at line 314 of file quadrature.h.

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

◆ 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 301 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: