libMesh::QGauss Class Reference

#include <quadrature_gauss.h>

Inheritance diagram for libMesh::QGauss:

Public Member Functions

 QGauss (const unsigned int _dim, const Order _order=INVALID_ORDER)
 
 ~QGauss ()
 
QuadratureType type () const
 
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
 
void init (const ElemType type=INVALID_ELEM, 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 AutoPtr< QBasebuild (const std::string &name, const unsigned int _dim, const Order _order=INVALID_ORDER)
 
static AutoPtr< 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 increment_constructor_count (const std::string &name)
 
void increment_destructor_count (const std::string &name)
 

Protected Attributes

libMesh::err<< "ERROR: Seems
as if this quadrature rule"
<< std::endl<< " is not
implemented for 2D."<< std::endl;libmesh_error();}#endif virtual void init_3D(const ElemType, unsigned int=0)#ifndef DEBUG{}#else{libMesh::err<< "ERROR: Seems as if this quadrature rule"<< std::endl<< " is not implemented for 3D."<< std::endl;libmesh_error();}#endif void tensor_product_quad(const QBase &q1D);void tensor_product_hex(const QBase &q1D);void tensor_product_prism(const QBase &q1D, const QBase &q2D);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

void init_1D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0)
 
void init_2D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0)
 
void init_3D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0)
 
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

This class implemenets specific orders of Gauss quadrature. Gauss quadrature rules of Order p have the property of integrating polynomials of degree p exactly.

Definition at line 43 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 113 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.

Definition at line 105 of file quadrature_gauss.h.

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

106  : QBase(d,o)
107 {
108  // explicitly call the init function in 1D since the
109  // other tensor-product rules require this one.
110  // note that EDGE will not be used internally, however
111  // if we called the function with INVALID_ELEM it would try to
112  // be smart and return, thinking it had already done the work.
113  if (_dim == 1)
114  init(EDGE2);
115 }
libMesh::QGauss::~QGauss ( )
inline

Destructor.

Definition at line 121 of file quadrature_gauss.h.

122 {
123 }

Member Function Documentation

AutoPtr< QBase > libMesh::QBase::build ( const std::string &  name,
const unsigned int  _dim,
const Order  _order = INVALID_ORDER 
)
staticinherited

Builds a specific quadrature rule, identified through the name string. An AutoPtr<QBase> is returned to prevent a memory leak. This way the user need not remember to delete the object. Enables run-time decision of the quadrature rule. The input parameter name must be mappable through the Utility::string_to_enum<>() function.

Definition at line 37 of file quadrature_build.C.

Referenced by libMesh::InfFE< friend_Dim, friend_T_radial, friend_T_map >::attach_quadrature_rule().

40 {
41  return QBase::build (Utility::string_to_enum<QuadratureType> (type),
42  _dim,
43  _order);
44 }
AutoPtr< QBase > libMesh::QBase::build ( const QuadratureType  _qt,
const unsigned int  _dim,
const Order  _order = INVALID_ORDER 
)
staticinherited

Builds a specific quadrature rule, identified through the QuadratureType. An AutoPtr<QBase> is returned to prevent a memory leak. This way the user need not remember to delete the object. Enables run-time decision of the quadrature rule.

Definition at line 48 of file quadrature_build.C.

References libMesh::err, libMeshEnums::FIRST, libMeshEnums::FORTYTHIRD, libMesh::out, libMeshEnums::QCLOUGH, libMeshEnums::QGAUSS, libMeshEnums::QJACOBI_1_0, libMeshEnums::QJACOBI_2_0, libMeshEnums::QSIMPSON, libMeshEnums::QTRAP, libMeshEnums::THIRD, and libMeshEnums::TWENTYTHIRD.

