libMesh
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Public Attributes | Protected Types | Protected Member Functions | Protected Attributes | Static Protected Attributes | Private Member Functions | List of all members
libMesh::QMonomial Class Referencefinal

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

#include <quadrature_monomial.h>

Inheritance diagram for libMesh::QMonomial:
[legend]

Public Member Functions

 QMonomial (unsigned int dim, Order order=INVALID_ORDER)
 Constructor.
 
 QMonomial (const QMonomial &)=default
 Copy/move ctor, copy/move assignment operator, and destructor are all explicitly defaulted for this simple class.
 
 QMonomial (QMonomial &&)=default
 
QMonomial & operator= (const QMonomial &)=default
 
QMonomial & operator= (QMonomial &&)=default
 
virtual ~QMonomial ()=default
 
virtual QuadratureType type () const override
 
virtual std::unique_ptr< QBase > clone () const override
 
ElemType get_elem_type () const
 
unsigned int get_p_level () const
 
unsigned int n_points () const
 
unsigned int size () const
 Alias for n_points() to enable use in index_range.
 
unsigned int get_dim () const
 
const std::vector< Point > & get_points () const
 
std::vector< Point > & get_points ()
 
const std::vector< Real > & get_weights () const
 
std::vector< Real > & get_weights ()
 
Point qp (const unsigned int i) const
 
Real w (const unsigned int i) const
 
virtual void init (const Elem &e, unsigned int p_level=invalid_uint)
 Initializes the data structures for a quadrature rule for the element e.
 
virtual void init (const ElemType type=INVALID_ELEM, unsigned int p_level=0, bool simple_type_only=false)
 Initializes the data structures for a quadrature rule for an element of type type.
 
virtual void init (const QBase &other_rule)
 Initializes the data structures for a quadrature rule based on the element, element type, and p_level settings of other_rule.
 
virtual void init (const Elem &elem, const std::vector< Real > &vertex_distance_func, unsigned int p_level=0)
 Initializes the data structures for an element potentially "cut" by a signed distance function.
 
Order get_order () const
 
Order get_base_order () const
 
void print_info (std::ostream &os=libMesh::out) const
 Prints information relevant to the quadrature rule, by default to libMesh::out.
 
void scale (std::pair< Real, Real > old_range, std::pair< Real, Real > new_range)
 Maps the points of a 1D quadrature rule defined by "old_range" to another 1D interval defined by "new_range" and scales the weights accordingly.
 
virtual bool shapes_need_reinit ()
 

Static Public Member Functions

static std::unique_ptr< QBase > build (std::string_view name, const unsigned int dim, const Order order=INVALID_ORDER)
 Builds a specific quadrature rule based on the name string.
 
static std::unique_ptr< QBase > build (const QuadratureType qt, const unsigned int dim, const Order order=INVALID_ORDER)
 Builds a specific quadrature rule based on the QuadratureType.
 
static void print_info (std::ostream &out_stream=libMesh::out)
 Prints the reference information, by default to libMesh::out.
 
static std::string get_info ()
 Gets a string containing the reference information.
 
static unsigned int n_objects ()
 Prints the number of outstanding (created, but not yet destroyed) objects.
 
static void enable_print_counter_info ()
 Methods to enable/disable the reference counter output from print_info().
 
static void disable_print_counter_info ()
 

Public Attributes

bool allow_rules_with_negative_weights
 Flag (default true) controlling the use of quadrature rules with negative weights.
 
bool allow_nodal_pyramid_quadrature
 The flag's value defaults to false so that one does not accidentally use a nodal quadrature rule on Pyramid elements, since evaluating the inverse element Jacobian (e.g.
 

Protected Types

typedef std::map< std::string, std::pair< unsigned int, unsigned int > > Counts
 Data structure to log the information.
 

Protected Member Functions

virtual void init_0D ()
 Initializes the 0D quadrature rule by filling the points and weights vectors with the appropriate values.
 
void tensor_product_quad (const QBase &q1D)
 Constructs a 2D rule from the tensor product of q1D with itself.
 
void tensor_product_hex (const QBase &q1D)
 Computes the tensor product quadrature rule [q1D x q1D x q1D] from the 1D rule q1D.
 
void tensor_product_prism (const QBase &q1D, const QBase &q2D)
 Computes the tensor product of a 1D quadrature rule and a 2D quadrature rule.
 
void increment_constructor_count (const std::string &name) noexcept
 Increments the construction counter.
 
void increment_destructor_count (const std::string &name) noexcept
 Increments the destruction counter.
 

Protected Attributes

unsigned int _dim
 The spatial dimension of the quadrature rule.
 
Order _order
 The polynomial order which the quadrature rule is capable of integrating exactly.
 
ElemType _type
 The type of element for which the current values have been computed.
 
const Elem * _elem
 The element for which the current values were computed, or nullptr if values were computed without a specific element.
 
unsigned int _p_level
 The p-level of the element for which the current values have been computed.
 
std::vector< Point > _points
 The locations of the quadrature points in reference element space.
 
std::vector< Real > _weights
 The quadrature weights.
 

Static Protected Attributes

static Counts _counts
 Actually holds the data.
 
static Threads::atomic< unsigned int > _n_objects
 The number of objects.
 
static Threads::spin_mutex _mutex
 Mutual exclusion object to enable thread-safe reference counting.
 
static bool _enable_print_counter = true
 Flag to control whether reference count information is printed when print_info is called.
 

Private Member Functions

virtual void init_1D () override
 Uses a Gauss rule in 1D.
 
virtual void init_2D () override
 Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate values.
 
virtual void init_3D () override
 Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate values.
 
void wissmann_rule (const Real rule_data[][3], const unsigned int n_pts)
 Wissmann published three interesting "partially symmetric" rules for integrating degree 4, 6, and 8 polynomials exactly on QUADs.
 
void stroud_rule (const Real rule_data[][3], const unsigned int *rule_symmetry, const unsigned int n_pts)
 Stroud's rules for quads and hexes can have one of several different types of symmetry.
 
void kim_rule (const Real rule_data[][4], const unsigned int *rule_id, const unsigned int n_pts)
 Rules from Kim and Song, Comm.
 

Detailed Description

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

While tensor product rules are optimal 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 provides quadrature rules which are more efficient than tensor product rules when they are available, and falls back on Gaussian quadrature rules otherwise.

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 when it is available.

Author
John W. Peterson
Date
2008

Implements quadrature rules for non-tensor polynomials.

Definition at line 58 of file quadrature_monomial.h.

Member Typedef Documentation

◆ Counts

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

Data structure to log the information.

The log is identified by the class name.

Definition at line 119 of file reference_counter.h.

Constructor & Destructor Documentation

◆ QMonomial() [1/3]

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

Constructor.

Declares the order of the quadrature rule.

Definition at line 65 of file quadrature_monomial.h.

66 :
67 QBase(dim,order)
68 {
69 if (dim == 1)
70 init(EDGE2);
71 }
dim
unsigned int dim
Definition adaptivity_ex3.C:114
libMesh::QBase::init
virtual void init(const Elem &e, unsigned int p_level=invalid_uint)
Initializes the data structures for a quadrature rule for the element e.
Definition quadrature.C:65
libMesh::QBase::QBase
QBase(unsigned int dim, Order order=INVALID_ORDER)
Constructor.
Definition quadrature.C:27

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

◆ QMonomial() [2/3]

libMesh::QMonomial::QMonomial ( const QMonomial &  )
default

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

◆ QMonomial() [3/3]

libMesh::QMonomial::QMonomial ( QMonomial &&  )
default

◆ ~QMonomial()

virtual libMesh::QMonomial::~QMonomial ( )
virtualdefault

Member Function Documentation

◆ build() [1/2]

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

Builds a specific quadrature rule based on the QuadratureType.

This enables selection of the quadrature rule at run-time.

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

Definition at line 54 of file quadrature_build.C.

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

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

◆ build() [2/2]

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

Builds a specific quadrature rule based on the name string.

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

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

Definition at line 43 of file quadrature_build.C.

46{
47 return QBase::build (Utility::string_to_enum<QuadratureType> (type),
48 _dim,
49 _order);
50}
libMesh::QBase::type
virtual QuadratureType type() const =0
libMesh::QBase::build
static std::unique_ptr< QBase > build(std::string_view name, const unsigned int dim, const Order order=INVALID_ORDER)
Builds a specific quadrature rule based on the name string.
Definition quadrature_build.C:43

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

Referenced by assemble_poisson(), libMesh::InfFE< Dim, T_radial, T_map >::attach_quadrature_rule(), libMesh::QBase::build(), libMesh::QBase::clone(), main(), libMesh::OverlapCoupling::OverlapCoupling(), libMesh::InfFE< Dim, T_radial, T_map >::reinit(), AllRBBTest::test_cylinder(), QuadratureTest::testBuild(), QuadratureTest::testJacobi(), and QuadratureTest::testPolynomials().

◆ clone()

std::unique_ptr< QBase > libMesh::QMonomial::clone ( ) const
overridevirtual
Returns
A copy of this quadrature rule wrapped in a smart pointer.

Reimplemented from libMesh::QBase.

Definition at line 37 of file quadrature_monomial.C.

38{
39 return std::make_unique<QMonomial>(*this);
40}

◆ disable_print_counter_info()

void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited

Definition at line 100 of file reference_counter.C.

101{
102 _enable_print_counter = false;
103 return;
104}
libMesh::ReferenceCounter::_enable_print_counter
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called.
Definition reference_counter.h:143

References libMesh::ReferenceCounter::_enable_print_counter.

◆ enable_print_counter_info()

void libMesh::ReferenceCounter::enable_print_counter_info ( )
staticinherited

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

Enabled by default.

Definition at line 94 of file reference_counter.C.

95{
96 _enable_print_counter = true;
97 return;
98}

References libMesh::ReferenceCounter::_enable_print_counter.

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

◆ get_base_order()

Order libMesh::QBase::get_base_order ( ) const
inlineinherited
Returns
The "base" order of the quadrature rule, independent of element.

This function should be used when comparing quadrature objects independently of their last initialization.

Definition at line 258 of file quadrature.h.

258{ return static_cast<Order>(_order); }

References libMesh::QBase::_order.

◆ get_dim()

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

Definition at line 150 of file quadrature.h.

150{ return _dim; }

References libMesh::QBase::_dim.

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

◆ get_elem_type()

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

Definition at line 121 of file quadrature.h.

121{ return _type; }
libMesh::QBase::_type
ElemType _type
The type of element for which the current values have been computed.
Definition quadrature.h:391

References libMesh::QBase::_type.

◆ get_info()

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

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

48{
49#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
50
51 std::ostringstream oss;
52
53 oss << '\n'
54 << " ---------------------------------------------------------------------------- \n"
55 << "| Reference count information |\n"
56 << " ---------------------------------------------------------------------------- \n";
57
58 for (const auto & [name, cd] : _counts)
59 oss << "| " << name << " reference count information:\n"
60 << "| Creations: " << cd.first << '\n'
61 << "| Destructions: " << cd.second << '\n';
62
63 oss << " ---------------------------------------------------------------------------- \n";
64
65 return oss.str();
66
67#else
68
69 return "";
70
71#endif
72}
libMesh::ReferenceCounter::_counts
static Counts _counts
Actually holds the data.
Definition reference_counter.h:124
libMesh::Quality::name
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition elem_quality.C:42

References libMesh::ReferenceCounter::_counts.

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

◆ get_order()

Order libMesh::QBase::get_order ( ) const
inlineinherited
Returns
The current "total" order of the quadrature rule which can vary element by element, depending on the Elem::p_level(), which gets passed to us during init().

Each additional power of p increases the quadrature order required to integrate the mass matrix by 2, hence the formula below.

Todo:
This function should also be used in all of the Order switch statements in the rules themselves.

Definition at line 249 of file quadrature.h.

249{ return static_cast<Order>(_order + 2 * _p_level); }
libMesh::QBase::_p_level
unsigned int _p_level
The p-level of the element for which the current values have been computed.
Definition quadrature.h:403

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

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::attach_quadrature_rule(), libMesh::QBase::clone(), libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), libMesh::RBEIMConstruction::enrich_eim_approximation_on_interiors(), libMesh::QConical::init_1D(), libMesh::QGauss::init_1D(), libMesh::QGaussLobatto::init_1D(), libMesh::QGrundmann_Moller::init_1D(), libMesh::QJacobi::init_1D(), libMesh::QGauss::init_2D(), libMesh::QGaussLobatto::init_2D(), libMesh::QGrundmann_Moller::init_2D(), init_2D(), libMesh::QGauss::init_3D(), libMesh::QGaussLobatto::init_3D(), libMesh::QGrundmann_Moller::init_3D(), init_3D(), and libMesh::RBParametrizedFunction::preevaluate_parametrized_function_on_mesh().

◆ get_p_level()

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

Definition at line 126 of file quadrature.h.

126{ return _p_level; }

References libMesh::QBase::_p_level.

◆ get_points() [1/2]

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

Definition at line 162 of file quadrature.h.

162{ return _points; }
libMesh::QBase::_points
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition quadrature.h:409

References libMesh::QBase::_points.

◆ get_points() [2/2]

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

Definition at line 156 of file quadrature.h.

156{ return _points; }

References libMesh::QBase::_points.

