libMesh::QMonomial Class Reference

#include <quadrature_monomial.h>

Inheritance diagram for libMesh::QMonomial:

Public Member Functions

 QMonomial (const unsigned int _dim, const Order _order=INVALID_ORDER)
 
 ~QMonomial ()
 
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, unsigned int=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 wissmann_rule (const Real rule_data[][3], const unsigned int n_pts)
 
void stroud_rule (const Real rule_data[][3], const unsigned int *rule_symmetry, const unsigned int n_pts)
 
void kim_rule (const Real rule_data[][4], const unsigned int *rule_id, const unsigned int n_pts)
 

Detailed Description

This class defines alternate quadrature rules on "tensor-product" elements (QUADs and HEXes) which can be useful when integrating monomial finite element bases.

While tensor product rules are ideal for integrating bi/tri-linear, bi/tri-quadratic, etc. (i.e. tensor product) bases (which consist of incomplete polynomials up to degree= dim*p) they are not optimal for the MONOMIAL or FEXYZ bases, which consist of complete polynomials of degree=p.

This class is implemented to provide quadrature rules which are more efficient than tensor product rules when they are available, and fall back on Gaussian quadrature rules when necessary.

A number of these rules have been helpfully collected in electronic form by:

Prof. Ronald Cools Katholieke Universiteit Leuven, Dept. Computerwetenschappen http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html

(A username and password to access the tables is available by request.)

We also provide the original reference for each rule, as available, in the source code file.

Author
John W. Peterson, 2008

Definition at line 65 of file quadrature_monomial.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::QMonomial::QMonomial ( const unsigned int  _dim,
const Order  _order = INVALID_ORDER 
)

Constructor. Declares the order of the quadrature rule.

Definition at line 33 of file quadrature_monomial.C.

34  : QBase(d,o)
35 {
36 }
libMesh::QMonomial::~QMonomial ( )

Destructor.

Definition at line 40 of file quadrature_monomial.C.

41 {
42 }

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::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(), libMesh::QGauss::init_2D(), init_2D(), libMesh::QGauss::init_3D(), and 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(), libMesh::QGauss::init_2D(), init_2D(), libMesh::QGauss::init_3D(), and 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(), libMesh::QGauss::init_2D(), libMesh::QSimpson::init_2D(), init_2D(), libMesh::QTrap::init_3D(), libMesh::QGrid::init_3D(), libMesh::QGauss::init_3D(), libMesh::QSimpson::init_3D(), init_3D(), libMesh::QGauss::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::QMonomial::init_1D ( const ElemType  type,
unsigned  p_level = 0 
)
inlineprivatevirtual

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 88 of file quadrature_monomial.h.

90  {
91  // See about making this non-pure virtual in the base class?
92  libmesh_error();
93  }
void libMesh::QMonomial::init_2D ( const ElemType  _type = INVALID_ELEM,
unsigned int  p_level = 0 
)
privatevirtual

More efficient rules for QUADs

Reimplemented from libMesh::QBase.

Definition at line 28 of file quadrature_monomial_2D.C.

References libMesh::QBase::_points, libMesh::QBase::_weights, data, libMeshEnums::EIGHTH, libMeshEnums::ELEVENTH, libMeshEnums::FIFTEENTH, libMeshEnums::FIFTH, libMeshEnums::FOURTEENTH, libMeshEnums::FOURTH, libMesh::QBase::get_points(), libMesh::QBase::get_weights(), libMesh::QBase::init(), libMeshEnums::NINTH, libMeshEnums::QUAD4, libMeshEnums::QUAD8, libMeshEnums::QUAD9, libMesh::Real, libMeshEnums::SECOND, libMeshEnums::SEVENTEENTH, libMeshEnums::SEVENTH, libMeshEnums::SIXTEENTH, libMeshEnums::SIXTH, stroud_rule(), libMeshEnums::TENTH, libMeshEnums::THIRTEENTH, libMeshEnums::TWELFTH, and wissmann_rule().

