libMesh::QGauss Class Reference

This class implements specific orders of Gauss quadrature. More...

#include <quadrature_gauss.h>

Inheritance diagram for libMesh::QGauss:

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

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

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

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

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

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

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

Private Member Functions
virtual void	init_1D () override
	Initializes the 1D quadrature rule by filling the points and weights vectors with the appropriate values.
virtual void	init_2D () override
	Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate values.
virtual void	init_3D () override
	Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate values.
void	dunavant_rule (const Real rule_data[][4], const unsigned int n_pts)
	The Dunavant rules are for triangles.
void	dunavant_rule2 (const Real *wts, const Real *a, const Real *b, const unsigned int *permutation_ids, const unsigned int n_wts)
void	keast_rule (const Real rule_data[][4], const unsigned int n_pts)
	The Keast rules are for tets.

Detailed Description

This class implements specific orders of Gauss quadrature.

The rules of Order p are capable of integrating polynomials of degree p exactly.

Author

Benjamin Kirk

John W. Peterson

Date

2003

Implements 1, 2, and 3D "Gaussian" quadrature rules.

Definition at line 39 of file quadrature_gauss.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

QGauss() [1/3]

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

inline

Constructor.

Declares the order of the quadrature rule. We explicitly call the init function in 1D since the other tensor-product rules require this one.

Note

The element type, EDGE2, will not be used internally, however if we called the function with INVALID_ELEM it would try to be smart and return, thinking it had already done the work.

Definition at line 52 of file quadrature_gauss.h.

                                     :
    QBase(dim, order)
  {
    if (dim == 1)
      init(EDGE2);
  }

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

QGauss() [2/3]

libMesh::QGauss::QGauss ( const QGauss &  )

default

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

QGauss() [3/3]

libMesh::QGauss::QGauss ( QGauss &&  )

default

~QGauss()

virtual libMesh::QGauss::~QGauss ( )

virtualdefault

Member Function Documentation

build() [1/2]

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

staticinherited

Builds a specific quadrature rule based on the QuadratureType.

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

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

Definition at line 54 of file quadrature_build.C.

{
  switch (_qt)
    {
 
    case QCLOUGH:
      {
#ifdef DEBUG
        if (_order > TWENTYTHIRD)
          {
            libMesh::out << "WARNING: Clough quadrature implemented" << std::endl
                         << " up to TWENTYTHIRD order." << std::endl;
          }
#endif
 
        return std::make_unique<QClough>(_dim, _order);
      }
 
    case QGAUSS:
      {
 
#ifdef DEBUG
        if (_order > FORTYTHIRD)
          {
            libMesh::out << "WARNING: Gauss quadrature implemented" << std::endl
                         << " up to FORTYTHIRD order." << std::endl;
          }
#endif
 
        return std::make_unique<QGauss>(_dim, _order);
      }
 
    case QJACOBI_1_0:
      {
 
#ifdef DEBUG
        if (_order > FORTYTHIRD)
          {
            libMesh::out << "WARNING: Jacobi(1,0) quadrature implemented" << std::endl
                         << " up to FORTYTHIRD order." << std::endl;
          }
 
        if (_dim > 1)
          {
            libMesh::out << "WARNING: Jacobi(1,0) quadrature implemented" << std::endl
                         << " in 1D only." << std::endl;
          }
#endif
 
        return std::make_unique<QJacobi>(_dim, _order, 1, 0);
      }
 
    case QJACOBI_2_0:
      {
 
#ifdef DEBUG
        if (_order > FORTYTHIRD)
          {
            libMesh::out << "WARNING: Jacobi(2,0) quadrature implemented" << std::endl
                         << " up to FORTYTHIRD order." << std::endl;
          }
 
        if (_dim > 1)
          {
            libMesh::out << "WARNING: Jacobi(2,0) quadrature implemented" << std::endl
                         << " in 1D only." << std::endl;
          }
#endif
 
        return std::make_unique<QJacobi>(_dim, _order, 2, 0);
      }
 
    case QSIMPSON:
      {
 
#ifdef DEBUG
        if (_order > THIRD)
          {
            libMesh::out << "WARNING: Simpson rule provides only" << std::endl
                         << " THIRD order!" << std::endl;
          }
#endif
 
        return std::make_unique<QSimpson>(_dim);
      }
 
    case QTRAP:
      {
 
#ifdef DEBUG
        if (_order > FIRST)
          {
            libMesh::out << "WARNING: Trapezoidal rule provides only" << std::endl
                         << " FIRST order!" << std::endl;
          }
#endif
 
        return std::make_unique<QTrap>(_dim);
      }
 
    case QGRID:
      return std::make_unique<QGrid>(_dim, _order);
 
    case QGRUNDMANN_MOLLER:
      return std::make_unique<QGrundmann_Moller>(_dim, _order);
 
    case QMONOMIAL:
      return std::make_unique<QMonomial>(_dim, _order);
 
    case QGAUSS_LOBATTO:
      return std::make_unique<QGaussLobatto>(_dim, _order);
 
    case QCONICAL:
      return std::make_unique<QConical>(_dim, _order);
 
    case QNODAL:
      return std::make_unique<QNodal>(_dim, _order);
 
    default:
      libmesh_error_msg("ERROR: Bad qt=" << _qt);
    }
}

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

build() [2/2]

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

staticinherited

Builds a specific quadrature rule based on the name string.

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

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

Definition at line 43 of file quadrature_build.C.

{
  return QBase::build (Utility::string_to_enum<QuadratureType> (type),
                       _dim,
                       _order);
}

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

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

clone()

std::unique_ptr< QBase > libMesh::QGauss::clone ( ) const

overridevirtual

Returns

A copy of this quadrature rule wrapped in a smart pointer.

Reimplemented from libMesh::QBase.

Definition at line 38 of file quadrature_gauss.C.

{
  return std::make_unique<QGauss>(*this);
}

disable_print_counter_info()

void libMesh::ReferenceCounter::disable_print_counter_info ( )

staticinherited

Definition at line 100 of file reference_counter.C.

{
  _enable_print_counter = false;
  return;
}

References libMesh::ReferenceCounter::_enable_print_counter.

dunavant_rule()

void libMesh::QGauss::dunavant_rule ( const Real  rule_data[][4], const unsigned int  n_pts  )

private

The Dunavant rules are for triangles.

This function takes permutation points and weights in a specific format as input and fills the _points and _weights vectors.

Definition at line 270 of file quadrature_gauss.C.

{
  // The input data array has 4 columns.  The first 3 are the permutation points.
  // The last column is the weights for a given set of permutation points.  A zero
  // in two of the first 3 columns implies the point is a 1-permutation (centroid).
  // A zero in one of the first 3 columns implies the point is a 3-permutation.
  // Otherwise each point is assumed to be a 6-permutation.
 
  // Always insert into the points & weights vector relative to the offset
  unsigned int offset=0;
 
 
  for (unsigned int p=0; p<n_pts; ++p)
    {
 
      // There must always be a non-zero entry to start the row
      libmesh_assert_not_equal_to ( rule_data[p][0], static_cast<Real>(0.0) );
 
      // A zero weight may imply you did not set up the raw data correctly
      libmesh_assert_not_equal_to ( rule_data[p][3], static_cast<Real>(0.0) );
 
      // What kind of point is this?
      // One non-zero entry in first 3 cols   ? 1-perm (centroid) point = 1
      // Two non-zero entries in first 3 cols ? 3-perm point            = 3
      // Three non-zero entries               ? 6-perm point            = 6
      unsigned int pointtype=1;
 
      if (rule_data[p][1] != static_cast<Real>(0.0))
        {
          if (rule_data[p][2] != static_cast<Real>(0.0))
            pointtype = 6;
          else
            pointtype = 3;
        }
 
      switch (pointtype)
        {
        case 1:
          {
            // Be sure we have enough space to insert this point
            libmesh_assert_less (offset + 0, _points.size());
 
            // The point has only a single permutation (the centroid!)
            _points[offset  + 0] = Point(rule_data[p][0], rule_data[p][0]);
 
            // The weight is always the last entry in the row.
            _weights[offset + 0] = rule_data[p][3];
 
            offset += 1;
            break;
          }
 
        case 3:
          {
            // Be sure we have enough space to insert these points
            libmesh_assert_less (offset + 2, _points.size());
 
            // Here it's understood the second entry is to be used twice, and
            // thus there are three possible permutations.
            _points[offset + 0] = Point(rule_data[p][0], rule_data[p][1]);
            _points[offset + 1] = Point(rule_data[p][1], rule_data[p][0]);
            _points[offset + 2] = Point(rule_data[p][1], rule_data[p][1]);
 
            for (unsigned int j=0; j<3; ++j)
              _weights[offset + j] = rule_data[p][3];
 
            offset += 3;
            break;
          }
 
        case 6:
          {
            // Be sure we have enough space to insert these points
            libmesh_assert_less (offset + 5, _points.size());
 
            // Three individual entries with six permutations.
            _points[offset + 0] = Point(rule_data[p][0], rule_data[p][1]);
            _points[offset + 1] = Point(rule_data[p][0], rule_data[p][2]);
            _points[offset + 2] = Point(rule_data[p][1], rule_data[p][0]);
            _points[offset + 3] = Point(rule_data[p][1], rule_data[p][2]);
            _points[offset + 4] = Point(rule_data[p][2], rule_data[p][0]);
            _points[offset + 5] = Point(rule_data[p][2], rule_data[p][1]);
 
            for (unsigned int j=0; j<6; ++j)
              _weights[offset + j] = rule_data[p][3];
 
            offset += 6;
            break;
          }
 
        default:
          libmesh_error_msg("Don't know what to do with this many permutation points!");
        }
    }
}

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

Referenced by init_2D().

dunavant_rule2()

void libMesh::QGauss::dunavant_rule2 ( const Real *  wts, const Real *  a, const Real *  b, const unsigned int *  permutation_ids, const unsigned int  n_wts  )

private

Definition at line 195 of file quadrature_gauss.C.

{
  // Figure out how many total points by summing up the entries
  // in the permutation_ids array, and resize the _points and _weights
  // vectors appropriately.
  unsigned int total_pts = 0;
  for (unsigned int p=0; p<n_wts; ++p)
    total_pts += permutation_ids[p];
 
  // Resize point and weight vectors appropriately.
  _points.resize(total_pts);
  _weights.resize(total_pts);
 
  // Always insert into the points & weights vector relative to the offset
  unsigned int offset=0;
 
  for (unsigned int p=0; p<n_wts; ++p)
    {
      switch (permutation_ids[p])
        {
        case 1:
          {
            // The point has only a single permutation (the centroid!)
            // So we don't even need to look in the a or b arrays.
            _points [offset  + 0] = Point(Real(1)/3, Real(1)/3);
            _weights[offset + 0] = wts[p];
 
            offset += 1;
            break;
          }
 
 
        case 3:
          {
            // For this type of rule, don't need to look in the b array.
            _points[offset + 0] = Point(a[p],         a[p]);         // (a,a)
            _points[offset + 1] = Point(a[p],         1-2*a[p]); // (a,1-2a)
            _points[offset + 2] = Point(1-2*a[p], a[p]);         // (1-2a,a)
 
            for (unsigned int j=0; j<3; ++j)
              _weights[offset + j] = wts[p];
 
            offset += 3;
            break;
          }
 
        case 6:
          {
            // This type of point uses all 3 arrays...
            _points[offset + 0] = Point(a[p], b[p]);
            _points[offset + 1] = Point(b[p], a[p]);
            _points[offset + 2] = Point(a[p], 1-a[p]-b[p]);
            _points[offset + 3] = Point(1-a[p]-b[p], a[p]);
            _points[offset + 4] = Point(b[p], 1-a[p]-b[p]);
            _points[offset + 5] = Point(1-a[p]-b[p], b[p]);
 
            for (unsigned int j=0; j<6; ++j)
              _weights[offset + j] = wts[p];
 
            offset += 6;
            break;
          }
 
        default:
          libmesh_error_msg("Unknown permutation id: " << permutation_ids[p] << "!");
        }
    }
 
}

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

Referenced by init_2D().

enable_print_counter_info()

void libMesh::ReferenceCounter::enable_print_counter_info ( )

staticinherited

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

Enabled by default.

Definition at line 94 of file reference_counter.C.

{
  _enable_print_counter = true;
  return;
}

References libMesh::ReferenceCounter::_enable_print_counter.

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

get_base_order()

Order libMesh::QBase::get_base_order ( ) const

inlineinherited

Returns

The "base" order of the quadrature rule, independent of element.

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

Definition at line 258 of file quadrature.h.

{ return static_cast<Order>(_order); }

References libMesh::QBase::_order.

get_dim()

unsigned int libMesh::QBase::get_dim ( ) const

inlineinherited

Returns

The spatial dimension of the quadrature rule.

Definition at line 150 of file quadrature.h.

{ return _dim; }

References libMesh::QBase::_dim.

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

get_elem_type()

ElemType libMesh::QBase::get_elem_type ( ) const

inlineinherited

Returns

The element type we're currently using.

Definition at line 121 of file quadrature.h.

{ return _type; }

References libMesh::QBase::_type.

get_info()

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

staticinherited

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
 
  std::ostringstream oss;
 
  oss << '\n'
      << " ---------------------------------------------------------------------------- \n"
      << "| Reference count information                                                |\n"
      << " ---------------------------------------------------------------------------- \n";
 
  for (const auto & [name, cd] : _counts)
    oss << "| " << name << " reference count information:\n"
        << "|  Creations:    " << cd.first    << '\n'
        << "|  Destructions: " << cd.second << '\n';
 
  oss << " ---------------------------------------------------------------------------- \n";
 
  return oss.str();
 
#else
 
  return "";
 
#endif
}

References libMesh::ReferenceCounter::_counts.

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

get_order()

Order libMesh::QBase::get_order ( ) const

inlineinherited

Returns

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

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

Todo:

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

Definition at line 249 of file quadrature.h.

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

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

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

get_p_level()

unsigned int libMesh::QBase::get_p_level ( ) const

inlineinherited

Returns

The p-refinement level we're currently using.

Definition at line 126 of file quadrature.h.

{ return _p_level; }

References libMesh::QBase::_p_level.

get_points() [1/2]

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

inlineinherited

Returns

A std::vector containing the quadrature point locations in reference element space as a writable reference.

Definition at line 162 of file quadrature.h.

{ return _points; }

References libMesh::QBase::_points.

get_points() [2/2]

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

inlineinherited

Returns

A std::vector containing the quadrature point locations in reference element space.

Definition at line 156 of file quadrature.h.

{ return _points; }

References libMesh::QBase::_points.

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

get_weights() [1/2]

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

inlineinherited

Returns

A writable references to a std::vector containing the quadrature weights.

Definition at line 174 of file quadrature.h.

{ return _weights; }

References libMesh::QBase::_weights.

get_weights() [2/2]

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

inlineinherited

Returns

A constant reference to a std::vector containing the quadrature weights.

Definition at line 168 of file quadrature.h.

{ return _weights; }

References libMesh::QBase::_weights.

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

increment_constructor_count()

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

inlineprotectednoexceptinherited

Increments the construction counter.

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

Definition at line 183 of file reference_counter.h.

{
  libmesh_try
  {
    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
    std::pair<unsigned int, unsigned int> & p = _counts[name];
    p.first++;
  }
  libmesh_catch (...)
  {
    auto stream = libMesh::err.get();
    stream->exceptions(stream->goodbit); // stream must not throw
    libMesh::err << "Encountered unrecoverable error while calling "
                 << "ReferenceCounter::increment_constructor_count() "
                 << "for a(n) " << name << " object." << std::endl;
    std::terminate();
  }
}

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

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

increment_destructor_count()

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

inlineprotectednoexceptinherited

Increments the destruction counter.

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

Definition at line 207 of file reference_counter.h.

{
  libmesh_try
  {
    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
    std::pair<unsigned int, unsigned int> & p = _counts[name];
    p.second++;
  }
  libmesh_catch (...)
  {
    auto stream = libMesh::err.get();
    stream->exceptions(stream->goodbit); // stream must not throw
    libMesh::err << "Encountered unrecoverable error while calling "
                 << "ReferenceCounter::increment_destructor_count() "
                 << "for a(n) " << name << " object." << std::endl;
    std::terminate();
  }
}

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

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

init() [1/4]

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

virtualinherited

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

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

Definition at line 65 of file quadrature.C.

{
  libmesh_assert_equal_to(elem.dim(), _dim);
 
  // Default to the element p_level() value
  if (p == invalid_uint)
    p = elem.p_level();
 
  ElemType t = elem.type();
 
  // check to see if we have already
  // done the work for this quadrature rule
  //
  // If we have something like a Polygon subclass then we're going to
  // need to recompute to be safe; even if we're using the same
  // element, it might have been distorted enough that its subtriangle
  // triangulation has been changed.
  if (t == _type && p == _p_level && !elem.runtime_topology())
    return;
  else
    {
      _elem = &elem;
      _type = t;
      _p_level = p;
    }
 
  switch(_elem->dim())
    {
    case 0:
      this->init_0D();
 
      return;
 
    case 1:
      this->init_1D();
 
      return;
 
    case 2:
      this->init_2D();
 
      return;
 
    case 3:
      this->init_3D();
 
      return;
 
    default:
      libmesh_error_msg("Invalid dimension _dim = " << _dim);
    }
}

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

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

init() [2/4]

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

virtualinherited

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

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

Definition at line 188 of file quadrature.C.

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

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

init() [3/4]

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

virtualinherited

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

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

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

Definition at line 121 of file quadrature.C.

