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< QBase >	build (const std::string &name, const unsigned int _dim, const Order _order=INVALID_ORDER)
static AutoPtr< QBase >	build (const QuadratureType _qt, const unsigned int _dim, const Order _order=INVALID_ORDER)
static void	print_info (std::ostream &out=libMesh::out)
static std::string	get_info ()
static unsigned int	n_objects ()
static void	enable_print_counter_info ()
static void	disable_print_counter_info ()
Public Attributes
bool	allow_rules_with_negative_weights
Protected Types
typedef std::map< std::string, std::pair< unsigned int, unsigned int > >	Counts
Protected Member Functions
virtual void	init_0D (const ElemType type=INVALID_ELEM, unsigned int p_level=0)
void	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)
Friends
std::ostream &	operator<< (std::ostream &os, const QBase &q)

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 [protected, inherited]

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.

00034                                     : QBase(d,o)
00035 {
00036 }

libMesh::QMonomial::~QMonomial ( )

Destructor.

Definition at line 40 of file quadrature_monomial.C.

00041 {
00042 }

Member Function Documentation

AutoPtr< QBase > libMesh::QBase::build	(	const QuadratureType	_qt,
		const unsigned int	_dim,
		const Order	_order = `INVALID_ORDER`
	)			`[static, inherited]`

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.

00051 {
00052   switch (_qt)
00053     {
00054 
00055     case QCLOUGH:
00056       {
00057 #ifdef DEBUG
00058         if (_order > TWENTYTHIRD)
00059           {
00060             libMesh::out << "WARNING: Clough quadrature implemented" << std::endl
00061                           << " up to TWENTYTHIRD order." << std::endl;
00062           }
00063 #endif
00064 
00065         AutoPtr<QBase> ap(new QClough(_dim, _order));
00066         return ap;
00067       }
00068 
00069     case QGAUSS:
00070       {
00071 
00072 #ifdef DEBUG
00073         if (_order > FORTYTHIRD)
00074           {
00075             libMesh::out << "WARNING: Gauss quadrature implemented" << std::endl
00076                           << " up to FORTYTHIRD order." << std::endl;
00077           }
00078 #endif
00079 
00080         AutoPtr<QBase> ap(new QGauss(_dim, _order));
00081         return ap;
00082       }
00083 
00084     case QJACOBI_1_0:
00085       {
00086 
00087 #ifdef DEBUG
00088         if (_order > TWENTYTHIRD)
00089           {
00090             libMesh::out << "WARNING: Jacobi(1,0) quadrature implemented" << std::endl
00091                           << " up to TWENTYTHIRD order." << std::endl;
00092           }
00093 
00094         if (_dim > 1)
00095           {
00096             libMesh::out << "WARNING: Jacobi(1,0) quadrature implemented" << std::endl
00097                           << " in 1D only." << std::endl;
00098           }
00099 #endif
00100 
00101         AutoPtr<QBase> ap(new QJacobi(_dim, _order, 1, 0));
00102         return ap;
00103       }
00104 
00105     case QJACOBI_2_0:
00106       {
00107 
00108 #ifdef DEBUG
00109         if (_order > TWENTYTHIRD)
00110           {
00111             libMesh::out << "WARNING: Jacobi(2,0) quadrature implemented" << std::endl
00112                           << " up to TWENTYTHIRD order." << std::endl;
00113           }
00114 
00115         if (_dim > 1)
00116           {
00117             libMesh::out << "WARNING: Jacobi(2,0) quadrature implemented" << std::endl
00118                           << " in 1D only." << std::endl;
00119           }
00120 #endif
00121 
00122         AutoPtr<QBase> ap(new QJacobi(_dim, _order, 2, 0));
00123         return ap;
00124       }
00125 
00126     case QSIMPSON:
00127       {
00128 
00129 #ifdef DEBUG
00130         if (_order > THIRD)
00131           {
00132             libMesh::out << "WARNING: Simpson rule provides only" << std::endl
00133                           << " THIRD order!" << std::endl;
00134           }
00135 #endif
00136 
00137         AutoPtr<QBase> ap(new QSimpson(_dim));
00138         return ap;
00139       }
00140 
00141     case QTRAP:
00142       {
00143 
00144 #ifdef DEBUG
00145         if (_order > FIRST)
00146           {
00147             libMesh::out << "WARNING: Trapezoidal rule provides only" << std::endl
00148                           << " FIRST order!" << std::endl;
00149           }
00150 #endif
00151 
00152         AutoPtr<QBase> ap(new QTrap(_dim));
00153         return ap;
00154       }
00155 
00156 
00157     default:
00158       {
00159         libMesh::err << "ERROR: Bad qt=" << _qt << std::endl;
00160         libmesh_error();
00161       }
00162     }
00163 
00164 
00165   libmesh_error();
00166   AutoPtr<QBase> ap(NULL);
00167   return ap;
00168 }

AutoPtr< QBase > libMesh::QBase::build	(	const std::string &	name,
		const unsigned int	_dim,
		const Order	_order = `INVALID_ORDER`
	)			`[static, inherited]`

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.

References libMesh::Utility::string_to_enum< QuadratureType >().

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

00040 {
00041   return QBase::build (Utility::string_to_enum<QuadratureType> (type),
00042                        _dim,
00043                        _order);
00044 }

void libMesh::ReferenceCounter::disable_print_counter_info ( ) [static, inherited]

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

00107 {
00108   _enable_print_counter = false;
00109   return;
00110 }

void libMesh::ReferenceCounter::enable_print_counter_info ( ) [static, inherited]

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.

00101 {
00102   _enable_print_counter = true;
00103   return;
00104 }

unsigned int libMesh::QBase::get_dim ( ) const [inline, inherited]

Returns:: the dimension of the quadrature rule.

Definition at line 123 of file quadrature.h.

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

00123 { return _dim;  }

ElemType libMesh::QBase::get_elem_type ( ) const [inline, inherited]

Returns:: the current element type we're set up for

Definition at line 104 of file quadrature.h.

00105     { return _type; }

std::string libMesh::ReferenceCounter::get_info ( ) [static, inherited]

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

00048 {
00049 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
00050 
00051   std::ostringstream oss;
00052 
00053   oss << '\n'
00054       << " ---------------------------------------------------------------------------- \n"
00055       << "| Reference count information                                                |\n"
00056       << " ---------------------------------------------------------------------------- \n";
00057 
00058   for (Counts::iterator it = _counts.begin();
00059        it != _counts.end(); ++it)
00060     {
00061       const std::string name(it->first);
00062       const unsigned int creations    = it->second.first;
00063       const unsigned int destructions = it->second.second;
00064 
00065       oss << "| " << name << " reference count information:\n"
00066           << "|  Creations:    " << creations    << '\n'
00067           << "|  Destructions: " << destructions << '\n';
00068     }
00069 
00070   oss << " ---------------------------------------------------------------------------- \n";
00071 
00072   return oss.str();
00073 
00074 #else
00075 
00076   return "";
00077 
00078 #endif
00079 }

Order libMesh::QBase::get_order ( ) const [inline, inherited]

Returns:: the order of the quadrature rule.

Definition at line 169 of file quadrature.h.

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

00169 { return static_cast<Order>(_order + _p_level); }

unsigned int libMesh::QBase::get_p_level ( ) const [inline, inherited]

Returns:: the current p refinement level we're initialized with

Definition at line 110 of file quadrature.h.

00111     { return _p_level; }

std::vector<Point>& libMesh::QBase::get_points ( ) [inline, inherited]

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.

00135 { return _points;  }

const std::vector<Point>& libMesh::QBase::get_points ( ) const [inline, inherited]

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::FE< Dim, T >::edge_reinit(), libMesh::QClough::init_1D(), init_2D(), libMesh::QGauss::init_2D(), libMesh::QClough::init_2D(), init_3D(), libMesh::QGauss::init_3D(), libMesh::InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), libMesh::InfFE< Dim, T_radial, T_map >::init_shape_functions(), libMesh::InfFE< Dim, T_radial, T_map >::reinit(), libMesh::FEXYZ< Dim >::reinit(), and libMesh::FE< Dim, T >::reinit().

00129 { return _points;  }

std::vector<Real>& libMesh::QBase::get_weights ( ) [inline, inherited]

Returns:: a std::vector containing the quadrature weights.

Definition at line 145 of file quadrature.h.

References libMesh::QBase::_weights.

00145 { return _weights; }

const std::vector<Real>& libMesh::QBase::get_weights ( ) const [inline, inherited]

Returns:: a std::vector containing the quadrature weights.

Definition at line 140 of file quadrature.h.

References libMesh::QBase::_weights.

Referenced by libMesh::FE< Dim, T >::edge_reinit(), libMesh::QClough::init_1D(), init_2D(), libMesh::QGauss::init_2D(), libMesh::QClough::init_2D(), init_3D(), libMesh::QGauss::init_3D(), libMesh::InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), libMesh::InfFE< Dim, T_radial, T_map >::init_shape_functions(), libMesh::FEXYZ< Dim >::reinit(), and libMesh::FE< Dim, T >::reinit().

00140 { return _weights; }

void libMesh::ReferenceCounter::increment_constructor_count ( const std::string & name ) [inline, protected, inherited]

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, and libMesh::Threads::spin_mtx.

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

00164 {
00165   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
00166   std::pair<unsigned int, unsigned int>& p = _counts[name];
00167 
00168   p.first++;
00169 }

void libMesh::ReferenceCounter::increment_destructor_count ( const std::string & name ) [inline, protected, inherited]

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, and libMesh::Threads::spin_mtx.

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

00177 {
00178   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
00179   std::pair<unsigned int, unsigned int>& p = _counts[name];
00180 
00181   p.second++;
00182 }

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::FE< Dim, T >::edge_reinit(), libMesh::QClough::init_1D(), libMesh::QTrap::init_2D(), libMesh::QSimpson::init_2D(), init_2D(), libMesh::QGrid::init_2D(), libMesh::QGauss::init_2D(), libMesh::QClough::init_2D(), libMesh::QTrap::init_3D(), libMesh::QSimpson::init_3D(), init_3D(), libMesh::QGrid::init_3D(), libMesh::QGauss::init_3D(), libMesh::InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), libMesh::QGauss::QGauss(), libMesh::QJacobi::QJacobi(), libMesh::QSimpson::QSimpson(), libMesh::QTrap::QTrap(), libMesh::InfFE< Dim, T_radial, T_map >::reinit(), libMesh::FEXYZ< Dim >::reinit(), and libMesh::FE< Dim, T >::reinit().

00029 {
00030   // check to see if we have already
00031   // done the work for this quadrature rule
00032   if (t == _type && p == _p_level)
00033     return;
00034   else
00035     {
00036       _type = t;
00037       _p_level = p;
00038     }
00039 
00040 
00041 
00042   switch(_dim)
00043     {
00044     case 0:
00045       this->init_0D(_type,_p_level);
00046 
00047       return;
00048 
00049     case 1:
00050       this->init_1D(_type,_p_level);
00051 
00052       return;
00053 
00054     case 2:
00055       this->init_2D(_type,_p_level);
00056 
00057       return;
00058 
00059     case 3:
00060       this->init_3D(_type,_p_level);
00061 
00062       return;
00063 
00064     default:
00065       libmesh_error();
00066     }
00067 }

void libMesh::QBase::init_0D	(	const ElemType	type = `INVALID_ELEM`,
		unsigned int	p_level = `0`
	)			`[protected, virtual, inherited]`

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

00073 {
00074   _points.resize(1);
00075   _weights.resize(1);
00076   _points[0] = Point(0.);
00077   _weights[0] = 1.0;
00078 }

void libMesh::QMonomial::init_1D	(	const	type,
		unsigned	p_level = `0`
	)			`[inline, private, virtual]`

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.

00090   {
00091     // See about making this non-pure virtual in the base class?
00092     libmesh_error();
00093   }

void libMesh::QMonomial::init_2D	(	const ElemType	_type = `INVALID_ELEM`,
		unsigned int	p_level = `0`
	)			`[private, virtual]`

More efficient rules for QUADs

Reimplemented from libMesh::QBase.

Definition at line 28 of file quadrature_monomial_2D.C.