30 {
31 
32  switch (type_in)
33  {
34  //---------------------------------------------
35  // Quadrilateral quadrature rules
36  case QUAD4:
37  case QUAD8:
38  case QUAD9:
39  {
40  switch(_order + 2*p)
41  {
42  case SECOND:
43  {
44  // A degree=2 rule for the QUAD with 3 points.
45  // A tensor product degree-2 Gauss would have 4 points.
46  // This rule (or a variation on it) is probably available in
47  //
48  // A.H. Stroud, Approximate calculation of multiple integrals,
49  // Prentice-Hall, Englewood Cliffs, N.J., 1971.
50  //
51  // though I have never actually seen a reference for it.
52  // Luckily it's fairly easy to derive, which is what I've done
53  // here [JWP].
54  const Real
55  s=std::sqrt(1./3.),
56  t=std::sqrt(2./3.);
57 
58  const Real data[2][3] =
59  {
60  {0.0, s, 2.0},
61  { t, -s, 1.0}
62  };
63 
64  _points.resize(3);
65  _weights.resize(3);
66 
67  wissmann_rule(data, 2);
68 
69  return;
70  } // end case SECOND
71 
72 
73 
74  // For third-order, fall through to default case, use 2x2 Gauss product rule.
75  // case THIRD:
76  // {
77  // } // end case THIRD
78 
79  case FOURTH:
80  {
81  // A pair of degree=4 rules for the QUAD "C2" due to
82  // Wissmann and Becker. These rules both have six points.
83  // A tensor product degree-4 Gauss would have 9 points.
84  //
85  // J. W. Wissmann and T. Becker, Partially symmetric cubature
86  // formulas for even degrees of exactness, SIAM J. Numer. Anal. 23
87  // (1986), 676--685.
88  const Real data[4][3] =
89  {
90  // First of 2 degree-4 rules given by Wissmann
91  {0.0000000000000000e+00, 0.0000000000000000e+00, 1.1428571428571428e+00},
92  {0.0000000000000000e+00, 9.6609178307929590e-01, 4.3956043956043956e-01},
93  {8.5191465330460049e-01, 4.5560372783619284e-01, 5.6607220700753210e-01},
94  {6.3091278897675402e-01, -7.3162995157313452e-01, 6.4271900178367668e-01}
95  //
96  // Second of 2 degree-4 rules given by Wissmann. These both
97  // yield 4th-order accurate rules, I just chose the one that
98  // happened to contain the origin.
99  // {0.000000000000000, -0.356822089773090, 1.286412084888852},
100  // {0.000000000000000, 0.934172358962716, 0.491365692888926},
101  // {0.774596669241483, 0.390885162530071, 0.761883709085613},
102  // {0.774596669241483, -0.852765377881771, 0.349227402025498}
103  };
104 
105  _points.resize(6);
106  _weights.resize(6);
107 
108  wissmann_rule(data, 4);
109 
110  return;
111  } // end case FOURTH
112 
113 
114 
115 
116  case FIFTH:
117  {
118  // A degree 5, 7-point rule due to Stroud.
119  //
120  // A.H. Stroud, Approximate calculation of multiple integrals,
121  // Prentice-Hall, Englewood Cliffs, N.J., 1971.
122  //
123  // This rule is provably minimal in the number of points.
124  // A tensor-product rule accurate for "bi-quintic" polynomials would have 9 points.
125  const Real data[3][3] =
126  {
127  { 0.L, 0.L, static_cast<Real>(8.L / 7.L)}, // 1
128  { 0.L, static_cast<Real>(std::sqrt(14.L/15.L)), static_cast<Real>(20.L / 63.L)}, // 2
129  {static_cast<Real>(std::sqrt(3.L/5.L)), static_cast<Real>(std::sqrt(1.L/3.L)), static_cast<Real>(20.L / 36.L)} // 4
130  };
131 
132  const unsigned int symmetry[3] = {
133  0, // Origin
134  7, // Central Symmetry
135  6 // Rectangular
136  };
137 
138  _points.resize (7);
139  _weights.resize(7);
140 
141  stroud_rule(data, symmetry, 3);
142 
143  return;
144  } // end case FIFTH
145 
146 
147 
148 
149  case SIXTH:
150  {
151  // A pair of degree=6 rules for the QUAD "C2" due to
152  // Wissmann and Becker. These rules both have 10 points.
153  // A tensor product degree-6 Gauss would have 16 points.
154  //
155  // J. W. Wissmann and T. Becker, Partially symmetric cubature
156  // formulas for even degrees of exactness, SIAM J. Numer. Anal. 23
157  // (1986), 676--685.
158  const Real data[6][3] =
159  {
160  // First of 2 degree-6, 10 point rules given by Wissmann
161  // {0.000000000000000, 0.836405633697626, 0.455343245714174},
162  // {0.000000000000000, -0.357460165391307, 0.827395973202966},
163  // {0.888764014654765, 0.872101531193131, 0.144000884599645},
164  // {0.604857639464685, 0.305985162155427, 0.668259104262665},
165  // {0.955447506641064, -0.410270899466658, 0.225474004890679},
166  // {0.565459993438754, -0.872869311156879, 0.320896396788441}
167  //
168  // Second of 2 degree-6, 10 point rules given by Wissmann.
169  // Either of these will work, I just chose the one with points
170  // slightly further into the element interior.
171  {0.0000000000000000e+00, 8.6983337525005900e-01, 3.9275059096434794e-01},
172  {0.0000000000000000e+00, -4.7940635161211124e-01, 7.5476288124261053e-01},
173  {8.6374282634615388e-01, 8.0283751620765670e-01, 2.0616605058827902e-01},
174  {5.1869052139258234e-01, 2.6214366550805818e-01, 6.8999213848986375e-01},
175  {9.3397254497284950e-01, -3.6309658314806653e-01, 2.6051748873231697e-01},
176  {6.0897753601635630e-01, -8.9660863276245265e-01, 2.6956758608606100e-01}
177  };
178 
179  _points.resize(10);
180  _weights.resize(10);
181 
182  wissmann_rule(data, 6);
183 
184  return;
185  } // end case SIXTH
186 
187 
188 
189 
190  case SEVENTH:
191  {
192  // A degree 7, 12-point rule due to Tyler, can be found in Stroud's book
193  //
194  // A.H. Stroud, Approximate calculation of multiple integrals,
195  // Prentice-Hall, Englewood Cliffs, N.J., 1971.
196  //
197  // This rule is fully-symmetric and provably minimal in the number of points.
198  // A tensor-product rule accurate for "bi-septic" polynomials would have 16 points.
199  const Real
200  r = std::sqrt(6.L/7.L),
201  s = std::sqrt( (114.L - 3.L*std::sqrt(583.L)) / 287.L ),
202  t = std::sqrt( (114.L + 3.L*std::sqrt(583.L)) / 287.L ),
203  B1 = 196.L / 810.L,
204  B2 = 4.L * (178981.L + 2769.L*std::sqrt(583.L)) / 1888920.L,
205  B3 = 4.L * (178981.L - 2769.L*std::sqrt(583.L)) / 1888920.L;
206 
207  const Real data[3][3] =
208  {
209  {r, 0.0, B1}, // 4
210  {s, 0.0, B2}, // 4
211  {t, 0.0, B3} // 4
212  };
213 
214  const unsigned int symmetry[3] = {
215  3, // Full Symmetry, (x,0)
216  2, // Full Symmetry, (x,x)
217  2 // Full Symmetry, (x,x)
218  };
219 
220  _points.resize (12);
221  _weights.resize(12);
222 
223  stroud_rule(data, symmetry, 3);
224 
225  return;
226  } // end case SEVENTH
227 
228 
229 
230 
231  case EIGHTH:
232  {
233  // A pair of degree=8 rules for the QUAD "C2" due to
234  // Wissmann and Becker. These rules both have 16 points.
235  // A tensor product degree-6 Gauss would have 25 points.
236  //
237  // J. W. Wissmann and T. Becker, Partially symmetric cubature
238  // formulas for even degrees of exactness, SIAM J. Numer. Anal. 23
239  // (1986), 676--685.
240  const Real data[10][3] =
241  {
242  // First of 2 degree-8, 16 point rules given by Wissmann
243  // {0.000000000000000, 0.000000000000000, 0.055364705621440},
244  // {0.000000000000000, 0.757629177660505, 0.404389368726076},
245  // {0.000000000000000, -0.236871842255702, 0.533546604952635},
246  // {0.000000000000000, -0.989717929044527, 0.117054188786739},
247  // {0.639091304900370, 0.950520955645667, 0.125614417613747},
248  // {0.937069076924990, 0.663882736885633, 0.136544584733588},
249  // {0.537083530541494, 0.304210681724104, 0.483408479211257},
250  // {0.887188506449625, -0.236496718536120, 0.252528506429544},
251  // {0.494698820670197, -0.698953476086564, 0.361262323882172},
252  // {0.897495818279768, -0.900390774211580, 0.085464254086247}
253  //
254  // Second of 2 degree-8, 16 point rules given by Wissmann.
255  // Either of these will work, I just chose the one with points
256  // further into the element interior.
257  {0.0000000000000000e+00, 6.5956013196034176e-01, 4.5027677630559029e-01},
258  {0.0000000000000000e+00, -9.4914292304312538e-01, 1.6657042677781274e-01},
259  {9.5250946607156228e-01, 7.6505181955768362e-01, 9.8869459933431422e-02},
260  {5.3232745407420624e-01, 9.3697598108841598e-01, 1.5369674714081197e-01},
261  {6.8473629795173504e-01, 3.3365671773574759e-01, 3.9668697607290278e-01},
262  {2.3314324080140552e-01, -7.9583272377396852e-02, 3.5201436794569501e-01},
263  {9.2768331930611748e-01, -2.7224008061253425e-01, 1.8958905457779799e-01},
264  {4.5312068740374942e-01, -6.1373535339802760e-01, 3.7510100114758727e-01},
265  {8.3750364042281223e-01, -8.8847765053597136e-01, 1.2561879164007201e-01}
266  };
267 
268  _points.resize(16);
269  _weights.resize(16);
270 
271  wissmann_rule(data, /*10*/ 9);
272 
273  return;
274  } // end case EIGHTH
275 
276 
277 
278 
279  case NINTH:
280  {
281  // A degree 9, 17-point rule due to Moller.
282  //
283  // H.M. Moller, Kubaturformeln mit minimaler Knotenzahl,
284  // Numer. Math. 25 (1976), 185--200.
285  //
286  // This rule is provably minimal in the number of points.
287  // A tensor-product rule accurate for "bi-ninth" degree polynomials would have 25 points.
288  const Real data[5][3] =
289  {
290  {0.0000000000000000e+00, 0.0000000000000000e+00, 5.2674897119341563e-01}, // 1
291  {6.3068011973166885e-01, 9.6884996636197772e-01, 8.8879378170198706e-02}, // 4
292  {9.2796164595956966e-01, 7.5027709997890053e-01, 1.1209960212959648e-01}, // 4
293  {4.5333982113564719e-01, 5.2373582021442933e-01, 3.9828243926207009e-01}, // 4
294  {8.5261572933366230e-01, 7.6208328192617173e-02, 2.6905133763978080e-01} // 4
295  };
296 
297  const unsigned int symmetry[5] = {
298  0, // Single point
299  4, // Rotational Invariant
300  4, // Rotational Invariant
301  4, // Rotational Invariant
302  4 // Rotational Invariant
303  };
304 
305  _points.resize (17);
306  _weights.resize(17);
307 
308  stroud_rule(data, symmetry, 5);
309 
310  return;
311  } // end case NINTH
312 
313 
314 
315 
316  case TENTH:
317  case ELEVENTH:
318  {
319  // A degree 11, 24-point rule due to Cools and Haegemans.
320  //
321  // R. Cools and A. Haegemans, Another step forward in searching for
322  // cubature formulae with a minimal number of knots for the square,
323  // Computing 40 (1988), 139--146.
324  //
325  // P. Verlinden and R. Cools, The algebraic construction of a minimal
326  // cubature formula of degree 11 for the square, Cubature Formulas
327  // and their Applications (Russian) (Krasnoyarsk) (M.V. Noskov, ed.),
328  // 1994, pp. 13--23.
329  //
330  // This rule is provably minimal in the number of points.
331  // A tensor-product rule accurate for "bi-tenth" or "bi-eleventh" degree polynomials would have 36 points.
332  const Real data[6][3] =
333  {
334  {6.9807610454956756e-01, 9.8263922354085547e-01, 4.8020763350723814e-02}, // 4
335  {9.3948638281673690e-01, 8.2577583590296393e-01, 6.6071329164550595e-02}, // 4
336  {9.5353952820153201e-01, 1.8858613871864195e-01, 9.7386777358668164e-02}, // 4
337  {3.1562343291525419e-01, 8.1252054830481310e-01, 2.1173634999894860e-01}, // 4
338  {7.1200191307533630e-01, 5.2532025036454776e-01, 2.2562606172886338e-01}, // 4
339  {4.2484724884866925e-01, 4.1658071912022368e-02, 3.5115871839824543e-01} // 4
340  };
341 
342  const unsigned int symmetry[6] = {
343  4, // Rotational Invariant
344  4, // Rotational Invariant
345  4, // Rotational Invariant
346  4, // Rotational Invariant
347  4, // Rotational Invariant
348  4 // Rotational Invariant
349  };
350 
351  _points.resize (24);
352  _weights.resize(24);
353 
354  stroud_rule(data, symmetry, 6);
355 
356  return;
357  } // end case TENTH,ELEVENTH
358 
359 
360 
361 
362  case TWELFTH:
363  case THIRTEENTH:
364  {
365  // A degree 13, 33-point rule due to Cools and Haegemans.
366  //
367  // R. Cools and A. Haegemans, Another step forward in searching for
368  // cubature formulae with a minimal number of knots for the square,
369  // Computing 40 (1988), 139--146.
370  //
371  // A tensor-product rule accurate for "bi-12" or "bi-13" degree polynomials would have 49 points.
372  const Real data[9][3] =
373  {
374  {0.0000000000000000e+00, 0.0000000000000000e+00, 3.0038211543122536e-01}, // 1
375  {9.8348668243987226e-01, 7.7880971155441942e-01, 2.9991838864499131e-02}, // 4
376  {8.5955600564163892e-01, 9.5729769978630736e-01, 3.8174421317083669e-02}, // 4
377  {9.5892517028753485e-01, 1.3818345986246535e-01, 6.0424923817749980e-02}, // 4
378  {3.9073621612946100e-01, 9.4132722587292523e-01, 7.7492738533105339e-02}, // 4
379  {8.5007667369974857e-01, 4.7580862521827590e-01, 1.1884466730059560e-01}, // 4
380  {6.4782163718701073e-01, 7.5580535657208143e-01, 1.2976355037000271e-01}, // 4
381  {7.0741508996444936e-02, 6.9625007849174941e-01, 2.1334158145718938e-01}, // 4
382  {4.0930456169403884e-01, 3.4271655604040678e-01, 2.5687074948196783e-01} // 4
383  };
384 
385  const unsigned int symmetry[9] = {
386  0, // Single point
387  4, // Rotational Invariant
388  4, // Rotational Invariant
389  4, // Rotational Invariant
390  4, // Rotational Invariant
391  4, // Rotational Invariant
392  4, // Rotational Invariant
393  4, // Rotational Invariant
394  4 // Rotational Invariant
395  };
396 
397  _points.resize (33);
398  _weights.resize(33);
399 
400  stroud_rule(data, symmetry, 9);
401 
402  return;
403  } // end case TWELFTH,THIRTEENTH
404 
405 
406 
407 
408  case FOURTEENTH:
409  case FIFTEENTH:
410  {
411  // A degree-15, 48 point rule originally due to Rabinowitz and Richter,
412  // can be found in Cools' 1971 book.
413  //
414  // A.H. Stroud, Approximate calculation of multiple integrals,
415  // Prentice-Hall, Englewood Cliffs, N.J., 1971.
416  //
417  // The product Gauss rule for this order has 8^2=64 points.
418  const Real data[9][3] =
419  {
420  {9.915377816777667e-01L, 0.0000000000000000e+00, 3.01245207981210e-02L}, // 4
421  {8.020163879230440e-01L, 0.0000000000000000e+00, 8.71146840209092e-02L}, // 4
422  {5.648674875232742e-01L, 0.0000000000000000e+00, 1.250080294351494e-01L}, // 4
423  {9.354392392539896e-01L, 0.0000000000000000e+00, 2.67651407861666e-02L}, // 4
424  {7.624563338825799e-01L, 0.0000000000000000e+00, 9.59651863624437e-02L}, // 4
425  {2.156164241427213e-01L, 0.0000000000000000e+00, 1.750832998343375e-01L}, // 4
426  {9.769662659711761e-01L, 6.684480048977932e-01L, 2.83136372033274e-02L}, // 4
427  {8.937128379503403e-01L, 3.735205277617582e-01L, 8.66414716025093e-02L}, // 4
428  {6.122485619312083e-01L, 4.078983303613935e-01L, 1.150144605755996e-01L} // 4
429  };
430 
431  const unsigned int symmetry[9] = {
432  3, // Full Symmetry, (x,0)
433  3, // Full Symmetry, (x,0)
434  3, // Full Symmetry, (x,0)
435  2, // Full Symmetry, (x,x)
436  2, // Full Symmetry, (x,x)
437  2, // Full Symmetry, (x,x)
438  1, // Full Symmetry, (x,y)
439  1, // Full Symmetry, (x,y)
440  1, // Full Symmetry, (x,y)
441  };
442 
443  _points.resize (48);
444  _weights.resize(48);
445 
446  stroud_rule(data, symmetry, 9);
447 
448  return;
449  } // case FOURTEENTH, FIFTEENTH:
450 
451 
452 
453 
454  case SIXTEENTH:
455  case SEVENTEENTH:
456  {
457  // A degree 17, 60-point rule due to Cools and Haegemans.
458  //
459  // R. Cools and A. Haegemans, Another step forward in searching for
460  // cubature formulae with a minimal number of knots for the square,
461  // Computing 40 (1988), 139--146.
462  //
463  // A tensor-product rule accurate for "bi-14" or "bi-15" degree polynomials would have 64 points.
464  // A tensor-product rule accurate for "bi-16" or "bi-17" degree polynomials would have 81 points.
465  const Real data[10][3] =
466  {
467  {9.8935307451260049e-01, 0.0000000000000000e+00, 2.0614915919990959e-02}, // 4
468  {3.7628520715797329e-01, 0.0000000000000000e+00, 1.2802571617990983e-01}, // 4
469  {9.7884827926223311e-01, 0.0000000000000000e+00, 5.5117395340318905e-03}, // 4
470  {8.8579472916411612e-01, 0.0000000000000000e+00, 3.9207712457141880e-02}, // 4
471  {1.7175612383834817e-01, 0.0000000000000000e+00, 7.6396945079863302e-02}, // 4
472  {5.9049927380600241e-01, 3.1950503663457394e-01, 1.4151372994997245e-01}, // 8
473  {7.9907913191686325e-01, 5.9797245192945738e-01, 8.3903279363797602e-02}, // 8
474  {8.0374396295874471e-01, 5.8344481776550529e-02, 6.0394163649684546e-02}, // 8
475  {9.3650627612749478e-01, 3.4738631616620267e-01, 5.7387752969212695e-02}, // 8
476  {9.8132117980545229e-01, 7.0600028779864611e-01, 2.1922559481863763e-02}, // 8
477  };
478 
479  const unsigned int symmetry[10] = {
480  3, // Fully symmetric (x,0)
481  3, // Fully symmetric (x,0)
482  2, // Fully symmetric (x,x)
483  2, // Fully symmetric (x,x)
484  2, // Fully symmetric (x,x)
485  1, // Fully symmetric (x,y)
486  1, // Fully symmetric (x,y)
487  1, // Fully symmetric (x,y)
488  1, // Fully symmetric (x,y)
489  1 // Fully symmetric (x,y)
490  };
491 
492  _points.resize (60);
493  _weights.resize(60);
494 
495  stroud_rule(data, symmetry, 10);
496 
497  return;
498  } // end case FOURTEENTH through SEVENTEENTH
499 
500 
501 
502  // By default: construct and use a Gauss quadrature rule
503  default:
504  {
505  // Break out and fall down into the default: case for the
506  // outer switch statement.
507  break;
508  }
509 
510  } // end switch(_order + 2*p)
511  } // end case QUAD4/8/9
512 
513 
514  // By default: construct and use a Gauss quadrature rule
515  default:
516  {
517  QGauss gauss_rule(2, _order);
518  gauss_rule.init(type_in, p);
519 
520  // Swap points and weights with the about-to-be destroyed rule.
521  _points.swap (gauss_rule.get_points() );
522  _weights.swap(gauss_rule.get_weights());
523 
524  return;
525  }
526  } // end switch (type_in)
527 }
void libMesh::QMonomial::init_3D ( const ElemType  _type = INVALID_ELEM,
unsigned int  p_level = 0 
)
private

