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 (const unsigned int _dim, const Order _order=INVALID_ORDER)
 Constructor. More...
 
 ~QGauss ()
 Destructor. More...
 
virtual QuadratureType type () const libmesh_override
 
ElemType get_elem_type () const
 
unsigned int get_p_level () const
 
unsigned int n_points () const
 
unsigned int get_dim () const
 
const std::vector< Point > & get_points () const
 
std::vector< Point > & get_points ()
 
const std::vector< Real > & get_weights () const
 
std::vector< Real > & get_weights ()
 
Point qp (const unsigned int i) const
 
Real w (const unsigned int i) const
 
virtual void init (const ElemType type=INVALID_ELEM, unsigned int p_level=0)
 Initializes the data structures for a quadrature rule for an element of type type. More...
 
virtual void init (const Elem &elem, const std::vector< Real > &vertex_distance_func, unsigned int p_level=0)
 Initializes the data structures for an element potentially "cut" by a signed distance function. More...
 
Order get_order () const
 
void print_info (std::ostream &os=libMesh::out) const
 Prints information relevant to the quadrature rule, by default to libMesh::out. More...
 
void scale (std::pair< Real, Real > old_range, std::pair< Real, Real > new_range)
 Maps the points of a 1D quadrature rule defined by "old_range" to another 1D interval defined by "new_range" and scales the weights accordingly. More...
 
virtual bool shapes_need_reinit ()
 

Static Public Member Functions

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

Public Attributes

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

Protected Types

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

Protected Member Functions

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

Protected Attributes

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

Static Protected Attributes

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

Private Member Functions

virtual void init_1D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0) libmesh_override
 Initializes the 1D quadrature rule by filling the points and weights vectors with the appropriate values. More...
 
virtual void init_2D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0) libmesh_override
 Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate values. More...
 
virtual void init_3D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0) libmesh_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

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

libMesh::QGauss::QGauss ( const unsigned int  _dim,
const Order  _order = INVALID_ORDER 
)

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 libMesh::EDGE2, and libMesh::QBase::init().

53  :
54  QBase(_dim, _order)
55  {
56  if (_dim == 1)
57  init(EDGE2);
58  }
virtual void init(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Initializes the data structures for a quadrature rule for an element of type type.
Definition: quadrature.C:28
const unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:311
QBase(const unsigned int _dim, const Order _order=INVALID_ORDER)
Constructor.
Definition: quadrature.h:350
const Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:317
libMesh::QGauss::~QGauss ( )

Destructor.

Definition at line 63 of file quadrature_gauss.h.

63 {}

Member Function Documentation

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

Builds a specific quadrature rule based on the name string.

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

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

Definition at line 42 of file quadrature_build.C.

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

45 {
46  return QBase::build (Utility::string_to_enum<QuadratureType> (type),
47  _dim,
48  _order);
49 }
const unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:311
virtual QuadratureType type() const =0
static UniquePtr< QBase > build(const std::string &name, const unsigned int dim, const Order order=INVALID_ORDER)
Builds a specific quadrature rule based on the name string.
const Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:317
UniquePtr< 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 UniquePtr<QBase> is returned so that the user does not accidentally leak it.

Definition at line 53 of file quadrature_build.C.

References 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::QSIMPSON, libMesh::QTRAP, libMesh::THIRD, and libMesh::TWENTYTHIRD.

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

Definition at line 107 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

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

108 {
109  _enable_print_counter = false;
110  return;
111 }
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called...
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 260 of file quadrature_gauss.C.

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

Referenced by init_2D(), and type().

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

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

Referenced by init_2D(), and type().

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

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

Definition at line 101 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

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

102 {
103  _enable_print_counter = true;
104  return;
105 }
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called...
unsigned int libMesh::QBase::get_dim ( ) const
inherited
Returns
The spatial dimension of the quadrature rule.

Definition at line 122 of file quadrature.h.

References libMesh::QBase::_dim.

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::attach_quadrature_rule(), and libMesh::RBEIMAssembly::evaluate_basis_function().

122 { return _dim; }
const unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:311
ElemType libMesh::QBase::get_elem_type ( ) const
inherited
Returns
The element type we're currently using.

Definition at line 103 of file quadrature.h.

References libMesh::QBase::_type.

103 { return _type; }
ElemType _type
The type of element for which the current values have been computed.
Definition: quadrature.h:323
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 (Counts::iterator it = _counts.begin();
59  it != _counts.end(); ++it)
60  {
61  const std::string name(it->first);
62  const unsigned int creations = it->second.first;
63  const unsigned int destructions = it->second.second;
64 
65  oss << "| " << name << " reference count information:\n"
66  << "| Creations: " << creations << '\n'
67  << "| Destructions: " << destructions << '\n';
68  }
69 
70  oss << " ---------------------------------------------------------------------------- \n";
71 
72  return oss.str();
73 
74 #else
75 
76  return "";
77 
78 #endif
79 }
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:39
static Counts _counts
Actually holds the data.
Order libMesh::QBase::get_order ( ) const
inherited
Returns
The order of the quadrature rule.

Definition at line 189 of file quadrature.h.

References libMesh::QBase::_order, libMesh::QBase::_p_level, libMesh::QBase::operator<<, libMesh::out, libMesh::QBase::print_info(), and libMesh::QBase::scale().

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::attach_quadrature_rule(), and libMesh::RBEIMAssembly::evaluate_basis_function().

189 { return static_cast<Order>(_order + _p_level); }
unsigned int _p_level
The p-level of the element for which the current values have been computed.
Definition: quadrature.h:329
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
const Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:317
unsigned int libMesh::QBase::get_p_level ( ) const
inherited
Returns
The p-refinement level we're currently using.

Definition at line 108 of file quadrature.h.

References libMesh::QBase::_p_level.

108 { return _p_level; }
unsigned int _p_level
The p-level of the element for which the current values have been computed.
Definition: quadrature.h:329
const std::vector<Point>& libMesh::QBase::get_points ( ) const
inherited
Returns
A std::vector containing the quadrature point locations in reference element space.

Definition at line 128 of file quadrature.h.

References libMesh::QBase::_points.

Referenced by assemble_ellipticdg(), and libMesh::FESubdivision::attach_quadrature_rule().

128 { return _points; }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:335
std::vector<Point>& libMesh::QBase::get_points ( )
inherited
Returns
A std::vector containing the quadrature point locations in reference element space as a writable reference.

Definition at line 134 of file quadrature.h.

References libMesh::QBase::_points.

134 { return _points; }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:335
const std::vector<Real>& libMesh::QBase::get_weights ( ) const
inherited
Returns
A constant reference to a std::vector containing the quadrature weights.

Definition at line 140 of file quadrature.h.

References libMesh::QBase::_weights.

Referenced by libMesh::FESubdivision::attach_quadrature_rule().

140 { return _weights; }
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:341
std::vector<Real>& libMesh::QBase::get_weights ( )
inherited
Returns
A writable references to a std::vector containing the quadrature weights.

Definition at line 146 of file quadrature.h.

References libMesh::QBase::_weights.

146 { return _weights; }
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:341
void libMesh::ReferenceCounter::increment_constructor_count ( const std::string &  name)
protectedinherited

Increments the construction counter.

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

Definition at line 185 of file reference_counter.h.

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

Referenced by libMesh::ReferenceCounter::n_objects(), and libMesh::ReferenceCountedObject< RBParametrized >::ReferenceCountedObject().

186 {
187  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
188  std::pair<unsigned int, unsigned int> & p = _counts[name];
189 
190  p.first++;
191 }
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:39
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:29
static Counts _counts
Actually holds the data.
void libMesh::ReferenceCounter::increment_destructor_count ( const std::string &  name)
protectedinherited

Increments the destruction counter.

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

Definition at line 198 of file reference_counter.h.

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

Referenced by libMesh::ReferenceCounter::n_objects(), and libMesh::ReferenceCountedObject< RBParametrized >::~ReferenceCountedObject().

199 {
200  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
201  std::pair<unsigned int, unsigned int> & p = _counts[name];
202 
203  p.second++;
204 }
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:39
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:29
static Counts _counts
Actually holds the data.
void libMesh::QBase::init ( const ElemType  type = INVALID_ELEM,
unsigned int  p_level = 0 
)
virtualinherited

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

Definition at line 28 of file quadrature.C.

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

Referenced by libMesh::FESubdivision::attach_quadrature_rule(), libMesh::QBase::init(), libMesh::QTrap::init_2D(), libMesh::QSimpson::init_2D(), init_2D(), libMesh::QTrap::init_3D(), init_3D(), libMesh::QSimpson::init_3D(), QGauss(), libMesh::QSimpson::QSimpson(), libMesh::QTrap::QTrap(), and libMesh::QBase::w().

