libMesh::QGauss Class Reference

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

#include <quadrature_gauss.h>

Inheritance diagram for libMesh::QGauss:

Public Member Functions

 QGauss (const unsigned int _dim, const Order _order=INVALID_ORDER)
 
 ~QGauss ()
 
virtual QuadratureType type () const libmesh_override
 
ElemType get_elem_type () const
 
unsigned int get_p_level () const
 
unsigned int n_points () const
 
unsigned int get_dim () const
 
const std::vector< Point > & get_points () const
 
std::vector< Point > & get_points ()
 
const std::vector< Real > & get_weights () const
 
std::vector< Real > & get_weights ()
 
Point qp (const unsigned int i) const
 
Real w (const unsigned int i) const
 
virtual void init (const ElemType type=INVALID_ELEM, unsigned int p_level=0)
 
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 UniquePtr< QBasebuild (const std::string &name, const unsigned int dim, const Order order=INVALID_ORDER)
 
static UniquePtr< QBasebuild (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

const unsigned int _dim
 
const 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) libmesh_override
 
virtual void init_2D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0) libmesh_override
 
virtual void init_3D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0) libmesh_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 implemenets 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

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

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

Definition at line 119 of file reference_counter.h.

Constructor & Destructor Documentation

libMesh::QGauss::QGauss ( const unsigned int  _dim,
const Order  _order = INVALID_ORDER 
)
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  }
virtual void init(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Definition: quadrature.C:28
const unsigned int _dim
Definition: quadrature.h:311
QBase(const unsigned int _dim, const Order _order=INVALID_ORDER)
Definition: quadrature.h:350
const Order _order
Definition: quadrature.h:317
libMesh::QGauss::~QGauss ( )
inline

Destructor.

Definition at line 63 of file quadrature_gauss.h.

63 {}

Member Function Documentation

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

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

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

Definition at line 42 of file quadrature_build.C.

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

45 {
46  return QBase::build (Utility::string_to_enum<QuadratureType> (type),
47  _dim,
48  _order);
49 }
const unsigned int _dim
Definition: quadrature.h:311
virtual QuadratureType type() const =0
static UniquePtr< QBase > build(const std::string &name, const unsigned int dim, const Order order=INVALID_ORDER)
const Order _order
Definition: quadrature.h:317
UniquePtr< QBase > libMesh::QBase::build ( const QuadratureType  qt,
const unsigned int  dim,
const Order  order = INVALID_ORDER 
)
staticinherited

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

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

Definition at line 53 of file quadrature_build.C.

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

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

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

Definition at line 260 of file quadrature_gauss.C.

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

Referenced by init_2D(), and type().

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

Definition at line 185 of file quadrature_gauss.C.

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

Referenced by init_2D(), and type().

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

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

Definition at line 101 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

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

102 {
103  _enable_print_counter = true;
104  return;
105 }
unsigned int libMesh::QBase::get_dim ( ) const
inlineinherited
Returns
The spatial dimension of the quadrature rule.

Definition at line 122 of file quadrature.h.

References libMesh::QBase::_dim.

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

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

Definition at line 103 of file quadrature.h.

References libMesh::QBase::_type.

103 { return _type; }
ElemType _type
Definition: quadrature.h:323
std::string libMesh::ReferenceCounter::get_info ( )
staticinherited

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

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

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

48 {
49 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
50 
51  std::ostringstream oss;
52 
53  oss << '\n'
54  << " ---------------------------------------------------------------------------- \n"
55  << "| Reference count information |\n"
56  << " ---------------------------------------------------------------------------- \n";
57 
58  for (Counts::iterator it = _counts.begin();
59  it != _counts.end(); ++it)
60  {
61  const std::string name(it->first);
62  const unsigned int creations = it->second.first;
63  const unsigned int destructions = it->second.second;
64 
65  oss << "| " << name << " reference count information:\n"
66  << "| Creations: " << creations << '\n'
67  << "| Destructions: " << destructions << '\n';
68  }
69 
70  oss << " ---------------------------------------------------------------------------- \n";
71 
72  return oss.str();
73 
74 #else
75 
76  return "";
77 
78 #endif
79 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
Order libMesh::QBase::get_order ( ) const
inlineinherited
Returns
The order of the quadrature rule.

Definition at line 189 of file quadrature.h.

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

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

189 { return static_cast<Order>(_order + _p_level); }
unsigned int _p_level
Definition: quadrature.h:329
const Order _order
Definition: quadrature.h:317
unsigned int libMesh::QBase::get_p_level ( ) const
inlineinherited
Returns
The p-refinement level we're currently using.

Definition at line 108 of file quadrature.h.

References libMesh::QBase::_p_level.

108 { return _p_level; }
unsigned int _p_level
Definition: quadrature.h:329
const std::vector<Point>& libMesh::QBase::get_points ( ) const
inlineinherited
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 134 of file quadrature.h.

References libMesh::QBase::_points.

134 { return _points; }
std::vector< Point > _points
Definition: quadrature.h:335
const std::vector<Real>& libMesh::QBase::get_weights ( ) const
inlineinherited
std::vector<Real>& libMesh::QBase::get_weights ( )
inlineinherited
Returns
A writable references to a std::vector containing the quadrature weights.

Definition at line 146 of file quadrature.h.

References libMesh::QBase::_weights.

146 { return _weights; }
std::vector< Real > _weights
Definition: quadrature.h:341
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 185 of file reference_counter.h.

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

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

186 {
187  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
188  std::pair<unsigned int, unsigned int> & p = _counts[name];
189 
190  p.first++;
191 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
spin_mutex spin_mtx
Definition: threads.C:29
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 198 of file reference_counter.h.

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

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

199 {
200  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
201  std::pair<unsigned int, unsigned int> & p = _counts[name];
202 
203  p.second++;
204 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
spin_mutex spin_mtx
Definition: threads.C:29
void libMesh::QBase::init ( const ElemType  type = INVALID_ELEM,
unsigned int  p_level = 0 
)
virtualinherited

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

Definition at line 28 of file quadrature.C.

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

Referenced by libMesh::FESubdivision::attach_quadrature_rule(), libMesh::QBase::init(), libMesh::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(), libMesh::QSimpson::init_3D(), init_3D(), libMesh::QMonomial::init_3D(), QGauss(), libMesh::QGaussLobatto::QGaussLobatto(), libMesh::QJacobi::QJacobi(), libMesh::QSimpson::QSimpson(), libMesh::QTrap::QTrap(), and libMesh::QBase::w().

30 {
31  // check to see if we have already
32  // done the work for this quadrature rule
33  if (t == _type && p == _p_level)
34  return;
35  else
36  {
37  _type = t;
38  _p_level = p;
39  }
40 
41 
42 
43  switch(_dim)
44  {
45  case 0:
46  this->init_0D(_type,_p_level);
47 
48  return;
49 
50  case 1:
51  this->init_1D(_type,_p_level);
52 
53  return;
54 
55  case 2:
56  this->init_2D(_type,_p_level);
57 
58  return;
59 
60  case 3:
61  this->init_3D(_type,_p_level);
62 
63  return;
64 
65  default:
66  libmesh_error_msg("Invalid dimension _dim = " << _dim);
67  }
68 }
const unsigned int _dim
Definition: quadrature.h:311
ElemType _type
Definition: quadrature.h:323
unsigned int _p_level
Definition: quadrature.h:329
virtual void init_2D(const ElemType, unsigned int=0)
Definition: quadrature.h:263
virtual void init_1D(const ElemType type=INVALID_ELEM, unsigned int p_level=0)=0
virtual void init_0D(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Definition: quadrature.C:82
virtual void init_3D(const ElemType, unsigned int=0)
Definition: quadrature.h:279
void libMesh::QBase::init ( const Elem elem,
const std::vector< Real > &  vertex_distance_func,
unsigned int  p_level = 0 
)
virtualinherited

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

Definition at line 72 of file quadrature.C.

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

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

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

Definition at line 82 of file quadrature.C.

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

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

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

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

Implements libMesh::QBase.

Definition at line 30 of file quadrature_gauss_1D.C.

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

Referenced by type().

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

Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate values. The order of the rule will be defined by the implementing class. Should not be pure virtual since a derived quadrature rule may only be defined in 1D. If not redefined, gives an error (when DEBUG is defined) when called.

Reimplemented from libMesh::QBase.

Definition at line 28 of file quadrature_gauss_2D.C.

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

Referenced by type().

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

Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate values. The order of the rule will be defined by the implementing class. Should not be pure virtual since a derived quadrature rule may only be defined in 1D. If not redefined, gives an error (when DEBUG is defined) when called.

Reimplemented from libMesh::QBase.

Definition at line 28 of file quadrature_gauss_3D.C.

References libMesh::QBase::_order, libMesh::QBase::_points, libMesh::QBase::_weights, libMesh::QBase::allow_rules_with_negative_weights, libMesh::CONSTANT, libMesh::EDGE2, libMesh::EIGHTH, libMesh::FIFTH, libMesh::FIRST, libMesh::FOURTH, libMesh::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.

Referenced by type().

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

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

Definition at line 33 of file quadrature_gauss.C.

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

Referenced by init_3D(), and type().

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

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

Definition at line 88 of file reference_counter.C.

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

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

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

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

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

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

Definition at line 151 of file quadrature.h.

References libMesh::QBase::_points.

Referenced by libMesh::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().

152  {
153  libmesh_assert_less (i, _points.size());
154  return _points[i];
155  }
std::vector< Point > _points
Definition: quadrature.h:335
void libMesh::QBase::scale ( std::pair< Real, Real old_range,
std::pair< Real, Real new_range 
)
inherited

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

Definition at line 93 of file quadrature.C.

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

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

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

217 { return false; }
void libMesh::QBase::tensor_product_hex ( const QBase q1D)
protectedinherited

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

Definition at line 154 of file quadrature.C.

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

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

155 {
156  const unsigned int np = q1D.n_points();
157 
158  _points.resize(np * np * np);
159 
160  _weights.resize(np * np * np);
161 
162  unsigned int q=0;
163 
164  for (unsigned int k=0; k<np; k++)
165  for (unsigned int j=0; j<np; j++)
166  for (unsigned int i=0; i<np; i++)
167  {
168  _points[q](0) = q1D.qp(i)(0);
169  _points[q](1) = q1D.qp(j)(0);
170  _points[q](2) = q1D.qp(k)(0);
171 
172  _weights[q] = q1D.w(i) * q1D.w(j) * q1D.w(k);
173 
174  q++;
175  }
176 }
std::vector< Point > _points
Definition: quadrature.h:335
std::vector< Real > _weights
Definition: quadrature.h:341
void libMesh::QBase::tensor_product_prism ( const QBase q1D,
const QBase q2D 
)
protectedinherited

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

Definition at line 181 of file quadrature.C.

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

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

182 {
183  const unsigned int n_points1D = q1D.n_points();
184  const unsigned int n_points2D = q2D.n_points();
185 
186  _points.resize (n_points1D * n_points2D);
187  _weights.resize (n_points1D * n_points2D);
188 
189  unsigned int q=0;
190 
191  for (unsigned int j=0; j<n_points1D; j++)
192  for (unsigned int i=0; i<n_points2D; i++)
193  {
194  _points[q](0) = q2D.qp(i)(0);
195  _points[q](1) = q2D.qp(i)(1);
196  _points[q](2) = q1D.qp(j)(0);
197 
198  _weights[q] = q2D.w(i) * q1D.w(j);
199 
200  q++;
201  }
202 
203 }
std::vector< Point > _points
Definition: quadrature.h:335
std::vector< Real > _weights
Definition: quadrature.h:341
void libMesh::QBase::tensor_product_quad ( const QBase q1D)
protectedinherited

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

Definition at line 127 of file quadrature.C.

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

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

128 {
129 
130  const unsigned int np = q1D.n_points();
131 
132  _points.resize(np * np);
133 
134  _weights.resize(np * np);
135 
136  unsigned int q=0;
137 
138  for (unsigned int j=0; j<np; j++)
139  for (unsigned int i=0; i<np; i++)
140  {
141  _points[q](0) = q1D.qp(i)(0);
142  _points[q](1) = q1D.qp(j)(0);
143 
144  _weights[q] = q1D.w(i)*q1D.w(j);
145 
146  q++;
147  }
148 }
std::vector< Point > _points
Definition: quadrature.h:335
std::vector< Real > _weights
Definition: quadrature.h:341
virtual QuadratureType libMesh::QGauss::type ( ) const
inlinevirtual
Real libMesh::QBase::w ( const unsigned int  i) const
inlineinherited

Member Data Documentation

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

The spatial dimension of the quadrature rule.

Definition at line 311 of file quadrature.h.

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

bool libMesh::ReferenceCounter::_enable_print_counter = true
staticprotectedinherited

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

Definition at line 143 of file reference_counter.h.

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

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

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 137 of file reference_counter.h.

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

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

Definition at line 132 of file reference_counter.h.

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

unsigned int libMesh::QBase::_p_level
protectedinherited

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

Definition at line 329 of file quadrature.h.

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

bool libMesh::QBase::allow_rules_with_negative_weights
inherited

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

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

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

Definition at line 232 of file quadrature.h.

Referenced by 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: