libMesh
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. More...
 
 QMonomial (const QMonomial &)=default
 Copy/move ctor, copy/move assignment operator, and destructor are all explicitly defaulted for this simple class. More...
 
 QMonomial (QMonomial &&)=default
 
QMonomialoperator= (const QMonomial &)=default
 
QMonomialoperator= (QMonomial &&)=default
 
virtual ~QMonomial ()=default
 
virtual QuadratureType type () const override
 
virtual std::unique_ptr< QBaseclone () 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. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
virtual bool shapes_need_reinit ()
 

Static Public Member Functions

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

Protected Types

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

Protected Member Functions

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

Protected Attributes

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

Static Protected Attributes

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

Private Member Functions

virtual void init_1D () override
 Uses a Gauss rule in 1D. More...
 
virtual void init_2D () override
 Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate values. More...
 
virtual void init_3D () override
 Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate values. More...
 
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. More...
 
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. More...
 
void kim_rule (const Real rule_data[][4], const unsigned int *rule_id, const unsigned int n_pts)
 Rules from Kim and Song, Comm. More...
 

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.

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

66  :
67  QBase(dim,order)
68  {
69  if (dim == 1)
70  init(EDGE2);
71  }
unsigned int dim
QBase(unsigned int dim, Order order=INVALID_ORDER)
Constructor.
Definition: quadrature.C:27
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

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

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

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

46 {
47  return QBase::build (Utility::string_to_enum<QuadratureType> (type),
48  _dim,
49  _order);
50 }
unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:379
virtual QuadratureType type() const =0
Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:385
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.

◆ build() [2/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.

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.

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 }
unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:379
Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:385
OStreamProxy out

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

References libMesh::ReferenceCounter::_enable_print_counter.

101 {
102  _enable_print_counter = false;
103  return;
104 }
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called...

◆ enable_print_counter_info()

void libMesh::ReferenceCounter::enable_print_counter_info ( )
staticinherited

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

Definition at line 94 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

95 {
96  _enable_print_counter = true;
97  return;
98 }
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called...

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

References libMesh::QBase::_order.

258 { return static_cast<Order>(_order); }
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:385

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

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

150 { return _dim; }
unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:379

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

References libMesh::QBase::_type.

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

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

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

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

48 {
49 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
50 
51  std::ostringstream oss;
52 
53  oss << '\n'
54  << " ---------------------------------------------------------------------------- \n"
55  << "| Reference count information |\n"
56  << " ---------------------------------------------------------------------------- \n";
57 
58  for (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 }
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
static Counts _counts
Actually holds the data.

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

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::QGaussLobatto::init_1D(), libMesh::QGauss::init_1D(), libMesh::QConical::init_1D(), libMesh::QJacobi::init_1D(), libMesh::QGrundmann_Moller::init_1D(), libMesh::QGaussLobatto::init_2D(), libMesh::QGauss::init_2D(), init_2D(), libMesh::QGrundmann_Moller::init_2D(), libMesh::QGaussLobatto::init_3D(), libMesh::QGauss::init_3D(), init_3D(), libMesh::QGrundmann_Moller::init_3D(), and libMesh::RBParametrizedFunction::preevaluate_parametrized_function_on_mesh().

249 { return static_cast<Order>(_order + 2 * _p_level); }
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
unsigned int _p_level
The p-level of the element for which the current values have been computed.
Definition: quadrature.h:403
Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:385

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

References libMesh::QBase::_p_level.

126 { return _p_level; }
unsigned int _p_level
The p-level of the element for which the current values have been computed.
Definition: quadrature.h:403

◆ get_points() [1/2]

const std::vector<Point>& libMesh::QBase::get_points ( ) const
inlineinherited

◆ get_points() [2/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.

References libMesh::QBase::_points.

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

◆ get_weights() [1/2]

const std::vector<Real>& libMesh::QBase::get_weights ( ) const
inlineinherited

◆ get_weights() [2/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.

References libMesh::QBase::_weights.

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

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

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

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

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 }
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
OStreamProxy err
static Counts _counts
Actually holds the data.
streamT * get()
Rather than implement every ostream/ios/ios_base function, we&#39;ll be lazy and make esoteric uses go th...
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:30

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

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

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

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 }
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
OStreamProxy err
static Counts _counts
Actually holds the data.
streamT * get()
Rather than implement every ostream/ios/ios_base function, we&#39;ll be lazy and make esoteric uses go th...
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:30

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

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::QClough::init_1D(), libMesh::QConical::init_1D(), libMesh::QNodal::init_1D(), init_1D(), libMesh::QGaussLobatto::init_2D(), libMesh::QClough::init_2D(), libMesh::QGauss::init_2D(), libMesh::QNodal::init_2D(), init_2D(), libMesh::QGaussLobatto::init_3D(), libMesh::QGauss::init_3D(), libMesh::QSimpson::init_3D(), libMesh::QTrap::init_3D(), libMesh::QGrid::init_3D(), libMesh::QNodal::init_3D(), 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, LAGRANGE_VEC >::reinit_default_dual_shape_coeffs().

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 }
ElemType
Defines an enum for geometric element types.
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:310
virtual void init_2D()
Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:208
ElemType _type
The type of element for which the current values have been computed.
Definition: quadrature.h:391
unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:379
const Elem * _elem
The element for which the current values were computed, or nullptr if values were computed without a ...
Definition: quadrature.h:397
unsigned int _p_level
The p-level of the element for which the current values have been computed.
Definition: quadrature.h:403
virtual void init_0D()
Initializes the 0D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:198
virtual void init_1D()=0
Initializes the 1D quadrature rule by filling the points and weights vectors with the appropriate val...
virtual unsigned short dim() const =0
virtual void init_3D()
Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:215