More efficient rules for HEXes

Definition at line 28 of file quadrature_monomial_3D.C.

References libMesh::QBase::_points, libMesh::QBase::_weights, A, libMesh::QBase::allow_rules_with_negative_weights, data, libMeshEnums::EIGHTH, libMeshEnums::FIFTH, libMeshEnums::FOURTH, libMesh::QBase::get_points(), libMesh::QBase::get_weights(), libMeshEnums::HEX20, libMeshEnums::HEX27, libMeshEnums::HEX8, libMesh::QBase::init(), kim_rule(), libMesh::Real, libMeshEnums::SECOND, libMeshEnums::SEVENTH, libMeshEnums::SIXTH, and libMeshEnums::THIRD.

30 {
31 
32  switch (type_in)
33  {
34  //---------------------------------------------
35  // Hex quadrature rules
36  case HEX8:
37  case HEX20:
38  case HEX27:
39  {
40  switch(_order + 2*p)
41  {
42 
43  // The CONSTANT/FIRST rule is the 1-point Gauss "product" rule...we fall
44  // through to the default case for this rule.
45 
46  case SECOND:
47  case THIRD:
48  {
49  // A degree 3, 6-point, "rotationally-symmetric" rule by
50  // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
51  //
52  // Warning: this rule contains points on the boundary of the reference
53  // element, and therefore may be unsuitable for some problems. The alternative
54  // would be a 2x2x2 Gauss product rule.
55  const Real data[1][4] =
56  {
57  {1.0L, 0.0L, 0.0L, static_cast<Real>(4.0L/3.0L)}
58  };
59 
60  const unsigned int rule_id[1] = {
61  1 // (x,0,0) -> 6 permutations
62  };
63 
64  _points.resize(6);
65  _weights.resize(6);
66 
67  kim_rule(data, rule_id, 1);
68  return;
69  } // end case SECOND,THIRD
70 
71  case FOURTH:
72  case FIFTH:
73  {
74  // A degree 5, 13-point rule by Stroud,
75  // AH Stroud, "Some Fifth Degree Integration Formulas for Symmetric Regions II.",
76  // Numerische Mathematik 9, pp. 460-468 (1967).
77  //
78  // This rule is provably minimal in the number of points. The equations given for
79  // the n-cube on pg. 466 of the paper for mu/gamma and gamma are wrong, at least for
80  // the n=3 case. The analytical values given here were computed by me [JWP] in Maple.
81 
82  // Convenient intermediate values.
83  const Real sqrt19 = std::sqrt(19.L);
84  const Real tp = std::sqrt(71440.L + 6802.L*sqrt19);
85 
86  // Point data for permutations.
87  const Real eta = 0.00000000000000000000000000000000e+00L;
88 
89  const Real lambda = std::sqrt(1919.L/3285.L - 148.L*sqrt19/3285.L + 4.L*tp/3285.L);
90  // 8.8030440669930978047737818209860e-01L;
91 
92  const Real xi = -std::sqrt(1121.L/3285.L + 74.L*sqrt19/3285.L - 2.L*tp/3285.L);
93  // -4.9584817142571115281421242364290e-01L;
94 
95  const Real mu = std::sqrt(1121.L/3285.L + 74.L*sqrt19/3285.L + 2.L*tp/3285.L);
96  // 7.9562142216409541542982482567580e-01L;
97 
98  const Real gamma = std::sqrt(1919.L/3285.L - 148.L*sqrt19/3285.L - 4.L*tp/3285.L);
99  // 2.5293711744842581347389255929324e-02L;
100 
101  // Weights: the centroid weight is given analytically. Weight B (resp C) goes
102  // with the {lambda,xi} (resp {gamma,mu}) permutation. The single-precision
103  // results reported by Stroud are given for reference.
104 
105  const Real A = 32.0L / 19.0L;
106  // Stroud: 0.21052632 * 8.0 = 1.684210560;
107 
108  const Real B = 1.L / ( 260072.L/133225.L - 1520*sqrt19/133225.L + (133.L - 37.L*sqrt19)*tp/133225.L );
109  // 5.4498735127757671684690782180890e-01L; // Stroud: 0.068123420 * 8.0 = 0.544987360;
110 
111  const Real C = 1.L / ( 260072.L/133225.L - 1520*sqrt19/133225.L - (133.L - 37.L*sqrt19)*tp/133225.L );
112  // 5.0764422766979170420572375713840e-01L; // Stroud: 0.063455527 * 8.0 = 0.507644216;
113 
114  _points.resize(13);
115  _weights.resize(13);
116 
117  unsigned int c=0;
118 
119  // Point with weight A (origin)
120  _points[c] = Point(eta, eta, eta);
121  _weights[c++] = A;
122 
123  // Points with weight B
124  _points[c] = Point(lambda, xi, xi);
125  _weights[c++] = B;
126  _points[c] = -_points[c-1];
127  _weights[c++] = B;
128 
129  _points[c] = Point(xi, lambda, xi);
130  _weights[c++] = B;
131  _points[c] = -_points[c-1];
132  _weights[c++] = B;
133 
134  _points[c] = Point(xi, xi, lambda);
135  _weights[c++] = B;
136  _points[c] = -_points[c-1];
137  _weights[c++] = B;
138 
139  // Points with weight C
140  _points[c] = Point(mu, mu, gamma);
141  _weights[c++] = C;
142  _points[c] = -_points[c-1];
143  _weights[c++] = C;
144 
145  _points[c] = Point(mu, gamma, mu);
146  _weights[c++] = C;
147  _points[c] = -_points[c-1];
148  _weights[c++] = C;
149 
150  _points[c] = Point(gamma, mu, mu);
151  _weights[c++] = C;
152  _points[c] = -_points[c-1];
153  _weights[c++] = C;
154 
155  return;
156 
157 
158 // // A degree 5, 14-point, "rotationally-symmetric" rule by
159 // // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
160 // // Was also reported in Stroud's 1971 book.
161 // const Real data[2][4] =
162 // {
163 // {7.95822425754221463264548820476135e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 8.86426592797783933518005540166204e-01L},
164 // {7.58786910639328146269034278112267e-01L, 7.58786910639328146269034278112267e-01L, 7.58786910639328146269034278112267e-01L, 3.35180055401662049861495844875346e-01L}
165 // };
166 
167 // const unsigned int rule_id[2] = {
168 // 1, // (x,0,0) -> 6 permutations
169 // 4 // (x,x,x) -> 8 permutations
170 // };
171 
172 // _points.resize(14);
173 // _weights.resize(14);
174 
175 // kim_rule(data, rule_id, 2);
176 // return;
177  } // end case FOURTH,FIFTH
178 
179  case SIXTH:
180  case SEVENTH:
181  {
183  {
184  // A degree 7, 31-point, "rotationally-symmetric" rule by
185  // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
186  // This rule contains a negative weight, so only use it if such type of
187  // rules are allowed.
188  const Real data[3][4] =
189  {
190  {0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, -1.27536231884057971014492753623188e+00L},
191  {5.85540043769119907612630781744060e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 8.71111111111111111111111111111111e-01L},
192  {6.94470135991704766602025803883310e-01L, 9.37161638568208038511047377665396e-01L, 4.15659267604065126239606672567031e-01L, 1.68695652173913043478260869565217e-01L}
193  };
194 
195  const unsigned int rule_id[3] = {
196  0, // (0,0,0) -> 1 permutation
197  1, // (x,0,0) -> 6 permutations
198  6 // (x,y,z) -> 24 permutations
199  };
200 
201  _points.resize(31);
202  _weights.resize(31);
203 
204  kim_rule(data, rule_id, 3);
205  return;
206  } // end if (allow_rules_with_negative_weights)
207 
208 
209  // A degree 7, 34-point, "fully-symmetric" rule, first published in
210  // P.C. Hammer and A.H. Stroud, "Numerical Evaluation of Multiple Integrals II",
211  // Mathmatical Tables and Other Aids to Computation, vol 12., no 64, 1958, pp. 272-280
212  //
213  // This rule happens to fall under the same general
214  // construction as the Kim rules, so we've re-used
215  // that code here. Stroud gives 16 digits for his rule,
216  // and this is the most accurate version I've found.
217  //
218  // For comparison, a SEVENTH-order Gauss product rule
219  // (which integrates tri-7th order polynomials) would
220  // have 4^3=64 points.
221  const Real
222  r = std::sqrt(6.L/7.L),
223  s = std::sqrt((960.L - 3.L*std::sqrt(28798.L)) / 2726.L),
224  t = std::sqrt((960.L + 3.L*std::sqrt(28798.L)) / 2726.L),
225  B1 = 8624.L / 29160.L,
226  B2 = 2744.L / 29160.L,
227  B3 = 8.L*(774.L*t*t - 230.L)/(9720.L*(t*t-s*s)),
228  B4 = 8.L*(230.L - 774.L*s*s)/(9720.L*(t*t-s*s));
229 
230  const Real data[4][4] =
231  {
232  {r, 0.L, 0.L, B1},
233  {r, r, 0.L, B2},
234  {s, s, s, B3},
235  {t, t, t, B4}
236  };
237 
238  const unsigned int rule_id[4] = {
239  1, // (x,0,0) -> 6 permutations
240  2, // (x,x,0) -> 12 permutations
241  4, // (x,x,x) -> 8 permutations
242  4 // (x,x,x) -> 8 permutations
243  };
244 
245  _points.resize(34);
246  _weights.resize(34);
247 
248  kim_rule(data, rule_id, 4);
249  return;
250 
251 
252 // // A degree 7, 38-point, "rotationally-symmetric" rule by
253 // // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
254 // //
255 // // This rule is obviously inferior to the 34-point rule above...
256 // const Real data[3][4] =
257 // {
258 // {9.01687807821291289082811566285950e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 2.95189738262622903181631100062774e-01L},
259 // {4.08372221499474674069588900002128e-01L, 4.08372221499474674069588900002128e-01L, 4.08372221499474674069588900002128e-01L, 4.04055417266200582425904380777126e-01L},
260 // {8.59523090201054193116477875786220e-01L, 8.59523090201054193116477875786220e-01L, 4.14735913727987720499709244748633e-01L, 1.24850759678944080062624098058597e-01L}
261 // };
262 //
263 // const unsigned int rule_id[3] = {
264 // 1, // (x,0,0) -> 6 permutations
265 // 4, // (x,x,x) -> 8 permutations
266 // 5 // (x,x,z) -> 24 permutations
267 // };
268 //
269 // _points.resize(38);
270 // _weights.resize(38);
271 //
272 // kim_rule(data, rule_id, 3);
273 // return;
274  } // end case SIXTH,SEVENTH
275 
276  case EIGHTH:
277  {
278  // A degree 8, 47-point, "rotationally-symmetric" rule by
279  // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
280  //
281  // A EIGHTH-order Gauss product rule (which integrates tri-8th order polynomials)
282  // would have 5^3=125 points.
283  const Real data[5][4] =
284  {
285  {0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 4.51903714875199690490763818699555e-01L},
286  {7.82460796435951590652813975429717e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 2.99379177352338919703385618576171e-01L},
287  {4.88094669706366480526729301468686e-01L, 4.88094669706366480526729301468686e-01L, 4.88094669706366480526729301468686e-01L, 3.00876159371240019939698689791164e-01L},
288  {8.62218927661481188856422891110042e-01L, 8.62218927661481188856422891110042e-01L, 8.62218927661481188856422891110042e-01L, 4.94843255877038125738173175714853e-02L},
289  {2.81113909408341856058098281846420e-01L, 9.44196578292008195318687494773744e-01L, 6.97574833707236996779391729948984e-01L, 1.22872389222467338799199767122592e-01L}
290  };
291 
292  const unsigned int rule_id[5] = {
293  0, // (0,0,0) -> 1 permutation
294  1, // (x,0,0) -> 6 permutations
295  4, // (x,x,x) -> 8 permutations
296  4, // (x,x,x) -> 8 permutations
297  6 // (x,y,z) -> 24 permutations
298  };
299 
300  _points.resize(47);
301  _weights.resize(47);
302 
303  kim_rule(data, rule_id, 5);
304  return;
305  } // end case EIGHTH
306 
307 
308  // By default: construct and use a Gauss quadrature rule
309  default:
310  {
311  // Break out and fall down into the default: case for the
312  // outer switch statement.
313  break;
314  }
315 
316  } // end switch(_order + 2*p)
317  } // end case HEX8/20/27
318 
319 
320  // By default: construct and use a Gauss quadrature rule
321  default:
322  {
323  QGauss gauss_rule(3, _order);
324  gauss_rule.init(type_in, p);
325 
326  // Swap points and weights with the about-to-be destroyed rule.
327  _points.swap (gauss_rule.get_points() );
328  _weights.swap(gauss_rule.get_weights());
329 
330  return;
331  }
332  } // end switch (type_in)
333 }
void libMesh::QMonomial::kim_rule ( const Real  rule_data[][4],
const unsigned int *  rule_id,
const unsigned int  n_pts 
)
private