Referenced by assemble_ellipticdg(), libMesh::VariationalSmootherSystem::get_target_to_reference_jacobian(), libMesh::QClough::init_1D(), libMesh::QConical::init_1D(), libMesh::QGrundmann_Moller::init_1D(), init_1D(), libMesh::QNodal::init_1D(), libMesh::QClough::init_2D(), libMesh::QGauss::init_2D(), libMesh::QGaussLobatto::init_2D(), init_2D(), libMesh::QNodal::init_2D(), libMesh::QGauss::init_3D(), libMesh::QGaussLobatto::init_3D(), init_3D(), libMesh::QNodal::init_3D(), libMesh::OverlapCoupling::operator()(), and libMesh::FE< Dim, T >::reinit_default_dual_shape_coeffs().

◆ get_weights() [1/2]

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

Definition at line 174 of file quadrature.h.

174{ return _weights; }
libMesh::QBase::_weights
std::vector< Real > _weights
The quadrature weights.
Definition quadrature.h:415

References libMesh::QBase::_weights.

◆ get_weights() [2/2]

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

Definition at line 168 of file quadrature.h.

168{ return _weights; }

References libMesh::QBase::_weights.

Referenced by libMesh::VariationalSmootherSystem::compute_mesh_quality_info(), libMesh::VariationalSmootherSystem::element_time_derivative(), libMesh::VariationalSmootherSystem::get_target_to_reference_jacobian(), libMesh::QClough::init_1D(), libMesh::QConical::init_1D(), libMesh::QGrundmann_Moller::init_1D(), init_1D(), libMesh::QNodal::init_1D(), libMesh::QClough::init_2D(), libMesh::QGauss::init_2D(), libMesh::QGaussLobatto::init_2D(), init_2D(), libMesh::QNodal::init_2D(), libMesh::QGauss::init_3D(), libMesh::QGaussLobatto::init_3D(), init_3D(), libMesh::QNodal::init_3D(), libMesh::OverlapCoupling::operator()(), and libMesh::FE< Dim, T >::reinit_default_dual_shape_coeffs().

◆ increment_constructor_count()

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

Increments the construction counter.

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

Definition at line 183 of file reference_counter.h.

184{
185 libmesh_try
186 {
187 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
188 std::pair<unsigned int, unsigned int> & p = _counts[name];
189 p.first++;
190 }
191 libmesh_catch (...)
192 {
193 auto stream = libMesh::err.get();
194 stream->exceptions(stream->goodbit); // stream must not throw
195 libMesh::err << "Encountered unrecoverable error while calling "
196 << "ReferenceCounter::increment_constructor_count() "
197 << "for a(n) " << name << " object." << std::endl;
198 std::terminate();
199 }
200}
libMesh::BasicOStreamProxy::get
streamT * get()
Rather than implement every ostream/ios/ios_base function, we'll be lazy and make esoteric uses go th...
Definition ostream_proxy.h:217
libMesh::Threads::spin_mtx
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition threads.C:30
libMesh::err
OStreamProxy err