{
  // Some element types require data from a specific element, so can
  // only be used with newer APIs.
  if (t == C0POLYGON || t == C0POLYHEDRON)
    libmesh_error_msg("Code (see stack trace) used an outdated quadrature function overload.\n"
                      "Quadrature rules on a C0Polygon are not defined by its ElemType alone.");
 
  // This API is dangerous to use on general meshes, which may include
  // element types where the desired quadrature depends on the
  // physical element, but we still want to be able to initialize
  // based on only a type for certain simple cases
  if (t != EDGE2 && !simple_type_only)
    libmesh_deprecated();
 
  // check to see if we have already
  // done the work for this quadrature rule
  if (t == _type && p == _p_level)
    return;
  else
    {
      _elem = nullptr;
      _type = t;
      _p_level = p;
    }
 
  switch(_dim)
    {
    case 0:
      this->init_0D();
 
      return;
 
    case 1:
      this->init_1D();
 
      return;
 
    case 2:
      this->init_2D();
 
      return;
 
    case 3:
      this->init_3D();
 
      return;
 
    default:
      libmesh_error_msg("Invalid dimension _dim = " << _dim);
    }
}

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

init() [4/4]

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

virtualinherited

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

Definition at line 178 of file quadrature.C.

{
  if (other_rule._elem)
    this->init(*other_rule._elem, other_rule._p_level);
  else
    this->init(other_rule._type, other_rule._p_level, true);
}

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

init_0D()

void libMesh::QBase::init_0D ( )

protectedvirtualinherited

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

Generally this is just one point with weight 1.

Definition at line 198 of file quadrature.C.

{
  _points.resize(1);
  _weights.resize(1);
  _points[0] = Point(0.);
  _weights[0] = 1.0;
}

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

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

init_1D()

void libMesh::QGauss::init_1D ( )

overrideprivatevirtual

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

The order of the rule will be defined by the implementing class. It is assumed that derived quadrature rules will at least define the init_1D function, therefore it is pure virtual.

Implements libMesh::QBase.

Definition at line 30 of file quadrature_gauss_1D.C.