51 {
52  switch (_qt)
53  {
54 
55  case QCLOUGH:
56  {
57 #ifdef DEBUG
58  if (_order > TWENTYTHIRD)
59  {
60  libMesh::out << "WARNING: Clough quadrature implemented" << std::endl
61  << " up to TWENTYTHIRD order." << std::endl;
62  }
63 #endif
64 
65  AutoPtr<QBase> ap(new QClough(_dim, _order));
66  return ap;
67  }
68 
69  case QGAUSS:
70  {
71 
72 #ifdef DEBUG
73  if (_order > FORTYTHIRD)
74  {
75  libMesh::out << "WARNING: Gauss quadrature implemented" << std::endl
76  << " up to FORTYTHIRD order." << std::endl;
77  }
78 #endif
79 
80  AutoPtr<QBase> ap(new QGauss(_dim, _order));
81  return ap;
82  }
83 
84  case QJACOBI_1_0:
85  {
86 
87 #ifdef DEBUG
88  if (_order > TWENTYTHIRD)
89  {
90  libMesh::out << "WARNING: Jacobi(1,0) quadrature implemented" << std::endl
91  << " up to TWENTYTHIRD order." << std::endl;
92  }
93 
94  if (_dim > 1)
95  {
96  libMesh::out << "WARNING: Jacobi(1,0) quadrature implemented" << std::endl
97  << " in 1D only." << std::endl;
98  }
99 #endif
100 
101  AutoPtr<QBase> ap(new QJacobi(_dim, _order, 1, 0));
102  return ap;
103  }
104 
105  case QJACOBI_2_0:
106  {
107 
108 #ifdef DEBUG
109  if (_order > TWENTYTHIRD)
110  {
111  libMesh::out << "WARNING: Jacobi(2,0) quadrature implemented" << std::endl
112  << " up to TWENTYTHIRD order." << std::endl;
113  }
114 
115  if (_dim > 1)
116  {
117  libMesh::out << "WARNING: Jacobi(2,0) quadrature implemented" << std::endl
118  << " in 1D only." << std::endl;
119  }
120 #endif
121 
122  AutoPtr<QBase> ap(new QJacobi(_dim, _order, 2, 0));
123  return ap;
124  }
125 
126  case QSIMPSON:
127  {
128 
129 #ifdef DEBUG
130  if (_order > THIRD)
131  {
132  libMesh::out << "WARNING: Simpson rule provides only" << std::endl
133  << " THIRD order!" << std::endl;
134  }
135 #endif
136 
137  AutoPtr<QBase> ap(new QSimpson(_dim));
138  return ap;
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  AutoPtr<QBase> ap(new QTrap(_dim));
153  return ap;
154  }
155 
156 
157  default:
158  {
159  libMesh::err << "ERROR: Bad qt=" << _qt << std::endl;
160  libmesh_error();
161  }
162  }
163 
164 
165  libmesh_error();
166  AutoPtr<QBase> ap(NULL);
167  return ap;
168 }
void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

107 {
108  _enable_print_counter = false;
109  return;
110 }
void libMesh::QGauss::dunavant_rule ( const Real  rule_data[][4],
const unsigned int  n_pts 
)
private

The Dunavant rule is for triangles. It takes permutation points and weights in a specific format as input and fills the pre-sized _points and _weights vectors. This function is only used internally by the TRI geometric elements.

Definition at line 266 of file quadrature_gauss.C.

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

Referenced by init_2D().

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

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

Referenced by init_2D().

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

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

Definition at line 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

101 {
102  _enable_print_counter = true;
103  return;
104 }
unsigned int libMesh::QBase::get_dim ( ) const
inlineinherited
ElemType libMesh::QBase::get_elem_type ( ) const
inlineinherited
Returns
the current element type we're set up for

Definition at line 104 of file quadrature.h.

105  { return _type; }
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 }
Order libMesh::QBase::get_order ( ) const
inlineinherited
Returns
the order of the quadrature rule.

Definition at line 169 of file quadrature.h.

Referenced by libMesh::InfFE< friend_Dim, friend_T_radial, friend_T_map >::attach_quadrature_rule().

169 { return static_cast<Order>(_order + _p_level); }
unsigned int libMesh::QBase::get_p_level ( ) const
inlineinherited
Returns
the current p refinement level we're initialized with

Definition at line 110 of file quadrature.h.

111  { return _p_level; }
const std::vector<Point>& libMesh::QBase::get_points ( ) const
inlineinherited
Returns
a std::vector containing the quadrature point locations on a reference object.

Definition at line 129 of file quadrature.h.

References libMesh::QBase::_points.

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

129 { return _points; }
std::vector<Point>& libMesh::QBase::get_points ( )
inlineinherited
Returns
a std::vector containing the quadrature point locations on a reference object as a writeable reference.

Definition at line 135 of file quadrature.h.

References libMesh::QBase::_points.

135 { return _points; }
const std::vector<Real>& libMesh::QBase::get_weights ( ) const
inlineinherited
Returns
a std::vector containing the quadrature weights.

Definition at line 140 of file quadrature.h.

References libMesh::QBase::_weights.

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

140 { return _weights; }
std::vector<Real>& libMesh::QBase::get_weights ( )
inlineinherited
Returns
a std::vector containing the quadrature weights.

