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

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

#include <quadrature_gauss.h>

Inheritance diagram for libMesh::QGauss:
Inheritance graph
[legend]

Public Member Functions

 QGauss (unsigned int dim, Order order=INVALID_ORDER)
 
 QGauss (const QGauss &)=default
 
 QGauss (QGauss &&)=default
 
QGauss & operator= (const QGauss &)=default
 
QGauss & operator= (QGauss &&)=default
 
virtual ~QGauss ()=default
 
virtual QuadratureType type () const override
 
ElemType get_elem_type () const
 
unsigned int get_p_level () const
 
unsigned int n_points () const
 
unsigned int get_dim () const
 
const std::vector< Point > & get_points () const
 
std::vector< Point > & get_points ()
 
const std::vector< Real > & get_weights () const
 
std::vector< Real > & get_weights ()
 
Point qp (const unsigned int i) const
 
Real w (const unsigned int i) const
 
virtual void init (const ElemType type=INVALID_ELEM, unsigned int p_level=0)
 
virtual void init (const Elem &elem, const std::vector< Real > &vertex_distance_func, unsigned int p_level=0)
 
Order get_order () const
 
void print_info (std::ostream &os=libMesh::out) const
 
void scale (std::pair< Real, Real > old_range, std::pair< Real, Real > new_range)
 
virtual bool shapes_need_reinit ()
 

Static Public Member Functions

static std::unique_ptr< QBase > build (const std::string &name, const unsigned int dim, const Order order=INVALID_ORDER)
 
static std::unique_ptr< QBase > build (const QuadratureType qt, const unsigned int dim, const Order order=INVALID_ORDER)
 
static void print_info (std::ostream &out=libMesh::out)
 
static std::string get_info ()
 
static unsigned int n_objects ()
 
static void enable_print_counter_info ()
 
static void disable_print_counter_info ()
 

Public Attributes

bool allow_rules_with_negative_weights
 

Protected Types

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

Protected Member Functions

virtual void init_0D (const ElemType type=INVALID_ELEM, unsigned int p_level=0)
 
void tensor_product_quad (const QBase &q1D)
 
void tensor_product_hex (const QBase &q1D)
 
void tensor_product_prism (const QBase &q1D, const QBase &q2D)
 
void increment_constructor_count (const std::string &name)
 
void increment_destructor_count (const std::string &name)
 

Protected Attributes

unsigned int _dim
 
Order _order
 
ElemType _type
 
unsigned int _p_level
 
std::vector< Point > _points
 
std::vector< Real > _weights
 

Static Protected Attributes

static Counts _counts
 
static Threads::atomic< unsigned int > _n_objects
 
static Threads::spin_mutex _mutex
 
static bool _enable_print_counter = true
 

Private Member Functions

virtual void init_1D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0) override
 
virtual void init_2D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0) override
 
virtual void init_3D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0) override
 
void dunavant_rule (const Real rule_data[][4], const unsigned int n_pts)
 
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)
 

Detailed Description

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

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

Definition at line 39 of file quadrature_gauss.h.

Member Typedef Documentation

◆ Counts

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

Data structure to log the information. The log is identified by the class name.

Definition at line 117 of file reference_counter.h.

Constructor & Destructor Documentation

◆ QGauss() [1/3]

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

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

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

Definition at line 52 of file quadrature_gauss.h.

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

53  :
54  QBase(dim, order)
55  {
56  if (dim == 1)
57  init(EDGE2);
58  }
libMesh::QBase::init
virtual void init(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Definition: quadrature.C:28
libMesh::QBase::QBase
QBase(unsigned int dim, Order order=INVALID_ORDER)
Definition: quadrature.h:364

◆ QGauss() [2/3]

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

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

◆ QGauss() [3/3]

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

◆ ~QGauss()

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

Member Function Documentation

◆ build() [1/2]

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

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

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

Definition at line 40 of file quadrature_build.C.

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

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

43 {
44  return QBase::build (Utility::string_to_enum<QuadratureType> (type),
45  _dim,
46  _order);
47 }
libMesh::QBase::_dim
unsigned int _dim
Definition: quadrature.h:325
libMesh::QBase::type
virtual QuadratureType type() const =0
libMesh::QBase::build
static std::unique_ptr< QBase > build(const std::string &name, const unsigned int dim, const Order order=INVALID_ORDER)
Definition: quadrature_build.C:40

◆ build() [2/2]

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

Builds a specific quadrature rule based on the QuadratureType. This enables selection of the quadrature rule at run-time.

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

Definition at line 51 of file quadrature_build.C.

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

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

◆ disable_print_counter_info()

void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

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

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

◆ dunavant_rule()

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

The Dunavant rules are for triangles. This function takes permutation points and weights in a specific format as input and fills the _points and _weights vectors.

Definition at line 265 of file quadrature_gauss.C.

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

Referenced by init_2D().

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

◆ dunavant_rule2()

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

Definition at line 190 of file quadrature_gauss.C.

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

Referenced by init_2D().

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

◆ enable_print_counter_info()

void libMesh::ReferenceCounter::enable_print_counter_info ( )
staticinherited

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

Definition at line 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

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

◆ get_dim()

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

Definition at line 136 of file quadrature.h.

References libMesh::QBase::_dim.

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

136 { return _dim; }
libMesh::QBase::_dim
unsigned int _dim
Definition: quadrature.h:325

◆ get_elem_type()

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

Definition at line 117 of file quadrature.h.

References libMesh::QBase::_type.

117 { return _type; }
libMesh::QBase::_type
ElemType _type
Definition: quadrature.h:337

◆ get_info()

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

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

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

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

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

◆ get_order()

Order libMesh::QBase::get_order ( ) const
inlineinherited
Returns
The order of the quadrature rule.

Definition at line 203 of file quadrature.h.

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

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

203 { return static_cast<Order>(_order + _p_level); }
libMesh::QBase::_p_level
unsigned int _p_level
Definition: quadrature.h:343

◆ get_p_level()

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

Definition at line 122 of file quadrature.h.

References libMesh::QBase::_p_level.

122 { return _p_level; }
libMesh::QBase::_p_level
unsigned int _p_level
Definition: quadrature.h:343

◆ get_points() [1/2]

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

Definition at line 142 of file quadrature.h.

References libMesh::QBase::_points.

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

142 { return _points; }
libMesh::QBase::_points
std::vector< Point > _points
Definition: quadrature.h:349

◆ get_points() [2/2]

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

Definition at line 148 of file quadrature.h.

References libMesh::QBase::_points.

148 { return _points; }
libMesh::QBase::_points
std::vector< Point > _points
Definition: quadrature.h:349

◆ get_weights() [1/2]

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

Definition at line 154 of file quadrature.h.

References libMesh::QBase::_weights.

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

154 { return _weights; }
libMesh::QBase::_weights
std::vector< Real > _weights
Definition: quadrature.h:355

◆ get_weights() [2/2]

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

Definition at line 160 of file quadrature.h.

References libMesh::QBase::_weights.

160 { return _weights; }
libMesh::QBase::_weights
std::vector< Real > _weights
Definition: quadrature.h:355

◆ increment_constructor_count()

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

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

Definition at line 181 of file reference_counter.h.

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

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

182 {
183  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
184  std::pair<unsigned int, unsigned int> & p = _counts[name];
185 
186  p.first++;
187 }
libMesh::Quality::name
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
libMesh::Threads::spin_mtx
spin_mutex spin_mtx
Definition: threads.C:29

◆ increment_destructor_count()

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

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

Definition at line 194 of file reference_counter.h.

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

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

195 {
196  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
197  std::pair<unsigned int, unsigned int> & p = _counts[name];
198 
199  p.second++;
200 }
libMesh::Quality::name
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
libMesh::Threads::spin_mtx
spin_mutex spin_mtx
Definition: threads.C:29

◆ init() [1/2]

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

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

Definition at line 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::QBase::init(), libMesh::QClough::init_1D(), libMesh::QMonomial::init_1D(), libMesh::QClough::init_2D(), libMesh::QGaussLobatto::init_2D(), libMesh::QGrid::init_2D(), libMesh::QTrap::init_2D(), libMesh::QSimpson::init_2D(), init_2D(), libMesh::QMonomial::init_2D(), libMesh::QGaussLobatto::init_3D(), libMesh::QTrap::init_3D(), libMesh::QGrid::init_3D(), init_3D(), libMesh::QSimpson::init_3D(), libMesh::QMonomial::init_3D(), QGauss(), libMesh::QGaussLobatto::QGaussLobatto(), libMesh::QJacobi::QJacobi(), libMesh::QSimpson::QSimpson(), and libMesh::QTrap::QTrap().

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 }
libMesh::QBase::_type
ElemType _type
Definition: quadrature.h:337
libMesh::QBase::_dim
unsigned int _dim
Definition: quadrature.h:325
libMesh::QBase::_p_level
unsigned int _p_level
Definition: quadrature.h:343
libMesh::QBase::init_2D
virtual void init_2D(const ElemType, unsigned int=0)
Definition: quadrature.h:277
libMesh::QBase::init_1D
virtual void init_1D(const ElemType type=INVALID_ELEM, unsigned int p_level=0)=0
libMesh::QBase::init_0D
virtual void init_0D(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Definition: quadrature.C:82
libMesh::QBase::init_3D
virtual void init_3D(const ElemType, unsigned int=0)
Definition: quadrature.h:293

◆ init() [2/2]

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

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

Definition at line 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 }
libMesh::QBase::init
virtual void init(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Definition: quadrature.C:28

◆ init_0D()

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

Initializes the 0D quadrature rule by filling the points and weights vectors with the appropriate values. Generally this is just one point with weight 1.

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 }
libMesh::QBase::_points
std::vector< Point > _points
Definition: quadrature.h:349
libMesh::QBase::_weights
std::vector< Real > _weights
Definition: quadrature.h:355