{
  //----------------------------------------------------------------------
  // 1D quadrature rules
  switch(get_order())
    {
    case CONSTANT:
    case FIRST:
      {
        _points.resize (1);
        _weights.resize(1);
 
        _points[0](0)  = 0.;
 
        _weights[0]    = 2.;
 
        return;
      }
    case SECOND:
    case THIRD:
      {
        _points.resize (2);
        _weights.resize(2);
 
        _points[0](0) = -std::sqrt(Real(3))/3;
        _points[1]    = -_points[0];
 
        _weights[0]   = 1.;
        _weights[1]   = _weights[0];
 
        return;
      }
    case FOURTH:
    case FIFTH:
      {
        _points.resize (3);
        _weights.resize(3);
 
        _points[ 0](0) = -7.7459666924148337703585307995648e-01_R;
        _points[ 1](0) = 0.;
        _points[ 2]    = -_points[0];
 
        _weights[ 0]   = Real(5)/9;
        _weights[ 1]   = Real(8)/9;
        _weights[ 2]   = _weights[0];
 
        return;
      }
    case SIXTH:
    case SEVENTH:
      {
        _points.resize (4);
        _weights.resize(4);
 
        _points[ 0](0) = -8.6113631159405257522394648889281e-01_R;
        _points[ 1](0) = -3.3998104358485626480266575910324e-01_R;
        _points[ 2]    = -_points[1];
        _points[ 3]    = -_points[0];
 
        _weights[ 0]   = 3.4785484513745385737306394922200e-01_R;
        _weights[ 1]   = 6.5214515486254614262693605077800e-01_R;
        _weights[ 2]   = _weights[1];
        _weights[ 3]   = _weights[0];
 
        return;
      }
    case EIGHTH:
    case NINTH:
      {
        _points.resize (5);
        _weights.resize(5);
 
        _points[ 0](0) = -9.0617984593866399279762687829939e-01_R;
        _points[ 1](0) = -5.3846931010568309103631442070021e-01_R;
        _points[ 2](0) = 0.;
        _points[ 3]    = -_points[1];
        _points[ 4]    = -_points[0];
 
        _weights[ 0]   = 2.3692688505618908751426404071992e-01_R;
        _weights[ 1]   = 4.7862867049936646804129151483564e-01_R;
        _weights[ 2]   = 5.6888888888888888888888888888889e-01_R;
        _weights[ 3]   = _weights[1];
        _weights[ 4]   = _weights[0];
 
        return;
      }
    case TENTH:
    case ELEVENTH:
      {
        _points.resize (6);
        _weights.resize(6);
 
        _points[ 0](0) = -9.3246951420315202781230155449399e-01_R;
        _points[ 1](0) = -6.6120938646626451366139959501991e-01_R;
        _points[ 2](0) = -2.3861918608319690863050172168071e-01_R;
        _points[ 3]    = -_points[2];
        _points[ 4]    = -_points[1];
        _points[ 5]    = -_points[0];
 
        _weights[ 0]   = 1.7132449237917034504029614217273e-01_R;
        _weights[ 1]   = 3.6076157304813860756983351383772e-01_R;
        _weights[ 2]   = 4.6791393457269104738987034398955e-01_R;
        _weights[ 3]   = _weights[2];
        _weights[ 4]   = _weights[1];
        _weights[ 5]   = _weights[0];
 
        return;
      }
    case TWELFTH:
    case THIRTEENTH:
      {
        _points.resize (7);
        _weights.resize(7);
 
        _points[ 0](0) = -9.4910791234275852452618968404785e-01_R;
        _points[ 1](0) = -7.4153118559939443986386477328079e-01_R;
        _points[ 2](0) = -4.0584515137739716690660641207696e-01_R;
        _points[ 3](0) = 0.;
        _points[ 4]    = -_points[2];
        _points[ 5]    = -_points[1];
        _points[ 6]    = -_points[0];
 
        _weights[ 0]   = 1.2948496616886969327061143267908e-01_R;
        _weights[ 1]   = 2.7970539148927666790146777142378e-01_R;
        _weights[ 2]   = 3.8183005050511894495036977548898e-01_R;
        _weights[ 3]   = 4.1795918367346938775510204081633e-01_R;
        _weights[ 4]   = _weights[2];
        _weights[ 5]   = _weights[1];
        _weights[ 6]   = _weights[0];
 
        return;
      }
    case FOURTEENTH:
    case FIFTEENTH:
      {
        _points.resize (8);
        _weights.resize(8);
 
        _points[ 0](0) = -9.6028985649753623168356086856947e-01_R;
        _points[ 1](0) = -7.9666647741362673959155393647583e-01_R;
        _points[ 2](0) = -5.2553240991632898581773904918925e-01_R;
        _points[ 3](0) = -1.8343464249564980493947614236018e-01_R;
        _points[ 4]    = -_points[3];
        _points[ 5]    = -_points[2];
        _points[ 6]    = -_points[1];
        _points[ 7]    = -_points[0];
 
        _weights[ 0]   = 1.0122853629037625915253135430996e-01_R;
        _weights[ 1]   = 2.2238103445337447054435599442624e-01_R;
        _weights[ 2]   = 3.1370664587788728733796220198660e-01_R;
        _weights[ 3]   = 3.6268378337836198296515044927720e-01_R;
        _weights[ 4]   = _weights[3];
        _weights[ 5]   = _weights[2];
        _weights[ 6]   = _weights[1];
        _weights[ 7]   = _weights[0];
 
        return;
      }
    case SIXTEENTH:
    case SEVENTEENTH:
      {
        _points.resize (9);
        _weights.resize(9);
 
        _points[ 0](0) = -9.6816023950762608983557620290367e-01_R;
        _points[ 1](0) = -8.3603110732663579429942978806973e-01_R;
        _points[ 2](0) = -6.1337143270059039730870203934147e-01_R;
        _points[ 3](0) = -3.2425342340380892903853801464334e-01_R;
        _points[ 4](0) = 0.;
        _points[ 5]    = -_points[3];
        _points[ 6]    = -_points[2];
        _points[ 7]    = -_points[1];
        _points[ 8]    = -_points[0];
 
        _weights[ 0]   = 8.1274388361574411971892158110524e-02_R;
        _weights[ 1]   = 1.8064816069485740405847203124291e-01_R;
        _weights[ 2]   = 2.6061069640293546231874286941863e-01_R;
        _weights[ 3]   = 3.1234707704000284006863040658444e-01_R;
        _weights[ 4]   = 3.3023935500125976316452506928697e-01_R;
        _weights[ 5]   = _weights[3];
        _weights[ 6]   = _weights[2];
        _weights[ 7]   = _weights[1];
        _weights[ 8]   = _weights[0];
 
        return;
      }
    case EIGHTTEENTH:
    case NINETEENTH:
      {
        _points.resize (10);
        _weights.resize(10);
 
        _points[ 0](0) = -9.7390652851717172007796401208445e-01_R;
        _points[ 1](0) = -8.6506336668898451073209668842349e-01_R;
        _points[ 2](0) = -6.7940956829902440623432736511487e-01_R;
        _points[ 3](0) = -4.3339539412924719079926594316578e-01_R;
        _points[ 4](0) = -1.4887433898163121088482600112972e-01_R;
        _points[ 5]    = -_points[4];
        _points[ 6]    = -_points[3];
        _points[ 7]    = -_points[2];
        _points[ 8]    = -_points[1];
        _points[ 9]    = -_points[0];
 
        _weights[ 0]   = 6.6671344308688137593568809893332e-02_R;
        _weights[ 1]   = 1.4945134915058059314577633965770e-01_R;
        _weights[ 2]   = 2.1908636251598204399553493422816e-01_R;
        _weights[ 3]   = 2.6926671930999635509122692156947e-01_R;
        _weights[ 4]   = 2.9552422471475287017389299465134e-01_R;
        _weights[ 5]   = _weights[4];
        _weights[ 6]   = _weights[3];
        _weights[ 7]   = _weights[2];
        _weights[ 8]   = _weights[1];
        _weights[ 9]   = _weights[0];
 
        return;
      }
 
    case TWENTIETH:
    case TWENTYFIRST:
      {
        _points.resize (11);
        _weights.resize(11);
 
        _points[ 0](0) = -9.7822865814605699280393800112286e-01_R;
        _points[ 1](0) = -8.8706259976809529907515776930393e-01_R;
        _points[ 2](0) = -7.3015200557404932409341625203115e-01_R;
        _points[ 3](0) = -5.1909612920681181592572566945861e-01_R;
        _points[ 4](0) = -2.6954315595234497233153198540086e-01_R;
        _points[ 5](0) = 0.;
        _points[ 6]    = -_points[4];
        _points[ 7]    = -_points[3];
        _points[ 8]    = -_points[2];
        _points[ 9]    = -_points[1];
        _points[10]    = -_points[0];
 
        _weights[ 0]   = 5.5668567116173666482753720442549e-02_R;
        _weights[ 1]   = 1.2558036946490462463469429922394e-01_R;
        _weights[ 2]   = 1.8629021092773425142609764143166e-01_R;
        _weights[ 3]   = 2.3319376459199047991852370484318e-01_R;
        _weights[ 4]   = 2.6280454451024666218068886989051e-01_R;
        _weights[ 5]   = 2.7292508677790063071448352833634e-01_R;
        _weights[ 6]   = _weights[4];
        _weights[ 7]   = _weights[3];
        _weights[ 8]   = _weights[2];
        _weights[ 9]   = _weights[1];
        _weights[10]   = _weights[0];
 
        return;
      }
 
    case TWENTYSECOND:
    case TWENTYTHIRD:
      {
        _points.resize (12);
        _weights.resize(12);
 
        _points[ 0](0) = -9.8156063424671925069054909014928e-01_R;
        _points[ 1](0) = -9.0411725637047485667846586611910e-01_R;
        _points[ 2](0) = -7.6990267419430468703689383321282e-01_R;
        _points[ 3](0) = -5.8731795428661744729670241894053e-01_R;
        _points[ 4](0) = -3.6783149899818019375269153664372e-01_R;
        _points[ 5](0) = -1.2523340851146891547244136946385e-01_R;
        _points[ 6]    = -_points[5];
        _points[ 7]    = -_points[4];
        _points[ 8]    = -_points[3];
        _points[ 9]    = -_points[2];
        _points[10]    = -_points[1];
        _points[11]    = -_points[0];
 
        _weights[ 0]   = 4.7175336386511827194615961485017e-02_R;
        _weights[ 1]   = 1.0693932599531843096025471819400e-01_R;
        _weights[ 2]   = 1.6007832854334622633465252954336e-01_R;
        _weights[ 3]   = 2.0316742672306592174906445580980e-01_R;
        _weights[ 4]   = 2.3349253653835480876084989892483e-01_R;
        _weights[ 5]   = 2.4914704581340278500056243604295e-01_R;
        _weights[ 6]   = _weights[5];
        _weights[ 7]   = _weights[4];
        _weights[ 8]   = _weights[3];
        _weights[ 9]   = _weights[2];
        _weights[10]   = _weights[1];
        _weights[11]   = _weights[0];
 
        return;
      }
 
    case TWENTYFOURTH:
    case TWENTYFIFTH:
      {
        _points.resize (13);
        _weights.resize(13);
 
        _points[ 0](0) = -9.8418305471858814947282944880711e-01_R;
        _points[ 1](0) = -9.1759839922297796520654783650072e-01_R;
        _points[ 2](0) = -8.0157809073330991279420648958286e-01_R;
        _points[ 3](0) = -6.4234933944034022064398460699552e-01_R;
        _points[ 4](0) = -4.4849275103644685287791285212764e-01_R;
        _points[ 5](0) = -2.3045831595513479406552812109799e-01_R;
        _points[ 6](0) = 0.;
        _points[ 7]    = -_points[5];
        _points[ 8]    = -_points[4];
        _points[ 9]    = -_points[3];
        _points[10]    = -_points[2];
        _points[11]    = -_points[1];
        _points[12]    = -_points[0];
 
        _weights[ 0]   = 4.0484004765315879520021592200986e-02_R;
        _weights[ 1]   = 9.2121499837728447914421775953797e-02_R;
        _weights[ 2]   = 1.3887351021978723846360177686887e-01_R;
        _weights[ 3]   = 1.7814598076194573828004669199610e-01_R;
        _weights[ 4]   = 2.0781604753688850231252321930605e-01_R;
        _weights[ 5]   = 2.2628318026289723841209018603978e-01_R;
        _weights[ 6]   = 2.3255155323087391019458951526884e-01_R;
        _weights[ 7]   = _weights[5];
        _weights[ 8]   = _weights[4];
        _weights[ 9]   = _weights[3];
        _weights[10]   = _weights[2];
        _weights[11]   = _weights[1];
        _weights[12]   = _weights[0];
 
        return;
      }
 
    case TWENTYSIXTH:
    case TWENTYSEVENTH:
      {
        _points.resize (14);
        _weights.resize(14);
 
        _points[ 0](0) = -9.8628380869681233884159726670405e-01_R;
        _points[ 1](0) = -9.2843488366357351733639113937787e-01_R;
        _points[ 2](0) = -8.2720131506976499318979474265039e-01_R;
        _points[ 3](0) = -6.8729290481168547014801980301933e-01_R;
        _points[ 4](0) = -5.1524863635815409196529071855119e-01_R;
        _points[ 5](0) = -3.1911236892788976043567182416848e-01_R;
        _points[ 6](0) = -1.0805494870734366206624465021983e-01_R;
        _points[ 7]    = -_points[6];
        _points[ 8]    = -_points[5];
        _points[ 9]    = -_points[4];
        _points[10]    = -_points[3];
        _points[11]    = -_points[2];
        _points[12]    = -_points[1];
        _points[13]    = -_points[0];
 
        _weights[ 0]   = 3.5119460331751863031832876138192e-02_R;
        _weights[ 1]   = 8.0158087159760209805633277062854e-02_R;
        _weights[ 2]   = 1.2151857068790318468941480907248e-01_R;
        _weights[ 3]   = 1.5720316715819353456960193862384e-01_R;
        _weights[ 4]   = 1.8553839747793781374171659012516e-01_R;
        _weights[ 5]   = 2.0519846372129560396592406566122e-01_R;
        _weights[ 6]   = 2.1526385346315779019587644331626e-01_R;
        _weights[ 7]   = _weights[6];
        _weights[ 8]   = _weights[5];
        _weights[ 9]   = _weights[4];
        _weights[10]   = _weights[3];
        _weights[11]   = _weights[2];
        _weights[12]   = _weights[1];
        _weights[13]   = _weights[0];
 
        return;
      }
 
    case TWENTYEIGHTH:
    case TWENTYNINTH:
      {
        _points.resize (15);
        _weights.resize(15);
 
        _points[ 0](0) = -9.8799251802048542848956571858661e-01_R;
        _points[ 1](0) = -9.3727339240070590430775894771021e-01_R;
        _points[ 2](0) = -8.4820658341042721620064832077422e-01_R;
        _points[ 3](0) = -7.2441773136017004741618605461394e-01_R;
        _points[ 4](0) = -5.7097217260853884753722673725391e-01_R;
        _points[ 5](0) = -3.9415134707756336989720737098105e-01_R;
        _points[ 6](0) = -2.0119409399743452230062830339460e-01_R;
        _points[ 7](0) = 0.;
        _points[ 8]    = -_points[6];
        _points[ 9]    = -_points[5];
        _points[10]    = -_points[4];
        _points[11]    = -_points[3];
        _points[12]    = -_points[2];
        _points[13]    = -_points[1];
        _points[14]    = -_points[0];
 
        _weights[ 0]   = 3.0753241996117268354628393577204e-02_R;
        _weights[ 1]   = 7.0366047488108124709267416450667e-02_R;
        _weights[ 2]   = 1.0715922046717193501186954668587e-01_R;
        _weights[ 3]   = 1.3957067792615431444780479451103e-01_R;
        _weights[ 4]   = 1.6626920581699393355320086048121e-01_R;
        _weights[ 5]   = 1.8616100001556221102680056186642e-01_R;
        _weights[ 6]   = 1.9843148532711157645611832644384e-01_R;
        _weights[ 7]   = 2.0257824192556127288062019996752e-01_R;
        _weights[ 8]   = _weights[6];
        _weights[ 9]   = _weights[5];
        _weights[10]   = _weights[4];
        _weights[11]   = _weights[3];
        _weights[12]   = _weights[2];
        _weights[13]   = _weights[1];
        _weights[14]   = _weights[0];
 
        return;
      }
 
    case THIRTIETH:
    case THIRTYFIRST:
      {
        _points.resize (16);
        _weights.resize(16);
 
        _points[ 0](0) = -9.8940093499164993259615417345033e-01_R;
        _points[ 1](0) = -9.4457502307323257607798841553461e-01_R;
        _points[ 2](0) = -8.6563120238783174388046789771239e-01_R;
        _points[ 3](0) = -7.5540440835500303389510119484744e-01_R;
        _points[ 4](0) = -6.1787624440264374844667176404879e-01_R;
        _points[ 5](0) = -4.5801677765722738634241944298358e-01_R;
        _points[ 6](0) = -2.8160355077925891323046050146050e-01_R;
        _points[ 7](0) = -9.5012509837637440185319335424958e-02_R;
        _points[ 8]    = -_points[7];
        _points[ 9]    = -_points[6];
        _points[10]    = -_points[5];
        _points[11]    = -_points[4];
        _points[12]    = -_points[3];
        _points[13]    = -_points[2];
        _points[14]    = -_points[1];
        _points[15]    = -_points[0];
 
        _weights[ 0]   = 2.7152459411754094851780572456018e-02_R;
        _weights[ 1]   = 6.2253523938647892862843836994378e-02_R;
        _weights[ 2]   = 9.5158511682492784809925107602246e-02_R;
        _weights[ 3]   = 1.2462897125553387205247628219202e-01_R;
        _weights[ 4]   = 1.4959598881657673208150173054748e-01_R;
        _weights[ 5]   = 1.6915651939500253818931207903033e-01_R;
        _weights[ 6]   = 1.8260341504492358886676366796922e-01_R;
        _weights[ 7]   = 1.8945061045506849628539672320828e-01_R;
        _weights[ 8]   = _weights[7];
        _weights[ 9]   = _weights[6];
        _weights[10]   = _weights[5];
        _weights[11]   = _weights[4];
        _weights[12]   = _weights[3];
        _weights[13]   = _weights[2];
        _weights[14]   = _weights[1];
        _weights[15]   = _weights[0];
 
        return;
      }
 
    case THIRTYSECOND:
    case THIRTYTHIRD:
      {
        _points.resize (17);
        _weights.resize(17);
 
        _points[ 0](0) = -9.9057547531441733567543401994067e-01_R;
        _points[ 1](0) = -9.5067552176876776122271695789580e-01_R;
        _points[ 2](0) = -8.8023915372698590212295569448816e-01_R;
        _points[ 3](0) = -7.8151400389680140692523005552048e-01_R;
        _points[ 4](0) = -6.5767115921669076585030221664300e-01_R;
        _points[ 5](0) = -5.1269053708647696788624656862955e-01_R;
        _points[ 6](0) = -3.5123176345387631529718551709535e-01_R;
        _points[ 7](0) = -1.7848418149584785585067749365407e-01_R;
        _points[ 8](0) = 0.;
        _points[ 9]    = -_points[7];
        _points[10]    = -_points[6];
        _points[11]    = -_points[5];
        _points[12]    = -_points[4];
        _points[13]    = -_points[3];
        _points[14]    = -_points[2];
        _points[15]    = -_points[1];
        _points[16]    = -_points[0];
 
        _weights[ 0]   = 2.4148302868547931960110026287565e-02_R;
        _weights[ 1]   = 5.5459529373987201129440165358245e-02_R;
        _weights[ 2]   = 8.5036148317179180883535370191062e-02_R;
        _weights[ 3]   = 1.1188384719340397109478838562636e-01_R;
        _weights[ 4]   = 1.3513636846852547328631998170235e-01_R;
        _weights[ 5]   = 1.5404576107681028808143159480196e-01_R;
        _weights[ 6]   = 1.6800410215645004450997066378832e-01_R;
        _weights[ 7]   = 1.7656270536699264632527099011320e-01_R;
        _weights[ 8]   = 1.7944647035620652545826564426189e-01_R;
        _weights[ 9]   = _weights[7];
        _weights[10]   = _weights[6];
        _weights[11]   = _weights[5];
        _weights[12]   = _weights[4];
        _weights[13]   = _weights[3];
        _weights[14]   = _weights[2];
        _weights[15]   = _weights[1];
        _weights[16]   = _weights[0];
 
        return;
      }
 
    case THIRTYFOURTH:
    case THIRTYFIFTH:
      {
        _points.resize (18);
        _weights.resize(18);
 
        _points[ 0](0) = -9.9156516842093094673001600470615e-01_R;
        _points[ 1](0) = -9.5582394957139775518119589292978e-01_R;
        _points[ 2](0) = -8.9260246649755573920606059112715e-01_R;
        _points[ 3](0) = -8.0370495897252311568241745501459e-01_R;
        _points[ 4](0) = -6.9168704306035320787489108128885e-01_R;
        _points[ 5](0) = -5.5977083107394753460787154852533e-01_R;
        _points[ 6](0) = -4.1175116146284264603593179383305e-01_R;
        _points[ 7](0) = -2.5188622569150550958897285487791e-01_R;
        _points[ 8](0) = -8.4775013041735301242261852935784e-02_R;
        _points[ 9]    = -_points[8];
        _points[10]    = -_points[7];
        _points[11]    = -_points[6];
        _points[12]    = -_points[5];
        _points[13]    = -_points[4];
        _points[14]    = -_points[3];
        _points[15]    = -_points[2];
        _points[16]    = -_points[1];
        _points[17]    = -_points[0];
 
        _weights[ 0]   = 2.1616013526483310313342710266452e-02_R;
        _weights[ 1]   = 4.9714548894969796453334946202639e-02_R;
        _weights[ 2]   = 7.6425730254889056529129677616637e-02_R;
        _weights[ 3]   = 1.0094204410628716556281398492483e-01_R;
        _weights[ 4]   = 1.2255520671147846018451912680020e-01_R;
        _weights[ 5]   = 1.4064291467065065120473130375195e-01_R;
        _weights[ 6]   = 1.5468467512626524492541800383637e-01_R;
        _weights[ 7]   = 1.6427648374583272298605377646593e-01_R;
        _weights[ 8]   = 1.6914238296314359184065647013499e-01_R;
        _weights[ 9]   = _weights[8];
        _weights[10]   = _weights[7];
        _weights[11]   = _weights[6];
        _weights[12]   = _weights[5];
        _weights[13]   = _weights[4];
        _weights[14]   = _weights[3];
        _weights[15]   = _weights[2];
        _weights[16]   = _weights[1];
        _weights[17]   = _weights[0];
 
        return;
      }
 
    case THIRTYSIXTH:
    case THIRTYSEVENTH:
      {
        _points.resize (19);
        _weights.resize(19);
 
        _points[ 0](0) = -9.9240684384358440318901767025326e-01_R;
        _points[ 1](0) = -9.6020815213483003085277884068765e-01_R;
        _points[ 2](0) = -9.0315590361481790164266092853231e-01_R;
        _points[ 3](0) = -8.2271465653714282497892248671271e-01_R;
        _points[ 4](0) = -7.2096617733522937861709586082378e-01_R;
        _points[ 5](0) = -6.0054530466168102346963816494624e-01_R;
        _points[ 6](0) = -4.6457074137596094571726714810410e-01_R;
        _points[ 7](0) = -3.1656409996362983199011732884984e-01_R;
        _points[ 8](0) = -1.6035864564022537586809611574074e-01_R;
        _points[ 9](0) = 0.;
        _points[10]    = -_points[8];
        _points[11]    = -_points[7];
        _points[12]    = -_points[6];
        _points[13]    = -_points[5];
        _points[14]    = -_points[4];
        _points[15]    = -_points[3];
        _points[16]    = -_points[2];
        _points[17]    = -_points[1];
        _points[18]    = -_points[0];
 
        _weights[ 0]   = 1.9461788229726477036312041464438e-02_R;
        _weights[ 1]   = 4.4814226765699600332838157401994e-02_R;
        _weights[ 2]   = 6.9044542737641226580708258006013e-02_R;
        _weights[ 3]   = 9.1490021622449999464462094123840e-02_R;
        _weights[ 4]   = 1.1156664554733399471602390168177e-01_R;
        _weights[ 5]   = 1.2875396253933622767551578485688e-01_R;
        _weights[ 6]   = 1.4260670217360661177574610944190e-01_R;
        _weights[ 7]   = 1.5276604206585966677885540089766e-01_R;
        _weights[ 8]   = 1.5896884339395434764995643946505e-01_R;
        _weights[ 9]   = 1.6105444984878369597916362532092e-01_R;
        _weights[10]   = _weights[8];
        _weights[11]   = _weights[7];
        _weights[12]   = _weights[6];
        _weights[13]   = _weights[5];
        _weights[14]   = _weights[4];
        _weights[15]   = _weights[3];
        _weights[16]   = _weights[2];
        _weights[17]   = _weights[1];
        _weights[18]   = _weights[0];
 
        return;
      }
 
    case THIRTYEIGHTH:
    case THIRTYNINTH:
      {
        _points.resize (20);
        _weights.resize(20);
 
        _points[ 0](0) = -9.9312859918509492478612238847132e-01_R;
        _points[ 1](0) = -9.6397192727791379126766613119728e-01_R;
        _points[ 2](0) = -9.1223442825132590586775244120330e-01_R;
        _points[ 3](0) = -8.3911697182221882339452906170152e-01_R;
        _points[ 4](0) = -7.4633190646015079261430507035564e-01_R;
        _points[ 5](0) = -6.3605368072651502545283669622629e-01_R;
        _points[ 6](0) = -5.1086700195082709800436405095525e-01_R;
        _points[ 7](0) = -3.7370608871541956067254817702493e-01_R;
        _points[ 8](0) = -2.2778585114164507808049619536857e-01_R;
        _points[ 9](0) = -7.6526521133497333754640409398838e-02_R;
        _points[10]    = -_points[9];
        _points[11]    = -_points[8];
        _points[12]    = -_points[7];
        _points[13]    = -_points[6];
        _points[14]    = -_points[5];
        _points[15]    = -_points[4];
        _points[16]    = -_points[3];
        _points[17]    = -_points[2];
        _points[18]    = -_points[1];
        _points[19]    = -_points[0];
 
        _weights[ 0]   = 1.7614007139152118311861962351853e-02_R;
        _weights[ 1]   = 4.0601429800386941331039952274932e-02_R;
        _weights[ 2]   = 6.2672048334109063569506535187042e-02_R;
        _weights[ 3]   = 8.3276741576704748724758143222046e-02_R;
        _weights[ 4]   = 1.0193011981724043503675013548035e-01_R;
        _weights[ 5]   = 1.1819453196151841731237737771138e-01_R;
        _weights[ 6]   = 1.3168863844917662689849449974816e-01_R;
        _weights[ 7]   = 1.4209610931838205132929832506716e-01_R;
        _weights[ 8]   = 1.4917298647260374678782873700197e-01_R;
        _weights[ 9]   = 1.5275338713072585069808433195510e-01_R;
        _weights[10]   = _weights[9];
        _weights[11]   = _weights[8];
        _weights[12]   = _weights[7];
        _weights[13]   = _weights[6];
        _weights[14]   = _weights[5];
        _weights[15]   = _weights[4];
        _weights[16]   = _weights[3];
        _weights[17]   = _weights[2];
        _weights[18]   = _weights[1];
        _weights[19]   = _weights[0];
 
        return;
      }
 
    case FORTIETH:
    case FORTYFIRST:
      {
        _points.resize (21);
        _weights.resize(21);
 
        _points[ 0](0) = -9.9375217062038950026024203593794e-01_R;
        _points[ 1](0) = -9.6722683856630629431662221490770e-01_R;
        _points[ 2](0) = -9.2009933415040082879018713371497e-01_R;
        _points[ 3](0) = -8.5336336458331728364725063858757e-01_R;
        _points[ 4](0) = -7.6843996347567790861587785130623e-01_R;
        _points[ 5](0) = -6.6713880419741231930596666999034e-01_R;
        _points[ 6](0) = -5.5161883588721980705901879672431e-01_R;
        _points[ 7](0) = -4.2434212020743878357366888854379e-01_R;
        _points[ 8](0) = -2.8802131680240109660079251606460e-01_R;
        _points[ 9](0) = -1.4556185416089509093703098233869e-01_R;
        _points[10](0) = 0.;
        _points[11]    = -_points[9];
        _points[12]    = -_points[8];
        _points[13]    = -_points[7];
        _points[14]    = -_points[6];
        _points[15]    = -_points[5];
        _points[16]    = -_points[4];
        _points[17]    = -_points[3];
        _points[18]    = -_points[2];
        _points[19]    = -_points[1];
        _points[20]    = -_points[0];
 
        _weights[ 0]   = 1.6017228257774333324224616858471e-02_R;
        _weights[ 1]   = 3.6953789770852493799950668299330e-02_R;
        _weights[ 2]   = 5.7134425426857208283635826472448e-02_R;
        _weights[ 3]   = 7.6100113628379302017051653300183e-02_R;
        _weights[ 4]   = 9.3444423456033861553289741113932e-02_R;
        _weights[ 5]   = 1.0879729916714837766347457807011e-01_R;
        _weights[ 6]   = 1.2183141605372853419536717712572e-01_R;
        _weights[ 7]   = 1.3226893863333746178105257449678e-01_R;
        _weights[ 8]   = 1.3988739479107315472213342386758e-01_R;
        _weights[ 9]   = 1.4452440398997005906382716655375e-01_R;
        _weights[10]   = 1.4608113364969042719198514768337e-01_R;
        _weights[11]   = _weights[9];
        _weights[12]   = _weights[8];
        _weights[13]   = _weights[7];
        _weights[14]   = _weights[6];
        _weights[15]   = _weights[5];
        _weights[16]   = _weights[4];
        _weights[17]   = _weights[3];
        _weights[18]   = _weights[2];
        _weights[19]   = _weights[1];
        _weights[20]   = _weights[0];
 
        return;
      }
 
    case FORTYSECOND:
    case FORTYTHIRD:
      {
        _points.resize (22);
        _weights.resize(22);
 
        _points[ 0](0) = -9.9429458548239929207303142116130e-01_R;
        _points[ 1](0) = -9.7006049783542872712395098676527e-01_R;
        _points[ 2](0) = -9.2695677218717400052069293925905e-01_R;
        _points[ 3](0) = -8.6581257772030013653642563701938e-01_R;
        _points[ 4](0) = -7.8781680597920816200427795540835e-01_R;
        _points[ 5](0) = -6.9448726318668278005068983576226e-01_R;
        _points[ 6](0) = -5.8764040350691159295887692763865e-01_R;
        _points[ 7](0) = -4.6935583798675702640633071096641e-01_R;
        _points[ 8](0) = -3.4193582089208422515814742042738e-01_R;
        _points[ 9](0) = -2.0786042668822128547884653391955e-01_R;
        _points[10](0) = -6.9739273319722221213841796118628e-02_R;
        _points[11]    = -_points[10];
        _points[12]    = -_points[9];
        _points[13]    = -_points[8];
        _points[14]    = -_points[7];
        _points[15]    = -_points[6];
        _points[16]    = -_points[5];
        _points[17]    = -_points[4];
        _points[18]    = -_points[3];
        _points[19]    = -_points[2];
        _points[20]    = -_points[1];
        _points[21]    = -_points[0];
 
        _weights[ 0]   = 1.4627995298272200684991098047185e-02_R;
        _weights[ 1]   = 3.3774901584814154793302246865913e-02_R;
        _weights[ 2]   = 5.2293335152683285940312051273211e-02_R;
        _weights[ 3]   = 6.9796468424520488094961418930218e-02_R;
        _weights[ 4]   = 8.5941606217067727414443681372703e-02_R;
        _weights[ 5]   = 1.0041414444288096493207883783054e-01_R;
        _weights[ 6]   = 1.1293229608053921839340060742175e-01_R;
        _weights[ 7]   = 1.2325237681051242428556098615481e-01_R;
        _weights[ 8]   = 1.3117350478706237073296499253031e-01_R;
        _weights[ 9]   = 1.3654149834601517135257383123152e-01_R;
        _weights[10]   = 1.3925187285563199337541024834181e-01_R;
        _weights[11]   = _weights[10];
        _weights[12]   = _weights[9];
        _weights[13]   = _weights[8];
        _weights[14]   = _weights[7];
        _weights[15]   = _weights[6];
        _weights[16]   = _weights[5];
        _weights[17]   = _weights[4];
        _weights[18]   = _weights[3];
        _weights[19]   = _weights[2];
        _weights[20]   = _weights[1];
        _weights[21]   = _weights[0];
 
        return;
      }
 
 
    default:
      libmesh_error_msg("Quadrature rule " << _order << " not supported!");
    }
}

References libMesh::QBase::_order, libMesh::QBase::_points, libMesh::QBase::_weights, libMesh::CONSTANT, libMesh::EIGHTH, libMesh::EIGHTTEENTH, libMesh::ELEVENTH, libMesh::FIFTEENTH, libMesh::FIFTH, libMesh::FIRST, libMesh::FORTIETH, libMesh::FORTYFIRST, libMesh::FORTYSECOND, libMesh::FORTYTHIRD, libMesh::FOURTEENTH, libMesh::FOURTH, libMesh::QBase::get_order(), libMesh::NINETEENTH, libMesh::NINTH, libMesh::Real, libMesh::SECOND, libMesh::SEVENTEENTH, libMesh::SEVENTH, libMesh::SIXTEENTH, libMesh::SIXTH, libMesh::TENTH, libMesh::THIRD, libMesh::THIRTEENTH, libMesh::THIRTIETH, libMesh::THIRTYEIGHTH, libMesh::THIRTYFIFTH, libMesh::THIRTYFIRST, libMesh::THIRTYFOURTH, libMesh::THIRTYNINTH, libMesh::THIRTYSECOND, libMesh::THIRTYSEVENTH, libMesh::THIRTYSIXTH, libMesh::THIRTYTHIRD, libMesh::TWELFTH, libMesh::TWENTIETH, libMesh::TWENTYEIGHTH, libMesh::TWENTYFIFTH, libMesh::TWENTYFIRST, libMesh::TWENTYFOURTH, libMesh::TWENTYNINTH, libMesh::TWENTYSECOND, libMesh::TWENTYSEVENTH, libMesh::TWENTYSIXTH, and libMesh::TWENTYTHIRD.