References libMesh::err, libMesh::BasicOStreamProxy< charT, traits >::get(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< T >::ReferenceCountedObject(), libMesh::ReferenceCountedObject< T >::ReferenceCountedObject(), and libMesh::ReferenceCountedObject< T >::ReferenceCountedObject().

◆ increment_destructor_count()

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

Increments the destruction counter.

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

Definition at line 207 of file reference_counter.h.

208{
209 libmesh_try
210 {
211 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
212 std::pair<unsigned int, unsigned int> & p = _counts[name];
213 p.second++;
214 }
215 libmesh_catch (...)
216 {
217 auto stream = libMesh::err.get();
218 stream->exceptions(stream->goodbit); // stream must not throw
219 libMesh::err << "Encountered unrecoverable error while calling "
220 << "ReferenceCounter::increment_destructor_count() "
221 << "for a(n) " << name << " object." << std::endl;
222 std::terminate();
223 }
224}

References libMesh::err, libMesh::BasicOStreamProxy< charT, traits >::get(), and libMesh::Threads::spin_mtx.

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

◆ init() [1/4]

void libMesh::QBase::init ( const Elem &  e,
unsigned int  p_level = invalid_uint 
)
virtualinherited

Initializes the data structures for a quadrature rule for the element e.

If p_level is specified it overrides the element p_level() elevation to use.

Definition at line 65 of file quadrature.C.

67{
68 libmesh_assert_equal_to(elem.dim(), _dim);
69
70 // Default to the element p_level() value
71 if (p == invalid_uint)
72 p = elem.p_level();
73
74 ElemType t = elem.type();
75
76 // check to see if we have already
77 // done the work for this quadrature rule
78 //
79 // If we have something like a Polygon subclass then we're going to
80 // need to recompute to be safe; even if we're using the same
81 // element, it might have been distorted enough that its subtriangle
82 // triangulation has been changed.
83 if (t == _type && p == _p_level && !elem.runtime_topology())
84 return;
85 else
86 {
87 _elem = &elem;
88 _type = t;
89 _p_level = p;
90 }
91
92 switch(_elem->dim())
93 {
94 case 0:
95 this->init_0D();
96
97 return;
98
99 case 1:
100 this->init_1D();
101
102 return;
103
104 case 2:
105 this->init_2D();
106
107 return;
108
109 case 3:
110 this->init_3D();
111
112 return;
113
114 default:
115 libmesh_error_msg("Invalid dimension _dim = " << _dim);
116 }
117}
libMesh::Elem::dim
virtual unsigned short dim() const =0
libMesh::QBase::init_2D
virtual void init_2D()
Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition quadrature.C:208
libMesh::QBase::init_0D
virtual void init_0D()
Initializes the 0D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition quadrature.C:198
libMesh::QBase::init_1D
virtual void init_1D()=0
Initializes the 1D quadrature rule by filling the points and weights vectors with the appropriate val...
libMesh::QBase::_elem
const Elem * _elem
The element for which the current values were computed, or nullptr if values were computed without a ...
Definition quadrature.h:397
libMesh::QBase::init_3D
virtual void init_3D()
Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition quadrature.C:215
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition libmesh.h:303

References libMesh::QBase::_dim, libMesh::QBase::_elem, libMesh::QBase::_p_level, libMesh::QBase::_type, libMesh::Elem::dim(), libMesh::QBase::init_0D(), libMesh::QBase::init_1D(), libMesh::QBase::init_2D(), libMesh::QBase::init_3D(), libMesh::invalid_uint, libMesh::Elem::p_level(), libMesh::Elem::runtime_topology(), and libMesh::Elem::type().

Referenced by libMesh::QBase::init(), libMesh::QBase::init(), libMesh::QClough::init_1D(), libMesh::QConical::init_1D(), init_1D(), libMesh::QNodal::init_1D(), libMesh::QClough::init_2D(), libMesh::QGauss::init_2D(), libMesh::QGaussLobatto::init_2D(), init_2D(), libMesh::QNodal::init_2D(), libMesh::QGauss::init_3D(), libMesh::QGaussLobatto::init_3D(), libMesh::QGrid::init_3D(), init_3D(), libMesh::QNodal::init_3D(), libMesh::QSimpson::init_3D(), libMesh::QTrap::init_3D(), libMesh::OverlapCoupling::operator()(), libMesh::QClough::QClough(), libMesh::QConical::QConical(), libMesh::QGauss::QGauss(), libMesh::QGaussLobatto::QGaussLobatto(), libMesh::QGrid::QGrid(), libMesh::QGrundmann_Moller::QGrundmann_Moller(), libMesh::QJacobi::QJacobi(), QMonomial(), libMesh::QNodal::QNodal(), libMesh::QSimpson::QSimpson(), libMesh::QTrap::QTrap(), and libMesh::FE< Dim, T >::reinit_default_dual_shape_coeffs().

◆ init() [2/4]

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

Initializes the data structures for an element potentially "cut" by a signed distance function.

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

Definition at line 188 of file quadrature.C.

191{
192 // dispatch generic implementation
193 this->init(elem.type(), p_level);
194}

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

◆ init() [3/4]

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

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

Some types, such as Polygon subclasses, might require more detailed element information and so might not be compatible with this API.

New code should use the Elem-based API for most use cases, but some code may initialize quadrature rules with simple ElemType values like triangles and edges for use in tensor product or conical product constructions; this code can set the simple_type_only flag to avoid being identified as deprecated.

Definition at line 121 of file quadrature.C.

124{
125 // Some element types require data from a specific element, so can
126 // only be used with newer APIs.
127 if (t == C0POLYGON || t == C0POLYHEDRON)
128 libmesh_error_msg("Code (see stack trace) used an outdated quadrature function overload.\n"
129 "Quadrature rules on a C0Polygon are not defined by its ElemType alone.");
130
131 // This API is dangerous to use on general meshes, which may include
132 // element types where the desired quadrature depends on the
133 // physical element, but we still want to be able to initialize
134 // based on only a type for certain simple cases
135 if (t != EDGE2 && !simple_type_only)
136 libmesh_deprecated();
137
138 // check to see if we have already
139 // done the work for this quadrature rule
140 if (t == _type && p == _p_level)
141 return;
142 else
143 {
144 _elem = nullptr;
145 _type = t;
146 _p_level = p;
147 }
148
149 switch(_dim)
150 {
151 case 0:
152 this->init_0D();
153
154 return;
155
156 case 1:
157 this->init_1D();
158
159 return;
160
161 case 2:
162 this->init_2D();
163
164 return;
165
166 case 3:
167 this->init_3D();
168
169 return;
170
171 default:
172 libmesh_error_msg("Invalid dimension _dim = " << _dim);
173 }
174}

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

◆ init() [4/4]

void libMesh::QBase::init ( const QBase &  other_rule)
virtualinherited

Initializes the data structures for a quadrature rule based on the element, element type, and p_level settings of other_rule.

Definition at line 178 of file quadrature.C.

179{
180 if (other_rule._elem)
181 this->init(*other_rule._elem, other_rule._p_level);
182 else
183 this->init(other_rule._type, other_rule._p_level, true);
184}

References libMesh::QBase::_elem, libMesh::QBase::_p_level, libMesh::QBase::_type, and libMesh::QBase::init().

◆ init_0D()

void libMesh::QBase::init_0D ( )
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 198 of file quadrature.C.

199{
200 _points.resize(1);
201 _weights.resize(1);
202 _points[0] = Point(0.);
203 _weights[0] = 1.0;
204}

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

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

◆ init_1D()

void libMesh::QMonomial::init_1D ( )
overrideprivatevirtual

Uses a Gauss rule in 1D.

More efficient rules for non tensor product bases on quadrilaterals and hexahedra.

Implements libMesh::QBase.

Definition at line 29 of file quadrature_monomial_1D.C.

30{
31 QGauss gauss_rule(1, _order);
32 gauss_rule.init(*this);
33
34 _points.swap(gauss_rule.get_points());
35 _weights.swap(gauss_rule.get_weights());
36}

References libMesh::QBase::_order, libMesh::QBase::_points, libMesh::QBase::_weights, libMesh::QBase::get_points(), libMesh::QBase::get_weights(), and libMesh::QBase::init().

◆ init_2D()

void libMesh::QMonomial::init_2D ( )
overrideprivatevirtual

Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate values.

The order of the rule will be defined by the implementing class. Should not be pure virtual since a derived quadrature rule may only be defined in 1D. If not overridden, throws an error.

Reimplemented from libMesh::QBase.

Definition at line 28 of file quadrature_monomial_2D.C.

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

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

◆ init_3D()

void libMesh::QMonomial::init_3D ( )
overrideprivatevirtual

Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate values.

The order of the rule will be defined by the implementing class. Should not be pure virtual since a derived quadrature rule may only be defined in 1D. If not overridden, throws an error.

Reimplemented from libMesh::QBase.

Definition at line 28 of file quadrature_monomial_3D.C.

29{
30
31 switch (_type)
32 {
33 //---------------------------------------------
34 // Hex quadrature rules
35 case HEX8:
36 case HEX20:
37 case HEX27:
38 {
39 switch(get_order())
40 {
41
42 // The CONSTANT/FIRST rule is the 1-point Gauss "product" rule...we fall
43 // through to the default case for this rule.
44
45 case SECOND:
46 case THIRD:
47 {
48 // A degree 3, 6-point, "rotationally-symmetric" rule by
49 // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
50 //
51 // Warning: this rule contains points on the boundary of the reference
52 // element, and therefore may be unsuitable for some problems. The alternative
53 // would be a 2x2x2 Gauss product rule.
54 const Real data[1][4] =
55 {
56 {1, 0, 0, Real(4)/3}
57 };
58
59 const unsigned int rule_id[1] = {
60 1 // (x,0,0) -> 6 permutations
61 };
62
63 _points.resize(6);
64 _weights.resize(6);
65
66 kim_rule(data, rule_id, 1);
67 return;
68 } // end case SECOND,THIRD
69
70 case FOURTH:
71 case FIFTH:
72 {
73 // A degree 5, 13-point rule by Stroud,
74 // AH Stroud, "Some Fifth Degree Integration Formulas for Symmetric Regions II.",
75 // Numerische Mathematik 9, pp. 460-468 (1967).
76 //
77 // This rule is provably minimal in the number of points. The equations given for
78 // the n-cube on pg. 466 of the paper for mu/gamma and gamma are wrong, at least for
79 // the n=3 case. The analytical values given here were computed by me [JWP] in Maple.
80
81 // Convenient intermediate values.
82 const Real sqrt19 = std::sqrt(Real(19)); // ~4.35889894354
83 const Real tp = std::sqrt(71440 + 6802*sqrt19); // ~317.945326454
84
85 // Point data for permutations.
86 const Real eta = 0;
87
88 const Real lambda = std::sqrt(Real(1919)/3285 - 148*sqrt19/3285 + 4*tp/3285);
89 // 8.8030440669930978047737818209860e-01_R;
90
91 const Real xi = -std::sqrt(Real(1121)/3285 + 74*sqrt19/3285 - 2*tp/3285);
92 // -4.9584817142571115281421242364290e-01_R;
93
94 const Real mu = std::sqrt(Real(1121)/3285 + 74*sqrt19/3285 + 2*tp/3285);
95 // 7.9562142216409541542982482567580e-01_R;
96
97 const Real gamma = std::sqrt(Real(1919)/3285 - 148*sqrt19/3285 - 4*tp/3285);
98 // 2.5293711744842581347389255929324e-02_R;
99
100 // Weights: the centroid weight is given analytically. Weight B (resp C) goes
101 // with the {lambda,xi} (resp {gamma,mu}) permutation. The single-precision
102 // results reported by Stroud are given for reference.
103
104 const Real A = Real(32)/19; // ~1.684210560
105 // Stroud: 0.21052632 * 8.0;
106
107 const Real B = Real(1) / (Real(260072)/133225 - 1520*sqrt19/133225 + (133 - 37*sqrt19)*tp/133225);
108 // 5.4498735127757671684690782180890e-01_R; // Stroud: 0.068123420 * 8.0 = 0.544987360;
109
110 const Real C = Real(1) / (Real(260072)/133225 - 1520*sqrt19/133225 - (133 - 37*sqrt19)*tp/133225);
111 // 5.0764422766979170420572375713840e-01_R; // Stroud: 0.063455527 * 8.0 = 0.507644216;
112
113 _points.resize(13);
114 _weights.resize(13);
115
116 unsigned int c=0;
117
118 // Point with weight A (origin)
119 _points[c] = Point(eta, eta, eta);
120 _weights[c++] = A;
121
122 // Points with weight B
123 _points[c] = Point(lambda, xi, xi);
124 _weights[c++] = B;
125 _points[c] = -_points[c-1];
126 _weights[c++] = B;
127
128 _points[c] = Point(xi, lambda, xi);
129 _weights[c++] = B;
130 _points[c] = -_points[c-1];
131 _weights[c++] = B;
132
133 _points[c] = Point(xi, xi, lambda);
134 _weights[c++] = B;
135 _points[c] = -_points[c-1];
136 _weights[c++] = B;
137
138 // Points with weight C
139 _points[c] = Point(mu, mu, gamma);
140 _weights[c++] = C;
141 _points[c] = -_points[c-1];
142 _weights[c++] = C;
143
144 _points[c] = Point(mu, gamma, mu);
145 _weights[c++] = C;
146 _points[c] = -_points[c-1];
147 _weights[c++] = C;
148
149 _points[c] = Point(gamma, mu, mu);
150 _weights[c++] = C;
151 _points[c] = -_points[c-1];
152 _weights[c++] = C;
153
154 return;
155
156
157 // // A degree 5, 14-point, "rotationally-symmetric" rule by
158 // // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
159 // // Was also reported in Stroud's 1971 book.
160 // const Real data[2][4] =
161 // {
162 // {7.95822425754221463264548820476135e-01_R, 0.00000000000000000000000000000000e+00_R, 0.00000000000000000000000000000000e+00_R, 8.86426592797783933518005540166204e-01_R},
163 // {7.58786910639328146269034278112267e-01_R, 7.58786910639328146269034278112267e-01_R, 7.58786910639328146269034278112267e-01_R, 3.35180055401662049861495844875346e-01_R}
164 // };
165
166 // const unsigned int rule_id[2] = {
167 // 1, // (x,0,0) -> 6 permutations
168 // 4 // (x,x,x) -> 8 permutations
169 // };
170
171 // _points.resize(14);
172 // _weights.resize(14);
173
174 // kim_rule(data, rule_id, 2);
175 // return;
176 } // end case FOURTH,FIFTH
177
178
179 case SIXTH:
180 case SEVENTH:
181 {
182 if (allow_rules_with_negative_weights)
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+00_R, 0.00000000000000000000000000000000e+00_R, 0.00000000000000000000000000000000e+00_R, -1.27536231884057971014492753623188e+00_R},
191 {5.85540043769119907612630781744060e-01_R, 0.00000000000000000000000000000000e+00_R, 0.00000000000000000000000000000000e+00_R, 8.71111111111111111111111111111111e-01_R},
192 {6.94470135991704766602025803883310e-01_R, 9.37161638568208038511047377665396e-01_R, 4.15659267604065126239606672567031e-01_R, 1.68695652173913043478260869565217e-01_R}
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 // Mathematical 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(Real(6)/7), // ~0.92582009977
223 s = std::sqrt((Real(960) - 3*std::sqrt(Real(28798))) / 2726), // ~0.40670318642
224 t = std::sqrt((Real(960) + 3*std::sqrt(Real(28798))) / 2726), // ~0.73411252875
225 B1 = Real(8624)/29160, // ~0.29574759945
226 B2 = Real(2744)/29160, // ~0.09410150891
227 B3 = 8*(774*t*t - 230)/(9720*(t*t-s*s)), // ~0.41233386227
228 B4 = 8*(230 - 774*s*s)/(9720*(t*t-s*s)); // ~0.22470317477
229
230 const Real data[4][4] =
231 {
232 {r, 0, 0, B1},
233 {r, r, 0, 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-01_R, 0.00000000000000000000000000000000e+00_R, 0.00000000000000000000000000000000e+00_R, 2.95189738262622903181631100062774e-01_R},
259 // {4.08372221499474674069588900002128e-01_R, 4.08372221499474674069588900002128e-01_R, 4.08372221499474674069588900002128e-01_R, 4.04055417266200582425904380777126e-01_R},
260 // {8.59523090201054193116477875786220e-01_R, 8.59523090201054193116477875786220e-01_R, 4.14735913727987720499709244748633e-01_R, 1.24850759678944080062624098058597e-01_R}
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+00_R, 0.00000000000000000000000000000000e+00_R, 0.00000000000000000000000000000000e+00_R, 4.51903714875199690490763818699555e-01_R},
286 {7.82460796435951590652813975429717e-01_R, 0.00000000000000000000000000000000e+00_R, 0.00000000000000000000000000000000e+00_R, 2.99379177352338919703385618576171e-01_R},
287 {4.88094669706366480526729301468686e-01_R, 4.88094669706366480526729301468686e-01_R, 4.88094669706366480526729301468686e-01_R, 3.00876159371240019939698689791164e-01_R},
288 {8.62218927661481188856422891110042e-01_R, 8.62218927661481188856422891110042e-01_R, 8.62218927661481188856422891110042e-01_R, 4.94843255877038125738173175714853e-02_R},
289 {2.81113909408341856058098281846420e-01_R, 9.44196578292008195318687494773744e-01_R, 6.97574833707236996779391729948984e-01_R, 1.22872389222467338799199767122592e-01_R}
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 libmesh_fallthrough();
320
321 // By default: construct and use a Gauss quadrature rule
322 default:
323 {
324 QGauss gauss_rule(3, _order);
325 gauss_rule.init(*this);
326
327 // Swap points and weights with the about-to-be destroyed rule.
328 _points.swap (gauss_rule.get_points() );
329 _weights.swap(gauss_rule.get_weights());
330
331 return;
332 }
333 } // end switch (_type)
334}
libMesh::QBase::allow_rules_with_negative_weights
bool allow_rules_with_negative_weights
Flag (default true) controlling the use of quadrature rules with negative weights.
Definition quadrature.h:301
libMesh::QMonomial::kim_rule
void kim_rule(const Real rule_data[][4], const unsigned int *rule_id, const unsigned int n_pts)
Rules from Kim and Song, Comm.
Definition quadrature_monomial.C:210
B
Definition assembly.h:39

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

◆ kim_rule()

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, quadrature 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 210 of file quadrature_monomial.C.

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

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

Referenced by init_3D().

◆ n_objects()

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

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

Definition at line 85 of file reference_counter.h.

86 { return _n_objects; }
libMesh::ReferenceCounter::_n_objects
static Threads::atomic< unsigned int > _n_objects
The number of objects.
Definition reference_counter.h:132

References libMesh::ReferenceCounter::_n_objects.

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

◆ n_points()

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

Definition at line 131 of file quadrature.h.

132 {
133 libmesh_assert (!_points.empty());
134 return cast_int<unsigned int>(_points.size());
135 }
libMesh::libmesh_assert
libmesh_assert(ctx)

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

Referenced by libMesh::ExactSolution::_compute_error(), alternative_fe_assembly(), LinearElasticity::assemble(), assemble(), assemble(), assemble_1D(), AssembleOptimization::assemble_A_and_F(), libMesh::ClawSystem::assemble_advection_matrices(), libMesh::ClawSystem::assemble_avg_coupling_matrices(), libMesh::ClawSystem::assemble_boundary_condition_matrices(), assemble_cd(), assemble_cd(), assemble_divgrad(), assemble_elasticity(), assemble_ellipticdg(), assemble_graddiv(), assemble_helmholtz(), libMesh::ClawSystem::assemble_jump_coupling_matrix(), assemble_mass(), libMesh::ClawSystem::assemble_mass_matrix(), assemble_matrices(), assemble_poisson(), assemble_poisson(), assemble_poisson(), assemble_shell(), assemble_stokes(), assemble_temperature_jump(), assemble_wave(), assembly_with_dg_fem_context(), AssemblyF0::boundary_assembly(), AssemblyF1::boundary_assembly(), AssemblyF2::boundary_assembly(), AssemblyA0::boundary_assembly(), AssemblyA1::boundary_assembly(), AssemblyA2::boundary_assembly(), A2::boundary_assembly(), A3::boundary_assembly(), F0::boundary_assembly(), Output0::boundary_assembly(), libMesh::System::calculate_norm(), compute_enriched_soln(), compute_jacobian(), libMesh::FEGenericBase< OutputType >::compute_periodic_constraints(), libMesh::FEGenericBase< OutputType >::compute_proj_constraints(), compute_residual(), libMesh::FirstOrderUnsteadySolver::compute_second_order_eqns(), LinearElasticity::compute_stresses(), LargeDeformationElasticity::compute_stresses(), LinearElasticityWithContact::compute_stresses(), compute_stresses(), libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), libMesh::HDGProblem::create_identity_jacobian(), libMesh::HDGProblem::create_identity_residual(), SecondOrderScalarSystemSecondOrderTimeSolverBase::damping_residual(), SecondOrderScalarSystemFirstOrderTimeSolverBase::damping_residual(), CoupledSystem::element_constraint(), NavierSystem::element_constraint(), LaplaceSystem::element_postprocess(), PoissonSystem::element_postprocess(), LaplaceQoI::element_qoi(), HeatSystem::element_qoi(), HeatSystem::element_qoi_derivative(), LaplaceSystem::element_qoi_derivative(), LaplaceQoI::element_qoi_derivative(), LaplaceSystem::element_time_derivative(), CoupledSystem::element_time_derivative(), HeatSystem::element_time_derivative(), PoissonSystem::element_time_derivative(), NavierSystem::element_time_derivative(), SolidSystem::element_time_derivative(), ElasticitySystem::element_time_derivative(), CurlCurlSystem::element_time_derivative(), SigmaPhysics::element_time_derivative(), FirstOrderScalarSystemBase::element_time_derivative(), SecondOrderScalarSystemFirstOrderTimeSolverBase::element_time_derivative(), HilbertSystem::element_time_derivative(), fe_assembly(), form_functionA(), form_functionB(), form_matrixA(), libMesh::VariationalSmootherSystem::get_target_to_reference_jacobian(), A0::interior_assembly(), A1::interior_assembly(), A2::interior_assembly(), F0::interior_assembly(), OutputAssembly::interior_assembly(), B::interior_assembly(), M0::interior_assembly(), EIM_IP_assembly::interior_assembly(), AssemblyA0::interior_assembly(), AssemblyA1::interior_assembly(), AssemblyA2::interior_assembly(), InnerProductAssembly::interior_assembly(), AssemblyF0::interior_assembly(), AssemblyF1::interior_assembly(), Ex6InnerProduct::interior_assembly(), AcousticsInnerProduct::interior_assembly(), LargeDeformationElasticity::jacobian(), LaplaceYoung::jacobian(), libMesh::FEMPhysics::mass_residual(), NavierSystem::mass_residual(), ElasticitySystem::mass_residual(), FirstOrderScalarSystemBase::mass_residual(), SecondOrderScalarSystemSecondOrderTimeSolverBase::mass_residual(), SecondOrderScalarSystemFirstOrderTimeSolverBase::mass_residual(), libMesh::FEAbstract::n_quadrature_points(), periodic_bc_test_poisson(), libMesh::QBase::print_info(), LargeDeformationElasticity::residual(), LaplaceYoung::residual(), LinearElasticityWithContact::residual_and_jacobian(), LaplaceSystem::side_constraint(), LaplaceSystem::side_postprocess(), CoupledSystemQoI::side_qoi(), LaplaceSystem::side_qoi_derivative(), CoupledSystemQoI::side_qoi_derivative(), ElasticitySystem::side_time_derivative(), CurlCurlSystem::side_time_derivative(), libMesh::QBase::size(), libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), libMesh::QBase::tensor_product_quad(), and QuadratureTest::testPolynomial().

◆ operator=() [1/2]

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

◆ operator=() [2/2]

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

◆ print_info() [1/2]

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

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

Definition at line 43 of file quadrature.C.

44{
45 libmesh_assert(!_points.empty());
46 libmesh_assert(!_weights.empty());
47
48 Real summed_weights=0;
49 os << "N_Q_Points=" << this->n_points() << std::endl << std::endl;
50 for (auto qpoint: index_range(_points))
51 {
52 os << " Point " << qpoint << ":\n"
53 << " "
54 << _points[qpoint]
55 << "\n Weight:\n "
56 << " w=" << _weights[qpoint] << "\n" << std::endl;
57
58 summed_weights += _weights[qpoint];
59 }
60 os << "Summed Weights: " << summed_weights << std::endl;
61}
libMesh::QBase::n_points
unsigned int n_points() const
Definition quadrature.h:131
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition int_range.h:153

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

◆ print_info() [2/2]

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

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

Definition at line 81 of file reference_counter.C.

82{
83 if (_enable_print_counter)
84 out_stream << ReferenceCounter::get_info();
85}
libMesh::ReferenceCounter::get_info
static std::string get_info()
Gets a string containing the reference information.
Definition reference_counter.C:47

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

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

◆ qp()

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

Definition at line 179 of file quadrature.h.

180 {
181 libmesh_assert_less (i, _points.size());
182 return _points[i];
183 }

References libMesh::QBase::_points.

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

◆ scale()

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

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

Definition at line 222 of file quadrature.C.

224{
225 // Make sure we are in 1D
226 libmesh_assert_equal_to (_dim, 1);
227
228 Real
229 h_new = new_range.second - new_range.first,
230 h_old = old_range.second - old_range.first;
231
232 // Make sure that we have sane ranges
233 libmesh_assert_greater (h_new, 0.);
234 libmesh_assert_greater (h_old, 0.);
235
236 // Make sure there are some points
237 libmesh_assert_greater (_points.size(), 0);
238
239 // Compute the scale factor
240 Real scfact = h_new/h_old;
241
242 // We're mapping from old_range -> new_range
243 for (auto i : index_range(_points))
244 {
245 _points[i](0) = new_range.first +
246 (_points[i](0) - old_range.first) * scfact;
247
248 // Scale the weights
249 _weights[i] *= scfact;
250 }
251}

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

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

◆ shapes_need_reinit()

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

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

Definition at line 286 of file quadrature.h.

286{ return false; }

◆ size()

unsigned int libMesh::QBase::size ( ) const
inlineinherited

Alias for n_points() to enable use in index_range.

Returns
The number of points associated with the quadrature rule.

Definition at line 142 of file quadrature.h.

143 {
144 return n_points();
145 }

References libMesh::QBase::n_points().

◆ stroud_rule()

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 62 of file quadrature_monomial.C.

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

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

Referenced by init_2D().

◆ tensor_product_hex()

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

Computes the tensor product quadrature rule [q1D x q1D x q1D] from the 1D rule q1D.

Used in the init_3D routines for hexahedral element types.

Definition at line 283 of file quadrature.C.

284{
285 const unsigned int np = q1D.n_points();
286
287 _points.resize(np * np * np);
288
289 _weights.resize(np * np * np);
290
291 unsigned int q=0;
292
293 for (unsigned int k=0; k<np; k++)
294 for (unsigned int j=0; j<np; j++)
295 for (unsigned int i=0; i<np; i++)
296 {
297 _points[q](0) = q1D.qp(i)(0);
298 _points[q](1) = q1D.qp(j)(0);
299 _points[q](2) = q1D.qp(k)(0);
300
301 _weights[q] = q1D.w(i) * q1D.w(j) * q1D.w(k);
302
303 q++;
304 }
305}

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

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

◆ tensor_product_prism()

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

Computes the tensor product of a 1D quadrature rule and a 2D quadrature rule.

Used in the init_3D routines for prismatic element types.

Definition at line 310 of file quadrature.C.

311{
312 const unsigned int n_points1D = q1D.n_points();
313 const unsigned int n_points2D = q2D.n_points();
314
315 _points.resize (n_points1D * n_points2D);
316 _weights.resize (n_points1D * n_points2D);
317
318 unsigned int q=0;
319
320 for (unsigned int j=0; j<n_points1D; j++)
321 for (unsigned int i=0; i<n_points2D; i++)
322 {
323 _points[q](0) = q2D.qp(i)(0);
324 _points[q](1) = q2D.qp(i)(1);
325 _points[q](2) = q1D.qp(j)(0);
326
327 _weights[q] = q2D.w(i) * q1D.w(j);
328
329 q++;
330 }
331
332}

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

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

◆ tensor_product_quad()

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

Constructs a 2D rule from the tensor product of q1D with itself.

Used in the init_2D() routines for quadrilateral element types.

Definition at line 256 of file quadrature.C.

257{
258
259 const unsigned int np = q1D.n_points();
260
261 _points.resize(np * np);
262
263 _weights.resize(np * np);
264
265 unsigned int q=0;
266
267 for (unsigned int j=0; j<np; j++)
268 for (unsigned int i=0; i<np; i++)
269 {
270 _points[q](0) = q1D.qp(i)(0);
271 _points[q](1) = q1D.qp(j)(0);
272
273 _weights[q] = q1D.w(i)*q1D.w(j);
274
275 q++;
276 }
277}

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

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

◆ type()

QuadratureType libMesh::QMonomial::type ( ) const
overridevirtual
Returns
QMONOMIAL.

Implements libMesh::QBase.

Definition at line 32 of file quadrature_monomial.C.

33{
34 return QMONOMIAL;
35}

References libMesh::QMONOMIAL.

◆ w()

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

Definition at line 188 of file quadrature.h.

189 {
190 libmesh_assert_less (i, _weights.size());
191 return _weights[i];
192 }

References libMesh::QBase::_weights.

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

◆ wissmann_rule()

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 42 of file quadrature_monomial.C.

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

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

Referenced by init_2D().

Member Data Documentation

◆ _counts

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

Actually holds the data.

Definition at line 124 of file reference_counter.h.

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

◆ _dim

unsigned int libMesh::QBase::_dim
protectedinherited

The spatial dimension of the quadrature rule.

Definition at line 379 of file quadrature.h.

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

◆ _elem

const Elem* libMesh::QBase::_elem
protectedinherited

The element for which the current values were computed, or nullptr if values were computed without a specific element.

Definition at line 397 of file quadrature.h.

Referenced by libMesh::QBase::init(), libMesh::QBase::init(), libMesh::QBase::init(), libMesh::QGauss::init_2D(), libMesh::QNodal::init_2D(), libMesh::QGauss::init_3D(), and libMesh::QNodal::init_3D().

◆ _enable_print_counter

bool libMesh::ReferenceCounter::_enable_print_counter = true
staticprotectedinherited

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

Definition at line 143 of file reference_counter.h.

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

◆ _mutex

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

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 137 of file reference_counter.h.

◆ _n_objects

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

The number of objects.

Print the reference count information when the number returns to 0.

Definition at line 132 of file reference_counter.h.

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

◆ _order

Order libMesh::QBase::_order
protectedinherited

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

Definition at line 385 of file quadrature.h.

Referenced by libMesh::QBase::build(), libMesh::QBase::build(), libMesh::QBase::get_base_order(), libMesh::QBase::get_order(), libMesh::QClough::init_1D(), libMesh::QGauss::init_1D(), libMesh::QGaussLobatto::init_1D(), libMesh::QGrid::init_1D(), libMesh::QJacobi::init_1D(), init_1D(), libMesh::QNodal::init_1D(), libMesh::QClough::init_2D(), libMesh::QGauss::init_2D(), libMesh::QGaussLobatto::init_2D(), libMesh::QGrid::init_2D(), init_2D(), libMesh::QNodal::init_2D(), libMesh::QGauss::init_3D(), libMesh::QGaussLobatto::init_3D(), libMesh::QGrid::init_3D(), init_3D(), and libMesh::QNodal::init_3D().

◆ _p_level

unsigned int libMesh::QBase::_p_level
protectedinherited

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

Definition at line 403 of file quadrature.h.

Referenced by libMesh::QBase::get_order(), libMesh::QBase::get_p_level(), libMesh::QBase::init(), libMesh::QBase::init(), libMesh::QBase::init(), libMesh::QClough::init_2D(), libMesh::QGauss::init_2D(), libMesh::QGauss::init_3D(), libMesh::QSimpson::init_3D(), and libMesh::QTrap::init_3D().

◆ _points

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

The locations of the quadrature points in reference element space.

Definition at line 409 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::QBase::get_points(), libMesh::QGrundmann_Moller::gm_rule(), libMesh::QBase::init_0D(), libMesh::QClough::init_1D(), libMesh::QConical::init_1D(), libMesh::QGauss::init_1D(), libMesh::QGaussLobatto::init_1D(), libMesh::QGrundmann_Moller::init_1D(), libMesh::QGrid::init_1D(), libMesh::QJacobi::init_1D(), init_1D(), libMesh::QNodal::init_1D(), libMesh::QSimpson::init_1D(), libMesh::QTrap::init_1D(), libMesh::QClough::init_2D(), libMesh::QGauss::init_2D(), libMesh::QGaussLobatto::init_2D(), libMesh::QGrid::init_2D(), init_2D(), libMesh::QNodal::init_2D(), libMesh::QSimpson::init_2D(), libMesh::QTrap::init_2D(), libMesh::QGauss::init_3D(), libMesh::QGaussLobatto::init_3D(), libMesh::QGrundmann_Moller::init_3D(), libMesh::QGrid::init_3D(), init_3D(), libMesh::QNodal::init_3D(), libMesh::QSimpson::init_3D(), libMesh::QTrap::init_3D(), libMesh::QGauss::keast_rule(), kim_rule(), libMesh::QBase::n_points(), libMesh::QBase::print_info(), libMesh::QBase::qp(), libMesh::QBase::scale(), stroud_rule(), libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), libMesh::QBase::tensor_product_quad(), and wissmann_rule().