◆ init_1D()

void libMesh::QGauss::init_1D ( const ElemType  type = INVALID_ELEM,
unsigned int  p_level = 0 
)
overrideprivatevirtual

Initializes the 1D quadrature rule by filling the points and weights vectors with the appropriate values. The order of the rule will be defined by the implementing class. It is assumed that derived quadrature rules will at least define the init_1D function, therefore it is pure virtual.

Implements libMesh::QBase.

Definition at line 30 of file quadrature_gauss_1D.C.

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

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) = Real(-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) = Real(-7.7459666924148337703585307995648e-01L);
70  _points[ 1](0) = 0.;
71  _points[ 2] = -_points[0];
72 
73  _weights[ 0] = Real(5.5555555555555555555555555555556e-01L);
74  _weights[ 1] = Real(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) = Real(-8.6113631159405257522394648889281e-01L);
86  _points[ 1](0) = Real(-3.3998104358485626480266575910324e-01L);
87  _points[ 2] = -_points[1];
88  _points[ 3] = -_points[0];
89 
90  _weights[ 0] = Real(3.4785484513745385737306394922200e-01L);
91  _weights[ 1] = Real(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) = Real(-9.0617984593866399279762687829939e-01L);
104  _points[ 1](0) = Real(-5.3846931010568309103631442070021e-01L);
105  _points[ 2](0) = 0.;
106  _points[ 3] = -_points[1];
107  _points[ 4] = -_points[0];
108 
109  _weights[ 0] = Real(2.3692688505618908751426404071992e-01L);
110  _weights[ 1] = Real(4.7862867049936646804129151483564e-01L);
111  _weights[ 2] = Real(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) = Real(-9.3246951420315202781230155449399e-01L);
124  _points[ 1](0) = Real(-6.6120938646626451366139959501991e-01L);
125  _points[ 2](0) = Real(-2.3861918608319690863050172168071e-01L);
126  _points[ 3] = -_points[2];
127  _points[ 4] = -_points[1];
128  _points[ 5] = -_points[0];
129 
130  _weights[ 0] = Real(1.7132449237917034504029614217273e-01L);
131  _weights[ 1] = Real(3.6076157304813860756983351383772e-01L);
132  _weights[ 2] = Real(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) = Real(-9.4910791234275852452618968404785e-01L);
146  _points[ 1](0) = Real(-7.4153118559939443986386477328079e-01L);
147  _points[ 2](0) = Real(-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] = Real(1.2948496616886969327061143267908e-01L);
154  _weights[ 1] = Real(2.7970539148927666790146777142378e-01L);
155  _weights[ 2] = Real(3.8183005050511894495036977548898e-01L);
156  _weights[ 3] = Real(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) = Real(-9.6028985649753623168356086856947e-01L);
170  _points[ 1](0) = Real(-7.9666647741362673959155393647583e-01L);
171  _points[ 2](0) = Real(-5.2553240991632898581773904918925e-01L);
172  _points[ 3](0) = Real(-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] = Real(1.0122853629037625915253135430996e-01L);
179  _weights[ 1] = Real(2.2238103445337447054435599442624e-01L);
180  _weights[ 2] = Real(3.1370664587788728733796220198660e-01L);
181  _weights[ 3] = Real(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) = Real(-9.6816023950762608983557620290367e-01L);
196  _points[ 1](0) = Real(-8.3603110732663579429942978806973e-01L);
197  _points[ 2](0) = Real(-6.1337143270059039730870203934147e-01L);
198  _points[ 3](0) = Real(-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] = Real(8.1274388361574411971892158110524e-02L);
206  _weights[ 1] = Real(1.8064816069485740405847203124291e-01L);
207  _weights[ 2] = Real(2.6061069640293546231874286941863e-01L);
208  _weights[ 3] = Real(3.1234707704000284006863040658444e-01L);
209  _weights[ 4] = Real(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) = Real(-9.7390652851717172007796401208445e-01L);
224  _points[ 1](0) = Real(-8.6506336668898451073209668842349e-01L);
225  _points[ 2](0) = Real(-6.7940956829902440623432736511487e-01L);
226  _points[ 3](0) = Real(-4.3339539412924719079926594316578e-01L);
227  _points[ 4](0) = Real(-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] = Real(6.6671344308688137593568809893332e-02L);
235  _weights[ 1] = Real(1.4945134915058059314577633965770e-01L);
236  _weights[ 2] = Real(2.1908636251598204399553493422816e-01L);
237  _weights[ 3] = Real(2.6926671930999635509122692156947e-01L);
238  _weights[ 4] = Real(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) = Real(-9.7822865814605699280393800112286e-01L);
255  _points[ 1](0) = Real(-8.8706259976809529907515776930393e-01L);
256  _points[ 2](0) = Real(-7.3015200557404932409341625203115e-01L);
257  _points[ 3](0) = Real(-5.1909612920681181592572566945861e-01L);
258  _points[ 4](0) = Real(-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] = Real(5.5668567116173666482753720442549e-02L);
267  _weights[ 1] = Real(1.2558036946490462463469429922394e-01L);
268  _weights[ 2] = Real(1.8629021092773425142609764143166e-01L);
269  _weights[ 3] = Real(2.3319376459199047991852370484318e-01L);
270  _weights[ 4] = Real(2.6280454451024666218068886989051e-01L);
271  _weights[ 5] = Real(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) = Real(-9.8156063424671925069054909014928e-01L);
288  _points[ 1](0) = Real(-9.0411725637047485667846586611910e-01L);
289  _points[ 2](0) = Real(-7.6990267419430468703689383321282e-01L);
290  _points[ 3](0) = Real(-5.8731795428661744729670241894053e-01L);
291  _points[ 4](0) = Real(-3.6783149899818019375269153664372e-01L);
292  _points[ 5](0) = Real(-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] = Real(4.7175336386511827194615961485017e-02L);
301  _weights[ 1] = Real(1.0693932599531843096025471819400e-01L);
302  _weights[ 2] = Real(1.6007832854334622633465252954336e-01L);
303  _weights[ 3] = Real(2.0316742672306592174906445580980e-01L);
304  _weights[ 4] = Real(2.3349253653835480876084989892483e-01L);
305  _weights[ 5] = Real(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) = Real(-9.8418305471858814947282944880711e-01L);
323  _points[ 1](0) = Real(-9.1759839922297796520654783650072e-01L);
324  _points[ 2](0) = Real(-8.0157809073330991279420648958286e-01L);
325  _points[ 3](0) = Real(-6.4234933944034022064398460699552e-01L);
326  _points[ 4](0) = Real(-4.4849275103644685287791285212764e-01L);
327  _points[ 5](0) = Real(-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] = Real(4.0484004765315879520021592200986e-02L);
337  _weights[ 1] = Real(9.2121499837728447914421775953797e-02L);
338  _weights[ 2] = Real(1.3887351021978723846360177686887e-01L);
339  _weights[ 3] = Real(1.7814598076194573828004669199610e-01L);
340  _weights[ 4] = Real(2.0781604753688850231252321930605e-01L);
341  _weights[ 5] = Real(2.2628318026289723841209018603978e-01L);
342  _weights[ 6] = Real(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) = Real(-9.8628380869681233884159726670405e-01L);
360  _points[ 1](0) = Real(-9.2843488366357351733639113937787e-01L);
361  _points[ 2](0) = Real(-8.2720131506976499318979474265039e-01L);
362  _points[ 3](0) = Real(-6.8729290481168547014801980301933e-01L);
363  _points[ 4](0) = Real(-5.1524863635815409196529071855119e-01L);
364  _points[ 5](0) = Real(-3.1911236892788976043567182416848e-01L);
365  _points[ 6](0) = Real(-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] = Real(3.5119460331751863031832876138192e-02L);
375  _weights[ 1] = Real(8.0158087159760209805633277062854e-02L);
376  _weights[ 2] = Real(1.2151857068790318468941480907248e-01L);
377  _weights[ 3] = Real(1.5720316715819353456960193862384e-01L);
378  _weights[ 4] = Real(1.8553839747793781374171659012516e-01L);
379  _weights[ 5] = Real(2.0519846372129560396592406566122e-01L);
380  _weights[ 6] = Real(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) = Real(-9.8799251802048542848956571858661e-01L);
399  _points[ 1](0) = Real(-9.3727339240070590430775894771021e-01L);
400  _points[ 2](0) = Real(-8.4820658341042721620064832077422e-01L);
401  _points[ 3](0) = Real(-7.2441773136017004741618605461394e-01L);
402  _points[ 4](0) = Real(-5.7097217260853884753722673725391e-01L);
403  _points[ 5](0) = Real(-3.9415134707756336989720737098105e-01L);
404  _points[ 6](0) = Real(-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] = Real(3.0753241996117268354628393577204e-02L);
415  _weights[ 1] = Real(7.0366047488108124709267416450667e-02L);
416  _weights[ 2] = Real(1.0715922046717193501186954668587e-01L);
417  _weights[ 3] = Real(1.3957067792615431444780479451103e-01L);
418  _weights[ 4] = Real(1.6626920581699393355320086048121e-01L);
419  _weights[ 5] = Real(1.8616100001556221102680056186642e-01L);
420  _weights[ 6] = Real(1.9843148532711157645611832644384e-01L);
421  _weights[ 7] = Real(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) = Real(-9.8940093499164993259615417345033e-01L);
440  _points[ 1](0) = Real(-9.4457502307323257607798841553461e-01L);
441  _points[ 2](0) = Real(-8.6563120238783174388046789771239e-01L);
442  _points[ 3](0) = Real(-7.5540440835500303389510119484744e-01L);
443  _points[ 4](0) = Real(-6.1787624440264374844667176404879e-01L);
444  _points[ 5](0) = Real(-4.5801677765722738634241944298358e-01L);
445  _points[ 6](0) = Real(-2.8160355077925891323046050146050e-01L);
446  _points[ 7](0) = Real(-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] = Real(2.7152459411754094851780572456018e-02L);
457  _weights[ 1] = Real(6.2253523938647892862843836994378e-02L);
458  _weights[ 2] = Real(9.5158511682492784809925107602246e-02L);
459  _weights[ 3] = Real(1.2462897125553387205247628219202e-01L);
460  _weights[ 4] = Real(1.4959598881657673208150173054748e-01L);
461  _weights[ 5] = Real(1.6915651939500253818931207903033e-01L);
462  _weights[ 6] = Real(1.8260341504492358886676366796922e-01L);
463  _weights[ 7] = Real(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) = Real(-9.9057547531441733567543401994067e-01L);
483  _points[ 1](0) = Real(-9.5067552176876776122271695789580e-01L);
484  _points[ 2](0) = Real(-8.8023915372698590212295569448816e-01L);
485  _points[ 3](0) = Real(-7.8151400389680140692523005552048e-01L);
486  _points[ 4](0) = Real(-6.5767115921669076585030221664300e-01L);
487  _points[ 5](0) = Real(-5.1269053708647696788624656862955e-01L);
488  _points[ 6](0) = Real(-3.5123176345387631529718551709535e-01L);
489  _points[ 7](0) = Real(-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] = Real(2.4148302868547931960110026287565e-02L);
501  _weights[ 1] = Real(5.5459529373987201129440165358245e-02L);
502  _weights[ 2] = Real(8.5036148317179180883535370191062e-02L);
503  _weights[ 3] = Real(1.1188384719340397109478838562636e-01L);
504  _weights[ 4] = Real(1.3513636846852547328631998170235e-01L);
505  _weights[ 5] = Real(1.5404576107681028808143159480196e-01L);
506  _weights[ 6] = Real(1.6800410215645004450997066378832e-01L);
507  _weights[ 7] = Real(1.7656270536699264632527099011320e-01L);
508  _weights[ 8] = Real(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) = Real(-9.9156516842093094673001600470615e-01L);
528  _points[ 1](0) = Real(-9.5582394957139775518119589292978e-01L);
529  _points[ 2](0) = Real(-8.9260246649755573920606059112715e-01L);
530  _points[ 3](0) = Real(-8.0370495897252311568241745501459e-01L);
531  _points[ 4](0) = Real(-6.9168704306035320787489108128885e-01L);
532  _points[ 5](0) = Real(-5.5977083107394753460787154852533e-01L);
533  _points[ 6](0) = Real(-4.1175116146284264603593179383305e-01L);
534  _points[ 7](0) = Real(-2.5188622569150550958897285487791e-01L);
535  _points[ 8](0) = Real(-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] = Real(2.1616013526483310313342710266452e-02L);
547  _weights[ 1] = Real(4.9714548894969796453334946202639e-02L);
548  _weights[ 2] = Real(7.6425730254889056529129677616637e-02L);
549  _weights[ 3] = Real(1.0094204410628716556281398492483e-01L);
550  _weights[ 4] = Real(1.2255520671147846018451912680020e-01L);
551  _weights[ 5] = Real(1.4064291467065065120473130375195e-01L);
552  _weights[ 6] = Real(1.5468467512626524492541800383637e-01L);
553  _weights[ 7] = Real(1.6427648374583272298605377646593e-01L);
554  _weights[ 8] = Real(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) = Real(-9.9240684384358440318901767025326e-01L);
575  _points[ 1](0) = Real(-9.6020815213483003085277884068765e-01L);
576  _points[ 2](0) = Real(-9.0315590361481790164266092853231e-01L);
577  _points[ 3](0) = Real(-8.2271465653714282497892248671271e-01L);
578  _points[ 4](0) = Real(-7.2096617733522937861709586082378e-01L);
579  _points[ 5](0) = Real(-6.0054530466168102346963816494624e-01L);
580  _points[ 6](0) = Real(-4.6457074137596094571726714810410e-01L);
581  _points[ 7](0) = Real(-3.1656409996362983199011732884984e-01L);
582  _points[ 8](0) = Real(-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] = Real(1.9461788229726477036312041464438e-02L);
595  _weights[ 1] = Real(4.4814226765699600332838157401994e-02L);
596  _weights[ 2] = Real(6.9044542737641226580708258006013e-02L);
597  _weights[ 3] = Real(9.1490021622449999464462094123840e-02L);
598  _weights[ 4] = Real(1.1156664554733399471602390168177e-01L);
599  _weights[ 5] = Real(1.2875396253933622767551578485688e-01L);
600  _weights[ 6] = Real(1.4260670217360661177574610944190e-01L);
601  _weights[ 7] = Real(1.5276604206585966677885540089766e-01L);
602  _weights[ 8] = Real(1.5896884339395434764995643946505e-01L);
603  _weights[ 9] = Real(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) = Real(-9.9312859918509492478612238847132e-01L);
624  _points[ 1](0) = Real(-9.6397192727791379126766613119728e-01L);
625  _points[ 2](0) = Real(-9.1223442825132590586775244120330e-01L);
626  _points[ 3](0) = Real(-8.3911697182221882339452906170152e-01L);
627  _points[ 4](0) = Real(-7.4633190646015079261430507035564e-01L);
628  _points[ 5](0) = Real(-6.3605368072651502545283669622629e-01L);
629  _points[ 6](0) = Real(-5.1086700195082709800436405095525e-01L);
630  _points[ 7](0) = Real(-3.7370608871541956067254817702493e-01L);
631  _points[ 8](0) = Real(-2.2778585114164507808049619536857e-01L);
632  _points[ 9](0) = Real(-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] = Real(1.7614007139152118311861962351853e-02L);
645  _weights[ 1] = Real(4.0601429800386941331039952274932e-02L);
646  _weights[ 2] = Real(6.2672048334109063569506535187042e-02L);
647  _weights[ 3] = Real(8.3276741576704748724758143222046e-02L);
648  _weights[ 4] = Real(1.0193011981724043503675013548035e-01L);
649  _weights[ 5] = Real(1.1819453196151841731237737771138e-01L);
650  _weights[ 6] = Real(1.3168863844917662689849449974816e-01L);
651  _weights[ 7] = Real(1.4209610931838205132929832506716e-01L);
652  _weights[ 8] = Real(1.4917298647260374678782873700197e-01L);
653  _weights[ 9] = Real(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) = Real(-9.9375217062038950026024203593794e-01L);
675  _points[ 1](0) = Real(-9.6722683856630629431662221490770e-01L);
676  _points[ 2](0) = Real(-9.2009933415040082879018713371497e-01L);
677  _points[ 3](0) = Real(-8.5336336458331728364725063858757e-01L);
678  _points[ 4](0) = Real(-7.6843996347567790861587785130623e-01L);
679  _points[ 5](0) = Real(-6.6713880419741231930596666999034e-01L);
680  _points[ 6](0) = Real(-5.5161883588721980705901879672431e-01L);
681  _points[ 7](0) = Real(-4.2434212020743878357366888854379e-01L);
682  _points[ 8](0) = Real(-2.8802131680240109660079251606460e-01L);
683  _points[ 9](0) = Real(-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] = Real(1.6017228257774333324224616858471e-02L);
697  _weights[ 1] = Real(3.6953789770852493799950668299330e-02L);
698  _weights[ 2] = Real(5.7134425426857208283635826472448e-02L);
699  _weights[ 3] = Real(7.6100113628379302017051653300183e-02L);
700  _weights[ 4] = Real(9.3444423456033861553289741113932e-02L);
701  _weights[ 5] = Real(1.0879729916714837766347457807011e-01L);
702  _weights[ 6] = Real(1.2183141605372853419536717712572e-01L);
703  _weights[ 7] = Real(1.3226893863333746178105257449678e-01L);
704  _weights[ 8] = Real(1.3988739479107315472213342386758e-01L);
705  _weights[ 9] = Real(1.4452440398997005906382716655375e-01L);
706  _weights[10] = Real(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) = Real(-9.9429458548239929207303142116130e-01L);
728  _points[ 1](0) = Real(-9.7006049783542872712395098676527e-01L);
729  _points[ 2](0) = Real(-9.2695677218717400052069293925905e-01L);
730  _points[ 3](0) = Real(-8.6581257772030013653642563701938e-01L);
731  _points[ 4](0) = Real(-7.8781680597920816200427795540835e-01L);
732  _points[ 5](0) = Real(-6.9448726318668278005068983576226e-01L);
733  _points[ 6](0) = Real(-5.8764040350691159295887692763865e-01L);
734  _points[ 7](0) = Real(-4.6935583798675702640633071096641e-01L);
735  _points[ 8](0) = Real(-3.4193582089208422515814742042738e-01L);
736  _points[ 9](0) = Real(-2.0786042668822128547884653391955e-01L);
737  _points[10](0) = Real(-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] = Real(1.4627995298272200684991098047185e-02L);
751  _weights[ 1] = Real(3.3774901584814154793302246865913e-02L);
752  _weights[ 2] = Real(5.2293335152683285940312051273211e-02L);
753  _weights[ 3] = Real(6.9796468424520488094961418930218e-02L);
754  _weights[ 4] = Real(8.5941606217067727414443681372703e-02L);
755  _weights[ 5] = Real(1.0041414444288096493207883783054e-01L);
756  _weights[ 6] = Real(1.1293229608053921839340060742175e-01L);
757  _weights[ 7] = Real(1.2325237681051242428556098615481e-01L);
758  _weights[ 8] = Real(1.3117350478706237073296499253031e-01L);
759  _weights[ 9] = Real(1.3654149834601517135257383123152e-01L);
760  _weights[10] = Real(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 }
libMesh::QBase::_points
std::vector< Point > _points
Definition: quadrature.h:349
libMesh::QBase::_weights
std::vector< Real > _weights
Definition: quadrature.h:355
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131