Definition at line 145 of file quadrature.h.

References libMesh::QBase::_weights.

145 { return _weights; }
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 163 of file reference_counter.h.

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

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

164 {
165  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
166  std::pair<unsigned int, unsigned int>& p = _counts[name];
167 
168  p.first++;
169 }
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 176 of file reference_counter.h.

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

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

177 {
178  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
179  std::pair<unsigned int, unsigned int>& p = _counts[name];
180 
181  p.second++;
182 }
void libMesh::QBase::init ( const ElemType  type = INVALID_ELEM,
unsigned int  p_level = 0 
)
inherited

Initializes the data structures to contain a quadrature rule for an object of type type.

Definition at line 27 of file quadrature.C.

References libMesh::QBase::init_0D(), libMesh::QBase::init_1D(), and libMesh::QBase::init_2D().

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

29 {
30  // check to see if we have already
31  // done the work for this quadrature rule
32  if (t == _type && p == _p_level)
33  return;
34  else
35  {
36  _type = t;
37  _p_level = p;
38  }
39 
40 
41 
42  switch(_dim)
43  {
44  case 0:
45  this->init_0D(_type,_p_level);
46 
47  return;
48 
49  case 1:
50  this->init_1D(_type,_p_level);
51 
52  return;
53 
54  case 2:
55  this->init_2D(_type,_p_level);
56 
57  return;
58 
59  case 3:
60  this->init_3D(_type,_p_level);
61 
62  return;
63 
64  default:
65  libmesh_error();
66  }
67 }
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 71 of file quadrature.C.

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

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

73 {
74  _points.resize(1);
75  _weights.resize(1);
76  _points[0] = Point(0.);
77  _weights[0] = 1.0;
78 }
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::_points, libMesh::QBase::_weights, libMeshEnums::CONSTANT, libMeshEnums::EIGHTH, libMeshEnums::EIGHTTEENTH, libMeshEnums::ELEVENTH, libMesh::err, libMeshEnums::FIFTEENTH, libMeshEnums::FIFTH, libMeshEnums::FIRST, libMeshEnums::FORTIETH, libMeshEnums::FORTYFIRST, libMeshEnums::FORTYSECOND, libMeshEnums::FORTYTHIRD, libMeshEnums::FOURTEENTH, libMeshEnums::FOURTH, libMeshEnums::NINTEENTH, libMeshEnums::NINTH, libMeshEnums::SECOND, libMeshEnums::SEVENTEENTH, libMeshEnums::SEVENTH, libMeshEnums::SIXTEENTH, libMeshEnums::SIXTH, libMeshEnums::TENTH, libMeshEnums::THIRD, libMeshEnums::THIRTEENTH, libMeshEnums::THIRTIETH, libMeshEnums::THIRTYEIGHTH, libMeshEnums::THIRTYFIFTH, libMeshEnums::THIRTYFIRST, libMeshEnums::THIRTYFOURTH, libMeshEnums::THIRTYNINTH, libMeshEnums::THIRTYSECOND, libMeshEnums::THIRTYSEVENTH, libMeshEnums::THIRTYSIXTH, libMeshEnums::THIRTYTHIRD, libMeshEnums::TWELFTH, libMeshEnums::TWENTIETH, libMeshEnums::TWENTYEIGHTH, libMeshEnums::TWENTYFIFTH, libMeshEnums::TWENTYFIRST, libMeshEnums::TWENTYFOURTH, libMeshEnums::TWENTYNINTH, libMeshEnums::TWENTYSECOND, libMeshEnums::TWENTYSEVENTH, libMeshEnums::TWENTYSIXTH, and libMeshEnums::TWENTYTHIRD.

32 {
33  //----------------------------------------------------------------------
34  // 1D quadrature rules
35  switch(_order + 2*p)
36  {
37  case CONSTANT:
38  case FIRST:
39  {
40  _points.resize (1);
41  _weights.resize(1);
42 
43  _points[0](0) = 0.;
44 
45  _weights[0] = 2.;
46 
47  return;
48  }
49  case SECOND:
50  case THIRD:
51  {
52  _points.resize (2);
53  _weights.resize(2);
54 
55  _points[0](0) = -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 NINTEENTH:
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  {
779  libMesh::err << "Quadrature rule " << _order
780  << " not supported!" << std::endl;
781 
782  libmesh_error();
783  }
784  }
785 
786 
787 
788  return;
789 }
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 defined) when called.

Reimplemented from libMesh::QBase.