Rules from Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931. The rules are obtained by considering the group G^{rot} of rotations of the reference hex, and the invariant polynomials of this group.

In Kim and Song's rules, quadrauture points are described by the following points and their unique permutations under the G^{rot} group:

0.) (0,0,0) ( 1 perm ) -> [0, 0, 0] 1.) (x,0,0) ( 6 perms) -> [x, 0, 0], [0, -x, 0], [-x, 0, 0], [0, x, 0], [0, 0, -x], [0, 0, x] 2.) (x,x,0) (12 perms) -> [x, x, 0], [x, -x, 0], [-x, -x, 0], [-x, x, 0], [x, 0, -x], [x, 0, x], [0, x, -x], [0, x, x], [0, -x, -x], [-x, 0, -x], [0, -x, x], [-x, 0, x] 3.) (x,y,0) (24 perms) -> [x, y, 0], [y, -x, 0], [-x, -y, 0], [-y, x, 0], [x, 0, -y], [x, -y, 0], [x, 0, y], [0, y, -x], [-x, y, 0], [0, y, x], [y, 0, -x], [0, -y, -x], [-y, 0, -x], [y, x, 0], [-y, -x, 0], [y, 0, x], [0, -y, x], [-y, 0, x], [-x, 0, y], [0, -x, -y], [0, -x, y], [-x, 0, -y], [0, x, y], [0, x, -y] 4.) (x,x,x) ( 8 perms) -> [x, x, x], [x, -x, x], [-x, -x, x], [-x, x, x], [x, x, -x], [x, -x, -x], [-x, x, -x], [-x, -x, -x] 5.) (x,x,z) (24 perms) -> [x, x, z], [x, -x, z], [-x, -x, z], [-x, x, z], [x, z, -x], [x, -x, -z], [x, -z, x], [z, x, -x], [-x, x, -z], [-z, x, x], [x, -z, -x], [-z, -x, -x], [-x, z, -x], [x, x, -z], [-x, -x, -z], [x, z, x], [z, -x, x], [-x, -z, x], [-x, z, x], [z, -x, -x], [-z, -x, x], [-x, -z, -x], [z, x, x], [-z, x, -x] 6.) (x,y,z) (24 perms) -> [x, y, z], [y, -x, z], [-x, -y, z], [-y, x, z], [x, z, -y], [x, -y, -z], [x, -z, y], [z, y, -x], [-x, y, -z], [-z, y, x], [y, -z, -x], [-z, -y, -x], [-y, z, -x], [y, x, -z], [-y, -x, -z], [y, z, x], [z, -y, x], [-y, -z, x], [-x, z, y], [z, -x, -y], [-z, -x, y], [-x, -z, -y], [z, x, y], [-z, x, -y]