◆ init_2D()

void libMesh::QGauss::init_2D ( const ElemType  = INVALID_ELEM,
unsigned int  = 0 
)
overrideprivatevirtual

Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate values. The order of the rule will be defined by the implementing class. Should not be pure virtual since a derived quadrature rule may only be defined in 1D. If not 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::get_points(), libMesh::QBase::get_weights(), libMesh::QBase::init(), libMesh::NINETEENTH, libMesh::NINTH, libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::QUADSHELL4, libMesh::QUADSHELL8, libMesh::Real, libMesh::SECOND, libMesh::SEVENTEENTH, libMesh::SEVENTH, libMesh::SIXTEENTH, libMesh::SIXTH, 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.

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) = Real(1.5505102572168219018027159252941e-01L);
139  _points[0](1) = Real(1.7855872826361642311703513337422e-01L);
140  _points[1](0) = Real(6.4494897427831780981972840747059e-01L);
141  _points[1](1) = Real(7.5031110222608118177475598324603e-02L);
142  _points[2](0) = Real(1.5505102572168219018027159252941e-01L);
143  _points[2](1) = Real(6.6639024601470138670269327409637e-01L);
144  _points[3](0) = Real(6.4494897427831780981972840747059e-01L);
145  _points[3](1) = Real(2.8001991549907407200279599420481e-01L);
146 
147  _weights[0] = Real(1.5902069087198858469718450103758e-01L);
148  _weights[1] = Real(9.0979309128011415302815498962418e-02L);
149  _weights[2] = Real(1.5902069087198858469718450103758e-01L);
150  _weights[3] = Real(9.0979309128011415302815498962418e-02L);
151 
152  return;
153 
154 
155  // The following third-order rule is quite commonly cited
156  // in the literature and most likely works fine. However,
157  // we generally prefer a rule with all positive weights
158  // and an equal number of points, when available.
159  //
160  // (allow_rules_with_negative_weights)
161  // {
162  // // Exact for cubics
163  // _points.resize(4);
164  // _weights.resize(4);
165  //
166  // _points[0](0) = .33333333333333333333333333333333;
167  // _points[0](1) = .33333333333333333333333333333333;
168  //
169  // _points[1](0) = .2;
170  // _points[1](1) = .6;
171  //
172  // _points[2](0) = .2;
173  // _points[2](1) = .2;
174  //
175  // _points[3](0) = .6;
176  // _points[3](1) = .2;
177  //
178  //
179  // _weights[0] = -27./96.;
180  // _weights[1] = 25./96.;
181  // _weights[2] = 25./96.;
182  // _weights[3] = 25./96.;
183  //
184  // return;
185  // } // end if (allow_rules_with_negative_weights)
186  // Note: if !allow_rules_with_negative_weights, fall through to next case.
187  }
188 
189 
190 
191  // A degree 4 rule with six points. This rule can be found in many places
192  // including:
193  //
194  // J.N. Lyness and D. Jespersen, Moderate degree symmetric
195  // quadrature rules for the triangle, J. Inst. Math. Appl. 15 (1975),
196  // 19--32.
197  //
198  // We used the code in:
199  // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
200  // on triangles and tetrahedra" Journal of Computational Mathematics,
201  // v. 27, no. 1, 2009, pp. 89-96.
202  // to generate additional precision.
203  case FOURTH:
204  {
205  const unsigned int n_wts = 2;
206  const Real wts[n_wts] =
207  {
208  Real(1.1169079483900573284750350421656140e-01L),
209  Real(5.4975871827660933819163162450105264e-02L)
210  };
211 
212  const Real a[n_wts] =
213  {
214  Real(4.4594849091596488631832925388305199e-01L),
215  Real(9.1576213509770743459571463402201508e-02L)
216  };
217 
218  const Real b[n_wts] = {0., 0.}; // not used
219  const unsigned int permutation_ids[n_wts] = {3, 3};
220 
221  dunavant_rule2(wts, a, b, permutation_ids, n_wts); // 6 total points
222 
223  return;
224  }
225 
226 
227 
228  // Exact for quintics
229  // Can be found in "Quadrature on Simplices of Arbitrary
230  // Dimension" by Walkington.
231  case FIFTH:
232  {
233  const unsigned int n_wts = 3;
234  const Real wts[n_wts] =
235  {
236  Real(9)/80,
237  Real(31)/480 + std::sqrt(Real(15.0L))/2400,
238  Real(31)/480 - std::sqrt(Real(15.0L))/2400
239  };
240 
241  const Real a[n_wts] =
242  {
243  0., // 'a' parameter not used for origin permutation
244  Real(2)/7 + std::sqrt(Real(15.0L))/21,
245  Real(2)/7 - std::sqrt(Real(15.0L))/21
246  };
247 
248  const Real b[n_wts] = {0., 0., 0.}; // not used
249  const unsigned int permutation_ids[n_wts] = {1, 3, 3};
250 
251  dunavant_rule2(wts, a, b, permutation_ids, n_wts); // 7 total points
252 
253  return;
254  }
255 
256 
257 
258  // A degree 6 rule with 12 points. This rule can be found in many places
259  // including:
260  //
261  // J.N. Lyness and D. Jespersen, Moderate degree symmetric
262  // quadrature rules for the triangle, J. Inst. Math. Appl. 15 (1975),
263  // 19--32.
264  //
265  // We used the code in:
266  // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
267  // on triangles and tetrahedra" Journal of Computational Mathematics,
268  // v. 27, no. 1, 2009, pp. 89-96.
269  // to generate additional precision.
270  //
271  // Note that the following 7th-order Ro3-invariant rule also has only 12 points,
272  // which technically makes it the superior rule. This one is here for completeness.
273  case SIXTH:
274  {
275  const unsigned int n_wts = 3;
276  const Real wts[n_wts] =
277  {
278  Real(5.8393137863189683012644805692789721e-02L),
279  Real(2.5422453185103408460468404553434492e-02L),
280  Real(4.1425537809186787596776728210221227e-02L)
281  };
282 
283  const Real a[n_wts] =
284  {
285  Real(2.4928674517091042129163855310701908e-01L),
286  Real(6.3089014491502228340331602870819157e-02L),
287  Real(3.1035245103378440541660773395655215e-01L)
288  };
289 
290  const Real b[n_wts] =
291  {
292  0.,
293  0.,
294  Real(6.3650249912139864723014259441204970e-01L)
295  };
296 
297  const unsigned int permutation_ids[n_wts] = {3, 3, 6}; // 12 total points
298 
299  dunavant_rule2(wts, a, b, permutation_ids, n_wts);
300 
301  return;
302  }
303 
304 
305  // A degree 7 rule with 12 points. This rule can be found in:
306  //
307  // K. Gatermann, The construction of symmetric cubature
308  // formulas for the square and the triangle, Computing 40
309  // (1988), 229--240.
310  //
311  // This rule, which is provably minimal in the number of
312  // integration points, is said to be 'Ro3 invariant' which
313  // means that a given set of barycentric coordinates
314  // (z1,z2,z3) implies the quadrature points (z1,z2),
315  // (z3,z1), (z2,z3) which are formed by taking the first
316  // two entries in cyclic permutations of the barycentric
317  // point. Barycentric coordinates are related in the
318  // sense that: z3 = 1 - z1 - z2.
319  //
320  // The 12-point sixth-order rule for triangles given in
321  // Flaherty's (http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps)
322  // lecture notes has been removed in favor of this rule
323  // which is higher-order (for the same number of
324  // quadrature points) and has a few more digits of
325  // precision in the points and weights. Some 10-point
326  // degree 6 rules exist for the triangle but they have
327  // quadrature points outside the region of integration.
328  case SEVENTH:
329  {
330  _points.resize (12);
331  _weights.resize(12);
332 
333  const unsigned int nrows=4;
334 
335  // In each of the rows below, the first two entries are (z1, z2) which imply
336  // z3. The third entry is the weight for each of the points in the cyclic permutation.
337  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  // Real(2.6538900895116205835977487499847719e-02L),
369  // Real(3.5426541846066783659206291623201826e-02L),
370  // Real(3.4637341039708446756138297960207647e-02L)
371  // };
372  //
373  // const Real a[n_wts] =
374  // {
375  // Real(6.4930513159164863078379776030396538e-02L),
376  // Real(2.8457558424917033519741605734978046e-01L),
377  // Real(3.1355918438493150795585190219862865e-01L)
378  // };
379  //
380  // const Real b[n_wts] =
381  // {
382  // 0.,
383  // Real(1.9838447668150671917987659863332941e-01L),
384  // Real(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  Real(7.2157803838893584125545555244532310e-02L),
413  Real(4.7545817133642312396948052194292159e-02L),
414  Real(5.1608685267359125140895775146064515e-02L),
415  Real(1.6229248811599040155462964170890299e-02L),
416  Real(1.3615157087217497132422345036954462e-02L)
417  };
418 
419  const Real a[n_wts] =
420  {
421  0.0, // 'a' parameter not used for origin permutation
422  Real(4.5929258829272315602881551449416932e-01L),
423  Real(1.7056930775176020662229350149146450e-01L),
424  Real(5.0547228317030975458423550596598947e-02L),
425  Real(2.6311282963463811342178578628464359e-01L),
426  };
427 
428  const Real b[n_wts] =
429  {
430  0.,
431  0.,
432  0.,
433  0.,
434  Real(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  Real(4.8567898141399416909620991253644315e-02L),
461  Real(1.5667350113569535268427415643604658e-02L),
462  Real(1.2788837829349015630839399279499912e-02L),
463  Real(3.8913770502387139658369678149701978e-02L),
464  Real(3.9823869463605126516445887132022637e-02L),
465  Real(2.1641769688644688644688644688644689e-02L)
466  };
467 
468  const Real a[n_wts] =
469  {
470  0.0, // 'a' parameter not used for origin permutation
471  Real(4.8968251919873762778370692483619280e-01L),
472  Real(4.4729513394452709865106589966276365e-02L),
473  Real(4.3708959149293663726993036443535497e-01L),
474  Real(1.8820353561903273024096128046733557e-01L),
475  Real(2.2196298916076569567510252769319107e-01L)
476  };
477 
478  const Real b[n_wts] =
479  {
480  0.,
481  0.,
482  0.,
483  0.,
484  0.,
485  Real(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  Real(4.5408995191376790047643297550014267e-02L),
511  Real(1.8362978878233352358503035945683300e-02L),
512  Real(2.2660529717763967391302822369298659e-02L),
513  Real(3.6378958422710054302157588309680344e-02L),
514  Real(1.4163621265528742418368530791049552e-02L),
515  Real(4.7108334818664117299637354834434138e-03L)
516  };
517 
518  const Real a[n_wts] =
519  {
520  0.0, // 'a' parameter not used for origin permutation
521  Real(4.8557763338365737736750753220812615e-01L),
522  Real(1.0948157548503705479545863134052284e-01L),
523  Real(3.0793983876412095016515502293063162e-01L),
524  Real(2.4667256063990269391727646541117681e-01L),
525  Real(6.6803251012200265773540212762024737e-02L)
526  };
527 
528  const Real b[n_wts] =
529  {
530  0.,
531  0.,
532  0.,
533  Real(5.5035294182099909507816172659300821e-01L),
534  Real(7.2832390459741092000873505358107866e-01L),
535  Real(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  Real(3.6089021198604635216985338480426484e-02L),
567  Real(2.1607717807680420303346736867931050e-02L),
568  Real(3.1144524293927978774861144478241807e-03L),
569  Real(2.9086855161081509446654185084988077e-02L),
570  Real(8.4879241614917017182977532679947624e-03L),
571  Real(1.3795732078224796530729242858347546e-02L)
572  };
573 
574  const Real a[n_wts] =
575  {
576  Real(3.9355079629947969884346551941969960e-01L),
577  Real(4.7979065808897448654107733982929214e-01L),
578  Real(5.1003445645828061436081405648347852e-03L),
579  Real(2.6597620190330158952732822450744488e-01L),
580  Real(2.8536418538696461608233522814483715e-01L),
581  Real(1.3723536747817085036455583801851025e-01L)
582  };
583 
584  const Real b[n_wts] =
585  {
586  0.,
587  0.,
588  Real(5.6817155788572446538150614865768991e-02L),
589  Real(1.2539956353662088473247489775203396e-01L),
590  Real(1.2409970153698532116262152247041742e-02L),
591  Real(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  Real(3.0831305257795086169332418926151771e-03L),
620  Real(3.1429112108942550177135256546441273e-02L),
621  Real(1.7398056465354471494664198647499687e-02L),
622  Real(2.1846272269019201067728631278737487e-02L),
623  Real(1.2865533220227667708895461535782215e-02L),
624  Real(1.1178386601151722855919538351159995e-02L),
625  Real(8.6581155543294461858210504055170332e-03L),
626  Real(2.0185778883190464758914349626118386e-02L)
627  };
628 
629  const Real a[n_wts] =
630  {
631  Real(2.1317350453210370246856975515728246e-02L),
632  Real(2.7121038501211592234595134039689474e-01L),
633  Real(1.2757614554158592467389632515428357e-01L),
634  Real(4.3972439229446027297973662348436108e-01L),
635  Real(4.8821738977380488256466206525881104e-01L),
636  Real(2.8132558098993954824813069297455275e-01L),
637  Real(1.1625191590759714124135414784260182e-01L),
638  Real(2.7571326968551419397479634607976398e-01L)
639  };
640 
641  const Real b[n_wts] =
642  {
643  0.,
644  0.,
645  0.,
646  0.,
647  0.,
648  Real(6.9583608678780342214163552323607254e-01L),
649  Real(8.5801403354407263059053661662617818e-01L),
650  Real(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  Real(3.3980018293415822140887212340442440e-02L),
677  Real(2.7800983765226664353628733005230734e-02L),
678  Real(2.9139242559599990702383541756669905e-02L),
679  Real(3.0261685517695859208964000161454122e-03L),
680  Real(1.1997200964447365386855399725479827e-02L),
681  Real(1.7320638070424185232993414255459110e-02L),
682  Real(7.4827005525828336316229285664517190e-03L),
683  Real(1.2089519905796909568722872786530380e-02L),
684  Real(4.7953405017716313612975450830554457e-03L)
685  };
686 
687  const Real a[n_wts] =
688  {
689  0., // 'a' parameter not used for origin permutation
690  Real(4.2694141425980040602081253503137421e-01L),
691  Real(2.2137228629183290065481255470507908e-01L),
692  Real(2.1509681108843183869291313534052083e-02L),
693  Real(4.8907694645253934990068971909020439e-01L),
694  Real(3.0844176089211777465847185254124531e-01L),
695  Real(1.1092204280346339541286954522167452e-01L),
696  Real(1.6359740106785048023388790171095725e-01L),
697  Real(2.7251581777342966618005046435408685e-01L)
698  };
699 
700  const Real b[n_wts] =
701  {
702  0.,
703  0.,
704  0.,
705  0.,
706  0.,
707  Real(6.2354599555367557081585435318623659e-01L),
708  Real(8.6470777029544277530254595089569318e-01L),
709  Real(7.4850711589995219517301859578870965e-01L),
710  Real(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  Real(1.0941790684714445320422472981662986e-02L),
737  Real(1.6394176772062675320655489369312672e-02L),
738  Real(2.5887052253645793157392455083198201e-02L),
739  Real(2.1081294368496508769115218662093065e-02L),
740  Real(7.2168498348883338008549607403266583e-03L),
741  Real(2.4617018012000408409130117545210774e-03L),
742  Real(1.2332876606281836981437622591818114e-02L),
743  Real(1.9285755393530341614244513905205430e-02L),
744  Real(7.2181540567669202480443459995079017e-03L),
745  Real(2.5051144192503358849300465412445582e-03L)
746  };
747 
748  const Real a[n_wts] =
749  {
750  Real(4.8896391036217863867737602045239024e-01L),
751  Real(4.1764471934045392250944082218564344e-01L),
752  Real(2.7347752830883865975494428326269856e-01L),
753  Real(1.7720553241254343695661069046505908e-01L),
754  Real(6.1799883090872601267478828436935788e-02L),
755  Real(1.9390961248701048178250095054529511e-02L),
756  Real(1.7226668782135557837528960161365733e-01L),
757  Real(3.3686145979634500174405519708892539e-01L),
758  Real(2.9837288213625775297083151805961273e-01L),
759  Real(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  Real(7.7060855477499648258903327416742796e-01L),
771  Real(5.7022229084668317349769621336235426e-01L),
772  Real(6.8698016780808783735862715402031306e-01L),
773  Real(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  Real(2.4777380743035579804788826970198951e-02L),
808  Real(9.2433943023307730591540642828347660e-03L),
809  Real(2.2485768962175402793245929133296627e-03L),
810  Real(6.7052581900064143760518398833360903e-03L),
811  Real(1.9011381726930579256700190357527956e-02L),
812  Real(1.4605445387471889398286155981802858e-02L),
813  Real(1.5087322572773133722829435011138258e-02L),
814  Real(1.5630213780078803020711746273129099e-02L),
815  Real(6.1808086085778203192616856133701233e-03L),
816  Real(3.2209366452594664857296985751120513e-03L),
817  Real(5.8747373242569702667677969985668817e-03L)
818  };
819 
820  const Real a[n_wts] =
821  {
822  0.0, // 'a' parameter not used for origin
823  Real(7.9031013655541635005816956762252155e-02L),
824  Real(1.8789501810770077611247984432284226e-02L),
825  Real(4.9250168823249670532514526605352905e-01L),
826  Real(4.0886316907744105975059040108092775e-01L),
827  Real(5.3877851064220142445952549348423733e-01L),
828  Real(2.0250549804829997692885033941362673e-01L),
829  Real(5.5349674918711643207148086558288110e-01L),
830  Real(7.8345022567320812359258882143250181e-01L),
831  Real(8.9514624528794883409864566727625002e-01L),
832  Real(3.2515745241110782862789881780746490e-01L)
833  };
834 
835  const Real b[n_wts] =
836  {
837  0.,
838  0.,
839  0.,
840  0.,
841  0.,
842  Real(1.9412620368774630292701241080996842e-01L),
843  Real(9.8765911355712115933807754318089099e-02L),
844  Real(7.7663767064308164090246588765178087e-02L),
845  Real(2.1594628433980258573654682690950798e-02L),
846  Real(1.2563596287784997705599005477153617e-02L),
847  Real(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  Real(2.2668082505910087151996321171534230e-02L),
882  Real(8.4043060714818596159798961899306135e-03L),
883  Real(1.0850949634049747713966288634484161e-03L),
884  Real(7.2252773375423638869298219383808751e-03L),
885  Real(1.2997715227338366024036316182572871e-02L),
886  Real(2.0054466616677715883228810959112227e-02L),
887  Real(9.7299841600417010281624372720122710e-03L),
888  Real(1.1651974438298104227427176444311766e-02L),
889  Real(9.1291185550484450744725847363097389e-03L),
890  Real(3.5568614040947150231712567900113671e-03L),
891  Real(5.8355861686234326181790822005304303e-03L),
892  Real(4.7411314396804228041879331486234396e-03L)
893  };
894 
895  const Real a[n_wts] =
896  {
897  0.0, // 'a' parameter not used for centroid weight
898  Real(8.5402539407933203673769900926355911e-02L),
899  Real(1.2425572001444092841183633409631260e-02L),
900  Real(4.9174838341891594024701017768490960e-01L),
901  Real(4.5669426695387464162068900231444462e-01L),
902  Real(4.8506759880447437974189793537259677e-01L),
903  Real(2.0622099278664205707909858461264083e-01L),
904  Real(3.2374950270039093446805340265853956e-01L),
905  Real(7.3834330556606586255186213302750029e-01L),
906  Real(9.1210673061680792565673823935174611e-01L),
907  Real(6.6129919222598721544966837350891531e-01L),
908  Real(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  Real(3.2315912848634384647700266402091638e-01L),
919  Real(1.5341553679414688425981898952416987e-01L),
920  Real(7.4295478991330687632977899141707872e-02L),
921  Real(7.1278762832147862035977841733532020e-02L),
922  Real(1.6623223223705792825395256602140459e-02L),
923  Real(1.4160772533794791868984026749196156e-02L),
924  Real(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  Real(1.7464603792572004485690588092246146e-02L),
963  Real(5.9429003555801725246549713984660076e-03L),
964  Real(1.2490753345169579649319736639588729e-02L),
965  Real(1.5386987188875607593083456905596468e-02L),
966  Real(1.1185807311917706362674684312990270e-02L),
967  Real(1.0301845740670206831327304917180007e-02L),
968  Real(1.1767783072977049696840016810370464e-02L),
969  Real(3.8045312849431209558329128678945240e-03L),
970  Real(4.5139302178876351271037137230354382e-03L),
971  Real(2.2178812517580586419412547665472893e-03L),
972  Real(5.2216271537483672304731416553063103e-03L),
973  Real(9.8381136389470256422419930926212114e-04L)
974  };
975 
976  const Real a[n_wts] =
977  {
978  Real(2.8796825754667362165337965123570514e-01L),
979  Real(4.9216175986208465345536805750663939e-01L),
980  Real(4.6252866763171173685916780827044612e-01L),
981  Real(1.6730292951631792248498303276090273e-01L),
982  Real(1.5816335500814652972296428532213019e-01L),
983  Real(1.6352252138387564873002458959679529e-01L),
984  Real(6.2447680488959768233910286168417367e-01L),
985  Real(8.7317249935244454285263604347964179e-01L),
986  Real(3.4428164322282694677972239461699271e-01L),
987  Real(9.1584484467813674010523309855340209e-02L),
988  Real(2.0172088013378989086826623852040632e-01L),
989  Real(9.6538762758254643474731509845084691e-01L)
990  };
991 
992  const Real b[n_wts] =
993  {
994  0.0,
995  0.0,
996  0.0,
997  Real(3.4429160695501713926320695771253348e-01L),
998  Real(2.2541623431550639817203145525444726e-01L),
999  Real(8.0670083153531811694942222940484991e-02L),
1000  Real(6.5967451375050925655738829747288190e-02L),
1001  Real(4.5677879890996762665044366994439565e-02L),
1002  Real(1.1528411723154215812386518751976084e-02L),
1003  Real(9.3057714323900610398389176844165892e-03L),
1004  Real(1.5916814107619812717966560404970160e-02L),
1005  Real(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 }
libMesh::QBase::_points
std::vector< Point > _points
Definition: quadrature.h:349
libMesh::QBase::_weights
std::vector< Real > _weights
Definition: quadrature.h:355
libMesh::QGauss::dunavant_rule2
void dunavant_rule2(const Real *wts, const Real *a, const Real *b, const unsigned int *permutation_ids, const unsigned int n_wts)
Definition: quadrature_gauss.C:190
libMesh::QGauss::dunavant_rule
void dunavant_rule(const Real rule_data[][4], const unsigned int n_pts)
Definition: quadrature_gauss.C:265
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::QGauss::QGauss
QGauss(unsigned int dim, Order order=INVALID_ORDER)
Definition: quadrature_gauss.h:52
libMesh::QBase::tensor_product_quad
void tensor_product_quad(const QBase &q1D)
Definition: quadrature.C:127

◆ init_3D()

void libMesh::QGauss::init_3D ( const ElemType  = INVALID_ELEM,
unsigned int  = 0 
)
overrideprivatevirtual

Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate values. The order of the rule will be defined by the implementing class. Should not be pure virtual since a derived quadrature rule may only be defined in 1D. If not 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::QBase::get_points(), libMesh::QBase::get_weights(), libMesh::HEX20, libMesh::HEX27, libMesh::HEX8, libMesh::QBase::init(), keast_rule(), libMesh::PRISM15, libMesh::PRISM18, libMesh::PRISM6, libMesh::PYRAMID13, libMesh::PYRAMID14, libMesh::PYRAMID5, libMesh::Real, libMesh::SECOND, libMesh::SEVENTH, libMesh::SIXTH, libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), libMesh::TET10, libMesh::TET4, libMesh::THIRD, and libMesh::TRI3.

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  {
131  if (allow_rules_with_negative_weights)
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  {
199  if (allow_rules_with_negative_weights)
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  {
431  if (allow_rules_with_negative_weights)
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 }
libMesh::QBase::allow_rules_with_negative_weights
bool allow_rules_with_negative_weights
Definition: quadrature.h:246
libMesh::QBase::_points
std::vector< Point > _points
Definition: quadrature.h:349
libMesh::QBase::_weights
std::vector< Real > _weights
Definition: quadrature.h:355
libMesh::QBase::tensor_product_prism
void tensor_product_prism(const QBase &q1D, const QBase &q2D)
Definition: quadrature.C:181
libMesh::QBase::tensor_product_hex
void tensor_product_hex(const QBase &q1D)
Definition: quadrature.C:154
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::QGauss::keast_rule
void keast_rule(const Real rule_data[][4], const unsigned int n_pts)
Definition: quadrature_gauss.C:38
libMesh::QGauss::QGauss
QGauss(unsigned int dim, Order order=INVALID_ORDER)
Definition: quadrature_gauss.h:52

◆ keast_rule()

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

The Keast rules are for tets. This function takes permutation points and weights in a specific format as input and fills the _points and _weights vectors.

Definition at line 38 of file quadrature_gauss.C.

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

Referenced by init_3D().

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

◆ n_objects()

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

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

Definition at line 83 of file reference_counter.h.

References libMesh::ReferenceCounter::_n_objects.

84  { return _n_objects; }
libMesh::ReferenceCounter::_n_objects
static Threads::atomic< unsigned int > _n_objects
Definition: reference_counter.h:130

◆ n_points()

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

Definition at line 127 of file quadrature.h.

References libMesh::QBase::_points.

Referenced by libMesh::ExactSolution::_compute_error(), libMesh::System::calculate_norm(), libMesh::FirstOrderUnsteadySolver::compute_second_order_eqns(), libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), libMesh::QBase::print_info(), libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), and libMesh::QBase::tensor_product_quad().

128  {
129  libmesh_assert (!_points.empty());
130  return cast_int<unsigned int>(_points.size());
131  }
libMesh::QBase::_points
std::vector< Point > _points
Definition: quadrature.h:349

◆ operator=() [1/2]

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

◆ operator=() [2/2]

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

◆ print_info() [1/2]

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

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

Definition at line 87 of file reference_counter.C.

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

88 {
89  if (_enable_print_counter)
90  out_stream << ReferenceCounter::get_info();
91 }
libMesh::ReferenceCounter::get_info
static std::string get_info()
Definition: reference_counter.C:47

◆ print_info() [2/2]

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

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

Definition at line 378 of file quadrature.h.

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

Referenced by libMesh::operator<<().

379 {
380  libmesh_assert(!_points.empty());
381  libmesh_assert(!_weights.empty());
382 
383  Real summed_weights=0;
384  os << "N_Q_Points=" << this->n_points() << std::endl << std::endl;
385  for (unsigned int qpoint=0; qpoint<this->n_points(); qpoint++)
386  {
387  os << " Point " << qpoint << ":\n"
388  << " "
389  << _points[qpoint]
390  << "\n Weight:\n "
391  << " w=" << _weights[qpoint] << "\n" << std::endl;
392 
393  summed_weights += _weights[qpoint];
394  }
395  os << "Summed Weights: " << summed_weights << std::endl;
396 }
libMesh::QBase::_points
std::vector< Point > _points
Definition: quadrature.h:349
libMesh::QBase::_weights
std::vector< Real > _weights
Definition: quadrature.h:355
libMesh::QBase::n_points
unsigned int n_points() const
Definition: quadrature.h:127
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131

◆ qp()

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

Definition at line 165 of file quadrature.h.

References libMesh::QBase::_points.

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

166  {
167  libmesh_assert_less (i, _points.size());
168  return _points[i];
169  }
libMesh::QBase::_points
std::vector< Point > _points
Definition: quadrature.h:349

◆ scale()

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

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

Definition at line 93 of file quadrature.C.

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

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

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 }
libMesh::QBase::_dim
unsigned int _dim
Definition: quadrature.h:325
libMesh::QBase::_points
std::vector< Point > _points
Definition: quadrature.h:349
libMesh::QBase::_weights
std::vector< Real > _weights
Definition: quadrature.h:355
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131

◆ shapes_need_reinit()

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

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

Definition at line 231 of file quadrature.h.

231 { return false; }

◆ tensor_product_hex()

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

Computes the tensor product quadrature rule [q1D x q1D x q1D] from the 1D rule q1D. Used in the init_3D routines for hexahedral element types.

Definition at line 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::QGaussLobatto::init_3D(), libMesh::QTrap::init_3D(), libMesh::QGrid::init_3D(), libMesh::QSimpson::init_3D(), and 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 }
libMesh::QBase::_points
std::vector< Point > _points
Definition: quadrature.h:349
libMesh::QBase::_weights
std::vector< Real > _weights
Definition: quadrature.h:355

◆ tensor_product_prism()

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

Computes the tensor product of a 1D quadrature rule and a 2D quadrature rule. Used in the init_3D routines for prismatic element types.

Definition at line 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::QGrid::init_3D(), libMesh::QTrap::init_3D(), libMesh::QSimpson::init_3D(), and 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 }
libMesh::QBase::_points
std::vector< Point > _points
Definition: quadrature.h:349
libMesh::QBase::_weights
std::vector< Real > _weights
Definition: quadrature.h:355

◆ tensor_product_quad()

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

Constructs a 2D rule from the tensor product of q1D with itself. Used in the init_2D() routines for quadrilateral element types.

Definition at line 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::QGaussLobatto::init_2D(), libMesh::QTrap::init_2D(), libMesh::QGrid::init_2D(), libMesh::QSimpson::init_2D(), and init_2D().

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 }
libMesh::QBase::_points
std::vector< Point > _points
Definition: quadrature.h:349
libMesh::QBase::_weights
std::vector< Real > _weights
Definition: quadrature.h:355

◆ type()

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

Implements libMesh::QBase.

Definition at line 33 of file quadrature_gauss.C.

References libMesh::QGAUSS.

34 {
35  return QGAUSS;
36 }

◆ w()

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

Definition at line 174 of file quadrature.h.

References libMesh::QBase::_weights.

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

175  {
176  libmesh_assert_less (i, _weights.size());
177  return _weights[i];
178  }
libMesh::QBase::_weights
std::vector< Real > _weights
Definition: quadrature.h:355

Member Data Documentation

◆ _counts

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

Actually holds the data.

Definition at line 122 of file reference_counter.h.

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

◆ _dim

unsigned int libMesh::QBase::_dim
protectedinherited

The spatial dimension of the quadrature rule.

Definition at line 325 of file quadrature.h.

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

◆ _enable_print_counter

bool libMesh::ReferenceCounter::_enable_print_counter = true
staticprotectedinherited

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

Definition at line 141 of file reference_counter.h.

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

◆ _mutex

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

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 135 of file reference_counter.h.

◆ _n_objects

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

The number of objects. Print the reference count information when the number returns to 0.

Definition at line 130 of file reference_counter.h.

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

◆ _order

Order libMesh::QBase::_order
protectedinherited

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

Definition at line 331 of file quadrature.h.

Referenced by libMesh::QBase::build(), libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), libMesh::QBase::get_order(), libMesh::QClough::init_1D(), libMesh::QGaussLobatto::init_1D(), libMesh::QConical::init_1D(), libMesh::QGrid::init_1D(), init_1D(), libMesh::QJacobi::init_1D(), libMesh::QMonomial::init_1D(), libMesh::QGrundmann_Moller::init_1D(), libMesh::QClough::init_2D(), libMesh::QGaussLobatto::init_2D(), libMesh::QGrid::init_2D(), init_2D(), libMesh::QMonomial::init_2D(), libMesh::QGrundmann_Moller::init_2D(), libMesh::QGaussLobatto::init_3D(), libMesh::QGrid::init_3D(), init_3D(), libMesh::QMonomial::init_3D(), and libMesh::QGrundmann_Moller::init_3D().

◆ _p_level

unsigned int libMesh::QBase::_p_level
protectedinherited

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

Definition at line 343 of file quadrature.h.

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

◆ _points

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

The locations of the quadrature points in reference element space.

Definition at line 349 of file quadrature.h.

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

◆ _type

ElemType libMesh::QBase::_type
protectedinherited

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

Definition at line 337 of file quadrature.h.

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

◆ _weights

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

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

Definition at line 355 of file quadrature.h.

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

◆ allow_rules_with_negative_weights

bool libMesh::QBase::allow_rules_with_negative_weights
inherited

Flag (default true) controlling the use of quadrature rules with negative weights. Set this to false to require rules with all positive weights.

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

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

Definition at line 246 of file quadrature.h.

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

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