Definition at line 28 of file quadrature_gauss_2D.C.

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

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

Definition at line 28 of file quadrature_gauss_3D.C.

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

30 {
31 #if LIBMESH_DIM == 3
32 
33  //-----------------------------------------------------------------------
34  // 3D quadrature rules
35  switch (type_in)
36  {
37  //---------------------------------------------
38  // Hex quadrature rules
39  case HEX8:
40  case HEX20:
41  case HEX27:
42  {
43  // We compute the 3D quadrature rule as a tensor
44  // product of the 1D quadrature rule.
45  QGauss q1D(1,_order);
46  q1D.init(EDGE2,p);
47 
48  tensor_product_hex( q1D );
49 
50  return;
51  }
52 
53 
54 
55  //---------------------------------------------
56  // Tetrahedral quadrature rules
57  case TET4:
58  case TET10:
59  {
60  switch(_order + 2*p)
61  {
62  // Taken from pg. 222 of "The finite element method," vol. 1
63  // ed. 5 by Zienkiewicz & Taylor
64  case CONSTANT:
65  case FIRST:
66  {
67  // Exact for linears
68  _points.resize(1);
69  _weights.resize(1);
70 
71 
72  _points[0](0) = .25;
73  _points[0](1) = .25;
74  _points[0](2) = .25;
75 
76  _weights[0] = .1666666666666666666666666666666666666666666667L;
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] = .0416666666666666666666666666666666666666666667;
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) = .16666666666666666666666666666666666666666667;
143  _points[1](2) = .16666666666666666666666666666666666666666667;
144 
145  _points[2](0) = .16666666666666666666666666666666666666666667;
146  _points[2](1) = .5;
147  _points[2](2) = .16666666666666666666666666666666666666666667;
148 
149  _points[3](0) = .16666666666666666666666666666666666666666667;
150  _points[3](1) = .16666666666666666666666666666666666666666667;
151  _points[3](2) = .5;
152 
153  _points[4](0) = .16666666666666666666666666666666666666666667;
154  _points[4](1) = .16666666666666666666666666666666666666666667;
155  _points[4](2) = .16666666666666666666666666666666666666666667;
156 
157 
158  _weights[0] = -.133333333333333333333333333333333333333333333;
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 = 0.333333333333333333333333333333333333333;
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 = 0.7272727272727272727272727272727272727272727272727272727;
320  // const Real b = 0.0909090909090909090909090909090909090909090909090909091;
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 PYRAMID14:
529  {
530  // We compute the Pyramid rule as a conical product of a
531  // Jacobi rule with alpha==2 on the interval [0,1] two 1D
532  // Gauss rules with suitably modified points. The idea comes
533  // from: Stroud, A.H. "Approximate Calculation of Multiple
534  // Integrals."
535  QConical conical_rule(3, _order);
536  conical_rule.init(type_in, p);
537 
538  // Swap points and weights with the about-to-be destroyed rule.
539  _points.swap (conical_rule.get_points() );
540  _weights.swap(conical_rule.get_weights());
541 
542  return;
543 
544  }
545 
546 
547 
548  //---------------------------------------------
549  // Unsupported type
550  default:
551  {
552  libMesh::err << "ERROR: Unsupported type: " << type_in << std::endl;
553  libmesh_error();
554  }
555  }
556 
557  libmesh_error();
558 
559  return;
560 
561 #endif
562 }
void libMesh::QGauss::keast_rule ( const Real  rule_data[][4],
const unsigned int  n_pts 
)
private

The Keast rule is for tets. It takes permutation points and weights in a specific format as input and fills the pre-sized _points and _weights vectors. This function is only used internally by the TET geometric elements.

Definition at line 33 of file quadrature_gauss.C.

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

Referenced by init_3D().

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  {
172  libMesh::err << "Don't know what to do with this many permutation points!" << std::endl;
173  libmesh_error();
174  }
175  }
176 
177  }
178 
179 }
static unsigned int libMesh::ReferenceCounter::n_objects ( )
inlinestaticinherited

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

Definition at line 79 of file reference_counter.h.

References libMesh::ReferenceCounter::_n_objects.

80  { return _n_objects; }
unsigned int libMesh::QBase::n_points ( ) const
inlineinherited
Returns
the number of points associated with the quadrature rule.

Definition at line 116 of file quadrature.h.

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

Referenced by libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), libMesh::ProjectFEMSolution::operator()(), and libMesh::QBase::print_info().

117  { libmesh_assert (!_points.empty());
118  return libmesh_cast_int<unsigned int>(_points.size()); }
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().