00031 00032 00033     { 00034 00035 00036 00037 00038 00039       { 00040 00041           { 00042 00043 00044 00045 00046 00047 00048 00049 00050 00051 00052 00053 00054 00055 00056 00057 00058 00059 00060 00061 00062 00063 00064 00065 00066 00067 00068 00069 00070 00071 00072 00073 00074 00075 00076 00077 00078 00079 00080 00081 00082 00083 00084 00085 00086 00087 00088 00089 00090 00091 00092 00093 00094 00095 00096 00097 00098 00099 00100 00101 00102 00103 00104 00105 00106 00107 00108 00109 00110 00111 00112 00113 00114 00115 00116 00117 00118 00119 00120 00121 00122 00123 00124 00125 00126 00127 00128 00129 00130 00131 00132 00133 00134 00135 00136 00137 00138 00139 00140 00141 00142 00143 00144 00145 00146 00147 00148 00149 00150 00151 00152 00153 00154 00155 00156 00157 00158 00159 00160 00161 00162 00163 00164 00165 00166 00167 00168 00169 00170 00171 00172 00173 00174 00175 00176 00177 00178 00179 00180 00181 00182 00183 00184 00185 00186 00187 00188 00189 00190 00191 00192 00193 00194 00195 00196 00197 00198 00199 00200 00201 00202 00203 00204 00205 00206 00207 00208 00209 00210 00211 00212 00213 00214 00215 00216 00217 00218 00219 00220 00221 00222 00223 00224 00225 00226 00227 00228 00229 00230 00231 00232 00233 00234 00235 00236 00237 00238 00239 00240 00241 00242 00243 00244 00245 00246 00247 00248 00249 00250 00251 00252 00253 00254 00255 00256 00257 00258 00259 00260 00261 00262 00263 00264 00265 00266 00267 00268 00269 00270 00271 00272 00273 00274 00275 00276 00277 00278 00279 00280 00281 00282 00283 00284 00285 00286 00287 00288 00289 00290 00291 00292 00293 00294 00295 00296 00297 00298 00299 00300 00301 00302 00303 00304 00305 00306 00307 00308 00309 00310 00311 00312 00313 00314 00315 00316 00317 00318 00319 00320 00321 00322 00323 00324 00325 00326 00327 00328 00329 00330 00331 00332 00333 00334 00335 00336 00337 00338 00339 00340 00341 00342 00343 00344 00345 00346 00347 00348 00349 00350 00351 00352 00353 00354 00355 00356 00357 00358 00359 00360 00361 00362 00363 00364 00365 00366 00367 00368 00369 00370 00371 00372 00373 00374 00375 00376 00377 00378 00379 00380 00381 00382 00383 00384 00385 00386 00387 00388 00389 00390 00391 00392 00393 00394 00395 00396 00397 00398 00399 00400 00401 00402 00403 00404 00405 00406 00407 00408 00409 00410 00411 00412 00413 00414 00415 00416 00417 00418 00419 00420 00421 00422 00423 00424 00425 00426 00427 00428 00429 00430 00431 00432 00433 00434 00435 00436 00437 00438 00439 00440 00441 00442 00443 00444 00445 00446 00447 00448 00449 00450 00451 00452 00453 00454 00455 00456 00457 00458 00459 00460 00461 00462 00463 00464 00465 00466 00467 00468 00469 00470 00471 00472 00473 00474 00475 00476 00477 00478 00479 00480 00481 00482 00483 00484 00485 00486 00487 00488 00489 00490 00491 00492 00493 00494 00495 00496 00497 00498 00499 00500 00501 00502 00503 00504 00505 00506 00507 00508 00509 00510 00511 00512 00513 00514 00515 00516       { 00517 00518 00519 00520 00521 00522 00523 00524 00525       } 00526 00527 }

class="fragment">00030 { switch (type_in) //--------------------------------------------- // Quadrilateral quadrature rules case QUAD4: case QUAD8: case QUAD9: switch(_order + 2*p) case SECOND: { // A degree=2 rule for the QUAD with 3 points. // A tensor product degree-2 Gauss would have 4 points. // This rule (or a variation on it) is probably available in // // A.H. Stroud, Approximate calculation of multiple integrals, // Prentice-Hall, Englewood Cliffs, N.J., 1971. // // though I have never actually seen a reference for it. // Luckily it's fairly easy to derive, which is what I've done // here [JWP]. const Real s=std::sqrt(1./3.), t=std::sqrt(2./3.); const Real data[2][3] = { {0.0, s, 2.0}, { t, -s, 1.0} }; _points.resize(3); _weights.resize(3); wissmann_rule(data, 2); return; } // end case SECOND // For third-order, fall through to default case, use 2x2 Gauss product rule. // case THIRD: // { // } // end case THIRD case FOURTH: { // A pair of degree=4 rules for the QUAD "C2" due to // Wissmann and Becker. These rules both have six points. // A tensor product degree-4 Gauss would have 9 points. // // J. W. Wissmann and T. Becker, Partially symmetric cubature // formulas for even degrees of exactness, SIAM J. Numer. Anal. 23 // (1986), 676--685. const Real data[4][3] = { // First of 2 degree-4 rules given by Wissmann {0.0000000000000000e+00, 0.0000000000000000e+00, 1.1428571428571428e+00}, {0.0000000000000000e+00, 9.6609178307929590e-01, 4.3956043956043956e-01}, {8.5191465330460049e-01, 4.5560372783619284e-01, 5.6607220700753210e-01}, {6.3091278897675402e-01, -7.3162995157313452e-01, 6.4271900178367668e-01} // // Second of 2 degree-4 rules given by Wissmann. These both // yield 4th-order accurate rules, I just chose the one that // happened to contain the origin. // {0.000000000000000, -0.356822089773090, 1.286412084888852}, // {0.000000000000000, 0.934172358962716, 0.491365692888926}, // {0.774596669241483, 0.390885162530071, 0.761883709085613}, // {0.774596669241483, -0.852765377881771, 0.349227402025498} }; _points.resize(6); _weights.resize(6); wissmann_rule(data, 4); return; } // end case FOURTH case FIFTH: { // A degree 5, 7-point rule due to Stroud. // // A.H. Stroud, Approximate calculation of multiple integrals, // Prentice-Hall, Englewood Cliffs, N.J., 1971. // // This rule is provably minimal in the number of points. // A tensor-product rule accurate for "bi-quintic" polynomials would have 9 points. const Real data[3][3] = { { 0.L, 0.L, 8.L / 7.L}, // 1 { 0.L, static_cast<Real>(std::sqrt(14.L/15.L)), 20.L / 63.L}, // 2 {static_cast<Real>(std::sqrt(3.L/5.L)), static_cast<Real>(std::sqrt(1.L/3.L)), 20.L / 36.L} // 4 }; const unsigned int symmetry[3] = { 0, // Origin 7, // Central Symmetry 6 // Rectangular }; _points.resize (7); _weights.resize(7); stroud_rule(data, symmetry, 3); return; } // end case FIFTH case SIXTH: { // A pair of degree=6 rules for the QUAD "C2" due to // Wissmann and Becker. These rules both have 10 points. // A tensor product degree-6 Gauss would have 16 points. // // J. W. Wissmann and T. Becker, Partially symmetric cubature // formulas for even degrees of exactness, SIAM J. Numer. Anal. 23 // (1986), 676--685. const Real data[6][3] = { // First of 2 degree-6, 10 point rules given by Wissmann // {0.000000000000000, 0.836405633697626, 0.455343245714174}, // {0.000000000000000, -0.357460165391307, 0.827395973202966}, // {0.888764014654765, 0.872101531193131, 0.144000884599645}, // {0.604857639464685, 0.305985162155427, 0.668259104262665}, // {0.955447506641064, -0.410270899466658, 0.225474004890679}, // {0.565459993438754, -0.872869311156879, 0.320896396788441} // // Second of 2 degree-6, 10 point rules given by Wissmann. // Either of these will work, I just chose the one with points // slightly further into the element interior. {0.0000000000000000e+00, 8.6983337525005900e-01, 3.9275059096434794e-01}, {0.0000000000000000e+00, -4.7940635161211124e-01, 7.5476288124261053e-01}, {8.6374282634615388e-01, 8.0283751620765670e-01, 2.0616605058827902e-01}, {5.1869052139258234e-01, 2.6214366550805818e-01, 6.8999213848986375e-01}, {9.3397254497284950e-01, -3.6309658314806653e-01, 2.6051748873231697e-01}, {6.0897753601635630e-01, -8.9660863276245265e-01, 2.6956758608606100e-01} }; _points.resize(10); _weights.resize(10); wissmann_rule(data, 6); return; } // end case SIXTH case SEVENTH: { // A degree 7, 12-point rule due to Tyler, can be found in Stroud's book // // A.H. Stroud, Approximate calculation of multiple integrals, // Prentice-Hall, Englewood Cliffs, N.J., 1971. // // This rule is fully-symmetric and provably minimal in the number of points. // A tensor-product rule accurate for "bi-septic" polynomials would have 16 points. const Real r = std::sqrt(6.L/7.L), s = std::sqrt( (114.L - 3.L*std::sqrt(583.L)) / 287.L ), t = std::sqrt( (114.L + 3.L*std::sqrt(583.L)) / 287.L ), B1 = 196.L / 810.L, B2 = 4.L * (178981.L + 2769.L*std::sqrt(583.L)) / 1888920.L, B3 = 4.L * (178981.L - 2769.L*std::sqrt(583.L)) / 1888920.L; const Real data[3][3] = { {r, 0.0, B1}, // 4 {s, 0.0, B2}, // 4 {t, 0.0, B3} // 4 }; const unsigned int symmetry[3] = { 3, // Full Symmetry, (x,0) 2, // Full Symmetry, (x,x) 2 // Full Symmetry, (x,x) }; _points.resize (12); _weights.resize(12); stroud_rule(data, symmetry, 3); return; } // end case SEVENTH case EIGHTH: { // A pair of degree=8 rules for the QUAD "C2" due to // Wissmann and Becker. These rules both have 16 points. // A tensor product degree-6 Gauss would have 25 points. // // J. W. Wissmann and T. Becker, Partially symmetric cubature // formulas for even degrees of exactness, SIAM J. Numer. Anal. 23 // (1986), 676--685. const Real data[10][3] = { // First of 2 degree-8, 16 point rules given by Wissmann // {0.000000000000000, 0.000000000000000, 0.055364705621440}, // {0.000000000000000, 0.757629177660505, 0.404389368726076}, // {0.000000000000000, -0.236871842255702, 0.533546604952635}, // {0.000000000000000, -0.989717929044527, 0.117054188786739}, // {0.639091304900370, 0.950520955645667, 0.125614417613747}, // {0.937069076924990, 0.663882736885633, 0.136544584733588}, // {0.537083530541494, 0.304210681724104, 0.483408479211257}, // {0.887188506449625, -0.236496718536120, 0.252528506429544}, // {0.494698820670197, -0.698953476086564, 0.361262323882172}, // {0.897495818279768, -0.900390774211580, 0.085464254086247} // // Second of 2 degree-8, 16 point rules given by Wissmann. // Either of these will work, I just chose the one with points // further into the element interior. {0.0000000000000000e+00, 6.5956013196034176e-01, 4.5027677630559029e-01}, {0.0000000000000000e+00, -9.4914292304312538e-01, 1.6657042677781274e-01}, {9.5250946607156228e-01, 7.6505181955768362e-01, 9.8869459933431422e-02}, {5.3232745407420624e-01, 9.3697598108841598e-01, 1.5369674714081197e-01}, {6.8473629795173504e-01, 3.3365671773574759e-01, 3.9668697607290278e-01}, {2.3314324080140552e-01, -7.9583272377396852e-02, 3.5201436794569501e-01}, {9.2768331930611748e-01, -2.7224008061253425e-01, 1.8958905457779799e-01}, {4.5312068740374942e-01, -6.1373535339802760e-01, 3.7510100114758727e-01}, {8.3750364042281223e-01, -8.8847765053597136e-01, 1.2561879164007201e-01} }; _points.resize(16); _weights.resize(16); wissmann_rule(data, /*10*/ 9); return; } // end case EIGHTH case NINTH: { // A degree 9, 17-point rule due to Moller. // // H.M. Moller, Kubaturformeln mit minimaler Knotenzahl, // Numer. Math. 25 (1976), 185--200. // // This rule is provably minimal in the number of points. // A tensor-product rule accurate for "bi-ninth" degree polynomials would have 25 points. const Real data[5][3] = { {0.0000000000000000e+00, 0.0000000000000000e+00, 5.2674897119341563e-01}, // 1 {6.3068011973166885e-01, 9.6884996636197772e-01, 8.8879378170198706e-02}, // 4 {9.2796164595956966e-01, 7.5027709997890053e-01, 1.1209960212959648e-01}, // 4 {4.5333982113564719e-01, 5.2373582021442933e-01, 3.9828243926207009e-01}, // 4 {8.5261572933366230e-01, 7.6208328192617173e-02, 2.6905133763978080e-01} // 4 }; const unsigned int symmetry[5] = { 0, // Single point 4, // Rotational Invariant 4, // Rotational Invariant 4, // Rotational Invariant 4 // Rotational Invariant }; _points.resize (17); _weights.resize(17); stroud_rule(data, symmetry, 5); return; } // end case NINTH case TENTH: case ELEVENTH: { // A degree 11, 24-point rule due to Cools and Haegemans. // // R. Cools and A. Haegemans, Another step forward in searching for // cubature formulae with a minimal number of knots for the square, // Computing 40 (1988), 139--146. // // P. Verlinden and R. Cools, The algebraic construction of a minimal // cubature formula of degree 11 for the square, Cubature Formulas // and their Applications (Russian) (Krasnoyarsk) (M.V. Noskov, ed.), // 1994, pp. 13--23. // // This rule is provably minimal in the number of points. // A tensor-product rule accurate for "bi-tenth" or "bi-eleventh" degree polynomials would have 36 points. const Real data[6][3] = { {6.9807610454956756e-01, 9.8263922354085547e-01, 4.8020763350723814e-02}, // 4 {9.3948638281673690e-01, 8.2577583590296393e-01, 6.6071329164550595e-02}, // 4 {9.5353952820153201e-01, 1.8858613871864195e-01, 9.7386777358668164e-02}, // 4 {3.1562343291525419e-01, 8.1252054830481310e-01, 2.1173634999894860e-01}, // 4 {7.1200191307533630e-01, 5.2532025036454776e-01, 2.2562606172886338e-01}, // 4 {4.2484724884866925e-01, 4.1658071912022368e-02, 3.5115871839824543e-01} // 4 }; const unsigned int symmetry[6] = { 4, // Rotational Invariant 4, // Rotational Invariant 4, // Rotational Invariant 4, // Rotational Invariant 4, // Rotational Invariant 4 // Rotational Invariant }; _points.resize (24); _weights.resize(24); stroud_rule(data, symmetry, 6); return; } // end case TENTH,ELEVENTH case TWELFTH: case THIRTEENTH: { // A degree 13, 33-point rule due to Cools and Haegemans. // // R. Cools and A. Haegemans, Another step forward in searching for // cubature formulae with a minimal number of knots for the square, // Computing 40 (1988), 139--146. // // A tensor-product rule accurate for "bi-12" or "bi-13" degree polynomials would have 49 points. const Real data[9][3] = { {0.0000000000000000e+00, 0.0000000000000000e+00, 3.0038211543122536e-01}, // 1 {9.8348668243987226e-01, 7.7880971155441942e-01, 2.9991838864499131e-02}, // 4 {8.5955600564163892e-01, 9.5729769978630736e-01, 3.8174421317083669e-02}, // 4 {9.5892517028753485e-01, 1.3818345986246535e-01, 6.0424923817749980e-02}, // 4 {3.9073621612946100e-01, 9.4132722587292523e-01, 7.7492738533105339e-02}, // 4 {8.5007667369974857e-01, 4.7580862521827590e-01, 1.1884466730059560e-01}, // 4 {6.4782163718701073e-01, 7.5580535657208143e-01, 1.2976355037000271e-01}, // 4 {7.0741508996444936e-02, 6.9625007849174941e-01, 2.1334158145718938e-01}, // 4 {4.0930456169403884e-01, 3.4271655604040678e-01, 2.5687074948196783e-01} // 4 }; const unsigned int symmetry[9] = { 0, // Single point 4, // Rotational Invariant 4, // Rotational Invariant 4, // Rotational Invariant 4, // Rotational Invariant 4, // Rotational Invariant 4, // Rotational Invariant 4, // Rotational Invariant 4 // Rotational Invariant }; _points.resize (33); _weights.resize(33); stroud_rule(data, symmetry, 9); return; } // end case TWELFTH,THIRTEENTH case FOURTEENTH: case FIFTEENTH: { // A degree-15, 48 point rule originally due to Rabinowitz and Richter, // can be found in Cools' 1971 book. // // A.H. Stroud, Approximate calculation of multiple integrals, // Prentice-Hall, Englewood Cliffs, N.J., 1971. // // The product Gauss rule for this order has 8^2=64 points. const Real data[9][3] = { {9.915377816777667e-01L, 0.0000000000000000e+00, 3.01245207981210e-02L}, // 4 {8.020163879230440e-01L, 0.0000000000000000e+00, 8.71146840209092e-02L}, // 4 {5.648674875232742e-01L, 0.0000000000000000e+00, 1.250080294351494e-01L}, // 4 {9.354392392539896e-01L, 0.0000000000000000e+00, 2.67651407861666e-02L}, // 4 {7.624563338825799e-01L, 0.0000000000000000e+00, 9.59651863624437e-02L}, // 4 {2.156164241427213e-01L, 0.0000000000000000e+00, 1.750832998343375e-01L}, // 4 {9.769662659711761e-01L, 6.684480048977932e-01L, 2.83136372033274e-02L}, // 4 {8.937128379503403e-01L, 3.735205277617582e-01L, 8.66414716025093e-02L}, // 4 {6.122485619312083e-01L, 4.078983303613935e-01L, 1.150144605755996e-01L} // 4 }; const unsigned int symmetry[9] = { 3, // Full Symmetry, (x,0) 3, // Full Symmetry, (x,0) 3, // Full Symmetry, (x,0) 2, // Full Symmetry, (x,x) 2, // Full Symmetry, (x,x) 2, // Full Symmetry, (x,x) 1, // Full Symmetry, (x,y) 1, // Full Symmetry, (x,y) 1, // Full Symmetry, (x,y) }; _points.resize (48); _weights.resize(48); stroud_rule(data, symmetry, 9); return; } // case FOURTEENTH, FIFTEENTH: case SIXTEENTH: case SEVENTEENTH: { // A degree 17, 60-point rule due to Cools and Haegemans. // // R. Cools and A. Haegemans, Another step forward in searching for // cubature formulae with a minimal number of knots for the square, // Computing 40 (1988), 139--146. // // A tensor-product rule accurate for "bi-14" or "bi-15" degree polynomials would have 64 points. // A tensor-product rule accurate for "bi-16" or "bi-17" degree polynomials would have 81 points. const Real data[10][3] = { {9.8935307451260049e-01, 0.0000000000000000e+00, 2.0614915919990959e-02}, // 4 {3.7628520715797329e-01, 0.0000000000000000e+00, 1.2802571617990983e-01}, // 4 {9.7884827926223311e-01, 0.0000000000000000e+00, 5.5117395340318905e-03}, // 4 {8.8579472916411612e-01, 0.0000000000000000e+00, 3.9207712457141880e-02}, // 4 {1.7175612383834817e-01, 0.0000000000000000e+00, 7.6396945079863302e-02}, // 4 {5.9049927380600241e-01, 3.1950503663457394e-01, 1.4151372994997245e-01}, // 8 {7.9907913191686325e-01, 5.9797245192945738e-01, 8.3903279363797602e-02}, // 8 {8.0374396295874471e-01, 5.8344481776550529e-02, 6.0394163649684546e-02}, // 8 {9.3650627612749478e-01, 3.4738631616620267e-01, 5.7387752969212695e-02}, // 8 {9.8132117980545229e-01, 7.0600028779864611e-01, 2.1922559481863763e-02}, // 8 }; const unsigned int symmetry[10] = { 3, // Fully symmetric (x,0) 3, // Fully symmetric (x,0) 2, // Fully symmetric (x,x) 2, // Fully symmetric (x,x) 2, // Fully symmetric (x,x) 1, // Fully symmetric (x,y) 1, // Fully symmetric (x,y) 1, // Fully symmetric (x,y) 1, // Fully symmetric (x,y) 1 // Fully symmetric (x,y) }; _points.resize (60); _weights.resize(60); stroud_rule(data, symmetry, 10); return; } // end case FOURTEENTH through SEVENTEENTH // By default: construct and use a Gauss quadrature rule default: { // Break out and fall down into the default: case for the // outer switch statement. break; } } // end switch(_order + 2*p) } // end case QUAD4/8/9 // By default: construct and use a Gauss quadrature rule default: QGauss gauss_rule(2, _order); gauss_rule.init(type_in, p); // Swap points and weights with the about-to-be destroyed rule. _points.swap (gauss_rule.get_points() ); _weights.swap(gauss_rule.get_weights()); return; } // end switch (type_in)

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, 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.

00030 {
00031 
00032   switch (type_in)
00033     {
00034       //---------------------------------------------
00035       // Hex quadrature rules
00036     case HEX8:
00037     case HEX20:
00038     case HEX27:
00039       {
00040         switch(_order + 2*p)
00041           {
00042 
00043             // The CONSTANT/FIRST rule is the 1-point Gauss "product" rule...we fall
00044             // through to the default case for this rule.
00045 
00046           case SECOND:
00047           case THIRD:
00048             {
00049               // A degree 3, 6-point, "rotationally-symmetric" rule by
00050               // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
00051               //
00052               // Warning: this rule contains points on the boundary of the reference
00053               // element, and therefore may be unsuitable for some problems.  The alternative
00054               // would be a 2x2x2 Gauss product rule.
00055               const Real data[1][4] =
00056                 {
00057                   {1.0L, 0.0L, 0.0L, 4.0L/3.0L}
00058                 };
00059 
00060               const unsigned int rule_id[1] = {
00061                 1 // (x,0,0) -> 6 permutations
00062               };
00063 
00064               _points.resize(6);
00065               _weights.resize(6);
00066 
00067               kim_rule(data, rule_id, 1);
00068               return;
00069             } // end case SECOND,THIRD
00070 
00071           case FOURTH:
00072           case FIFTH:
00073             {
00074               // A degree 5, 13-point rule by Stroud,
00075               // AH Stroud, "Some Fifth Degree Integration Formulas for Symmetric Regions II.",
00076               // Numerische Mathematik 9, pp. 460-468 (1967).
00077               //
00078               // This rule is provably minimal in the number of points.  The equations given for
00079               // the n-cube on pg. 466 of the paper for mu/gamma and gamma are wrong, at least for
00080               // the n=3 case.  The analytical values given here were computed by me [JWP] in Maple.
00081 
00082               // Convenient intermediate values.
00083               const Real sqrt19 = std::sqrt(19.L);
00084               const Real tp     = std::sqrt(71440.L + 6802.L*sqrt19);
00085 
00086               // Point data for permutations.
00087               const Real eta    =  0.00000000000000000000000000000000e+00L;
00088 
00089               const Real lambda =  std::sqrt(1919.L/3285.L - 148.L*sqrt19/3285.L + 4.L*tp/3285.L);
00090               // 8.8030440669930978047737818209860e-01L;
00091 
00092               const Real xi     = -std::sqrt(1121.L/3285.L +  74.L*sqrt19/3285.L - 2.L*tp/3285.L);
00093               // -4.9584817142571115281421242364290e-01L;
00094 
00095               const Real mu     =  std::sqrt(1121.L/3285.L +  74.L*sqrt19/3285.L + 2.L*tp/3285.L);
00096               // 7.9562142216409541542982482567580e-01L;
00097 
00098               const Real gamma  =  std::sqrt(1919.L/3285.L - 148.L*sqrt19/3285.L - 4.L*tp/3285.L);
00099               // 2.5293711744842581347389255929324e-02L;
00100 
00101               // Weights: the centroid weight is given analytically.  Weight B (resp C) goes
00102               // with the {lambda,xi} (resp {gamma,mu}) permutation.  The single-precision
00103               // results reported by Stroud are given for reference.
00104 
00105               const Real A      = 32.0L / 19.0L;
00106               // Stroud: 0.21052632  * 8.0 = 1.684210560;
00107 
00108               const Real B      = 1.L / ( 260072.L/133225.L  - 1520*sqrt19/133225.L + (133.L - 37.L*sqrt19)*tp/133225.L );
00109               // 5.4498735127757671684690782180890e-01L; // Stroud: 0.068123420 * 8.0 = 0.544987360;
00110 
00111               const Real C      = 1.L / ( 260072.L/133225.L  - 1520*sqrt19/133225.L - (133.L - 37.L*sqrt19)*tp/133225.L );
00112               // 5.0764422766979170420572375713840e-01L; // Stroud: 0.063455527 * 8.0 = 0.507644216;
00113 
00114               _points.resize(13);
00115               _weights.resize(13);
00116 
00117               unsigned int c=0;
00118 
00119               // Point with weight A (origin)
00120               _points[c] = Point(eta, eta, eta);
00121               _weights[c++] = A;
00122 
00123               // Points with weight B
00124               _points[c] = Point(lambda, xi, xi);
00125               _weights[c++] = B;
00126               _points[c] = -_points[c-1];
00127               _weights[c++] = B;
00128 
00129               _points[c] = Point(xi, lambda, xi);
00130               _weights[c++] = B;
00131               _points[c] = -_points[c-1];
00132               _weights[c++] = B;
00133 
00134               _points[c] = Point(xi, xi, lambda);
00135               _weights[c++] = B;
00136               _points[c] = -_points[c-1];
00137               _weights[c++] = B;
00138 
00139               // Points with weight C
00140               _points[c] = Point(mu, mu, gamma);
00141               _weights[c++] = C;
00142               _points[c] = -_points[c-1];
00143               _weights[c++] = C;
00144 
00145               _points[c] = Point(mu, gamma, mu);
00146               _weights[c++] = C;
00147               _points[c] = -_points[c-1];
00148               _weights[c++] = C;
00149 
00150               _points[c] = Point(gamma, mu, mu);
00151               _weights[c++] = C;
00152               _points[c] = -_points[c-1];
00153               _weights[c++] = C;
00154 
00155               return;
00156 
00157 
00158 //            // A degree 5, 14-point, "rotationally-symmetric" rule by
00159 //            // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
00160 //            // Was also reported in Stroud's 1971 book.
00161 //            const Real data[2][4] =
00162 //              {
00163 //                {7.95822425754221463264548820476135e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 8.86426592797783933518005540166204e-01L},
00164 //                {7.58786910639328146269034278112267e-01L, 7.58786910639328146269034278112267e-01L, 7.58786910639328146269034278112267e-01L, 3.35180055401662049861495844875346e-01L}
00165 //              };
00166 
00167 //            const unsigned int rule_id[2] = {
00168 //              1, // (x,0,0) -> 6 permutations
00169 //              4  // (x,x,x) -> 8 permutations
00170 //            };
00171 
00172 //            _points.resize(14);
00173 //            _weights.resize(14);
00174 
00175 //            kim_rule(data, rule_id, 2);
00176 //            return;
00177             } // end case FOURTH,FIFTH
00178 
00179           case SIXTH:
00180           case SEVENTH:
00181             {
00182               if (allow_rules_with_negative_weights)
00183                 {
00184                   // A degree 7, 31-point, "rotationally-symmetric" rule by
00185                   // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
00186                   // This rule contains a negative weight, so only use it if such type of
00187                   // rules are allowed.
00188                   const Real data[3][4] =
00189                     {
00190                       {0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, -1.27536231884057971014492753623188e+00L},
00191                       {5.85540043769119907612630781744060e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L,  8.71111111111111111111111111111111e-01L},
00192                       {6.94470135991704766602025803883310e-01L, 9.37161638568208038511047377665396e-01L, 4.15659267604065126239606672567031e-01L,  1.68695652173913043478260869565217e-01L}
00193                     };
00194 
00195                   const unsigned int rule_id[3] = {
00196                     0, // (0,0,0) -> 1 permutation
00197                     1, // (x,0,0) -> 6 permutations
00198                     6  // (x,y,z) -> 24 permutations
00199                   };
00200 
00201                   _points.resize(31);
00202                   _weights.resize(31);
00203 
00204                   kim_rule(data, rule_id, 3);
00205                   return;
00206                 } // end if (allow_rules_with_negative_weights)
00207 
00208 
00209               // A degree 7, 34-point, "fully-symmetric" rule, first published in
00210               // P.C. Hammer and A.H. Stroud, "Numerical Evaluation of Multiple Integrals II",
00211               // Mathmatical Tables and Other Aids to Computation, vol 12., no 64, 1958, pp. 272-280
00212               //
00213               // This rule happens to fall under the same general
00214               // construction as the Kim rules, so we've re-used
00215               // that code here.  Stroud gives 16 digits for his rule,
00216               // and this is the most accurate version I've found.
00217               //
00218               // For comparison, a SEVENTH-order Gauss product rule
00219               // (which integrates tri-7th order polynomials) would
00220               // have 4^3=64 points.
00221               const Real
00222                 r  = std::sqrt(6.L/7.L),
00223                 s  = std::sqrt((960.L - 3.L*std::sqrt(28798.L)) / 2726.L),
00224                 t  = std::sqrt((960.L + 3.L*std::sqrt(28798.L)) / 2726.L),
00225                 B1 = 8624.L / 29160.L,
00226                 B2 = 2744.L / 29160.L,
00227                 B3 = 8.L*(774.L*t*t - 230.L)/(9720.L*(t*t-s*s)),
00228                 B4 = 8.L*(230.L - 774.L*s*s)/(9720.L*(t*t-s*s));
00229 
00230               const Real data[4][4] =
00231                 {
00232                   {r, 0.L, 0.L, B1},
00233                   {r,   r, 0.L, B2},
00234                   {s,   s,   s, B3},
00235                   {t,   t,   t, B4}
00236                 };
00237 
00238               const unsigned int rule_id[4] = {
00239                 1, // (x,0,0) -> 6 permutations
00240                 2, // (x,x,0) -> 12 permutations
00241                 4, // (x,x,x) -> 8 permutations
00242                 4  // (x,x,x) -> 8 permutations
00243                   };
00244 
00245               _points.resize(34);
00246               _weights.resize(34);
00247 
00248               kim_rule(data, rule_id, 4);
00249               return;
00250 
00251 
00252 //            // A degree 7, 38-point, "rotationally-symmetric" rule by
00253 //            // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
00254 //            //
00255 //            // This rule is obviously inferior to the 34-point rule above...
00256 //            const Real data[3][4] =
00257 //              {
00258 //                {9.01687807821291289082811566285950e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 2.95189738262622903181631100062774e-01L},
00259 //                {4.08372221499474674069588900002128e-01L, 4.08372221499474674069588900002128e-01L, 4.08372221499474674069588900002128e-01L, 4.04055417266200582425904380777126e-01L},
00260 //                {8.59523090201054193116477875786220e-01L, 8.59523090201054193116477875786220e-01L, 4.14735913727987720499709244748633e-01L, 1.24850759678944080062624098058597e-01L}
00261 //              };
00262 //
00263 //            const unsigned int rule_id[3] = {
00264 //              1, // (x,0,0) -> 6 permutations
00265 //              4, // (x,x,x) -> 8 permutations
00266 //              5  // (x,x,z) -> 24 permutations
00267 //            };
00268 //
00269 //            _points.resize(38);
00270 //            _weights.resize(38);
00271 //
00272 //            kim_rule(data, rule_id, 3);
00273 //            return;
00274             } // end case SIXTH,SEVENTH
00275 
00276           case EIGHTH:
00277             {
00278               // A degree 8, 47-point, "rotationally-symmetric" rule by
00279               // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
00280               //
00281               // A EIGHTH-order Gauss product rule (which integrates tri-8th order polynomials)
00282               // would have 5^3=125 points.
00283               const Real data[5][4] =
00284                 {
00285                   {0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 4.51903714875199690490763818699555e-01L},
00286                   {7.82460796435951590652813975429717e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 2.99379177352338919703385618576171e-01L},
00287                   {4.88094669706366480526729301468686e-01L, 4.88094669706366480526729301468686e-01L, 4.88094669706366480526729301468686e-01L, 3.00876159371240019939698689791164e-01L},
00288                   {8.62218927661481188856422891110042e-01L, 8.62218927661481188856422891110042e-01L, 8.62218927661481188856422891110042e-01L, 4.94843255877038125738173175714853e-02L},
00289                   {2.81113909408341856058098281846420e-01L, 9.44196578292008195318687494773744e-01L, 6.97574833707236996779391729948984e-01L, 1.22872389222467338799199767122592e-01L}
00290                 };
00291 
00292               const unsigned int rule_id[5] = {
00293                 0, // (0,0,0) -> 1 permutation
00294                 1, // (x,0,0) -> 6 permutations
00295                 4, // (x,x,x) -> 8 permutations
00296                 4, // (x,x,x) -> 8 permutations
00297                 6  // (x,y,z) -> 24 permutations
00298               };
00299 
00300               _points.resize(47);
00301               _weights.resize(47);
00302 
00303               kim_rule(data, rule_id, 5);
00304               return;
00305             } // end case EIGHTH
00306 
00307 
00308             // By default: construct and use a Gauss quadrature rule
00309           default:
00310             {
00311               // Break out and fall down into the default: case for the
00312               // outer switch statement.
00313               break;
00314             }
00315 
00316           } // end switch(_order + 2*p)
00317       } // end case HEX8/20/27
00318 
00319 
00320       // By default: construct and use a Gauss quadrature rule
00321     default:
00322       {
00323         QGauss gauss_rule(3, _order);
00324         gauss_rule.init(type_in, p);
00325 
00326         // Swap points and weights with the about-to-be destroyed rule.
00327         _points.swap (gauss_rule.get_points() );
00328         _weights.swap(gauss_rule.get_weights());
00329 
00330         return;
00331       }
00332     } // end switch (type_in)
00333 }

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

00218 {
00219   for (unsigned int i=0, c=0; i<n_pts; ++i)
00220     {
00221       const Real
00222         x=rule_data[i][0],
00223         y=rule_data[i][1],
00224         z=rule_data[i][2],
00225         wt=rule_data[i][3];
00226 
00227       switch(rule_id[i])
00228         {
00229         case 0: // (0,0,0) 1 permutation
00230           {
00231             _points[c]  = Point( x, y, z);          _weights[c++] = wt;
00232 
00233             break;
00234           }
00235         case 1: //  (x,0,0) 6 permutations
00236           {
00237             _points[c] = Point( x, 0., 0.);         _weights[c++] = wt;
00238             _points[c] = Point(0., -x, 0.);         _weights[c++] = wt;
00239             _points[c] = Point(-x, 0., 0.);         _weights[c++] = wt;
00240             _points[c] = Point(0.,  x, 0.);         _weights[c++] = wt;
00241             _points[c] = Point(0., 0., -x);         _weights[c++] = wt;
00242             _points[c] = Point(0., 0.,  x);         _weights[c++] = wt;
00243 
00244             break;
00245           }
00246         case 2: // (x,x,0) 12 permutations
00247           {
00248             _points[c] = Point( x,  x, 0.);         _weights[c++] = wt;
00249             _points[c] = Point( x, -x, 0.);         _weights[c++] = wt;
00250             _points[c] = Point(-x, -x, 0.);         _weights[c++] = wt;
00251             _points[c] = Point(-x,  x, 0.);         _weights[c++] = wt;
00252             _points[c] = Point( x, 0., -x);         _weights[c++] = wt;
00253             _points[c] = Point( x, 0.,  x);         _weights[c++] = wt;
00254             _points[c] = Point(0.,  x, -x);         _weights[c++] = wt;
00255             _points[c] = Point(0.,  x,  x);         _weights[c++] = wt;
00256             _points[c] = Point(0., -x, -x);         _weights[c++] = wt;
00257             _points[c] = Point(-x, 0., -x);         _weights[c++] = wt;
00258             _points[c] = Point(0., -x,  x);         _weights[c++] = wt;
00259             _points[c] = Point(-x, 0.,  x);         _weights[c++] = wt;
00260 
00261             break;
00262           }
00263         case 3: // (x,y,0) 24 permutations
00264           {
00265             _points[c] = Point( x,  y, 0.);         _weights[c++] = wt;
00266             _points[c] = Point( y, -x, 0.);         _weights[c++] = wt;
00267             _points[c] = Point(-x, -y, 0.);         _weights[c++] = wt;
00268             _points[c] = Point(-y,  x, 0.);         _weights[c++] = wt;
00269             _points[c] = Point( x, 0., -y);         _weights[c++] = wt;
00270             _points[c] = Point( x, -y, 0.);         _weights[c++] = wt;
00271             _points[c] = Point( x, 0.,  y);         _weights[c++] = wt;
00272             _points[c] = Point(0.,  y, -x);         _weights[c++] = wt;
00273             _points[c] = Point(-x,  y, 0.);         _weights[c++] = wt;
00274             _points[c] = Point(0.,  y,  x);         _weights[c++] = wt;
00275             _points[c] = Point( y, 0., -x);         _weights[c++] = wt;
00276             _points[c] = Point(0., -y, -x);         _weights[c++] = wt;
00277             _points[c] = Point(-y, 0., -x);         _weights[c++] = wt;
00278             _points[c] = Point( y,  x, 0.);         _weights[c++] = wt;
00279             _points[c] = Point(-y, -x, 0.);         _weights[c++] = wt;
00280             _points[c] = Point( y, 0.,  x);         _weights[c++] = wt;
00281             _points[c] = Point(0., -y,  x);         _weights[c++] = wt;
00282             _points[c] = Point(-y, 0.,  x);         _weights[c++] = wt;
00283             _points[c] = Point(-x, 0.,  y);         _weights[c++] = wt;
00284             _points[c] = Point(0., -x, -y);         _weights[c++] = wt;
00285             _points[c] = Point(0., -x,  y);         _weights[c++] = wt;
00286             _points[c] = Point(-x, 0., -y);         _weights[c++] = wt;
00287             _points[c] = Point(0.,  x,  y);         _weights[c++] = wt;
00288             _points[c] = Point(0.,  x, -y);         _weights[c++] = wt;
00289 
00290             break;
00291           }
00292         case 4: // (x,x,x) 8 permutations
00293           {
00294             _points[c] = Point( x,  x,  x);         _weights[c++] = wt;
00295             _points[c] = Point( x, -x,  x);         _weights[c++] = wt;
00296             _points[c] = Point(-x, -x,  x);         _weights[c++] = wt;
00297             _points[c] = Point(-x,  x,  x);         _weights[c++] = wt;
00298             _points[c] = Point( x,  x, -x);         _weights[c++] = wt;
00299             _points[c] = Point( x, -x, -x);         _weights[c++] = wt;
00300             _points[c] = Point(-x,  x, -x);         _weights[c++] = wt;
00301             _points[c] = Point(-x, -x, -x);         _weights[c++] = wt;
00302 
00303             break;
00304           }
00305         case 5: // (x,x,z) 24 permutations
00306           {
00307             _points[c] = Point( x,  x,  z);         _weights[c++] = wt;
00308             _points[c] = Point( x, -x,  z);         _weights[c++] = wt;
00309             _points[c] = Point(-x, -x,  z);         _weights[c++] = wt;
00310             _points[c] = Point(-x,  x,  z);         _weights[c++] = wt;
00311             _points[c] = Point( x,  z, -x);         _weights[c++] = wt;
00312             _points[c] = Point( x, -x, -z);         _weights[c++] = wt;
00313             _points[c] = Point( x, -z,  x);         _weights[c++] = wt;
00314             _points[c] = Point( z,  x, -x);         _weights[c++] = wt;
00315             _points[c] = Point(-x,  x, -z);         _weights[c++] = wt;
00316             _points[c] = Point(-z,  x,  x);         _weights[c++] = wt;
00317             _points[c] = Point( x, -z, -x);         _weights[c++] = wt;
00318             _points[c] = Point(-z, -x, -x);         _weights[c++] = wt;
00319             _points[c] = Point(-x,  z, -x);         _weights[c++] = wt;
00320             _points[c] = Point( x,  x, -z);         _weights[c++] = wt;
00321             _points[c] = Point(-x, -x, -z);         _weights[c++] = wt;
00322             _points[c] = Point( x,  z,  x);         _weights[c++] = wt;
00323             _points[c] = Point( z, -x,  x);         _weights[c++] = wt;
00324             _points[c] = Point(-x, -z,  x);         _weights[c++] = wt;
00325             _points[c] = Point(-x,  z,  x);         _weights[c++] = wt;
00326             _points[c] = Point( z, -x, -x);         _weights[c++] = wt;
00327             _points[c] = Point(-z, -x,  x);         _weights[c++] = wt;
00328             _points[c] = Point(-x, -z, -x);         _weights[c++] = wt;
00329             _points[c] = Point( z,  x,  x);         _weights[c++] = wt;
00330             _points[c] = Point(-z,  x, -x);         _weights[c++] = wt;
00331 
00332             break;
00333           }
00334         case 6: // (x,y,z) 24 permutations
00335           {
00336             _points[c] = Point( x,  y,  z);         _weights[c++] = wt;
00337             _points[c] = Point( y, -x,  z);         _weights[c++] = wt;
00338             _points[c] = Point(-x, -y,  z);         _weights[c++] = wt;
00339             _points[c] = Point(-y,  x,  z);         _weights[c++] = wt;
00340             _points[c] = Point( x,  z, -y);         _weights[c++] = wt;
00341             _points[c] = Point( x, -y, -z);         _weights[c++] = wt;
00342             _points[c] = Point( x, -z,  y);         _weights[c++] = wt;
00343             _points[c] = Point( z,  y, -x);         _weights[c++] = wt;
00344             _points[c] = Point(-x,  y, -z);         _weights[c++] = wt;
00345             _points[c] = Point(-z,  y,  x);         _weights[c++] = wt;
00346             _points[c] = Point( y, -z, -x);         _weights[c++] = wt;
00347             _points[c] = Point(-z, -y, -x);         _weights[c++] = wt;
00348             _points[c] = Point(-y,  z, -x);         _weights[c++] = wt;
00349             _points[c] = Point( y,  x, -z);         _weights[c++] = wt;
00350             _points[c] = Point(-y, -x, -z);         _weights[c++] = wt;
00351             _points[c] = Point( y,  z,  x);         _weights[c++] = wt;
00352             _points[c] = Point( z, -y,  x);         _weights[c++] = wt;
00353             _points[c] = Point(-y, -z,  x);         _weights[c++] = wt;
00354             _points[c] = Point(-x,  z,  y);         _weights[c++] = wt;
00355             _points[c] = Point( z, -x, -y);         _weights[c++] = wt;
00356             _points[c] = Point(-z, -x,  y);         _weights[c++] = wt;
00357             _points[c] = Point(-x, -z, -y);         _weights[c++] = wt;
00358             _points[c] = Point( z,  x,  y);         _weights[c++] = wt;
00359             _points[c] = Point(-z,  x, -y);         _weights[c++] = wt;
00360 
00361             break;
00362           }
00363         default:
00364           {
00365             libMesh::err << "Unknown rule ID: " << rule_id[i] << "!" << std::endl;
00366             libmesh_error();
00367           }
00368         } // end switch(rule_id[i])
00369     }
00370 }

static unsigned int libMesh::ReferenceCounter::n_objects ( ) [inline, static, inherited]

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.

00080   { return _n_objects; }

unsigned int libMesh::QBase::n_points ( ) const [inline, inherited]

Returns:: the number of points associated with the quadrature rule.

Definition at line 116 of file quadrature.h.

References libMesh::QBase::_points.

Referenced by libMesh::FEGenericBase< OutputType >::coarsened_dof_values(), libMesh::InfFE< Dim, T_radial, T_map >::combine_base_radial(), libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), libMesh::InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), libMesh::InfFE< Dim, T_radial, T_map >::init_shape_functions(), libMesh::FE< Dim, T >::n_quadrature_points(), libMesh::ProjectFEMSolution::operator()(), and libMesh::QBase::print_info().

00117     { libmesh_assert (!_points.empty());
00118       return libmesh_cast_int<unsigned int>(_points.size()); }

void libMesh::ReferenceCounter::print_info ( std::ostream & out = libMesh::out ) [static, inherited]

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

00089 {
00090   if( _enable_print_counter ) out_stream << ReferenceCounter::get_info();
00091 }

void libMesh::QBase::print_info ( std::ostream & os = libMesh::out ) const [inline, inherited]

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

Referenced by libMesh::operator<<().

00363 {
00364   libmesh_assert(!_points.empty());
00365   libmesh_assert(!_weights.empty());
00366 
00367   os << "N_Q_Points=" << this->n_points() << std::endl << std::endl;
00368   for (unsigned int qpoint=0; qpoint<this->n_points(); qpoint++)
00369     {
00370       os << " Point " << qpoint << ":\n"
00371          << "  "
00372          << _points[qpoint]
00373          << " Weight:\n "
00374          << "  w=" << _weights[qpoint] << "\n" << std::endl;
00375     }
00376 }

Point libMesh::QBase::qp ( const unsigned int i ) const [inline, inherited]

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

00151     { 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, and libMesh::Real.

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

00084 {
00085   // Make sure we are in 1D
00086   libmesh_assert_equal_to (_dim, 1);
00087 
00088   // Make sure that we have sane ranges
00089   libmesh_assert_greater (new_range.second, new_range.first);
00090   libmesh_assert_greater (old_range.second, old_range.first);
00091 
00092   // Make sure there are some points
00093   libmesh_assert_greater (_points.size(), 0);
00094 
00095   // We're mapping from old_range -> new_range
00096   for (unsigned int i=0; i<_points.size(); i++)
00097     {
00098       _points[i](0) =
00099         (_points[i](0) - old_range.first) *
00100         (new_range.second - new_range.first) /
00101         (old_range.second - old_range.first) +
00102         new_range.first;
00103     }
00104 
00105   // Compute the scale factor and scale the weights
00106   const Real scfact = (new_range.second - new_range.first) /
00107                       (old_range.second - old_range.first);
00108 
00109   for (unsigned int i=0; i<_points.size(); i++)
00110     _weights[i] *= scfact;
00111 }

virtual bool libMesh::QBase::shapes_need_reinit ( ) [inline, virtual, inherited]

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.

Referenced by libMesh::FE< Dim, T >::edge_reinit(), and libMesh::FE< Dim, T >::reinit().

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

00067 {
00068   for (unsigned int i=0, c=0; i<n_pts; ++i)
00069     {
00070       const Real
00071         x=rule_data[i][0],
00072         y=rule_data[i][1],
00073         wt=rule_data[i][2];
00074 
00075       switch(rule_symmetry[i])
00076         {
00077         case 0: // Single point (no symmetry)
00078           {
00079             _points[c]  = Point( x, y);
00080             _weights[c++] = wt;
00081 
00082             break;
00083           }
00084         case 1: // Fully-symmetric (x,y)
00085           {
00086             _points[c]    = Point( x, y);
00087             _weights[c++] = wt;
00088 
00089             _points[c]    = Point(-x, y);
00090             _weights[c++] = wt;
00091 
00092             _points[c]    = Point( x,-y);
00093             _weights[c++] = wt;
00094 
00095             _points[c]    = Point(-x,-y);
00096             _weights[c++] = wt;
00097 
00098             _points[c]    = Point( y, x);
00099             _weights[c++] = wt;
00100 
00101             _points[c]    = Point(-y, x);
00102             _weights[c++] = wt;
00103 
00104             _points[c]    = Point( y,-x);
00105             _weights[c++] = wt;
00106 
00107             _points[c]    = Point(-y,-x);
00108             _weights[c++] = wt;
00109 
00110             break;
00111           }
00112         case 2: // Fully-symmetric (x,x)
00113           {
00114             _points[c]    = Point( x, x);
00115             _weights[c++] = wt;
00116 
00117             _points[c]    = Point(-x, x);
00118             _weights[c++] = wt;
00119 
00120             _points[c]    = Point( x,-x);
00121             _weights[c++] = wt;
00122 
00123             _points[c]    = Point(-x,-x);
00124             _weights[c++] = wt;
00125 
00126             break;
00127           }
00128         case 3: // Fully-symmetric (x,0)
00129           {
00130             libmesh_assert_equal_to (y, 0.0);
00131 
00132             _points[c]    = Point( x,0.);
00133             _weights[c++] = wt;
00134 
00135             _points[c]    = Point(-x,0.);
00136             _weights[c++] = wt;
00137 
00138             _points[c]    = Point(0., x);
00139             _weights[c++] = wt;
00140 
00141             _points[c]    = Point(0.,-x);
00142             _weights[c++] = wt;
00143 
00144             break;
00145           }
00146         case 4: // Rotational invariant
00147           {
00148             _points[c]    = Point( x, y);
00149             _weights[c++] = wt;
00150 
00151             _points[c]    = Point(-x,-y);
00152             _weights[c++] = wt;
00153 
00154             _points[c]    = Point(-y, x);
00155             _weights[c++] = wt;
00156 
00157             _points[c]    = Point( y,-x);
00158             _weights[c++] = wt;
00159 
00160             break;
00161           }
00162         case 5: // Partial symmetry (Wissman's rules)
00163           {
00164             libmesh_assert_not_equal_to (x, 0.0);
00165 
00166             _points[c]    = Point( x, y);
00167             _weights[c++] = wt;
00168 
00169             _points[c]    = Point(-x, y);
00170             _weights[c++] = wt;
00171 
00172             break;
00173           }
00174         case 6: // Rectangular symmetry
00175           {
00176             _points[c]    = Point( x, y);
00177             _weights[c++] = wt;
00178 
00179             _points[c]    = Point(-x, y);
00180             _weights[c++] = wt;
00181 
00182             _points[c]    = Point(-x,-y);
00183             _weights[c++] = wt;
00184 
00185             _points[c]    = Point( x,-y);
00186             _weights[c++] = wt;
00187 
00188             break;
00189           }
00190         case 7: // Central symmetry
00191           {
00192             libmesh_assert_equal_to (x, 0.0);
00193             libmesh_assert_not_equal_to (y, 0.0);
00194 
00195             _points[c]    = Point(0., y);
00196             _weights[c++] = wt;
00197 
00198             _points[c]    = Point(0.,-y);
00199             _weights[c++] = wt;
00200 
00201             break;
00202           }
00203         default:
00204           {
00205             libMesh::err << "Unknown symmetry!" << std::endl;
00206             libmesh_error();
00207           }
00208         } // end switch(rule_symmetry[i])
00209     }
00210 }

QuadratureType libMesh::QMonomial::type ( ) const [inline, virtual]

Returns:: QMONOMIAL

Implements libMesh::QBase.

Definition at line 83 of file quadrature_monomial.h.

References libMeshEnums::QMONOMIAL.

00083 { return QMONOMIAL; }

Real libMesh::QBase::w ( const unsigned int i ) const [inline, inherited]

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

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

00046 {
00047   for (unsigned int i=0, c=0; i<n_pts; ++i)
00048     {
00049       _points[c]  = Point( rule_data[i][0], rule_data[i][1] );
00050       _weights[c++] = rule_data[i][2];
00051 
00052       // This may be an (x1,x2) -> (-x1,x2) point, in which case
00053       // we will also generate the mirror point using the same weight.
00054       if (rule_data[i][0] != static_cast<Real>(0.0))
00055         {
00056           _points[c]  = Point( -rule_data[i][0], rule_data[i][1] );
00057           _weights[c++] = rule_data[i][2];
00058         }
00059     }
00060 }

Friends And Related Function Documentation

std::ostream& operator<<	(	std::ostream &	os,
		const QBase &	q
	)			`[friend, inherited]`

Same as above, but allows you to use the stream syntax.

Member Data Documentation

ReferenceCounter::Counts libMesh::ReferenceCounter::_counts [static, protected, inherited]

Actually holds the data.

Definition at line 118 of file reference_counter.h.

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

bool libMesh::ReferenceCounter::_enable_print_counter = true [static, protected, inherited]

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 [static, protected, inherited]

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 [static, protected, inherited]

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 [protected, inherited]

Definition at line 332 of file quadrature.h.

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

std::vector<Real> libMesh::QBase::_weights [protected, inherited]

The value of the quadrature weights.

Definition at line 337 of file quadrature.h.

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

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

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

libMesh::QMonomial Class Reference

Public Member Functions

Static Public Member Functions

Public Attributes

Protected Types

Protected Member Functions

Protected Attributes

Static Protected Attributes

Private Member Functions

Friends

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Friends And Related Function Documentation

Member Data Documentation