◆ init() [2/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.

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

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 }
virtual void init_2D()
Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:208
ElemType _type
The type of element for which the current values have been computed.
Definition: quadrature.h:391
unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:379
const Elem * _elem
The element for which the current values were computed, or nullptr if values were computed without a ...
Definition: quadrature.h:397
unsigned int _p_level
The p-level of the element for which the current values have been computed.
Definition: quadrature.h:403
virtual void init_0D()
Initializes the 0D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:198
virtual void init_1D()=0
Initializes the 1D quadrature rule by filling the points and weights vectors with the appropriate val...
virtual void init_3D()
Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition: quadrature.C:215

◆ init() [3/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.

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

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 }
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

◆ init() [4/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.

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

191 {
192  // dispatch generic implementation
193  this->init(elem.type(), p_level);
194 }
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

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

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

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

199 {
200  _points.resize(1);
201  _weights.resize(1);
202  _points[0] = Point(0.);
203  _weights[0] = 1.0;
204 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415

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

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

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 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415
Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:385

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

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

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 }
void wissmann_rule(const Real rule_data[][3], const unsigned int n_pts)
Wissmann published three interesting "partially symmetric" rules for integrating degree 4...
ElemType _type
The type of element for which the current values have been computed.
Definition: quadrature.h:391
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415
Order get_order() const
Definition: quadrature.h:249
void stroud_rule(const Real rule_data[][3], const unsigned int *rule_symmetry, const unsigned int n_pts)
Stroud&#39;s rules for quads and hexes can have one of several different types of symmetry.
Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:385
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

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

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.

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  {
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 }
bool allow_rules_with_negative_weights
Flag (default true) controlling the use of quadrature rules with negative weights.
Definition: quadrature.h:301
ElemType _type
The type of element for which the current values have been computed.
Definition: quadrature.h:391
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415
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: assembly.h:38
Order get_order() const
Definition: quadrature.h:249
Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:385
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

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

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

Referenced by init_3D().

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 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

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

References libMesh::ReferenceCounter::_n_objects.

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

86  { return _n_objects; }
static Threads::atomic< unsigned int > _n_objects
The number of objects.

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

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

Referenced by libMesh::ExactSolution::_compute_error(), assemble(), assemble_1D(), assemble_ellipticdg(), assemble_helmholtz(), assemble_poisson(), assemble_shell(), assemble_wave(), assembly_with_dg_fem_context(), AssemblyA0::boundary_assembly(), AssemblyA1::boundary_assembly(), AssemblyF0::boundary_assembly(), AssemblyF1::boundary_assembly(), AssemblyA2::boundary_assembly(), AssemblyF2::boundary_assembly(), A2::boundary_assembly(), A3::boundary_assembly(), F0::boundary_assembly(), Output0::boundary_assembly(), compute_enriched_soln(), libMesh::FirstOrderUnsteadySolver::compute_second_order_eqns(), 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(), NavierSystem::element_constraint(), CoupledSystem::element_constraint(), PoissonSystem::element_postprocess(), LaplaceSystem::element_postprocess(), LaplaceQoI::element_qoi(), HeatSystem::element_qoi(), LaplaceQoI::element_qoi_derivative(), LaplaceSystem::element_qoi_derivative(), HeatSystem::element_qoi_derivative(), NavierSystem::element_time_derivative(), SolidSystem::element_time_derivative(), PoissonSystem::element_time_derivative(), LaplaceSystem::element_time_derivative(), CurlCurlSystem::element_time_derivative(), ElasticitySystem::element_time_derivative(), CoupledSystem::element_time_derivative(), HilbertSystem::element_time_derivative(), HeatSystem::element_time_derivative(), FirstOrderScalarSystemBase::element_time_derivative(), SecondOrderScalarSystemFirstOrderTimeSolverBase::element_time_derivative(), fe_assembly(), A0::interior_assembly(), B::interior_assembly(), M0::interior_assembly(), A1::interior_assembly(), AssemblyA0::interior_assembly(), AcousticsInnerProduct::interior_assembly(), AssemblyA1::interior_assembly(), A2::interior_assembly(), AssemblyA2::interior_assembly(), F0::interior_assembly(), OutputAssembly::interior_assembly(), EIM_IP_assembly::interior_assembly(), EIM_F::interior_assembly(), AssemblyEIM::interior_assembly(), InnerProductAssembly::interior_assembly(), AssemblyF0::interior_assembly(), AssemblyF1::interior_assembly(), Ex6InnerProduct::interior_assembly(), LaplaceYoung::jacobian(), NavierSystem::mass_residual(), ElasticitySystem::mass_residual(), libMesh::FEMPhysics::mass_residual(), FirstOrderScalarSystemBase::mass_residual(), SecondOrderScalarSystemSecondOrderTimeSolverBase::mass_residual(), SecondOrderScalarSystemFirstOrderTimeSolverBase::mass_residual(), libMesh::FEAbstract::n_quadrature_points(), libMesh::QBase::print_info(), LaplaceYoung::residual(), LaplaceSystem::side_constraint(), LaplaceSystem::side_postprocess(), CoupledSystemQoI::side_qoi(), CoupledSystemQoI::side_qoi_derivative(), LaplaceSystem::side_qoi_derivative(), SolidSystem::side_time_derivative(), CurlCurlSystem::side_time_derivative(), ElasticitySystem::side_time_derivative(), libMesh::QBase::size(), libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), libMesh::QBase::tensor_product_quad(), and QuadratureTest::testPolynomial().