Only two of Kim and Song's rules are particularly useful for FEM calculations: the degree 7, 38-point rule and their degree 8, 47-point rule. The others either contain negative weights or points outside the reference interval. The points and weights, to 32 digits, were obtained from: Ronald Cools' website (http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html) and the unique permutations of G^{rot} were computed by me [JWP] using Maple.

Definition at line 215 of file quadrature_monomial.C.

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

Referenced by init_3D().

218 {
219  for (unsigned int i=0, c=0; i<n_pts; ++i)
220  {
221  const Real
222  x=rule_data[i][0],
223  y=rule_data[i][1],
224  z=rule_data[i][2],
225  wt=rule_data[i][3];
226 
227  switch(rule_id[i])
228  {
229  case 0: // (0,0,0) 1 permutation
230  {
231  _points[c] = Point( x, y, z); _weights[c++] = wt;
232 
233  break;
234  }
235  case 1: // (x,0,0) 6 permutations
236  {
237  _points[c] = Point( x, 0., 0.); _weights[c++] = wt;
238  _points[c] = Point(0., -x, 0.); _weights[c++] = wt;
239  _points[c] = Point(-x, 0., 0.); _weights[c++] = wt;
240  _points[c] = Point(0., x, 0.); _weights[c++] = wt;
241  _points[c] = Point(0., 0., -x); _weights[c++] = wt;
242  _points[c] = Point(0., 0., x); _weights[c++] = wt;
243 
244  break;
245  }
246  case 2: // (x,x,0) 12 permutations
247  {
248  _points[c] = Point( x, x, 0.); _weights[c++] = wt;
249  _points[c] = Point( x, -x, 0.); _weights[c++] = wt;
250  _points[c] = Point(-x, -x, 0.); _weights[c++] = wt;
251  _points[c] = Point(-x, x, 0.); _weights[c++] = wt;
252  _points[c] = Point( x, 0., -x); _weights[c++] = wt;
253  _points[c] = Point( x, 0., x); _weights[c++] = wt;
254  _points[c] = Point(0., x, -x); _weights[c++] = wt;
255  _points[c] = Point(0., x, x); _weights[c++] = wt;
256  _points[c] = Point(0., -x, -x); _weights[c++] = wt;
257  _points[c] = Point(-x, 0., -x); _weights[c++] = wt;
258  _points[c] = Point(0., -x, x); _weights[c++] = wt;
259  _points[c] = Point(-x, 0., x); _weights[c++] = wt;
260 
261  break;
262  }
263  case 3: // (x,y,0) 24 permutations
264  {
265  _points[c] = Point( x, y, 0.); _weights[c++] = wt;
266  _points[c] = Point( y, -x, 0.); _weights[c++] = wt;
267  _points[c] = Point(-x, -y, 0.); _weights[c++] = wt;
268  _points[c] = Point(-y, x, 0.); _weights[c++] = wt;
269  _points[c] = Point( x, 0., -y); _weights[c++] = wt;
270  _points[c] = Point( x, -y, 0.); _weights[c++] = wt;
271  _points[c] = Point( x, 0., y); _weights[c++] = wt;
272  _points[c] = Point(0., y, -x); _weights[c++] = wt;
273  _points[c] = Point(-x, y, 0.); _weights[c++] = wt;
274  _points[c] = Point(0., y, x); _weights[c++] = wt;
275  _points[c] = Point( y, 0., -x); _weights[c++] = wt;
276  _points[c] = Point(0., -y, -x); _weights[c++] = wt;
277  _points[c] = Point(-y, 0., -x); _weights[c++] = wt;
278  _points[c] = Point( y, x, 0.); _weights[c++] = wt;
279  _points[c] = Point(-y, -x, 0.); _weights[c++] = wt;
280  _points[c] = Point( y, 0., x); _weights[c++] = wt;
281  _points[c] = Point(0., -y, x); _weights[c++] = wt;
282  _points[c] = Point(-y, 0., x); _weights[c++] = wt;
283  _points[c] = Point(-x, 0., y); _weights[c++] = wt;
284  _points[c] = Point(0., -x, -y); _weights[c++] = wt;
285  _points[c] = Point(0., -x, y); _weights[c++] = wt;
286  _points[c] = Point(-x, 0., -y); _weights[c++] = wt;
287  _points[c] = Point(0., x, y); _weights[c++] = wt;
288  _points[c] = Point(0., x, -y); _weights[c++] = wt;
289 
290  break;
291  }
292  case 4: // (x,x,x) 8 permutations
293  {
294  _points[c] = Point( x, x, x); _weights[c++] = wt;
295  _points[c] = Point( x, -x, x); _weights[c++] = wt;
296  _points[c] = Point(-x, -x, x); _weights[c++] = wt;
297  _points[c] = Point(-x, x, x); _weights[c++] = wt;
298  _points[c] = Point( x, x, -x); _weights[c++] = wt;
299  _points[c] = Point( x, -x, -x); _weights[c++] = wt;
300  _points[c] = Point(-x, x, -x); _weights[c++] = wt;
301  _points[c] = Point(-x, -x, -x); _weights[c++] = wt;
302 
303  break;
304  }
305  case 5: // (x,x,z) 24 permutations
306  {
307  _points[c] = Point( x, x, z); _weights[c++] = wt;
308  _points[c] = Point( x, -x, z); _weights[c++] = wt;
309  _points[c] = Point(-x, -x, z); _weights[c++] = wt;
310  _points[c] = Point(-x, x, z); _weights[c++] = wt;
311  _points[c] = Point( x, z, -x); _weights[c++] = wt;
312  _points[c] = Point( x, -x, -z); _weights[c++] = wt;
313  _points[c] = Point( x, -z, x); _weights[c++] = wt;
314  _points[c] = Point( z, x, -x); _weights[c++] = wt;
315  _points[c] = Point(-x, x, -z); _weights[c++] = wt;
316  _points[c] = Point(-z, x, x); _weights[c++] = wt;
317  _points[c] = Point( x, -z, -x); _weights[c++] = wt;
318  _points[c] = Point(-z, -x, -x); _weights[c++] = wt;
319  _points[c] = Point(-x, z, -x); _weights[c++] = wt;
320  _points[c] = Point( x, x, -z); _weights[c++] = wt;
321  _points[c] = Point(-x, -x, -z); _weights[c++] = wt;
322  _points[c] = Point( x, z, x); _weights[c++] = wt;
323  _points[c] = Point( z, -x, x); _weights[c++] = wt;
324  _points[c] = Point(-x, -z, x); _weights[c++] = wt;
325  _points[c] = Point(-x, z, x); _weights[c++] = wt;
326  _points[c] = Point( z, -x, -x); _weights[c++] = wt;
327  _points[c] = Point(-z, -x, x); _weights[c++] = wt;
328  _points[c] = Point(-x, -z, -x); _weights[c++] = wt;
329  _points[c] = Point( z, x, x); _weights[c++] = wt;
330  _points[c] = Point(-z, x, -x); _weights[c++] = wt;
331 
332  break;
333  }
334  case 6: // (x,y,z) 24 permutations
335  {
336  _points[c] = Point( x, y, z); _weights[c++] = wt;
337  _points[c] = Point( y, -x, z); _weights[c++] = wt;
338  _points[c] = Point(-x, -y, z); _weights[c++] = wt;
339  _points[c] = Point(-y, x, z); _weights[c++] = wt;
340  _points[c] = Point( x, z, -y); _weights[c++] = wt;
341  _points[c] = Point( x, -y, -z); _weights[c++] = wt;
342  _points[c] = Point( x, -z, y); _weights[c++] = wt;
343  _points[c] = Point( z, y, -x); _weights[c++] = wt;
344  _points[c] = Point(-x, y, -z); _weights[c++] = wt;
345  _points[c] = Point(-z, y, x); _weights[c++] = wt;
346  _points[c] = Point( y, -z, -x); _weights[c++] = wt;
347  _points[c] = Point(-z, -y, -x); _weights[c++] = wt;
348  _points[c] = Point(-y, z, -x); _weights[c++] = wt;
349  _points[c] = Point( y, x, -z); _weights[c++] = wt;
350  _points[c] = Point(-y, -x, -z); _weights[c++] = wt;
351  _points[c] = Point( y, z, x); _weights[c++] = wt;
352  _points[c] = Point( z, -y, x); _weights[c++] = wt;
353  _points[c] = Point(-y, -z, x); _weights[c++] = wt;
354  _points[c] = Point(-x, z, y); _weights[c++] = wt;
355  _points[c] = Point( z, -x, -y); _weights[c++] = wt;
356  _points[c] = Point(-z, -x, y); _weights[c++] = wt;
357  _points[c] = Point(-x, -z, -y); _weights[c++] = wt;
358  _points[c] = Point( z, x, y); _weights[c++] = wt;
359  _points[c] = Point(-z, x, -y); _weights[c++] = wt;
360 
361  break;
362  }
363  default:
364  {
365  libMesh::err << "Unknown rule ID: " << rule_id[i] << "!" << std::endl;
366  libmesh_error();
367  }
368  } // end switch(rule_id[i])
369  }
370 }
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; }
void libMesh::QMonomial::stroud_rule ( const Real  rule_data[][3],
const unsigned int *  rule_symmetry,
const unsigned int  n_pts 
)
private