89 {
91 }
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 362 of file quadrature.h.

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

Referenced by libMesh::operator<<().

363 {
364  libmesh_assert(!_points.empty());
365  libmesh_assert(!_weights.empty());
366 
367  os << "N_Q_Points=" << this->n_points() << std::endl << std::endl;
368  for (unsigned int qpoint=0; qpoint<this->n_points(); qpoint++)
369  {
370  os << " Point " << qpoint << ":\n"
371  << " "
372  << _points[qpoint]
373  << " Weight:\n "
374  << " w=" << _weights[qpoint] << "\n" << std::endl;
375  }
376 }
Point libMesh::QBase::qp ( const unsigned int  i) const
inlineinherited
Returns
the $ i^{th} $ quadrature point on the reference object.

Definition at line 150 of file quadrature.h.

References libMesh::QBase::_points.

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

151  { libmesh_assert_less (i, _points.size()); return _points[i]; }
void libMesh::QBase::scale ( std::pair< Real, Real old_range,
std::pair< Real, Real new_range 
)
inherited

Maps the points of a 1D interval quadrature rule (typically [-1,1]) to any other 1D interval (typically [0,1]) and scales the weights accordingly. The quadrature rule will be mapped from the entries of old_range to the entries of new_range.

Definition at line 82 of file quadrature.C.

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

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

84 {
85  // Make sure we are in 1D
86  libmesh_assert_equal_to (_dim, 1);
87 
88  // Make sure that we have sane ranges
89  libmesh_assert_greater (new_range.second, new_range.first);
90  libmesh_assert_greater (old_range.second, old_range.first);
91 
92  // Make sure there are some points
93  libmesh_assert_greater (_points.size(), 0);
94 
95  // We're mapping from old_range -> new_range
96  for (unsigned int i=0; i<_points.size(); i++)
97  {
98  _points[i](0) =
99  (_points[i](0) - old_range.first) *
100  (new_range.second - new_range.first) /
101  (old_range.second - old_range.first) +
102  new_range.first;
103  }
104 
105  // Compute the scale factor and scale the weights
106  const Real scfact = (new_range.second - new_range.first) /
107  (old_range.second - old_range.first);
108 
109  for (unsigned int i=0; i<_points.size(); i++)
110  _weights[i] *= scfact;
111 }
virtual bool libMesh::QBase::shapes_need_reinit ( )
inlinevirtualinherited

Returns true if the shape functions need to be recalculated.

This can happen if the number of points or their positions change.

By default this will return false.

Definition at line 198 of file quadrature.h.

198 { return false; }
QuadratureType libMesh::QGauss::type ( ) const
inlinevirtual
Returns
QGAUSS

Implements libMesh::QBase.

Definition at line 61 of file quadrature_gauss.h.

References libMeshEnums::QGAUSS.

61 { return QGAUSS; }
Real libMesh::QBase::w ( const unsigned int  i) const
inlineinherited
Returns
the $ i^{th} $ quadrature weight.

Definition at line 156 of file quadrature.h.

References libMesh::QBase::_weights.

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

157  { libmesh_assert_less (i, _weights.size()); return _weights[i]; }

Member Data Documentation

ReferenceCounter::Counts libMesh::ReferenceCounter::_counts
staticprotectedinherited
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 137 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 131 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 126 of file reference_counter.h.

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

libMesh::err<< "ERROR: Seems as if this quadrature rule" << std::endl << " is not implemented for 2D." << std::endl; libmesh_error(); }#endif virtual void init_3D (const ElemType, unsigned int =0)#ifndef DEBUG {}#else { libMesh::err << "ERROR: Seems as if this quadrature rule" << std::endl << " is not implemented for 3D." << std::endl; libmesh_error(); }#endif void tensor_product_quad (const QBase& q1D); void tensor_product_hex (const QBase& q1D); void tensor_product_prism (const QBase& q1D, const QBase& q2D); const unsigned int _dim; const Order _order; ElemType _type; unsigned int _p_level; std::vector<Point> libMesh::QBase::_points
protectedinherited
bool libMesh::QBase::allow_rules_with_negative_weights
inherited

Flag (default true) controlling the use of quadrature rules with negative weights. Set this to false to ONLY use (potentially) safer but more expensive rules with all positive weights.

Negative weights typically appear in Gaussian quadrature rules over three-dimensional elements. 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 215 of file quadrature.h.

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


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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:58:01 UTC

Hosted By:
SourceForge.net Logo