init_2D()

void libMesh::QGauss::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 31 of file quadrature_gauss_2D.C.

{
#if LIBMESH_DIM > 1
 
  //-----------------------------------------------------------------------
  // 2D quadrature rules
  switch (_type)
    {
 
 
      //---------------------------------------------
      // Quadrilateral quadrature rules
    case QUAD4:
    case QUADSHELL4:
    case QUAD8:
    case QUADSHELL8:
    case QUAD9:
    case QUADSHELL9:
      {
        // We compute the 2D quadrature rule as a tensor
        // product of the 1D quadrature rule.
        //
        // For QUADs, a quadrature rule of order 'p' must be able to integrate
        // bilinear (p=1), biquadratic (p=2), bicubic (p=3), etc. polynomials of the form
        //
        // (x^p + x^{p-1} + ... + 1) * (y^p + y^{p-1} + ... + 1)
        //
        // These polynomials have terms *up to* degree 2p but they are *not* complete
        // polynomials of degree 2p. For example, when p=2 we have
        //        1
        //     x      y
        // x^2    xy     y^2
        //    yx^2   xy^2
        //       x^2y^2
        QGauss q1D(1,get_order());
        tensor_product_quad( q1D );
        return;
      }
 
 
      //---------------------------------------------
      // Triangle quadrature rules
    case TRI3:
    case TRISHELL3:
    case TRI3SUBDIVISION:
    case TRI6:
    case TRI7:
      {
        switch(get_order())
          {
          case CONSTANT:
          case FIRST:
            {
              // Exact for linears
              _points.resize(1);
              _weights.resize(1);
 
              _points[0](0) = Real(1)/3;
              _points[0](1) = Real(1)/3;
 
              _weights[0] = 0.5;
 
              return;
            }
          case SECOND:
            {
              // Exact for quadratics
              _points.resize(3);
              _weights.resize(3);
 
              // Alternate rule with points on ref. elt. boundaries.
              // Not ideal for problems with material coefficient discontinuities
              // aligned along element boundaries.
              // _points[0](0) = .5;
              // _points[0](1) = .5;
              // _points[1](0) = 0.;
              // _points[1](1) = .5;
              // _points[2](0) = .5;
              // _points[2](1) = .0;
 
              _points[0](0) = Real(2)/3;
              _points[0](1) = Real(1)/6;
 
              _points[1](0) = Real(1)/6;
              _points[1](1) = Real(2)/3;
 
              _points[2](0) = Real(1)/6;
              _points[2](1) = Real(1)/6;
 
 
              _weights[0] = Real(1)/6;
              _weights[1] = Real(1)/6;
              _weights[2] = Real(1)/6;
 
              return;
            }
          case THIRD:
            {
              // Exact for cubics
              _points.resize(4);
              _weights.resize(4);
 
              // This rule is formed from a tensor product of
              // appropriately-scaled Gauss and Jacobi rules.  (See
              // also the QConical quadrature class, this is a
              // hard-coded version of one of those rules.)  For high
              // orders these rules generally have too many points,
              // but at extremely low order they are competitive and
              // have the additional benefit of having all positive
              // weights.
              _points[0](0) = 1.5505102572168219018027159252941e-01_R;
              _points[0](1) = 1.7855872826361642311703513337422e-01_R;
              _points[1](0) = 6.4494897427831780981972840747059e-01_R;
              _points[1](1) = 7.5031110222608118177475598324603e-02_R;
              _points[2](0) = 1.5505102572168219018027159252941e-01_R;
              _points[2](1) = 6.6639024601470138670269327409637e-01_R;
              _points[3](0) = 6.4494897427831780981972840747059e-01_R;
              _points[3](1) = 2.8001991549907407200279599420481e-01_R;
 
              _weights[0] = 1.5902069087198858469718450103758e-01_R;
              _weights[1] = 9.0979309128011415302815498962418e-02_R;
              _weights[2] = 1.5902069087198858469718450103758e-01_R;
              _weights[3] = 9.0979309128011415302815498962418e-02_R;
 
              return;
 
 
              // The following third-order rule is quite commonly cited
              // in the literature and most likely works fine.  However,
              // we generally prefer a rule with all positive weights
              // and an equal number of points, when available.
              //
              //  (allow_rules_with_negative_weights)
              // {
              //   // Exact for cubics
              //   _points.resize(4);
              //   _weights.resize(4);
              //
              //   _points[0](0) = .33333333333333333333333333333333;
              //   _points[0](1) = .33333333333333333333333333333333;
              //
              //   _points[1](0) = .2;
              //   _points[1](1) = .6;
              //
              //   _points[2](0) = .2;
              //   _points[2](1) = .2;
              //
              //   _points[3](0) = .6;
              //   _points[3](1) = .2;
              //
              //
              //   _weights[0] = -27./96.;
              //   _weights[1] =  25./96.;
              //   _weights[2] =  25./96.;
              //   _weights[3] =  25./96.;
              //
              //   return;
              // } // end if (allow_rules_with_negative_weights)
              // Note: if !allow_rules_with_negative_weights, fall through to next case.
            }
 
 
 
            // A degree 4 rule with six points.  This rule can be found in many places
            // including:
            //
            // J.N. Lyness and D. Jespersen, Moderate degree symmetric
            // quadrature rules for the triangle, J. Inst. Math. Appl.  15 (1975),
            // 19--32.
            //
            // We used the code in:
            // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
            // on triangles and tetrahedra"  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
            // to generate additional precision.
          case FOURTH:
            {
              const unsigned int n_wts = 2;
              const Real wts[n_wts] =
                {
                  1.1169079483900573284750350421656140e-01_R,
                  5.4975871827660933819163162450105264e-02_R
                };
 
              const Real a[n_wts] =
                {
                  4.4594849091596488631832925388305199e-01_R,
                  9.1576213509770743459571463402201508e-02_R
                };
 
              const Real b[n_wts] = {0., 0.}; // not used
              const unsigned int permutation_ids[n_wts] = {3, 3};
 
              dunavant_rule2(wts, a, b, permutation_ids, n_wts); // 6 total points
 
              return;
            }
 
 
 
            // Exact for quintics
            // Can be found in "Quadrature on Simplices of Arbitrary
            // Dimension" by Walkington.
          case FIFTH:
            {
              const unsigned int n_wts = 3;
              const Real wts[n_wts] =
                {
                  Real(9)/80,                              // ~0.1125
                  Real(31)/480 + std::sqrt(Real(15))/2400, // ~0.06619707639
                  Real(31)/480 - std::sqrt(Real(15))/2400  // ~0.06296959027
                };
 
              const Real a[n_wts] =
                {
                  0., // 'a' parameter not used for origin permutation
                  Real(2)/7 + std::sqrt(Real(15))/21, // ~0.4701420641
                  Real(2)/7 - std::sqrt(Real(15))/21  // ~0.1012865073
                };
 
              const Real b[n_wts] = {0., 0., 0.}; // not used
              const unsigned int permutation_ids[n_wts] = {1, 3, 3};
 
              dunavant_rule2(wts, a, b, permutation_ids, n_wts); // 7 total points
 
              return;
            }
 
 
 
            // A degree 6 rule with 12 points.  This rule can be found in many places
            // including:
            //
            // J.N. Lyness and D. Jespersen, Moderate degree symmetric
            // quadrature rules for the triangle, J. Inst. Math. Appl.  15 (1975),
            // 19--32.
            //
            // We used the code in:
            // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
            // on triangles and tetrahedra"  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
            // to generate additional precision.
            //
            // Note that the following 7th-order Ro3-invariant rule also has only 12 points,
            // which technically makes it the superior rule.  This one is here for completeness.
          case SIXTH:
            {
              const unsigned int n_wts = 3;
              const Real wts[n_wts] =
                {
                  5.8393137863189683012644805692789721e-02_R,
                  2.5422453185103408460468404553434492e-02_R,
                  4.1425537809186787596776728210221227e-02_R
                };
 
              const Real a[n_wts] =
                {
                  2.4928674517091042129163855310701908e-01_R,
                  6.3089014491502228340331602870819157e-02_R,
                  3.1035245103378440541660773395655215e-01_R
                };
 
              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  6.3650249912139864723014259441204970e-01_R
                };
 
              const unsigned int permutation_ids[n_wts] = {3, 3, 6}; // 12 total points
 
              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
 
              return;
            }
 
 
            // A degree 7 rule with 12 points.  This rule can be found in:
            //
            // K. Gatermann, The construction of symmetric cubature
            // formulas for the square and the triangle, Computing 40
            // (1988), 229--240.
            //
            // This rule, which is provably minimal in the number of
            // integration points, is said to be 'Ro3 invariant' which
            // means that a given set of barycentric coordinates
            // (z1,z2,z3) implies the quadrature points (z1,z2),
            // (z3,z1), (z2,z3) which are formed by taking the first
            // two entries in cyclic permutations of the barycentric
            // point.  Barycentric coordinates are related in the
            // sense that: z3 = 1 - z1 - z2.
            //
            // The 12-point sixth-order rule for triangles given in
            // Flaherty's (http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps)
            // lecture notes has been removed in favor of this rule
            // which is higher-order (for the same number of
            // quadrature points) and has a few more digits of
            // precision in the points and weights.  Some 10-point
            // degree 6 rules exist for the triangle but they have
            // quadrature points outside the region of integration.
          case SEVENTH:
            {
              _points.resize (12);
              _weights.resize(12);
 
              const unsigned int nrows=4;
 
              // In each of the rows below, the first two entries are (z1, z2) which imply
              // z3.  The third entry is the weight for each of the points in the cyclic permutation.
              // The original publication tabulated about 16 decimal digits for each point and weight
              // parameter. The additional digits shown here were obtained using a code in the
              // mp-quadrature library, https://github.com/jwpeterson/mp-quadrature
              const Real rule_data[nrows][3] = {
                {6.2382265094402118173683000996350e-02_R, 6.7517867073916085442557131050869e-02_R, 2.6517028157436251428754180460739e-02_R}, // group A
                {5.5225456656926611737479190275645e-02_R, 3.2150249385198182266630784919920e-01_R, 4.3881408714446055036769903139288e-02_R}, // group B
                {3.4324302945097146469630642483938e-02_R, 6.6094919618673565761198031019780e-01_R, 2.8775042784981585738445496900219e-02_R}, // group C
                {5.1584233435359177925746338682643e-01_R, 2.7771616697639178256958187139372e-01_R, 6.7493187009802774462697086166421e-02_R}  // group D
              };
 
              for (unsigned int i=0, offset=0; i<nrows; ++i)
                {
                  _points[offset + 0] = Point(rule_data[i][0],                    rule_data[i][1]); // (z1,z2)
                  _points[offset + 1] = Point(1.-rule_data[i][0]-rule_data[i][1], rule_data[i][0]); // (z3,z1)
                  _points[offset + 2] = Point(rule_data[i][1], 1.-rule_data[i][0]-rule_data[i][1]); // (z2,z3)
 
                  // All these points get the same weight
                  _weights[offset + 0] = rule_data[i][2];
                  _weights[offset + 1] = rule_data[i][2];
                  _weights[offset + 2] = rule_data[i][2];
 
                  // Increment offset
                  offset += 3;
                }
 
              return;
 
 
              //       // The following is an inferior 7th-order Lyness-style rule with 15 points.
              //       // It's here only for completeness and the Ro3-invariant rule above should
              //       // be used instead!
              //       const unsigned int n_wts = 3;
              //       const Real wts[n_wts] =
              // {
              //   2.6538900895116205835977487499847719e-02_R,
              //   3.5426541846066783659206291623201826e-02_R,
              //   3.4637341039708446756138297960207647e-02_R
              // };
              //
              //       const Real a[n_wts] =
              // {
              //   6.4930513159164863078379776030396538e-02_R,
              //   2.8457558424917033519741605734978046e-01_R,
              //   3.1355918438493150795585190219862865e-01_R
              // };
              //
              //       const Real b[n_wts] =
              // {
              //   0.,
              //   1.9838447668150671917987659863332941e-01_R,
              //   4.3863471792372471511798695971295936e-02_R
              // };
              //
              //       const unsigned int permutation_ids[n_wts] = {3, 6, 6}; // 15 total points
              //
              //       dunavant_rule2(wts, a, b, permutation_ids, n_wts);
              //
              //       return;
            }
 
 
 
 
            // Another Dunavant rule.  This one has all positive weights.  This rule has
            // 16 points while a comparable conical product rule would have 5*5=25.
            //
            // It was copied 23rd June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
            //
            // Additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
            // on triangles and tetrahedra"  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
          case EIGHTH:
            {
              const unsigned int n_wts = 5;
              const Real wts[n_wts] =
                {
                  7.2157803838893584125545555244532310e-02_R,
                  4.7545817133642312396948052194292159e-02_R,
                  5.1608685267359125140895775146064515e-02_R,
                  1.6229248811599040155462964170890299e-02_R,
                  1.3615157087217497132422345036954462e-02_R
                };
 
              const Real a[n_wts] =
                {
                  0.0, // 'a' parameter not used for origin permutation
                  4.5929258829272315602881551449416932e-01_R,
                  1.7056930775176020662229350149146450e-01_R,
                  5.0547228317030975458423550596598947e-02_R,
                  2.6311282963463811342178578628464359e-01_R,
                };
 
              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  0.,
                  0.,
                  7.2849239295540428124100037917606196e-01_R
                };
 
              const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 6}; // 16 total points
 
              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
 
              return;
            }
 
 
 
            // Another Dunavant rule.  This one has all positive weights.  This rule has 19
            // points. The comparable conical product rule would have 25.
            // It was copied 23rd June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
            //
            // Additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
            // on triangles and tetrahedra"  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
          case NINTH:
            {
              const unsigned int n_wts = 6;
              const Real wts[n_wts] =
                {
                  4.8567898141399416909620991253644315e-02_R,
                  1.5667350113569535268427415643604658e-02_R,
                  1.2788837829349015630839399279499912e-02_R,
                  3.8913770502387139658369678149701978e-02_R,
                  3.9823869463605126516445887132022637e-02_R,
                  2.1641769688644688644688644688644689e-02_R
                };
 
              const Real a[n_wts] =
                {
                  0.0, // 'a' parameter not used for origin permutation
                  4.8968251919873762778370692483619280e-01_R,
                  4.4729513394452709865106589966276365e-02_R,
                  4.3708959149293663726993036443535497e-01_R,
                  1.8820353561903273024096128046733557e-01_R,
                  2.2196298916076569567510252769319107e-01_R
                };
 
              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  0.,
                  0.,
                  0.,
                  7.4119859878449802069007987352342383e-01_R
                };
 
              const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 3, 6}; // 19 total points
 
              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
 
              return;
            }
 
 
            // Another Dunavant rule with all positive weights.  This rule has 25
            // points. The comparable conical product rule would have 36.
            // It was copied 23rd June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
            //
            // Additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
            // on triangles and tetrahedra"  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
          case TENTH:
            {
              const unsigned int n_wts = 6;
              const Real wts[n_wts] =
                {
                  4.5408995191376790047643297550014267e-02_R,
                  1.8362978878233352358503035945683300e-02_R,
                  2.2660529717763967391302822369298659e-02_R,
                  3.6378958422710054302157588309680344e-02_R,
                  1.4163621265528742418368530791049552e-02_R,
                  4.7108334818664117299637354834434138e-03_R
                };
 
              const Real a[n_wts] =
                {
                  0.0, // 'a' parameter not used for origin permutation
                  4.8557763338365737736750753220812615e-01_R,
                  1.0948157548503705479545863134052284e-01_R,
                  3.0793983876412095016515502293063162e-01_R,
                  2.4667256063990269391727646541117681e-01_R,
                  6.6803251012200265773540212762024737e-02_R
                };
 
              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  0.,
                  5.5035294182099909507816172659300821e-01_R,
                  7.2832390459741092000873505358107866e-01_R,
                  9.2365593358750027664630697761508843e-01_R
                };
 
              const unsigned int permutation_ids[n_wts] = {1, 3, 3, 6, 6, 6}; // 25 total points
 
              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
 
              return;
            }
 
 
            // Dunavant's 11th-order rule contains points outside the region of
            // integration, and is thus unacceptable for our FEM calculations.
            //
            // This 30-point, 11th-order rule was obtained by me [JWP] using the code in
            //
            // Additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
            // on triangles and tetrahedra"  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
            //
            // Note: the 28-point 11th-order rule obtained by Zhang in the paper above
            // does not appear to be unique.  It is a solution in the sense that it
            // minimizes the error in the least-squares minimization problem, but
            // it involves too many unknowns and the Jacobian is therefore singular
            // when attempting to improve the solution via Newton's method.
          case ELEVENTH:
            {
              const unsigned int n_wts = 6;
              const Real wts[n_wts] =
                {
                  3.6089021198604635216985338480426484e-02_R,
                  2.1607717807680420303346736867931050e-02_R,
                  3.1144524293927978774861144478241807e-03_R,
                  2.9086855161081509446654185084988077e-02_R,
                  8.4879241614917017182977532679947624e-03_R,
                  1.3795732078224796530729242858347546e-02_R
                };
 
              const Real a[n_wts] =
                {
                  3.9355079629947969884346551941969960e-01_R,
                  4.7979065808897448654107733982929214e-01_R,
                  5.1003445645828061436081405648347852e-03_R,
                  2.6597620190330158952732822450744488e-01_R,
                  2.8536418538696461608233522814483715e-01_R,
                  1.3723536747817085036455583801851025e-01_R
                };
 
              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  5.6817155788572446538150614865768991e-02_R,
                  1.2539956353662088473247489775203396e-01_R,
                  1.2409970153698532116262152247041742e-02_R,
                  5.2792057988217708934207928630851643e-02_R
                };
 
              const unsigned int permutation_ids[n_wts] = {3, 3, 6, 6, 6, 6}; // 30 total points
 
              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
 
              return;
            }
 
 
 
 
            // Another Dunavant rule with all positive weights.  This rule has 33
            // points. The comparable conical product rule would have 36 (ELEVENTH) or 49 (TWELFTH).
            //
            // It was copied 23rd June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
            //
            // Additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
            // on triangles and tetrahedra"  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
          case TWELFTH:
            {
              const unsigned int n_wts = 8;
              const Real wts[n_wts] =
                {
                  3.0831305257795086169332418926151771e-03_R,
                  3.1429112108942550177135256546441273e-02_R,
                  1.7398056465354471494664198647499687e-02_R,
                  2.1846272269019201067728631278737487e-02_R,
                  1.2865533220227667708895461535782215e-02_R,
                  1.1178386601151722855919538351159995e-02_R,
                  8.6581155543294461858210504055170332e-03_R,
                  2.0185778883190464758914349626118386e-02_R
                };
 
              const Real a[n_wts] =
                {
                  2.1317350453210370246856975515728246e-02_R,
                  2.7121038501211592234595134039689474e-01_R,
                  1.2757614554158592467389632515428357e-01_R,
                  4.3972439229446027297973662348436108e-01_R,
                  4.8821738977380488256466206525881104e-01_R,
                  2.8132558098993954824813069297455275e-01_R,
                  1.1625191590759714124135414784260182e-01_R,
                  2.7571326968551419397479634607976398e-01_R
                };
 
              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  0.,
                  0.,
                  0.,
                  6.9583608678780342214163552323607254e-01_R,
                  8.5801403354407263059053661662617818e-01_R,
                  6.0894323577978780685619243776371007e-01_R
                };
 
              const unsigned int permutation_ids[n_wts] = {3, 3, 3, 3, 3, 6, 6, 6}; // 33 total points
 
              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
 
              return;
            }
 
 
            // Another Dunavant rule with all positive weights.  This rule has 37
            // points. The comparable conical product rule would have 49 points.
            //
            // It was copied 23rd June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
            //
            // A second rule with additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
            // on triangles and tetrahedra"  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
          case THIRTEENTH:
            {
              const unsigned int n_wts = 9;
              const Real wts[n_wts] =
                {
                  3.3980018293415822140887212340442440e-02_R,
                  2.7800983765226664353628733005230734e-02_R,
                  2.9139242559599990702383541756669905e-02_R,
                  3.0261685517695859208964000161454122e-03_R,
                  1.1997200964447365386855399725479827e-02_R,
                  1.7320638070424185232993414255459110e-02_R,
                  7.4827005525828336316229285664517190e-03_R,
                  1.2089519905796909568722872786530380e-02_R,
                  4.7953405017716313612975450830554457e-03_R
                };
 
              const Real a[n_wts] =
                {
                  0., // 'a' parameter not used for origin permutation
                  4.2694141425980040602081253503137421e-01_R,
                  2.2137228629183290065481255470507908e-01_R,
                  2.1509681108843183869291313534052083e-02_R,
                  4.8907694645253934990068971909020439e-01_R,
                  3.0844176089211777465847185254124531e-01_R,
                  1.1092204280346339541286954522167452e-01_R,
                  1.6359740106785048023388790171095725e-01_R,
                  2.7251581777342966618005046435408685e-01_R
                };
 
              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  0.,
                  0.,
                  0.,
                  6.2354599555367557081585435318623659e-01_R,
                  8.6470777029544277530254595089569318e-01_R,
                  7.4850711589995219517301859578870965e-01_R,
                  7.2235779312418796526062013230478405e-01_R
                };
 
              const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 3, 6, 6, 6, 6}; // 37 total points
 
              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
 
              return;
            }
 
 
            // Another Dunavant rule.  This rule has 42 points, while
            // a comparable conical product rule would have 64.
            //
            // It was copied 23rd June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
            //
            // Additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
            // on triangles and tetrahedra"  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
          case FOURTEENTH:
            {
              const unsigned int n_wts = 10;
              const Real wts[n_wts] =
                {
                  1.0941790684714445320422472981662986e-02_R,
                  1.6394176772062675320655489369312672e-02_R,
                  2.5887052253645793157392455083198201e-02_R,
                  2.1081294368496508769115218662093065e-02_R,
                  7.2168498348883338008549607403266583e-03_R,
                  2.4617018012000408409130117545210774e-03_R,
                  1.2332876606281836981437622591818114e-02_R,
                  1.9285755393530341614244513905205430e-02_R,
                  7.2181540567669202480443459995079017e-03_R,
                  2.5051144192503358849300465412445582e-03_R
                };
 
              const Real a[n_wts] =
                {
                  4.8896391036217863867737602045239024e-01_R,
                  4.1764471934045392250944082218564344e-01_R,
                  2.7347752830883865975494428326269856e-01_R,
                  1.7720553241254343695661069046505908e-01_R,
                  6.1799883090872601267478828436935788e-02_R,
                  1.9390961248701048178250095054529511e-02_R,
                  1.7226668782135557837528960161365733e-01_R,
                  3.3686145979634500174405519708892539e-01_R,
                  2.9837288213625775297083151805961273e-01_R,
                  1.1897449769695684539818196192990548e-01_R
                };
 
              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  0.,
                  0.,
                  0.,
                  0.,
                  7.7060855477499648258903327416742796e-01_R,
                  5.7022229084668317349769621336235426e-01_R,
                  6.8698016780808783735862715402031306e-01_R,
                  8.7975717137017112951457163697460183e-01_R
                };
 
              const unsigned int permutation_ids[n_wts]
                = {3, 3, 3, 3, 3, 3, 6, 6, 6, 6}; // 42 total points
 
              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
 
              return;
            }
 
 
            // This 49-point rule was found by me [JWP] using the code in:
            //
            // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
            // on triangles and tetrahedra"  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
            //
            // A 54-point, 15th-order rule is reported by
            //
            // Stephen Wandzura, Hong Xiao,
            // Symmetric Quadrature Rules on a Triangle,
            // Computers and Mathematics with Applications,
            // Volume 45, Number 12, June 2003, pages 1829-1840.
            //
            // can be found here:
            // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
            //
            // but this 49-point rule is superior.
          case FIFTEENTH:
            {
              const unsigned int n_wts = 11;
              const Real wts[n_wts] =
                {
                  2.4777380743035579804788826970198951e-02_R,
                  9.2433943023307730591540642828347660e-03_R,
                  2.2485768962175402793245929133296627e-03_R,
                  6.7052581900064143760518398833360903e-03_R,
                  1.9011381726930579256700190357527956e-02_R,
                  1.4605445387471889398286155981802858e-02_R,
                  1.5087322572773133722829435011138258e-02_R,
                  1.5630213780078803020711746273129099e-02_R,
                  6.1808086085778203192616856133701233e-03_R,
                  3.2209366452594664857296985751120513e-03_R,
                  5.8747373242569702667677969985668817e-03_R
                };
 
              const Real a[n_wts] =
                {
                  0.0, // 'a' parameter not used for origin
                  7.9031013655541635005816956762252155e-02_R,
                  1.8789501810770077611247984432284226e-02_R,
                  4.9250168823249670532514526605352905e-01_R,
                  4.0886316907744105975059040108092775e-01_R,
                  5.3877851064220142445952549348423733e-01_R,
                  2.0250549804829997692885033941362673e-01_R,
                  5.5349674918711643207148086558288110e-01_R,
                  7.8345022567320812359258882143250181e-01_R,
                  8.9514624528794883409864566727625002e-01_R,
                  3.2515745241110782862789881780746490e-01_R
                };
 
              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  0.,
                  0.,
                  0.,
                  1.9412620368774630292701241080996842e-01_R,
                  9.8765911355712115933807754318089099e-02_R,
                  7.7663767064308164090246588765178087e-02_R,
                  2.1594628433980258573654682690950798e-02_R,
                  1.2563596287784997705599005477153617e-02_R,
                  1.5082654870922784345283124845552190e-02_R
                };
 
              const unsigned int permutation_ids[n_wts]
                = {1, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6}; // 49 total points
 
              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
 
              return;
            }
 
 
 
 
            // Dunavant's 16th-order rule contains points outside the region of
            // integration, and is thus unacceptable for our FEM calculations.
            //
            // This 55-point, 16th-order rule was obtained by me [JWP] using the code in
            //
            // Additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
            // on triangles and tetrahedra"  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
            //
            // Note: the 55-point 16th-order rule obtained by Zhang in the paper above
            // does not appear to be unique.  It is a solution in the sense that it
            // minimizes the error in the least-squares minimization problem, but
            // it involves too many unknowns and the Jacobian is therefore singular
            // when attempting to improve the solution via Newton's method.
          case SIXTEENTH:
            {
              const unsigned int n_wts = 12;
              const Real wts[n_wts] =
                {
                  2.2668082505910087151996321171534230e-02_R,
                  8.4043060714818596159798961899306135e-03_R,
                  1.0850949634049747713966288634484161e-03_R,
                  7.2252773375423638869298219383808751e-03_R,
                  1.2997715227338366024036316182572871e-02_R,
                  2.0054466616677715883228810959112227e-02_R,
                  9.7299841600417010281624372720122710e-03_R,
                  1.1651974438298104227427176444311766e-02_R,
                  9.1291185550484450744725847363097389e-03_R,
                  3.5568614040947150231712567900113671e-03_R,
                  5.8355861686234326181790822005304303e-03_R,
                  4.7411314396804228041879331486234396e-03_R
                };
 
              const Real a[n_wts] =
                {
                  0.0, // 'a' parameter not used for centroid weight
                  8.5402539407933203673769900926355911e-02_R,
                  1.2425572001444092841183633409631260e-02_R,
                  4.9174838341891594024701017768490960e-01_R,
                  4.5669426695387464162068900231444462e-01_R,
                  4.8506759880447437974189793537259677e-01_R,
                  2.0622099278664205707909858461264083e-01_R,
                  3.2374950270039093446805340265853956e-01_R,
                  7.3834330556606586255186213302750029e-01_R,
                  9.1210673061680792565673823935174611e-01_R,
                  6.6129919222598721544966837350891531e-01_R,
                  1.7807138906021476039088828811346122e-01_R
                };
 
              const Real b[n_wts] =
                {
                  0.0,
                  0.0,
                  0.0,
                  0.0,
                  0.0,
                  3.2315912848634384647700266402091638e-01_R,
                  1.5341553679414688425981898952416987e-01_R,
                  7.4295478991330687632977899141707872e-02_R,
                  7.1278762832147862035977841733532020e-02_R,
                  1.6623223223705792825395256602140459e-02_R,
                  1.4160772533794791868984026749196156e-02_R,
                  1.4539694958941854654807449467759690e-02_R
                };
 
              const unsigned int permutation_ids[n_wts]
                = {1, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6, 6}; // 55 total points
 
              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
 
              return;
            }
 
 
            // Dunavant's 17th-order rule has 61 points, while a
            // comparable conical product rule would have 81 (16th and 17th orders).
            //
            // It can be found here:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
            //
            // Zhang reports an identical rule in:
            // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
            // on triangles and tetrahedra"  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
            //
            // Note: the 61-point 17th-order rule obtained by Dunavant and Zhang
            // does not appear to be unique.  It is a solution in the sense that it
            // minimizes the error in the least-squares minimization problem, but
            // it involves too many unknowns and the Jacobian is therefore singular
            // when attempting to improve the solution via Newton's method.
            //
            // Therefore, we prefer the following 63-point rule which
            // I [JWP] found.  It appears to be more accurate than the
            // rule reported by Dunavant and Zhang, even though it has
            // a few more points.
          case SEVENTEENTH:
            {
              const unsigned int n_wts = 12;
              const Real wts[n_wts] =
                {
                  1.7464603792572004485690588092246146e-02_R,
                  5.9429003555801725246549713984660076e-03_R,
                  1.2490753345169579649319736639588729e-02_R,
                  1.5386987188875607593083456905596468e-02_R,
                  1.1185807311917706362674684312990270e-02_R,
                  1.0301845740670206831327304917180007e-02_R,
                  1.1767783072977049696840016810370464e-02_R,
                  3.8045312849431209558329128678945240e-03_R,
                  4.5139302178876351271037137230354382e-03_R,
                  2.2178812517580586419412547665472893e-03_R,
                  5.2216271537483672304731416553063103e-03_R,
                  9.8381136389470256422419930926212114e-04_R
                };
 
              const Real a[n_wts] =
                {
                  2.8796825754667362165337965123570514e-01_R,
                  4.9216175986208465345536805750663939e-01_R,
                  4.6252866763171173685916780827044612e-01_R,
                  1.6730292951631792248498303276090273e-01_R,
                  1.5816335500814652972296428532213019e-01_R,
                  1.6352252138387564873002458959679529e-01_R,
                  6.2447680488959768233910286168417367e-01_R,
                  8.7317249935244454285263604347964179e-01_R,
                  3.4428164322282694677972239461699271e-01_R,
                  9.1584484467813674010523309855340209e-02_R,
                  2.0172088013378989086826623852040632e-01_R,
                  9.6538762758254643474731509845084691e-01_R
                };
 
              const Real b[n_wts] =
                {
                  0.0,
                  0.0,
                  0.0,
                  3.4429160695501713926320695771253348e-01_R,
                  2.2541623431550639817203145525444726e-01_R,
                  8.0670083153531811694942222940484991e-02_R,
                  6.5967451375050925655738829747288190e-02_R,
                  4.5677879890996762665044366994439565e-02_R,
                  1.1528411723154215812386518751976084e-02_R,
                  9.3057714323900610398389176844165892e-03_R,
                  1.5916814107619812717966560404970160e-02_R,
                  1.0734733163764032541125434215228937e-02_R
                };
 
              const unsigned int permutation_ids[n_wts]
                = {3, 3, 3, 6, 6, 6, 6, 6, 6, 6, 6, 6}; // 63 total points
 
              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
 
              return;
 
              //       _points.resize (61);
              //       _weights.resize(61);
 
              //       // The raw data for the quadrature rule.
              //       const Real p[15][4] = {
              // {                1./3.,                    0.,                    0., 0.033437199290803e+00 / 2.0}, // 1-perm
              // {0.005658918886452e+00, 0.497170540556774e+00,                    0., 0.005093415440507e+00 / 2.0}, // 3-perm
              // {0.035647354750751e+00, 0.482176322624625e+00,                    0., 0.014670864527638e+00 / 2.0}, // 3-perm
              // {0.099520061958437e+00, 0.450239969020782e+00,                    0., 0.024350878353672e+00 / 2.0}, // 3-perm
              // {0.199467521245206e+00, 0.400266239377397e+00,                    0., 0.031107550868969e+00 / 2.0}, // 3-perm
              // {0.495717464058095e+00, 0.252141267970953e+00,                    0., 0.031257111218620e+00 / 2.0}, // 3-perm
              // {0.675905990683077e+00, 0.162047004658461e+00,                    0., 0.024815654339665e+00 / 2.0}, // 3-perm
              // {0.848248235478508e+00, 0.075875882260746e+00,                    0., 0.014056073070557e+00 / 2.0}, // 3-perm
              // {0.968690546064356e+00, 0.015654726967822e+00,                    0., 0.003194676173779e+00 / 2.0}, // 3-perm
              // {0.010186928826919e+00, 0.334319867363658e+00, 0.655493203809423e+00, 0.008119655318993e+00 / 2.0}, // 6-perm
              // {0.135440871671036e+00, 0.292221537796944e+00, 0.572337590532020e+00, 0.026805742283163e+00 / 2.0}, // 6-perm
              // {0.054423924290583e+00, 0.319574885423190e+00, 0.626001190286228e+00, 0.018459993210822e+00 / 2.0}, // 6-perm
              // {0.012868560833637e+00, 0.190704224192292e+00, 0.796427214974071e+00, 0.008476868534328e+00 / 2.0}, // 6-perm
              // {0.067165782413524e+00, 0.180483211648746e+00, 0.752351005937729e+00, 0.018292796770025e+00 / 2.0}, // 6-perm
              // {0.014663182224828e+00, 0.080711313679564e+00, 0.904625504095608e+00, 0.006665632004165e+00 / 2.0}  // 6-perm
              //       };
 
 
              //       // Now call the dunavant routine to generate _points and _weights
              //       dunavant_rule(p, 15);
 
              //       return;
            }
 
 
 
            // Dunavant's 18th-order rule contains points outside the region and is therefore unsuitable
            // for our FEM calculations.  His 19th-order rule has 73 points, compared with 100 points for
            // a comparable-order conical product rule.
            //
            // It was copied 23rd June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
          case EIGHTTEENTH:
          case NINETEENTH:
            {
              _points.resize (73);
              _weights.resize(73);
 
              // The raw data for the quadrature rule.
              const Real rule_data[17][4] = {
                {                1./3.,                    0.,                    0., 0.032906331388919e+00 / 2.0}, // 1-perm
                {0.020780025853987e+00, 0.489609987073006e+00,                    0., 0.010330731891272e+00 / 2.0}, // 3-perm
                {0.090926214604215e+00, 0.454536892697893e+00,                    0., 0.022387247263016e+00 / 2.0}, // 3-perm
                {0.197166638701138e+00, 0.401416680649431e+00,                    0., 0.030266125869468e+00 / 2.0}, // 3-perm
                {0.488896691193805e+00, 0.255551654403098e+00,                    0., 0.030490967802198e+00 / 2.0}, // 3-perm
                {0.645844115695741e+00, 0.177077942152130e+00,                    0., 0.024159212741641e+00 / 2.0}, // 3-perm
                {0.779877893544096e+00, 0.110061053227952e+00,                    0., 0.016050803586801e+00 / 2.0}, // 3-perm
                {0.888942751496321e+00, 0.055528624251840e+00,                    0., 0.008084580261784e+00 / 2.0}, // 3-perm
                {0.974756272445543e+00, 0.012621863777229e+00,                    0., 0.002079362027485e+00 / 2.0}, // 3-perm
                {0.003611417848412e+00, 0.395754787356943e+00, 0.600633794794645e+00, 0.003884876904981e+00 / 2.0}, // 6-perm
                {0.134466754530780e+00, 0.307929983880436e+00, 0.557603261588784e+00, 0.025574160612022e+00 / 2.0}, // 6-perm
                {0.014446025776115e+00, 0.264566948406520e+00, 0.720987025817365e+00, 0.008880903573338e+00 / 2.0}, // 6-perm
                {0.046933578838178e+00, 0.358539352205951e+00, 0.594527068955871e+00, 0.016124546761731e+00 / 2.0}, // 6-perm
                {0.002861120350567e+00, 0.157807405968595e+00, 0.839331473680839e+00, 0.002491941817491e+00 / 2.0}, // 6-perm
                {0.223861424097916e+00, 0.075050596975911e+00, 0.701087978926173e+00, 0.018242840118951e+00 / 2.0}, // 6-perm
                {0.034647074816760e+00, 0.142421601113383e+00, 0.822931324069857e+00, 0.010258563736199e+00 / 2.0}, // 6-perm
                {0.010161119296278e+00, 0.065494628082938e+00, 0.924344252620784e+00, 0.003799928855302e+00 / 2.0}  // 6-perm
              };
 
 
              // Now call the dunavant routine to generate _points and _weights
              dunavant_rule(rule_data, 17);
 
              return;
            }
 
 
            // 20th-order rule by Wandzura.
            //
            // Stephen Wandzura, Hong Xiao,
            // Symmetric Quadrature Rules on a Triangle,
            // Computers and Mathematics with Applications,
            // Volume 45, Number 12, June 2003, pages 1829-1840.
            //
            // Wandzura's work extends the work of Dunavant by providing degree
            // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
            //
            // Copied on 3rd July 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
          case TWENTIETH:
            {
              // The equivalent conical product rule would have 121 points
              _points.resize (85);
              _weights.resize(85);
 
              // The raw data for the quadrature rule.
              const Real rule_data[19][4] = {
                {0.33333333333333e+00,                  0.0,                  0.0, 0.2761042699769952e-01 / 2.0}, // 1-perm
                {0.00150064932443e+00, 0.49924967533779e+00,                  0.0, 0.1779029547326740e-02 / 2.0}, // 3-perm
                {0.09413975193895e+00, 0.45293012403052e+00,                  0.0, 0.2011239811396117e-01 / 2.0}, // 3-perm
                {0.20447212408953e+00, 0.39776393795524e+00,                  0.0, 0.2681784725933157e-01 / 2.0}, // 3-perm
                {0.47099959493443e+00, 0.26450020253279e+00,                  0.0, 0.2452313380150201e-01 / 2.0}, // 3-perm
                {0.57796207181585e+00, 0.21101896409208e+00,                  0.0, 0.1639457841069539e-01 / 2.0}, // 3-perm
                {0.78452878565746e+00, 0.10773560717127e+00,                  0.0, 0.1479590739864960e-01 / 2.0}, // 3-perm
                {0.92186182432439e+00, 0.03906908783780e+00,                  0.0, 0.4579282277704251e-02 / 2.0}, // 3-perm
                {0.97765124054134e+00, 0.01117437972933e+00,                  0.0, 0.1651826515576217e-02 / 2.0}, // 3-perm
                {0.00534961818734e+00, 0.06354966590835e+00, 0.93110071590431e+00, 0.2349170908575584e-02 / 2.0}, // 6-perm
                {0.00795481706620e+00, 0.15710691894071e+00, 0.83493826399309e+00, 0.4465925754181793e-02 / 2.0}, // 6-perm
                {0.01042239828126e+00, 0.39564211436437e+00, 0.59393548735436e+00, 0.6099566807907972e-02 / 2.0}, // 6-perm
                {0.01096441479612e+00, 0.27316757071291e+00, 0.71586801449097e+00, 0.6891081327188203e-02 / 2.0}, // 6-perm
                {0.03856671208546e+00, 0.10178538248502e+00, 0.85964790542952e+00, 0.7997475072478163e-02 / 2.0}, // 6-perm
                {0.03558050781722e+00, 0.44665854917641e+00, 0.51776094300637e+00, 0.7386134285336024e-02 / 2.0}, // 6-perm
                {0.04967081636276e+00, 0.19901079414950e+00, 0.75131838948773e+00, 0.1279933187864826e-01 / 2.0}, // 6-perm
                {0.05851972508433e+00, 0.32426118369228e+00, 0.61721909122339e+00, 0.1725807117569655e-01 / 2.0}, // 6-perm
                {0.12149778700439e+00, 0.20853136321013e+00, 0.66997084978547e+00, 0.1867294590293547e-01 / 2.0}, // 6-perm
                {0.14071084494394e+00, 0.32317056653626e+00, 0.53611858851980e+00, 0.2281822405839526e-01 / 2.0}  // 6-perm
              };
 
 
              // Now call the dunavant routine to generate _points and _weights
              dunavant_rule(rule_data, 19);
 
              return;
            }
 
 
 
            // 25th-order rule by Wandzura.
            //
            // Stephen Wandzura, Hong Xiao,
            // Symmetric Quadrature Rules on a Triangle,
            // Computers and Mathematics with Applications,
            // Volume 45, Number 12, June 2003, pages 1829-1840.
            //
            // Wandzura's work extends the work of Dunavant by providing degree
            // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
            //
            // Copied on 3rd July 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
            // case TWENTYFIRST: // fall through to 121 point conical product rule below
          case TWENTYSECOND:
          case TWENTYTHIRD:
          case TWENTYFOURTH:
          case TWENTYFIFTH:
            {
              // The equivalent conical product rule would have 169 points
              _points.resize (126);
              _weights.resize(126);
 
              // The raw data for the quadrature rule.
              const Real rule_data[26][4] = {
                {0.02794648307317e+00, 0.48602675846341e+00,                  0.0, 0.8005581880020417e-02 / 2.0},  // 3-perm
                {0.13117860132765e+00, 0.43441069933617e+00,                  0.0, 0.1594707683239050e-01 / 2.0},  // 3-perm
                {0.22022172951207e+00, 0.38988913524396e+00,                  0.0, 0.1310914123079553e-01 / 2.0},  // 3-perm
                {0.40311353196039e+00, 0.29844323401980e+00,                  0.0, 0.1958300096563562e-01 / 2.0},  // 3-perm
                {0.53191165532526e+00, 0.23404417233737e+00,                  0.0, 0.1647088544153727e-01 / 2.0},  // 3-perm
                {0.69706333078196e+00, 0.15146833460902e+00,                  0.0, 0.8547279074092100e-02 / 2.0},  // 3-perm
                {0.77453221290801e+00, 0.11273389354599e+00,                  0.0, 0.8161885857226492e-02 / 2.0},  // 3-perm
                {0.84456861581695e+00, 0.07771569209153e+00,                  0.0, 0.6121146539983779e-02 / 2.0},  // 3-perm
                {0.93021381277141e+00, 0.03489309361430e+00,                  0.0, 0.2908498264936665e-02 / 2.0},  // 3-perm
                {0.98548363075813e+00, 0.00725818462093e+00,                  0.0, 0.6922752456619963e-03 / 2.0},  // 3-perm
                {0.00129235270444e+00, 0.22721445215336e+00, 0.77149319514219e+00, 0.1248289199277397e-02 / 2.0},  // 6-perm
                {0.00539970127212e+00, 0.43501055485357e+00, 0.55958974387431e+00, 0.3404752908803022e-02 / 2.0},  // 6-perm
                {0.00638400303398e+00, 0.32030959927220e+00, 0.67330639769382e+00, 0.3359654326064051e-02 / 2.0},  // 6-perm
                {0.00502821150199e+00, 0.09175032228001e+00, 0.90322146621800e+00, 0.1716156539496754e-02 / 2.0},  // 6-perm
                {0.00682675862178e+00, 0.03801083585872e+00, 0.95516240551949e+00, 0.1480856316715606e-02 / 2.0},  // 6-perm
                {0.01001619963993e+00, 0.15742521848531e+00, 0.83255858187476e+00, 0.3511312610728685e-02 / 2.0},  // 6-perm
                {0.02575781317339e+00, 0.23988965977853e+00, 0.73435252704808e+00, 0.7393550149706484e-02 / 2.0},  // 6-perm
                {0.03022789811992e+00, 0.36194311812606e+00, 0.60782898375402e+00, 0.7983087477376558e-02 / 2.0},  // 6-perm
                {0.03050499010716e+00, 0.08355196095483e+00, 0.88594304893801e+00, 0.4355962613158041e-02 / 2.0},  // 6-perm
                {0.04595654736257e+00, 0.14844322073242e+00, 0.80560023190501e+00, 0.7365056701417832e-02 / 2.0},  // 6-perm
                {0.06744280054028e+00, 0.28373970872753e+00, 0.64881749073219e+00, 0.1096357284641955e-01 / 2.0},  // 6-perm
                {0.07004509141591e+00, 0.40689937511879e+00, 0.52305553346530e+00, 0.1174996174354112e-01 / 2.0},  // 6-perm
                {0.08391152464012e+00, 0.19411398702489e+00, 0.72197448833499e+00, 0.1001560071379857e-01 / 2.0},  // 6-perm
                {0.12037553567715e+00, 0.32413434700070e+00, 0.55549011732214e+00, 0.1330964078762868e-01 / 2.0},  // 6-perm
                {0.14806689915737e+00, 0.22927748355598e+00, 0.62265561728665e+00, 0.1415444650522614e-01 / 2.0},  // 6-perm
                {0.19177186586733e+00, 0.32561812259598e+00, 0.48261001153669e+00, 0.1488137956116801e-01 / 2.0}   // 6-perm
              };
 
 
              // Now call the dunavant routine to generate _points and _weights
              dunavant_rule(rule_data, 26);
 
              return;
            }
 
 
 
            // 30th-order rule by Wandzura.
            //
            // Stephen Wandzura, Hong Xiao,
            // Symmetric Quadrature Rules on a Triangle,
            // Computers and Mathematics with Applications,
            // Volume 45, Number 12, June 2003, pages 1829-1840.
            //
            // Wandzura's work extends the work of Dunavant by providing degree
            // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
            //
            // Copied on 3rd July 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
          case TWENTYSIXTH:
          case TWENTYSEVENTH:
          case TWENTYEIGHTH:
          case TWENTYNINTH:
          case THIRTIETH:
            {
              // The equivalent conical product rule would have 256 points
              _points.resize (175);
              _weights.resize(175);
 
              // The raw data for the quadrature rule.
              const Real rule_data[36][4] = {
                {0.33333333333333e+00,                  0.0,                  0.0, 0.1557996020289920e-01 / 2.0}, // 1-perm
                {0.00733011643277e+00, 0.49633494178362e+00,                  0.0, 0.3177233700534134e-02 / 2.0}, // 3-perm
                {0.08299567580296e+00, 0.45850216209852e+00,                  0.0, 0.1048342663573077e-01 / 2.0}, // 3-perm
                {0.15098095612541e+00, 0.42450952193729e+00,                  0.0, 0.1320945957774363e-01 / 2.0}, // 3-perm
                {0.23590585989217e+00, 0.38204707005392e+00,                  0.0, 0.1497500696627150e-01 / 2.0}, // 3-perm
                {0.43802430840785e+00, 0.28098784579608e+00,                  0.0, 0.1498790444338419e-01 / 2.0}, // 3-perm
                {0.54530204829193e+00, 0.22734897585403e+00,                  0.0, 0.1333886474102166e-01 / 2.0}, // 3-perm
                {0.65088177698254e+00, 0.17455911150873e+00,                  0.0, 0.1088917111390201e-01 / 2.0}, // 3-perm
                {0.75348314559713e+00, 0.12325842720144e+00,                  0.0, 0.8189440660893461e-02 / 2.0}, // 3-perm
                {0.83983154221561e+00, 0.08008422889220e+00,                  0.0, 0.5575387588607785e-02 / 2.0}, // 3-perm
                {0.90445106518420e+00, 0.04777446740790e+00,                  0.0, 0.3191216473411976e-02 / 2.0}, // 3-perm
                {0.95655897063972e+00, 0.02172051468014e+00,                  0.0, 0.1296715144327045e-02 / 2.0}, // 3-perm
                {0.99047064476913e+00, 0.00476467761544e+00,                  0.0, 0.2982628261349172e-03 / 2.0}, // 3-perm
                {0.00092537119335e+00, 0.41529527091331e+00, 0.58377935789334e+00, 0.9989056850788964e-03 / 2.0}, // 6-perm
                {0.00138592585556e+00, 0.06118990978535e+00, 0.93742416435909e+00, 0.4628508491732533e-03 / 2.0}, // 6-perm
                {0.00368241545591e+00, 0.16490869013691e+00, 0.83140889440718e+00, 0.1234451336382413e-02 / 2.0}, // 6-perm
                {0.00390322342416e+00, 0.02503506223200e+00, 0.97106171434384e+00, 0.5707198522432062e-03 / 2.0}, // 6-perm
                {0.00323324815501e+00, 0.30606446515110e+00, 0.69070228669389e+00, 0.1126946125877624e-02 / 2.0}, // 6-perm
                {0.00646743211224e+00, 0.10707328373022e+00, 0.88645928415754e+00, 0.1747866949407337e-02 / 2.0}, // 6-perm
                {0.00324747549133e+00, 0.22995754934558e+00, 0.76679497516308e+00, 0.1182818815031657e-02 / 2.0}, // 6-perm
                {0.00867509080675e+00, 0.33703663330578e+00, 0.65428827588746e+00, 0.1990839294675034e-02 / 2.0}, // 6-perm
                {0.01559702646731e+00, 0.05625657618206e+00, 0.92814639735063e+00, 0.1900412795035980e-02 / 2.0}, // 6-perm
                {0.01797672125369e+00, 0.40245137521240e+00, 0.57957190353391e+00, 0.4498365808817451e-02 / 2.0}, // 6-perm
                {0.01712424535389e+00, 0.24365470201083e+00, 0.73922105263528e+00, 0.3478719460274719e-02 / 2.0}, // 6-perm
                {0.02288340534658e+00, 0.16538958561453e+00, 0.81172700903888e+00, 0.4102399036723953e-02 / 2.0}, // 6-perm
                {0.03273759728777e+00, 0.09930187449585e+00, 0.86796052821639e+00, 0.4021761549744162e-02 / 2.0}, // 6-perm
                {0.03382101234234e+00, 0.30847833306905e+00, 0.65770065458860e+00, 0.6033164660795066e-02 / 2.0}, // 6-perm
                {0.03554761446002e+00, 0.46066831859211e+00, 0.50378406694787e+00, 0.3946290302129598e-02 / 2.0}, // 6-perm
                {0.05053979030687e+00, 0.21881529945393e+00, 0.73064491023920e+00, 0.6644044537680268e-02 / 2.0}, // 6-perm
                {0.05701471491573e+00, 0.37920955156027e+00, 0.56377573352399e+00, 0.8254305856078458e-02 / 2.0}, // 6-perm
                {0.06415280642120e+00, 0.14296081941819e+00, 0.79288637416061e+00, 0.6496056633406411e-02 / 2.0}, // 6-perm
                {0.08050114828763e+00, 0.28373128210592e+00, 0.63576756960645e+00, 0.9252778144146602e-02 / 2.0}, // 6-perm
                {0.10436706813453e+00, 0.19673744100444e+00, 0.69889549086103e+00, 0.9164920726294280e-02 / 2.0}, // 6-perm
                {0.11384489442875e+00, 0.35588914121166e+00, 0.53026596435959e+00, 0.1156952462809767e-01 / 2.0}, // 6-perm
                {0.14536348771552e+00, 0.25981868535191e+00, 0.59481782693256e+00, 0.1176111646760917e-01 / 2.0}, // 6-perm
                {0.18994565282198e+00, 0.32192318123130e+00, 0.48813116594672e+00, 0.1382470218216540e-01 / 2.0}  // 6-perm
              };
 
 
              // Now call the dunavant routine to generate _points and _weights
              dunavant_rule(rule_data, 36);
 
              return;
            }
 
 
            // By default, we fall back on the conical product rules.  If the user
            // requests an order higher than what is currently available in the 1D
            // rules, an error will be thrown from the respective 1D code.
          default:
            {
              // The following quadrature rules are generated as
              // conical products.  These tend to be non-optimal
              // (use too many points, cluster points in certain
              // regions of the domain) but they are quite easy to
              // automatically generate using a 1D Gauss rule on
              // [0,1] and two 1D Jacobi-Gauss rules on [0,1].
              QConical conical_rule(2, _order);
              conical_rule.init(*this);
 
              // Swap points and weights with the about-to-be destroyed rule.
              _points.swap (conical_rule.get_points() );
              _weights.swap(conical_rule.get_weights());
 
              return;
            }
          }
      }
 
 
      //---------------------------------------------
      // Arbitrary polygon quadrature rules
    case C0POLYGON:
      {
        QGauss tri_rule(2, _order);
        tri_rule.init(TRI3, _p_level, true);
 
        std::vector<Point> & tripoints = tri_rule.get_points();
        std::vector<Real> & triweights = tri_rule.get_weights();
 
        std::size_t numtripts = tripoints.size();
 
        // C0Polygon requires the newer Quadrature API
        if (!_elem)
          libmesh_error();
 
        libmesh_assert(_elem->type() == C0POLYGON);
 
        const C0Polygon & poly = *cast_ptr<const C0Polygon *>(_elem);
 
        std::size_t numtris = poly.n_subtriangles();
        _points.resize(numtripts*numtris);
        _weights.resize(numtripts*numtris);
        for (std::size_t t = 0; t != numtris; ++t)
          {
            auto master_points = poly.master_subtriangle(t);
 
            // The factor of one half from the triweights cancels out
            // the factor of two here, so we don't need to do so
            // ourselves.
            const Real twice_master_tri_area =
              (- master_points[1](1) * master_points[2](0)
               - master_points[0](1) * master_points[1](0)
               + master_points[0](1) * master_points[2](0)
               + master_points[0](0) * master_points[1](1)
               - master_points[0](0) * master_points[2](1)
               + master_points[1](0) * master_points[2](1));
 
            const Point v01 = master_points[1] - master_points[0];
            const Point v02 = master_points[2] - master_points[0];
 
            for (std::size_t i = 0; i != numtripts; ++i)
              {
                _points[numtripts*t+i] =
                  master_points[0] +
                  v01 * tripoints[i](0) +
                  v02 * tripoints[i](1);
                _weights[numtripts*t+i] = triweights[i] *
                                          twice_master_tri_area;
              }
          }
        return;
      }
 
      //---------------------------------------------
      // Unsupported type
    default:
      libmesh_error_msg("Element type not supported:" << Utility::enum_to_string(_type));
    }
#endif
}

References libMesh::QBase::_elem, libMesh::QBase::_order, libMesh::QBase::_p_level, libMesh::QBase::_points, libMesh::QBase::_type, libMesh::QBase::_weights, b, libMesh::C0POLYGON, libMesh::CONSTANT, dunavant_rule(), dunavant_rule2(), libMesh::EIGHTH, libMesh::EIGHTTEENTH, libMesh::ELEVENTH, libMesh::Utility::enum_to_string(), libMesh::FIFTEENTH, libMesh::FIFTH, libMesh::FIRST, libMesh::FOURTEENTH, libMesh::FOURTH, libMesh::QBase::get_order(), libMesh::QBase::get_points(), libMesh::QBase::get_weights(), libMesh::QBase::init(), libMesh::libmesh_assert(), libMesh::Polygon::master_subtriangle(), libMesh::Polygon::n_subtriangles(), libMesh::NINETEENTH, 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, libMesh::QBase::tensor_product_quad(), libMesh::TENTH, libMesh::THIRD, libMesh::THIRTEENTH, libMesh::THIRTIETH, libMesh::TRI3, libMesh::TRI3SUBDIVISION, libMesh::TRI6, libMesh::TRI7, libMesh::TRISHELL3, libMesh::TWELFTH, libMesh::TWENTIETH, libMesh::TWENTYEIGHTH, libMesh::TWENTYFIFTH, libMesh::TWENTYFOURTH, libMesh::TWENTYNINTH, libMesh::TWENTYSECOND, libMesh::TWENTYSEVENTH, libMesh::TWENTYSIXTH, libMesh::TWENTYTHIRD, and libMesh::Elem::type().