◆ _type

ElemType libMesh::QBase::_type
protectedinherited

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

Definition at line 391 of file quadrature.h.

Referenced by libMesh::QBase::get_elem_type(), libMesh::QBase::init(), libMesh::QBase::init(), libMesh::QBase::init(), libMesh::QNodal::init_1D(), libMesh::QClough::init_2D(), libMesh::QConical::init_2D(), libMesh::QGauss::init_2D(), libMesh::QGaussLobatto::init_2D(), libMesh::QGrundmann_Moller::init_2D(), libMesh::QGrid::init_2D(), init_2D(), libMesh::QNodal::init_2D(), libMesh::QSimpson::init_2D(), libMesh::QTrap::init_2D(), libMesh::QClough::init_3D(), libMesh::QConical::init_3D(), libMesh::QGauss::init_3D(), libMesh::QGaussLobatto::init_3D(), libMesh::QGrundmann_Moller::init_3D(), libMesh::QGrid::init_3D(), init_3D(), libMesh::QNodal::init_3D(), libMesh::QSimpson::init_3D(), and libMesh::QTrap::init_3D().

◆ _weights

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

The quadrature weights.

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

Definition at line 415 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::QBase::get_weights(), libMesh::QGrundmann_Moller::gm_rule(), libMesh::QBase::init_0D(), libMesh::QClough::init_1D(), libMesh::QConical::init_1D(), libMesh::QGauss::init_1D(), libMesh::QGaussLobatto::init_1D(), libMesh::QGrundmann_Moller::init_1D(), libMesh::QGrid::init_1D(), libMesh::QJacobi::init_1D(), init_1D(), libMesh::QNodal::init_1D(), libMesh::QSimpson::init_1D(), libMesh::QTrap::init_1D(), libMesh::QClough::init_2D(), libMesh::QGauss::init_2D(), libMesh::QGaussLobatto::init_2D(), libMesh::QGrid::init_2D(), init_2D(), libMesh::QNodal::init_2D(), libMesh::QSimpson::init_2D(), libMesh::QTrap::init_2D(), libMesh::QGauss::init_3D(), libMesh::QGaussLobatto::init_3D(), libMesh::QGrundmann_Moller::init_3D(), libMesh::QGrid::init_3D(), init_3D(), libMesh::QNodal::init_3D(), libMesh::QSimpson::init_3D(), libMesh::QTrap::init_3D(), libMesh::QGauss::keast_rule(), kim_rule(), libMesh::QBase::print_info(), libMesh::QBase::scale(), stroud_rule(), libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), libMesh::QBase::tensor_product_quad(), libMesh::QBase::w(), and wissmann_rule().

◆ allow_nodal_pyramid_quadrature

bool libMesh::QBase::allow_nodal_pyramid_quadrature
inherited

The flag's value defaults to false so that one does not accidentally use a nodal quadrature rule on Pyramid elements, since evaluating the inverse element Jacobian (e.g.

dphi) is not well-defined at the Pyramid apex because the element Jacobian is zero there.

We do not want to completely prevent someone from using a nodal quadrature rule on Pyramids, however, since there are legitimate use cases (lumped mass matrix) so the flag can be set to true to override this behavior.

Definition at line 314 of file quadrature.h.

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

◆ allow_rules_with_negative_weights

bool libMesh::QBase::allow_rules_with_negative_weights
inherited

Flag (default true) controlling the use of quadrature rules with negative weights.

Set this to false to require rules with all positive weights.

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

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

Definition at line 301 of file quadrature.h.

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

The documentation for this class was generated from the following files:
Generated on Mon Aug 24 2026 19:01:08 for libMesh by doxygen 1.9.8