132  {
133  libmesh_assert (!_points.empty());
134  return cast_int<unsigned int>(_points.size());
135  }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
libmesh_assert(ctx)

◆ 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::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.

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

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

82 {
84  out_stream << ReferenceCounter::get_info();
85 }
static std::string get_info()
Gets a string containing the reference information.
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called...

◆ print_info() [2/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.

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

Referenced by libMesh::operator<<().

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 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415
libmesh_assert(ctx)
unsigned int n_points() const
Definition: quadrature.h:131
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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:117

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

References libMesh::QBase::_points.

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

180  {
181  libmesh_assert_less (i, _points.size());
182  return _points[i];
183  }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409

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

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

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 }
unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:379
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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:117

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

References libMesh::QBase::n_points().

Referenced by libMesh::Elem::true_centroid().

143  {
144  return n_points();
145  }
unsigned int n_points() const
Definition: quadrature.h:131

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

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

Referenced by init_2D().

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 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

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

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

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

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 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415

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

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::QTrap::init_3D(), libMesh::QSimpson::init_3D(), and libMesh::QGrid::init_3D().

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 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415

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

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

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

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 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415

◆ type()

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

Implements libMesh::QBase.

Definition at line 32 of file quadrature_monomial.C.

References libMesh::QMONOMIAL.

33 {
34  return QMONOMIAL;
35 }

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

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

189  {
190  libmesh_assert_less (i, _weights.size());
191  return _weights[i];
192  }
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415

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

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

Referenced by init_2D().

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 }
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:409
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:415

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

◆ _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::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(), and libMesh::ReferenceCounter::~ReferenceCounter().

◆ _order

Order libMesh::QBase::_order
protectedinherited

◆ _p_level

unsigned int libMesh::QBase::_p_level
protectedinherited

◆ _points

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

◆ _type

ElemType libMesh::QBase::_type
protectedinherited

◆ _weights

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

◆ 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::QSimpson::init_3D(), libMesh::QTrap::init_3D(), and libMesh::QNodal::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(), init_3D(), and libMesh::QGrundmann_Moller::init_3D().


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