init_3D()

void libMesh::QGauss::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 30 of file quadrature_gauss_3D.C.

{
#if LIBMESH_DIM == 3
 
  //-----------------------------------------------------------------------
  // 3D quadrature rules
  switch (_type)
    {
      //---------------------------------------------
      // Hex quadrature rules
    case HEX8:
    case HEX20:
    case HEX27:
      {
        // We compute the 3D quadrature rule as a tensor
        // product of the 1D quadrature rule.
        QGauss q1D(1, get_order());
        tensor_product_hex( q1D );
        return;
      }
 
 
 
      //---------------------------------------------
      // Tetrahedral quadrature rules
    case TET4:
    case TET10:
    case TET14:
      {
        switch(get_order())
          {
            // Taken from pg. 222 of "The finite element method," vol. 1
            // ed. 5 by Zienkiewicz & Taylor
          case CONSTANT:
          case FIRST:
            {
              // Exact for linears
              _points.resize(1);
              _weights.resize(1);
 
 
              _points[0](0) = .25;
              _points[0](1) = .25;
              _points[0](2) = .25;
 
              _weights[0] = Real(1)/6;
 
              return;
            }
          case SECOND:
            {
              // Exact for quadratics
              _points.resize(4);
              _weights.resize(4);
 
 
              // Can't be constexpr with my version of Boost quad
              // precision
              const Real b = 0.25*(1-std::sqrt(Real(5))/5);
              const Real a = 1-3*b;
 
              _points[0](0) = a;
              _points[0](1) = b;
              _points[0](2) = b;
 
              _points[1](0) = b;
              _points[1](1) = a;
              _points[1](2) = b;
 
              _points[2](0) = b;
              _points[2](1) = b;
              _points[2](2) = a;
 
              _points[3](0) = b;
              _points[3](1) = b;
              _points[3](2) = b;
 
 
 
              _weights[0] = Real(1)/24;
              _weights[1] = _weights[0];
              _weights[2] = _weights[0];
              _weights[3] = _weights[0];
 
              return;
            }
 
 
 
            // Can be found in the class notes
            // http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps
            // by Flaherty.
            //
            // Caution: this rule has a negative weight and may be
            // unsuitable for some problems.
            // Exact for cubics.
            //
            // Note: Keast (see ref. elsewhere in this file) also gives
            // a third-order rule with positive weights, but it contains points
            // on the ref. elt. boundary, making it less suitable for FEM calculations.
          case THIRD:
            {
              if (allow_rules_with_negative_weights)
                {
                  _points.resize(5);
                  _weights.resize(5);
 
 
                  _points[0](0) = .25;
                  _points[0](1) = .25;
                  _points[0](2) = .25;
 
                  _points[1](0) = .5;
                  _points[1](1) = Real(1)/6;
                  _points[1](2) = Real(1)/6;
 
                  _points[2](0) = Real(1)/6;
                  _points[2](1) = .5;
                  _points[2](2) = Real(1)/6;
 
                  _points[3](0) = Real(1)/6;
                  _points[3](1) = Real(1)/6;
                  _points[3](2) = .5;
 
                  _points[4](0) = Real(1)/6;
                  _points[4](1) = Real(1)/6;
                  _points[4](2) = Real(1)/6;
 
 
                  _weights[0] = Real(-2)/15;
                  _weights[1] = .075;
                  _weights[2] = _weights[1];
                  _weights[3] = _weights[1];
                  _weights[4] = _weights[1];
 
                  return;
                } // end if (allow_rules_with_negative_weights)
              else
                {
                  // If a rule with positive weights is required, the 2x2x2 conical
                  // product rule is third-order accurate and has less points than
                  // the next-available positive-weight rule at FIFTH order.
                  QConical conical_rule(3, _order);
                  conical_rule.init(*this);
 
                  // Swap points and weights with the about-to-be destroyed rule.
                  _points.swap (conical_rule.get_points() );
                  _weights.swap(conical_rule.get_weights());
 
                  return;
                }
              // Note: if !allow_rules_with_negative_weights, fall through to next case.
            }
 
 
 
            // Originally a Keast rule,
            //    Patrick Keast,
            //    Moderate Degree Tetrahedral Quadrature Formulas,
            //    Computer Methods in Applied Mechanics and Engineering,
            //    Volume 55, Number 3, May 1986, pages 339-348.
            //
            // Can also be found the class notes
            // http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps
            // by Flaherty.
            //
            // Caution: this rule has a negative weight and may be
            // unsuitable for some problems.
          case FOURTH:
            {
              if (allow_rules_with_negative_weights)
                {
                  _points.resize(11);
                  _weights.resize(11);
 
                  // The raw data for the quadrature rule.
                  const Real rule_data[3][4] = {
                    {0.250000000000000000e+00_R,                         0._R,                            0._R,  -0.131555555555555556e-01_R},  // 1
                    {0.785714285714285714e+00_R,   0.714285714285714285e-01_R,                            0._R,   0.762222222222222222e-02_R},  // 4
                    {0.399403576166799219e+00_R,                         0._R,      0.100596423833200785e+00_R,   0.248888888888888889e-01_R}   // 6
                  };
 
 
                  // Now call the keast routine to generate _points and _weights
                  keast_rule(rule_data, 3);
 
                  return;
                } // end if (allow_rules_with_negative_weights)
              // Note: if !allow_rules_with_negative_weights, fall through to next case.
            }
 
            libmesh_fallthrough();
 
 
            // Walkington's fifth-order 14-point rule from
            // "Quadrature on Simplices of Arbitrary Dimension"
            //
            // We originally had a Keast rule here, but this rule had
            // more points than an equivalent rule by Walkington and
            // also contained points on the boundary of the ref. elt,
            // making it less suitable for FEM calculations.
          case FIFTH:
            {
              _points.resize(14);
              _weights.resize(14);
 
              // permutations of these points and suitably-modified versions of
              // these points are the quadrature point locations
              const Real a[3] = {0.31088591926330060980_R,    // a1 from the paper
                                 0.092735250310891226402_R,   // a2 from the paper
                                 0.045503704125649649492_R};  // a3 from the paper
 
              // weights.  a[] and wt[] are the only floating-point inputs required
              // for this rule.
              const Real wt[3] = {0.018781320953002641800_R,    // w1 from the paper
                                  0.012248840519393658257_R,    // w2 from the paper
                                  0.0070910034628469110730_R};  // w3 from the paper
 
              // The first two sets of 4 points are formed in a similar manner
              for (unsigned int i=0; i<2; ++i)
                {
                  // Where we will insert values into _points and _weights
                  const unsigned int offset=4*i;
 
                  // Stuff points and weights values into their arrays
                  const Real b = 1. - 3.*a[i];
 
                  // Here are the permutations.  Order of these is not important,
                  // all have the same weight
                  _points[offset + 0] = Point(a[i], a[i], a[i]);
                  _points[offset + 1] = Point(a[i],    b, a[i]);
                  _points[offset + 2] = Point(   b, a[i], a[i]);
                  _points[offset + 3] = Point(a[i], a[i],    b);
 
                  // These 4 points all have the same weights
                  for (unsigned int j=0; j<4; ++j)
                    _weights[offset + j] = wt[i];
                } // end for
 
 
              {
                // The third set contains 6 points and is formed a little differently
                const unsigned int offset = 8;
                const Real b = 0.5*(1. - 2.*a[2]);
 
                // Here are the permutations.  Order of these is not important,
                // all have the same weight
                _points[offset + 0] = Point(b   ,    b, a[2]);
                _points[offset + 1] = Point(b   , a[2], a[2]);
                _points[offset + 2] = Point(a[2], a[2],    b);
                _points[offset + 3] = Point(a[2],    b, a[2]);
                _points[offset + 4] = Point(   b, a[2],    b);
                _points[offset + 5] = Point(a[2],    b,    b);
 
                // These 6 points all have the same weights
                for (unsigned int j=0; j<6; ++j)
                  _weights[offset + j] = wt[2];
              }
 
 
              // Original rule by Keast, unsuitable because it has points on the
              // reference element boundary!
              //       _points.resize(15);
              //       _weights.resize(15);
 
              //       _points[0](0) = 0.25;
              //       _points[0](1) = 0.25;
              //       _points[0](2) = 0.25;
 
              //       {
              // const Real a = 0.;
              // const Real b = Real(1)/3;
 
              // _points[1](0) = a;
              // _points[1](1) = b;
              // _points[1](2) = b;
 
              // _points[2](0) = b;
              // _points[2](1) = a;
              // _points[2](2) = b;
 
              // _points[3](0) = b;
              // _points[3](1) = b;
              // _points[3](2) = a;
 
              // _points[4](0) = b;
              // _points[4](1) = b;
              // _points[4](2) = b;
              //       }
              //       {
              // const Real a = Real(8)/11;
              // const Real b = Real(1)/11;
 
              // _points[5](0) = a;
              // _points[5](1) = b;
              // _points[5](2) = b;
 
              // _points[6](0) = b;
              // _points[6](1) = a;
              // _points[6](2) = b;
 
              // _points[7](0) = b;
              // _points[7](1) = b;
              // _points[7](2) = a;
 
              // _points[8](0) = b;
              // _points[8](1) = b;
              // _points[8](2) = b;
              //       }
              //       {
              // const Real a = 0.066550153573664;
              // const Real b = 0.433449846426336;
 
              // _points[9](0) = b;
              // _points[9](1) = a;
              // _points[9](2) = a;
 
              // _points[10](0) = a;
              // _points[10](1) = a;
              // _points[10](2) = b;
 
              // _points[11](0) = a;
              // _points[11](1) = b;
              // _points[11](2) = b;
 
              // _points[12](0) = b;
              // _points[12](1) = a;
              // _points[12](2) = b;
 
              // _points[13](0) = b;
              // _points[13](1) = b;
              // _points[13](2) = a;
 
              // _points[14](0) = a;
              // _points[14](1) = b;
              // _points[14](2) = a;
              //       }
 
              //       _weights[0]  = 0.030283678097089;
              //       _weights[1]  = 0.006026785714286;
              //       _weights[2]  = _weights[1];
              //       _weights[3]  = _weights[1];
              //       _weights[4]  = _weights[1];
              //       _weights[5]  = 0.011645249086029;
              //       _weights[6]  = _weights[5];
              //       _weights[7]  = _weights[5];
              //       _weights[8]  = _weights[5];
              //       _weights[9]  = 0.010949141561386;
              //       _weights[10] = _weights[9];
              //       _weights[11] = _weights[9];
              //       _weights[12] = _weights[9];
              //       _weights[13] = _weights[9];
              //       _weights[14] = _weights[9];
 
              return;
            }
 
 
 
 
            // This rule is originally from Keast:
            //    Patrick Keast,
            //    Moderate Degree Tetrahedral Quadrature Formulas,
            //    Computer Methods in Applied Mechanics and Engineering,
            //    Volume 55, Number 3, May 1986, pages 339-348.
            //
            // It is accurate on 6th-degree polynomials and has 24 points
            // vs. 64 for the comparable conical product rule.
            //
            // Values copied 24th June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/keast/keast.f90
          case SIXTH:
            {
              _points.resize (24);
              _weights.resize(24);
 
              // The raw data for the quadrature rule.
              const Real rule_data[4][4] = {
                {0.356191386222544953e+00_R , 0.214602871259151684e+00_R ,                       0._R, 0.00665379170969464506e+00_R}, // 4
                {0.877978124396165982e+00_R , 0.0406739585346113397e+00_R,                       0._R, 0.00167953517588677620e+00_R}, // 4
                {0.0329863295731730594e+00_R, 0.322337890142275646e+00_R ,                       0._R, 0.00922619692394239843e+00_R}, // 4
                {0.0636610018750175299e+00_R, 0.269672331458315867e+00_R , 0.603005664791649076e+00_R, 0.00803571428571428248e+00_R}  // 12
              };
 
 
              // Now call the keast routine to generate _points and _weights
              keast_rule(rule_data, 4);
 
              return;
            }
 
 
            // Keast's 31 point, 7th-order rule contains points on the reference
            // element boundary, so we've decided not to include it here.
            //
            // Keast's 8th-order rule has 45 points.  and a negative
            // weight, so if you've explicitly disallowed such rules
            // you will fall through to the conical product rules
            // below.
          case SEVENTH:
          case EIGHTH:
            {
              if (allow_rules_with_negative_weights)
                {
                  _points.resize (45);
                  _weights.resize(45);
 
                  // The raw data for the quadrature rule.
                  const Real rule_data[7][4] = {
                    {0.250000000000000000e+00_R,                        0._R,                        0._R,   -0.393270066412926145e-01_R},  // 1
                    {0.617587190300082967e+00_R,  0.127470936566639015e+00_R,                        0._R,    0.408131605934270525e-02_R},  // 4
                    {0.903763508822103123e+00_R,  0.320788303926322960e-01_R,                        0._R,    0.658086773304341943e-03_R},  // 4
                    {0.497770956432810185e-01_R,                        0._R,  0.450222904356718978e+00_R,    0.438425882512284693e-02_R},  // 6
                    {0.183730447398549945e+00_R,                        0._R,  0.316269552601450060e+00_R,    0.138300638425098166e-01_R},  // 6
                    {0.231901089397150906e+00_R,  0.229177878448171174e-01_R,  0.513280033360881072e+00_R,    0.424043742468372453e-02_R},  // 12
                    {0.379700484718286102e-01_R,  0.730313427807538396e+00_R,  0.193746475248804382e+00_R,    0.223873973961420164e-02_R}   // 12
                  };
 
 
                  // Now call the keast routine to generate _points and _weights
                  keast_rule(rule_data, 7);
 
                  return;
                } // end if (allow_rules_with_negative_weights)
              // Note: if !allow_rules_with_negative_weights, fall through to next case.
            }
 
            libmesh_fallthrough();
 
 
            // Fall back on Grundmann-Moller or Conical Product rules at high orders.
          default:
            {
              if ((allow_rules_with_negative_weights) && (get_order() < 34))
                {
                  // The Grundmann-Moller rules are defined to arbitrary order and
                  // can have significantly fewer evaluation points than conical product
                  // rules.  If you allow rules with negative weights, the GM rules
                  // will be more efficient for degree up to 33 (for degree 34 and
                  // higher, CP is more efficient!) but may be more susceptible
                  // to round-off error.  Safest is to disallow rules with negative
                  // weights, but this decision should be made on a case-by-case basis.
                  QGrundmann_Moller gm_rule(3, _order);
                  gm_rule.init(*this);
 
                  // Swap points and weights with the about-to-be destroyed rule.
                  _points.swap (gm_rule.get_points() );
                  _weights.swap(gm_rule.get_weights());
 
                  return;
                }
 
              else
                {
                  // The following quadrature rules are generated as
                  // conical products.  These tend to be non-optimal
                  // (use too many points, cluster points in certain
                  // regions of the domain) but they are quite easy to
                  // automatically generate using a 1D Gauss rule on
                  // [0,1] and two 1D Jacobi-Gauss rules on [0,1].
                  QConical conical_rule(3, _order);
                  conical_rule.init(*this);
 
                  // Swap points and weights with the about-to-be destroyed rule.
                  _points.swap (conical_rule.get_points() );
                  _weights.swap(conical_rule.get_weights());
 
                  return;
                }
            }
          }
      } // end case TET
 
 
 
      //---------------------------------------------
      // Prism quadrature rules
    case PRISM6:
    case PRISM15:
    case PRISM18:
    case PRISM20:
    case PRISM21:
      {
        // We compute the 3D quadrature rule as a tensor
        // product of the 1D quadrature rule and a 2D
        // triangle quadrature rule
 
        QGauss q1D(1,get_order());
        QGauss q2D(2,_order);
 
        // Initialize the 2D rule (1D is pre-initialized)
        q2D.init(TRI3, _p_level, /*simple_type_only=*/true);
 
        tensor_product_prism(q1D, q2D);
 
        return;
      }
 
 
 
      //---------------------------------------------
      // Pyramid
    case PYRAMID5:
    case PYRAMID13:
    case PYRAMID14:
    case PYRAMID18:
      {
        // We compute the Pyramid rule as a conical product of a
        // Jacobi rule with alpha==2 on the interval [0,1] two 1D
        // Gauss rules with suitably modified points.  The idea comes
        // from: Stroud, A.H. "Approximate Calculation of Multiple
        // Integrals."
        QConical conical_rule(3, _order);
        conical_rule.init(*this);
 
        // Swap points and weights with the about-to-be destroyed rule.
        _points.swap (conical_rule.get_points() );
        _weights.swap(conical_rule.get_weights());
 
        return;
 
      }
 
 
      //---------------------------------------------
      // Arbitrary polyhedron quadrature rules
    case C0POLYHEDRON:
      {
        QGauss tet_rule(3, _order);
        tet_rule.init(TET4, _p_level, true);
 
        std::vector<Point> & tetpoints = tet_rule.get_points();
        std::vector<Real> & tetweights = tet_rule.get_weights();
 
        std::size_t numtetpts = tetpoints.size();
 
        // C0Polyhedron requires the newer Quadrature API
        if (!_elem)
          libmesh_error();
 
        libmesh_assert(_elem->type() == C0POLYHEDRON);
 
        const C0Polyhedron & poly = *cast_ptr<const C0Polyhedron *>(_elem);
 
        std::size_t numtets = poly.n_subelements();
        _points.resize(numtetpts*numtets);
        _weights.resize(numtetpts*numtets);
 
        for (std::size_t t = 0; t != numtets; ++t)
          {
            auto master_points = poly.master_subelement(t);
 
            const Point v01 = master_points[1] - master_points[0];
            const Point v02 = master_points[2] - master_points[0];
            const Point v03 = master_points[3] - master_points[0];
 
            // The factor of one sixth from the tetweights cancels out
            // the factor of six here, so we don't need to do so
            // ourselves.
            const Real six_master_tet_vol =
              triple_product(v01, v02, v03);
 
            for (std::size_t i = 0; i != numtetpts; ++i)
              {
                _points[numtetpts*t+i] =
                  master_points[0] +
                    v01 * tetpoints[i](0) +
                    v02 * tetpoints[i](1) +
                    v03 * tetpoints[i](2);
                _weights[numtetpts*t+i] = tetweights[i] *
                                          six_master_tet_vol;
              }
          }
        return;
      }
 
      //---------------------------------------------
      // Unsupported type
    default:
      libmesh_error_msg("ERROR: Unsupported type: " << Utility::enum_to_string(_type));
    }
#endif
}

References libMesh::QBase::_elem, libMesh::QBase::_order, libMesh::QBase::_p_level, libMesh::QBase::_points, libMesh::QBase::_type, libMesh::QBase::_weights, libMesh::QBase::allow_rules_with_negative_weights, b, libMesh::C0POLYHEDRON, libMesh::CONSTANT, libMesh::EIGHTH, libMesh::Utility::enum_to_string(), libMesh::FIFTH, libMesh::FIRST, libMesh::FOURTH, libMesh::QBase::get_order(), libMesh::QBase::get_points(), libMesh::QBase::get_weights(), libMesh::HEX20, libMesh::HEX27, libMesh::HEX8, libMesh::QBase::init(), keast_rule(), libMesh::libmesh_assert(), libMesh::Polyhedron::master_subelement(), libMesh::Polyhedron::n_subelements(), libMesh::PRISM15, libMesh::PRISM18, libMesh::PRISM20, libMesh::PRISM21, libMesh::PRISM6, libMesh::PYRAMID13, libMesh::PYRAMID14, libMesh::PYRAMID18, libMesh::PYRAMID5, libMesh::Real, libMesh::SECOND, libMesh::SEVENTH, libMesh::SIXTH, libMesh::QBase::tensor_product_hex(), libMesh::QBase::tensor_product_prism(), libMesh::TET10, libMesh::TET14, libMesh::TET4, libMesh::THIRD, libMesh::TRI3, libMesh::triple_product(), and libMesh::Elem::type().