30 {
31  // check to see if we have already
32  // done the work for this quadrature rule
33  if (t == _type && p == _p_level)
34  return;
35  else
36  {
37  _type = t;
38  _p_level = p;
39  }
40 
41 
42 
43  switch(_dim)
44  {
45  case 0:
46  this->init_0D(_type,_p_level);
47 
48  return;
49 
50  case 1:
51  this->init_1D(_type,_p_level);
52 
53  return;
54 
55  case 2:
56  this->init_2D(_type,_p_level);
57 
58  return;
59 
60  case 3:
61  this->init_3D(_type,_p_level);
62 
63  return;
64 
65  default:
66  libmesh_error_msg("Invalid dimension _dim = " << _dim);
67  }
68 }
const unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:311
ElemType _type
The type of element for which the current values have been computed.
Definition: quadrature.h:323
unsigned int _p_level
The p-level of the element for which the current values have been computed.
Definition: quadrature.h:329
virtual void init_2D(const ElemType, unsigned int=0)
Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.h:263
virtual void init_1D(const ElemType type=INVALID_ELEM, unsigned int p_level=0)=0
Initializes the 1D quadrature rule by filling the points and weights vectors with the appropriate val...
virtual void init_0D(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Initializes the 0D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:82
virtual void init_3D(const ElemType, unsigned int=0)
Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.h:279
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.

Reimplemented in libMesh::libmesh_final< T >.

Definition at line 72 of file quadrature.C.

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

75 {
76  // dispatch generic implementation
77  this->init(elem.type(), p_level);
78 }
virtual void init(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Initializes the data structures for a quadrature rule for an element of type type.
Definition: quadrature.C:28
void libMesh::QBase::init_0D ( const ElemType  type = INVALID_ELEM,
unsigned int  p_level = 0 
)
protectedvirtualinherited

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

Generally this is just one point with weight 1.

Definition at line 82 of file quadrature.C.

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

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

84 {
85  _points.resize(1);
86  _weights.resize(1);
87  _points[0] = Point(0.);
88  _weights[0] = 1.0;
89 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:335
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:341
void libMesh::QGauss::init_1D ( const ElemType  type = INVALID_ELEM,
unsigned int  p_level = 0 
)
privatevirtual

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::NINETEENTH, libMesh::NINTH, 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.

Referenced by type().

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

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 redefined, gives an error (when DEBUG is defined) when called.

Reimplemented from libMesh::QBase.

Definition at line 28 of file quadrature_gauss_2D.C.

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

Referenced by type().

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

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 redefined, gives an error (when DEBUG is defined) when called.

Reimplemented from libMesh::QBase.

Definition at line 28 of file quadrature_gauss_3D.C.

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

Referenced by type().

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

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

Referenced by init_3D(), and type().

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

Definition at line 113 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(), AssemblyA0::boundary_assembly(), AssemblyF0::boundary_assembly(), AssemblyA1::boundary_assembly(), AssemblyF1::boundary_assembly(), AssemblyF2::boundary_assembly(), A2::boundary_assembly(), AssemblyA2::boundary_assembly(), A3::boundary_assembly(), F0::boundary_assembly(), Output0::boundary_assembly(), libMesh::System::calculate_norm(), libMesh::FirstOrderUnsteadySolver::compute_second_order_eqns(), SecondOrderScalarSystemSecondOrderTimeSolverBase::damping_residual(), SecondOrderScalarSystemFirstOrderTimeSolverBase::damping_residual(), NavierSystem::element_constraint(), CoupledSystem::element_constraint(), PoissonSystem::element_postprocess(), LaplaceSystem::element_postprocess(), LaplaceQoI::element_qoi(), LaplaceQoI::element_qoi_derivative(), LaplaceSystem::element_qoi_derivative(), HeatSystem::element_qoi_derivative(), NavierSystem::element_time_derivative(), SolidSystem::element_time_derivative(), L2System::element_time_derivative(), PoissonSystem::element_time_derivative(), LaplaceSystem::element_time_derivative(), CurlCurlSystem::element_time_derivative(), ElasticitySystem::element_time_derivative(), CoupledSystem::element_time_derivative(), FirstOrderScalarSystemBase::element_time_derivative(), SecondOrderScalarSystemFirstOrderTimeSolverBase::element_time_derivative(), libMesh::RBEIMConstruction::enrich_RB_space(), LaplaceSystem::init_dirichlet_bcs(), B::interior_assembly(), A0::interior_assembly(), M0::interior_assembly(), A1::interior_assembly(), AssemblyA0::interior_assembly(), EIM_IP_assembly::interior_assembly(), AcousticsInnerProduct::interior_assembly(), AssemblyA1::interior_assembly(), A2::interior_assembly(), AssemblyA2::interior_assembly(), EIM_F::interior_assembly(), F0::interior_assembly(), OutputAssembly::interior_assembly(), InnerProductAssembly::interior_assembly(), AssemblyEIM::interior_assembly(), AssemblyF0::interior_assembly(), AssemblyF1::interior_assembly(), Ex6InnerProduct::interior_assembly(), Ex6EIMInnerProduct::interior_assembly(), LaplaceYoung::jacobian(), NavierSystem::mass_residual(), ElasticitySystem::mass_residual(), libMesh::FEMPhysics::mass_residual(), FirstOrderScalarSystemBase::mass_residual(), SecondOrderScalarSystemSecondOrderTimeSolverBase::mass_residual(), SecondOrderScalarSystemFirstOrderTimeSolverBase::mass_residual(), libMesh::QBase::print_info(), LaplaceYoung::residual(), LaplaceSystem::side_constraint(), LaplaceSystem::side_postprocess(), CoupledSystemQoI::side_qoi(), CoupledSystemQoI::side_qoi_derivative(), LaplaceSystem::side_qoi_derivative(), SolidSystem::side_time_derivative(), CurlCurlSystem::side_time_derivative(), ElasticitySystem::side_time_derivative(), libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), libMesh::QBase::tensor_product_quad(), and libMesh::RBEIMConstruction::truth_solve().

114  {
115  libmesh_assert (!_points.empty());
116  return cast_int<unsigned int>(_points.size());
117  }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:335
libmesh_assert(j)
void libMesh::ReferenceCounter::print_info ( std::ostream &  out = libMesh::out)
staticinherited

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

Definition at line 88 of file reference_counter.C.

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

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

89 {
91  out_stream << ReferenceCounter::get_info();
92 }
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...
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 364 of file quadrature.h.

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

Referenced by libMesh::QBase::get_order(), and libMesh::operator<<().

365 {
366  libmesh_assert(!_points.empty());
367  libmesh_assert(!_weights.empty());
368 
369  Real summed_weights=0;
370  os << "N_Q_Points=" << this->n_points() << std::endl << std::endl;
371  for (unsigned int qpoint=0; qpoint<this->n_points(); qpoint++)
372  {
373  os << " Point " << qpoint << ":\n"
374  << " "
375  << _points[qpoint]
376  << "\n Weight:\n "
377  << " w=" << _weights[qpoint] << "\n" << std::endl;
378 
379  summed_weights += _weights[qpoint];
380  }
381  os << "Summed Weights: " << summed_weights << std::endl;
382 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:335
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:341
libmesh_assert(j)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int n_points() const
Definition: quadrature.h:113
Point libMesh::QBase::qp ( const unsigned int  i) const
inherited
Returns
The $ i^{th} $ quadrature point in reference element space.

Definition at line 151 of file quadrature.h.

References libMesh::QBase::_points.

Referenced by libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), and libMesh::QBase::tensor_product_quad().

152  {
153  libmesh_assert_less (i, _points.size());
154  return _points[i];
155  }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:335
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 93 of file quadrature.C.

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

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

95 {
96  // Make sure we are in 1D
97  libmesh_assert_equal_to (_dim, 1);
98 
99  Real
100  h_new = new_range.second - new_range.first,
101  h_old = old_range.second - old_range.first;
102 
103  // Make sure that we have sane ranges
104  libmesh_assert_greater (h_new, 0.);
105  libmesh_assert_greater (h_old, 0.);
106 
107  // Make sure there are some points
108  libmesh_assert_greater (_points.size(), 0);
109 
110  // Compute the scale factor
111  Real scfact = h_new/h_old;
112 
113  // We're mapping from old_range -> new_range
114  for (std::size_t i=0; i<_points.size(); i++)
115  {
116  _points[i](0) = new_range.first +
117  (_points[i](0) - old_range.first) * scfact;
118 
119  // Scale the weights
120  _weights[i] *= scfact;
121  }
122 }
const unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:311
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:335
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:341
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual bool libMesh::QBase::shapes_need_reinit ( )
virtualinherited
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 217 of file quadrature.h.

217 { return false; }
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 154 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::QTrap::init_3D(), init_3D(), libMesh::QSimpson::init_3D(), and libMesh::QBase::init_3D().

155 {
156  const unsigned int np = q1D.n_points();
157 
158  _points.resize(np * np * np);
159 
160  _weights.resize(np * np * np);
161 
162  unsigned int q=0;
163 
164  for (unsigned int k=0; k<np; k++)
165  for (unsigned int j=0; j<np; j++)
166  for (unsigned int i=0; i<np; i++)
167  {
168  _points[q](0) = q1D.qp(i)(0);
169  _points[q](1) = q1D.qp(j)(0);
170  _points[q](2) = q1D.qp(k)(0);
171 
172  _weights[q] = q1D.w(i) * q1D.w(j) * q1D.w(k);
173 
174  q++;
175  }
176 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:335
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:341
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 181 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::QTrap::init_3D(), init_3D(), libMesh::QSimpson::init_3D(), and libMesh::QBase::init_3D().

182 {
183  const unsigned int n_points1D = q1D.n_points();
184  const unsigned int n_points2D = q2D.n_points();
185 
186  _points.resize (n_points1D * n_points2D);
187  _weights.resize (n_points1D * n_points2D);
188 
189  unsigned int q=0;
190 
191  for (unsigned int j=0; j<n_points1D; j++)
192  for (unsigned int i=0; i<n_points2D; i++)
193  {
194  _points[q](0) = q2D.qp(i)(0);
195  _points[q](1) = q2D.qp(i)(1);
196  _points[q](2) = q1D.qp(j)(0);
197 
198  _weights[q] = q2D.w(i) * q1D.w(j);
199 
200  q++;
201  }
202 
203 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:335
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:341
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 127 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::QTrap::init_2D(), init_2D(), libMesh::QSimpson::init_2D(), and libMesh::QBase::init_3D().

128 {
129 
130  const unsigned int np = q1D.n_points();
131 
132  _points.resize(np * np);
133 
134  _weights.resize(np * np);
135 
136  unsigned int q=0;
137 
138  for (unsigned int j=0; j<np; j++)
139  for (unsigned int i=0; i<np; i++)
140  {
141  _points[q](0) = q1D.qp(i)(0);
142  _points[q](1) = q1D.qp(j)(0);
143 
144  _weights[q] = q1D.w(i)*q1D.w(j);
145 
146  q++;
147  }
148 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:335
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:341
virtual QuadratureType libMesh::QGauss::type ( ) const
virtual
Real libMesh::QBase::w ( const unsigned int  i) const
inherited
Returns
The $ i^{th} $ quadrature weight.

Definition at line 160 of file quadrature.h.

References libMesh::QBase::_weights, libMesh::QBase::init(), libMesh::INVALID_ELEM, and libMesh::QBase::type().

Referenced by libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), and libMesh::QBase::tensor_product_quad().

161  {
162  libmesh_assert_less (i, _weights.size());
163  return _weights[i];
164  }
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:341

Member Data Documentation

ReferenceCounter::Counts libMesh::ReferenceCounter::_counts
staticprotectedinherited
const unsigned int libMesh::QBase::_dim
protectedinherited

The spatial dimension of the quadrature rule.

Definition at line 311 of file quadrature.h.

Referenced by libMesh::QBase::get_dim(), libMesh::QBase::init(), and libMesh::QBase::scale().

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

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.

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

const Order libMesh::QBase::_order
protectedinherited

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

Definition at line 317 of file quadrature.h.

Referenced by libMesh::QBase::get_order(), init_1D(), init_2D(), and init_3D().

unsigned int libMesh::QBase::_p_level
protectedinherited

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

Definition at line 329 of file quadrature.h.

Referenced by libMesh::QBase::get_order(), libMesh::QBase::get_p_level(), and libMesh::QBase::init().

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

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

Definition at line 323 of file quadrature.h.

Referenced by libMesh::QBase::get_elem_type(), libMesh::QBase::init(), libMesh::QTrap::type(), libMesh::QSimpson::type(), and type().

std::vector<Real> libMesh::QBase::_weights
protectedinherited
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 232 of file quadrature.h.

Referenced by init_3D().


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