Stroud's rules for QUADs and HEXes can have one of several different types of symmetry. The rule_symmetry array describes how the different lines of the rule_data array are to be applied. The different rule_symmetry possibilities are: 0) Origin or single-point: (x,y) Fully-symmetric, 3 cases: 1) (x,y) -> (x,y), (-x,y), (x,-y), (-x,-y) (y,x), (-y,x), (y,-x), (-y,-x) 2) (x,x) -> (x,x), (-x,x), (x,-x), (-x,-x) 3) (x,0) -> (x,0), (-x,0), (0, x), ( 0,-x) 4) Rotational Invariant, (x,y) -> (x,y), (-x,-y), (-y, x), (y,-x) 5) Partial Symmetry, (x,y) -> (x,y), (-x, y) [x!=0] 6) Rectangular Symmetry, (x,y) -> (x,y), (-x, y), (-x,-y), (x,-y) 7) Central Symmetry, (0,y) -> (0,y), ( 0,-y)

Not all rules with these symmetries are due to Stroud, however, his book is probably the most frequently-cited compendium of quadrature rules and later authors certainly built upon his work.

Definition at line 64 of file quadrature_monomial.C.

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

Referenced by init_2D().

67 {
68  for (unsigned int i=0, c=0; i<n_pts; ++i)
69  {
70  const Real
71  x=rule_data[i][0],
72  y=rule_data[i][1],
73  wt=rule_data[i][2];
74 
75  switch(rule_symmetry[i])
76  {
77  case 0: // Single point (no symmetry)
78  {
79  _points[c] = Point( x, y);
80  _weights[c++] = wt;
81 
82  break;
83  }
84  case 1: // Fully-symmetric (x,y)
85  {
86  _points[c] = Point( x, y);
87  _weights[c++] = wt;
88 
89  _points[c] = Point(-x, y);
90  _weights[c++] = wt;
91 
92  _points[c] = Point( x,-y);
93  _weights[c++] = wt;
94 
95  _points[c] = Point(-x,-y);
96  _weights[c++] = wt;
97 
98  _points[c] = Point( y, x);
99  _weights[c++] = wt;
100 
101  _points[c] = Point(-y, x);
102  _weights[c++] = wt;
103 
104  _points[c] = Point( y,-x);
105  _weights[c++] = wt;
106 
107  _points[c] = Point(-y,-x);
108  _weights[c++] = wt;
109 
110  break;
111  }
112  case 2: // Fully-symmetric (x,x)
113  {
114  _points[c] = Point( x, x);
115  _weights[c++] = wt;
116 
117  _points[c] = Point(-x, x);
118  _weights[c++] = wt;
119 
120  _points[c] = Point( x,-x);
121  _weights[c++] = wt;
122 
123  _points[c] = Point(-x,-x);
124  _weights[c++] = wt;
125 
126  break;
127  }
128  case 3: // Fully-symmetric (x,0)
129  {
130  libmesh_assert_equal_to (y, 0.0);
131 
132  _points[c] = Point( x,0.);
133  _weights[c++] = wt;
134 
135  _points[c] = Point(-x,0.);
136  _weights[c++] = wt;
137 
138  _points[c] = Point(0., x);
139  _weights[c++] = wt;
140 
141  _points[c] = Point(0.,-x);
142  _weights[c++] = wt;
143 
144  break;
145  }
146  case 4: // Rotational invariant
147  {
148  _points[c] = Point( x, y);
149  _weights[c++] = wt;
150 
151  _points[c] = Point(-x,-y);
152  _weights[c++] = wt;
153 
154  _points[c] = Point(-y, x);
155  _weights[c++] = wt;
156 
157  _points[c] = Point( y,-x);
158  _weights[c++] = wt;
159 
160  break;
161  }
162  case 5: // Partial symmetry (Wissman's rules)
163  {
164  libmesh_assert_not_equal_to (x, 0.0);
165 
166  _points[c] = Point( x, y);
167  _weights[c++] = wt;
168 
169  _points[c] = Point(-x, y);
170  _weights[c++] = wt;
171 
172  break;
173  }
174  case 6: // Rectangular symmetry
175  {
176  _points[c] = Point( x, y);
177  _weights[c++] = wt;
178 
179  _points[c] = Point(-x, y);
180  _weights[c++] = wt;
181 
182  _points[c] = Point(-x,-y);
183  _weights[c++] = wt;
184 
185  _points[c] = Point( x,-y);
186  _weights[c++] = wt;
187 
188  break;
189  }
190  case 7: // Central symmetry
191  {
192  libmesh_assert_equal_to (x, 0.0);
193  libmesh_assert_not_equal_to (y, 0.0);
194 
195  _points[c] = Point(0., y);
196  _weights[c++] = wt;
197 
198  _points[c] = Point(0.,-y);
199  _weights[c++] = wt;
200 
201  break;
202  }
203  default:
204  {
205  libMesh::err << "Unknown symmetry!" << std::endl;
206  libmesh_error();
207  }
208  } // end switch(rule_symmetry[i])
209  }
210 }
QuadratureType libMesh::QMonomial::type ( ) const
inlinevirtual
Returns
QMONOMIAL