keast_rule()

void libMesh::QGauss::keast_rule ( const Real  rule_data[][4], const unsigned int  n_pts  )

private

The Keast rules are for tets.

This function takes permutation points and weights in a specific format as input and fills the _points and _weights vectors.

Definition at line 43 of file quadrature_gauss.C.

{
  // Like the Dunavant rule, the input data should have 4 columns.  These columns
  // have the following format and implied permutations (w=weight).
  // {a, 0, 0, w} = 1-permutation  (a,a,a)
  // {a, b, 0, w} = 4-permutation  (a,b,b), (b,a,b), (b,b,a), (b,b,b)
  // {a, 0, b, w} = 6-permutation  (a,a,b), (a,b,b), (b,b,a), (b,a,b), (b,a,a), (a,b,a)
  // {a, b, c, w} = 12-permutation (a,a,b), (a,a,c), (b,a,a), (c,a,a), (a,b,a), (a,c,a)
  //                               (a,b,c), (a,c,b), (b,a,c), (b,c,a), (c,a,b), (c,b,a)
 
  // Always insert into the points & weights vector relative to the offset
  unsigned int offset=0;
 
 
  for (unsigned int p=0; p<n_pts; ++p)
    {
 
      // There must always be a non-zero entry to start the row
      libmesh_assert_not_equal_to (rule_data[p][0], static_cast<Real>(0.0));
 
      // A zero weight may imply you did not set up the raw data correctly
      libmesh_assert_not_equal_to (rule_data[p][3], static_cast<Real>(0.0));
 
      // What kind of point is this?
      // One non-zero entry in first 3 cols   ? 1-perm (centroid) point = 1
      // Two non-zero entries in first 3 cols ? 3-perm point            = 3
      // Three non-zero entries               ? 6-perm point            = 6
      unsigned int pointtype=1;
 
      if (rule_data[p][1] != static_cast<Real>(0.0))
        {
          if (rule_data[p][2] != static_cast<Real>(0.0))
            pointtype = 12;
          else
            pointtype = 4;
        }
      else
        {
          // The second entry is zero.  What about the third?
          if (rule_data[p][2] != static_cast<Real>(0.0))
            pointtype = 6;
        }
 
 
      switch (pointtype)
        {
        case 1:
          {
            // Be sure we have enough space to insert this point
            libmesh_assert_less (offset + 0, _points.size());
 
            const Real a = rule_data[p][0];
 
            // The point has only a single permutation (the centroid!)
            _points[offset  + 0] = Point(a,a,a);
 
            // The weight is always the last entry in the row.
            _weights[offset + 0] = rule_data[p][3];
 
            offset += pointtype;
            break;
          }
 
        case 4:
          {
            // Be sure we have enough space to insert these points
            libmesh_assert_less (offset + 3, _points.size());
 
            const Real a  = rule_data[p][0];
            const Real b  = rule_data[p][1];
            const Real wt = rule_data[p][3];
 
            // Here it's understood the second entry is to be used twice, and
            // thus there are three possible permutations.
            _points[offset + 0] = Point(a,b,b);
            _points[offset + 1] = Point(b,a,b);
            _points[offset + 2] = Point(b,b,a);
            _points[offset + 3] = Point(b,b,b);
 
            for (unsigned int j=0; j<pointtype; ++j)
              _weights[offset + j] = wt;
 
            offset += pointtype;
            break;
          }
 
        case 6:
          {
            // Be sure we have enough space to insert these points
            libmesh_assert_less (offset + 5, _points.size());
 
            const Real a  = rule_data[p][0];
            const Real b  = rule_data[p][2];
            const Real wt = rule_data[p][3];
 
            // Three individual entries with six permutations.
            _points[offset + 0] = Point(a,a,b);
            _points[offset + 1] = Point(a,b,b);
            _points[offset + 2] = Point(b,b,a);
            _points[offset + 3] = Point(b,a,b);
            _points[offset + 4] = Point(b,a,a);
            _points[offset + 5] = Point(a,b,a);
 
            for (unsigned int j=0; j<pointtype; ++j)
              _weights[offset + j] = wt;
 
            offset += pointtype;
            break;
          }
 
 
        case 12:
          {
            // Be sure we have enough space to insert these points
            libmesh_assert_less (offset + 11, _points.size());
 
            const Real a  = rule_data[p][0];
            const Real b  = rule_data[p][1];
            const Real c  = rule_data[p][2];
            const Real wt = rule_data[p][3];
 
            // Three individual entries with six permutations.
            _points[offset + 0] = Point(a,a,b);  _points[offset + 6]  = Point(a,b,c);
            _points[offset + 1] = Point(a,a,c); _points[offset + 7]  = Point(a,c,b);
            _points[offset + 2] = Point(b,a,a); _points[offset + 8]  = Point(b,a,c);
            _points[offset + 3] = Point(c,a,a); _points[offset + 9]  = Point(b,c,a);
            _points[offset + 4] = Point(a,b,a); _points[offset + 10] = Point(c,a,b);
            _points[offset + 5] = Point(a,c,a); _points[offset + 11] = Point(c,b,a);
 
            for (unsigned int j=0; j<pointtype; ++j)
              _weights[offset + j] = wt;
 
            offset += pointtype;
            break;
          }
 
        default:
          libmesh_error_msg("Don't know what to do with this many permutation points!");
        }
 
    }
 
}

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

Referenced by init_3D().

n_objects()

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

inlinestaticinherited

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

Definition at line 85 of file reference_counter.h.

  { return _n_objects; }

References libMesh::ReferenceCounter::_n_objects.

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

n_points()

unsigned int libMesh::QBase::n_points ( ) const

inlineinherited

Returns

The number of points associated with the quadrature rule.

Definition at line 131 of file quadrature.h.

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

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

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

operator=() [1/2]

QGauss & libMesh::QGauss::operator= ( const QGauss &  )

default

operator=() [2/2]

QGauss & libMesh::QGauss::operator= ( QGauss &&  )

default

print_info() [1/2]

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

inherited

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

Definition at line 43 of file quadrature.C.

{
  libmesh_assert(!_points.empty());
  libmesh_assert(!_weights.empty());
 
  Real summed_weights=0;
  os << "N_Q_Points=" << this->n_points() << std::endl << std::endl;
  for (auto qpoint: index_range(_points))
    {
      os << " Point " << qpoint << ":\n"
         << "  "
         << _points[qpoint]
         << "\n Weight:\n "
         << "  w=" << _weights[qpoint] << "\n" << std::endl;
 
      summed_weights += _weights[qpoint];
    }
  os << "Summed Weights: " << summed_weights << std::endl;
}

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

print_info() [2/2]

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

staticinherited

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

Definition at line 81 of file reference_counter.C.

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

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

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

qp()

Point libMesh::QBase::qp ( const unsigned int  i ) const

inlineinherited

Returns

The \( i^{th} \) quadrature point in reference element space.

Definition at line 179 of file quadrature.h.

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

References libMesh::QBase::_points.

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

scale()

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

inherited

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

Definition at line 222 of file quadrature.C.

{
  // Make sure we are in 1D
  libmesh_assert_equal_to (_dim, 1);
 
  Real
    h_new = new_range.second - new_range.first,
    h_old = old_range.second - old_range.first;
 
  // Make sure that we have sane ranges
  libmesh_assert_greater (h_new, 0.);
  libmesh_assert_greater (h_old, 0.);
 
  // Make sure there are some points
  libmesh_assert_greater (_points.size(), 0);
 
  // Compute the scale factor
  Real scfact = h_new/h_old;
 
  // We're mapping from old_range -> new_range
  for (auto i : index_range(_points))
    {
      _points[i](0) = new_range.first +
        (_points[i](0) - old_range.first) * scfact;
 
      // Scale the weights
      _weights[i] *= scfact;
    }
}

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

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

shapes_need_reinit()

virtual bool libMesh::QBase::shapes_need_reinit ( )

inlinevirtualinherited

Returns

true if the shape functions need to be recalculated, false otherwise.

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

Definition at line 286 of file quadrature.h.

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

  {
    return n_points();
  }

References libMesh::QBase::n_points().

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.

{
  const unsigned int np = q1D.n_points();
 
  _points.resize(np * np * np);
 
  _weights.resize(np * np * np);
 
  unsigned int q=0;
 
  for (unsigned int k=0; k<np; k++)
    for (unsigned int j=0; j<np; j++)
      for (unsigned int i=0; i<np; i++)
        {
          _points[q](0) = q1D.qp(i)(0);
          _points[q](1) = q1D.qp(j)(0);
          _points[q](2) = q1D.qp(k)(0);
 
          _weights[q] = q1D.w(i) * q1D.w(j) * q1D.w(k);
 
          q++;
        }
}

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

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

tensor_product_prism()

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

protectedinherited

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

Used in the init_3D routines for prismatic element types.

Definition at line 310 of file quadrature.C.

{
  const unsigned int n_points1D = q1D.n_points();
  const unsigned int n_points2D = q2D.n_points();
 
  _points.resize  (n_points1D * n_points2D);
  _weights.resize (n_points1D * n_points2D);
 
  unsigned int q=0;
 
  for (unsigned int j=0; j<n_points1D; j++)
    for (unsigned int i=0; i<n_points2D; i++)
      {
        _points[q](0) = q2D.qp(i)(0);
        _points[q](1) = q2D.qp(i)(1);
        _points[q](2) = q1D.qp(j)(0);
 
        _weights[q] = q2D.w(i) * q1D.w(j);
 
        q++;
      }
 
}

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

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

tensor_product_quad()

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

protectedinherited

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

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

Definition at line 256 of file quadrature.C.

{
 
  const unsigned int np = q1D.n_points();
 
  _points.resize(np * np);
 
  _weights.resize(np * np);
 
  unsigned int q=0;
 
  for (unsigned int j=0; j<np; j++)
    for (unsigned int i=0; i<np; i++)
      {
        _points[q](0) = q1D.qp(i)(0);
        _points[q](1) = q1D.qp(j)(0);
 
        _weights[q] = q1D.w(i)*q1D.w(j);
 
        q++;
      }
}

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

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

type()

QuadratureType libMesh::QGauss::type ( ) const

overridevirtual

Returns

QGAUSS.

Implements libMesh::QBase.

Definition at line 33 of file quadrature_gauss.C.

{
  return QGAUSS;
}

References libMesh::QGAUSS.

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.

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

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

Member Data Documentation

_counts

ReferenceCounter::Counts libMesh::ReferenceCounter::_counts

staticprotectedinherited

Actually holds the data.

Definition at line 124 of file reference_counter.h.

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

_dim

unsigned int libMesh::QBase::_dim

protectedinherited

The spatial dimension of the quadrature rule.

Definition at line 379 of file quadrature.h.

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

_elem

const Elem* libMesh::QBase::_elem

protectedinherited

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

Definition at line 397 of file quadrature.h.

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

_enable_print_counter

bool libMesh::ReferenceCounter::_enable_print_counter = true

staticprotectedinherited

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

Definition at line 143 of file reference_counter.h.

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

_mutex

Threads::spin_mutex libMesh::ReferenceCounter::_mutex

staticprotectedinherited

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 137 of file reference_counter.h.

_n_objects

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

staticprotectedinherited

The number of objects.

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

Definition at line 132 of file reference_counter.h.

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

_order

Order libMesh::QBase::_order

protectedinherited

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

Definition at line 385 of file quadrature.h.

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

_p_level

unsigned int libMesh::QBase::_p_level

protectedinherited

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

Definition at line 403 of file quadrature.h.

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

_points

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

protectedinherited

The locations of the quadrature points in reference element space.

Definition at line 409 of file quadrature.h.

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

_type

ElemType libMesh::QBase::_type

protectedinherited

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

Definition at line 391 of file quadrature.h.

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

_weights

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

protectedinherited

The quadrature weights.

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

Definition at line 415 of file quadrature.h.

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

allow_nodal_pyramid_quadrature

bool libMesh::QBase::allow_nodal_pyramid_quadrature

inherited

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

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

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

Definition at line 314 of file quadrature.h.

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

allow_rules_with_negative_weights

bool libMesh::QBase::allow_rules_with_negative_weights

inherited

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

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

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

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

Definition at line 301 of file quadrature.h.

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

---

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

• include/quadrature/quadrature_gauss.h
• src/quadrature/quadrature_gauss.C
• src/quadrature/quadrature_gauss_1D.C
• src/quadrature/quadrature_gauss_2D.C
• src/quadrature/quadrature_gauss_3D.C