Implements libMesh::QBase.

Definition at line 83 of file quadrature_monomial.h.

References libMeshEnums::QMONOMIAL.

83 { return QMONOMIAL; }
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]; }
void libMesh::QMonomial::wissmann_rule ( const Real  rule_data[][3],
const unsigned int  n_pts 
)
private

Wissmann published three interesting "partially symmetric" rules for integrating degree 4, 6, and 8 polynomials exactly on QUADs. These rules have all positive weights, all points inside the reference element, and have fewer points than tensor-product rules of equivalent order, making them superior to those rules for monomial bases.

J. W. Wissman and T. Becker, Partially symmetric cubature formulas for even degrees of exactness, SIAM J. Numer. Anal. 23 (1986), 676–685.

Definition at line 44 of file quadrature_monomial.C.

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

Referenced by init_2D().

46 {
47  for (unsigned int i=0, c=0; i<n_pts; ++i)
48  {
49  _points[c] = Point( rule_data[i][0], rule_data[i][1] );
50  _weights[c++] = rule_data[i][2];
51 
52  // This may be an (x1,x2) -> (-x1,x2) point, in which case
53  // we will also generate the mirror point using the same weight.
54  if (rule_data[i][0] != static_cast<Real>(0.0))
55  {
56  _points[c] = Point( -rule_data[i][0], rule_data[i][1] );
57  _weights[c++] = rule_data[i][2];
58  }
59  }
60 }

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