#include <quadrature_monomial.h>

Inheritance diagram for libMesh::QMonomial:

Public Member Functions
	QMonomial (const unsigned int _dim, const Order _order=INVALID_ORDER)
	~QMonomial ()
QuadratureType	type () const
ElemType	get_elem_type () const
unsigned int	get_p_level () const
unsigned int	n_points () const
unsigned int	get_dim () const
const std::vector< Point > &	get_points () const
std::vector< Point > &	get_points ()
const std::vector< Real > &	get_weights () const
std::vector< Real > &	get_weights ()
Point	qp (const unsigned int i) const
Real	w (const unsigned int i) const
virtual void	init (const ElemType type=INVALID_ELEM, unsigned int p_level=0)
virtual void	init (const Elem &elem, const std::vector< Real > &vertex_distance_func, unsigned int p_level=0)
Order	get_order () const
void	print_info (std::ostream &os=libMesh::out) const
void	scale (std::pair< Real, Real > old_range, std::pair< Real, Real > new_range)
virtual bool	shapes_need_reinit ()
Static Public Member Functions
static UniquePtr< QBase >	build (const std::string &name, const unsigned int _dim, const Order _order=INVALID_ORDER)
static UniquePtr< QBase >	build (const QuadratureType _qt, const unsigned int _dim, const Order _order=INVALID_ORDER)
static void	print_info (std::ostream &out=libMesh::out)
static std::string	get_info ()
static unsigned int	n_objects ()
static void	enable_print_counter_info ()
static void	disable_print_counter_info ()
Public Attributes
bool	allow_rules_with_negative_weights
Protected Types
typedef std::map< std::string, std::pair< unsigned int, unsigned int > >	Counts
Protected Member Functions
virtual void	init_0D (const ElemType type=INVALID_ELEM, unsigned int p_level=0)
	libmesh_error_msg ("ERROR: Seems as if this quadrature rule \nis not implemented for 2D.")
void	tensor_product_hex (const QBase &q1D)
void	tensor_product_prism (const QBase &q1D, const QBase &q2D)
void	increment_constructor_count (const std::string &name)
void	increment_destructor_count (const std::string &name)
Protected Attributes
const unsigned int	_dim
const Order	_order
ElemType	_type
unsigned int	_p_level
std::vector< Point >	_points
std::vector< Real >	_weights
Static Protected Attributes
static Counts	_counts
static Threads::atomic < unsigned int >	_n_objects
static Threads::spin_mutex	_mutex
static bool	_enable_print_counter = true
Private Member Functions
void	init_1D (const ElemType, unsigned int=0)
void	init_2D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0)
void	init_3D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0)
void	wissmann_rule (const Real rule_data[][3], const unsigned int n_pts)
void	stroud_rule (const Real rule_data[][3], const unsigned int *rule_symmetry, const unsigned int n_pts)
void	kim_rule (const Real rule_data[][4], const unsigned int *rule_id, const unsigned int n_pts)
Friends
std::ostream &	operator<< (std::ostream &os, const QBase &q)

Detailed Description

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

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

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

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

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

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

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

Author:: John W. Peterson, 2008

Definition at line 65 of file quadrature_monomial.h.

Member Typedef Documentation

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > libMesh::ReferenceCounter::Counts [protected, inherited]

Data structure to log the information. The log is identified by the class name.

Definition at line 113 of file reference_counter.h.

Constructor & Destructor Documentation

libMesh::QMonomial::QMonomial	(	const unsigned int	_dim,
		const Order	_order = `INVALID_ORDER`
	)

Constructor. Declares the order of the quadrature rule.

Definition at line 33 of file quadrature_monomial.C.

                                    : QBase(d,o)
{
}

libMesh::QMonomial::~QMonomial ( )

Destructor.

Definition at line 40 of file quadrature_monomial.C.

{
}

Member Function Documentation

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

Builds a specific quadrature rule, identified through the name string. An UniquePtr<QBase> is returned to prevent a memory leak. This way the user need not remember to delete the object. Enables run-time decision of the quadrature rule. The input parameter name must be mappable through the Utility::string_to_enum<>() function.

Definition at line 42 of file quadrature_build.C.

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

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

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

Builds a specific quadrature rule, identified through the QuadratureType. An UniquePtr<QBase> is returned to prevent a memory leak. This way the user need not remember to delete the object. Enables run-time decision of the quadrature rule.

Definition at line 53 of file quadrature_build.C.

References libMesh::FIRST, libMesh::FORTYTHIRD, libMesh::QBase::libmesh_error_msg(), 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::QSIMPSON, libMesh::QTRAP, libMesh::THIRD, and libMesh::TWENTYTHIRD.

{
  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 UniquePtr<QBase>(new 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 UniquePtr<QBase>(new 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 UniquePtr<QBase>(new 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 UniquePtr<QBase>(new 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 UniquePtr<QBase>(new QSimpson(_dim));
      }

    case QTRAP:
      {

#ifdef DEBUG
        if (_order > FIRST)
          {
            libMesh::out << "WARNING: Trapezoidal rule provides only" << std::endl
                         << " FIRST order!" << std::endl;
          }
#endif

        return UniquePtr<QBase>(new QTrap(_dim));
      }

    case QGRID:
      return UniquePtr<QBase>(new QGrid(_dim, _order));

    case QGRUNDMANN_MOLLER:
      return UniquePtr<QBase>(new QGrundmann_Moller(_dim, _order));

    case QMONOMIAL:
      return UniquePtr<QBase>(new QMonomial(_dim, _order));

    case QGAUSS_LOBATTO:
      return UniquePtr<QBase>(new QGaussLobatto(_dim, _order));

    case QCONICAL:
      return UniquePtr<QBase>(new QConical(_dim, _order));

    default:
      libmesh_error_msg("ERROR: Bad qt=" << _qt);
    }


  libmesh_error_msg("We'll never get here!");
  return UniquePtr<QBase>();
}

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

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

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

{
  _enable_print_counter = false;
  return;
}

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

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

Definition at line 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

{
  _enable_print_counter = true;
  return;
}

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

Returns:: the dimension of the quadrature rule.

Definition at line 125 of file quadrature.h.

References libMesh::QBase::_dim.

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

{ return _dim;  }

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

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

Definition at line 106 of file quadrature.h.

References libMesh::QBase::_type.

  { return _type; }

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

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

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

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

{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)

  std::ostringstream oss;

  oss << '\n'
      << " ---------------------------------------------------------------------------- \n"
      << "| Reference count information                                                |\n"
      << " ---------------------------------------------------------------------------- \n";

  for (Counts::iterator it = _counts.begin();
       it != _counts.end(); ++it)
    {
      const std::string name(it->first);
      const unsigned int creations    = it->second.first;
      const unsigned int destructions = it->second.second;

      oss << "| " << name << " reference count information:\n"
          << "|  Creations:    " << creations    << '\n'
          << "|  Destructions: " << destructions << '\n';
    }

  oss << " ---------------------------------------------------------------------------- \n";

  return oss.str();

#else

  return "";

#endif
}

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

Returns:: the order of the quadrature rule.

Definition at line 183 of file quadrature.h.

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

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

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

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

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

Definition at line 112 of file quadrature.h.

References libMesh::QBase::_p_level.

  { return _p_level; }

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

Returns:: a std::vector containing the quadrature point locations on a reference object.

Definition at line 131 of file quadrature.h.

References libMesh::QBase::_points.

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

{ return _points;  }

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

Returns:: a std::vector containing the quadrature point locations on a reference object as a writeable reference.

Definition at line 137 of file quadrature.h.

References libMesh::QBase::_points.

{ return _points;  }

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

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

Definition at line 142 of file quadrature.h.

References libMesh::QBase::_weights.

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

{ return _weights; }

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

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

Definition at line 147 of file quadrature.h.

References libMesh::QBase::_weights.

{ return _weights; }

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

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

Definition at line 163 of file reference_counter.h.

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

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

{
  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
  std::pair<unsigned int, unsigned int>& p = _counts[name];

  p.first++;
}

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

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

Definition at line 176 of file reference_counter.h.

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

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

{
  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
  std::pair<unsigned int, unsigned int>& p = _counts[name];

  p.second++;
}

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

Initializes the data structures to contain a quadrature rule for an object of type type.

Definition at line 28 of file quadrature.C.

References libMesh::QBase::_dim, libMesh::QBase::_p_level, libMesh::QBase::_type, libMesh::QBase::init_0D(), libMesh::QBase::init_1D(), libMesh::QBase::init_2D(), libMesh::QBase::init_3D(), and libMesh::QBase::libmesh_error_msg().

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

{
  // check to see if we have already
  // done the work for this quadrature rule
  if (t == _type && p == _p_level)
    return;
  else
    {
      _type = t;
      _p_level = p;
    }



  switch(_dim)
    {
    case 0:
      this->init_0D(_type,_p_level);

      return;

    case 1:
      this->init_1D(_type,_p_level);

      return;

    case 2:
      this->init_2D(_type,_p_level);

      return;

    case 3:
      this->init_3D(_type,_p_level);

      return;

    default:
      libmesh_error_msg("Invalid dimension _dim = " << _dim);
    }
}

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

Initializes the data structures for a specific, potentially cut element. The array vertex_distance_func contains vertex values of a signed distance function that cuts the element. This interface is indended to be extended by derived classes that can cut the element into subelements, for example, and constuct a composite quadrature rule for the cut element.

Definition at line 72 of file quadrature.C.

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

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

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

Initializes the 0D quadrature rule by filling the points and weights vectors with the appropriate values. Generally this is just one point with weight 1.

Definition at line 82 of file quadrature.C.

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

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

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

void libMesh::QMonomial::init_1D	(	const ElemType	_elemtype,
		unsigned int	p = `0`
	)		`[private, virtual]`

Just uses a Gauss rule in 1D.

Implements libMesh::QBase.

Definition at line 29 of file quadrature_monomial_1D.C.

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

{
  QGauss gauss_rule(1, _order);
  gauss_rule.init(_elemtype, p);

  _points.swap(gauss_rule.get_points());
  _weights.swap(gauss_rule.get_weights());
}

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

More efficient rules for QUADs

Reimplemented from libMesh::QBase.

Definition at line 28 of file quadrature_monomial_2D.C.

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

{

  switch (type_in)
    {
      //---------------------------------------------
      // Quadrilateral quadrature rules
    case QUAD4:
    case QUAD8:
    case QUAD9:
      {
        switch(_order + 2*p)
          {
          case SECOND:
            {
              // A degree=2 rule for the QUAD with 3 points.
              // A tensor product degree-2 Gauss would have 4 points.
              // This rule (or a variation on it) is probably available in
              //
              // A.H. Stroud, Approximate calculation of multiple integrals,
              // Prentice-Hall, Englewood Cliffs, N.J., 1971.
              //
              // though I have never actually seen a reference for it.
              // Luckily it's fairly easy to derive, which is what I've done
              // here [JWP].
              const Real
                s=std::sqrt(1./3.),
                t=std::sqrt(2./3.);

              const Real data[2][3] =
                {
                  {0.0,  s,  2.0},
                  {  t, -s,  1.0}
                };

              _points.resize(3);
              _weights.resize(3);

              wissmann_rule(data, 2);

              return;
            } // end case SECOND



            // For third-order, fall through to default case, use 2x2 Gauss product rule.
            // case THIRD:
            //   {
            //   }  // end case THIRD

          case FOURTH:
            {
              // A pair of degree=4 rules for the QUAD "C2" due to
              // Wissmann and Becker. These rules both have six points.
              // A tensor product degree-4 Gauss would have 9 points.
              //
              // J. W. Wissmann and T. Becker, Partially symmetric cubature
              // formulas for even degrees of exactness, SIAM J. Numer. Anal.  23
              // (1986), 676--685.
              const Real data[4][3] =
                {
                  // First of 2 degree-4 rules given by Wissmann
                  {0.0000000000000000e+00,  0.0000000000000000e+00,  1.1428571428571428e+00},
                  {0.0000000000000000e+00,  9.6609178307929590e-01,  4.3956043956043956e-01},
                  {8.5191465330460049e-01,  4.5560372783619284e-01,  5.6607220700753210e-01},
                  {6.3091278897675402e-01, -7.3162995157313452e-01,  6.4271900178367668e-01}
                  //
                  // Second of 2 degree-4 rules given by Wissmann.  These both
                  // yield 4th-order accurate rules, I just chose the one that
                  // happened to contain the origin.
                  // {0.000000000000000, -0.356822089773090,  1.286412084888852},
                  // {0.000000000000000,  0.934172358962716,  0.491365692888926},
                  // {0.774596669241483,  0.390885162530071,  0.761883709085613},
                  // {0.774596669241483, -0.852765377881771,  0.349227402025498}
                };

              _points.resize(6);
              _weights.resize(6);

              wissmann_rule(data, 4);

              return;
            } // end case FOURTH




          case FIFTH:
            {
              // A degree 5, 7-point rule due to Stroud.
              //
              // A.H. Stroud, Approximate calculation of multiple integrals,
              // Prentice-Hall, Englewood Cliffs, N.J., 1971.
              //
              // This rule is provably minimal in the number of points.
              // A tensor-product rule accurate for "bi-quintic" polynomials would have 9 points.
              const Real data[3][3] =
                {
                  {                                  0.L,                                     0.L, static_cast<Real>(8.L  /  7.L)}, // 1
                  {                                  0.L, static_cast<Real>(std::sqrt(14.L/15.L)), static_cast<Real>(20.L / 63.L)}, // 2
                  {static_cast<Real>(std::sqrt(3.L/5.L)),   static_cast<Real>(std::sqrt(1.L/3.L)), static_cast<Real>(20.L / 36.L)}  // 4
                };

              const unsigned int symmetry[3] = {
                0, // Origin
                7, // Central Symmetry
                6  // Rectangular
              };

              _points.resize (7);
              _weights.resize(7);

              stroud_rule(data, symmetry, 3);

              return;
            } // end case FIFTH




          case SIXTH:
            {
              // A pair of degree=6 rules for the QUAD "C2" due to
              // Wissmann and Becker. These rules both have 10 points.
              // A tensor product degree-6 Gauss would have 16 points.
              //
              // J. W. Wissmann and T. Becker, Partially symmetric cubature
              // formulas for even degrees of exactness, SIAM J. Numer. Anal.  23
              // (1986), 676--685.
              const Real data[6][3] =
                {
                  // First of 2 degree-6, 10 point rules given by Wissmann
                  // {0.000000000000000,  0.836405633697626,  0.455343245714174},
                  // {0.000000000000000, -0.357460165391307,  0.827395973202966},
                  // {0.888764014654765,  0.872101531193131,  0.144000884599645},
                  // {0.604857639464685,  0.305985162155427,  0.668259104262665},
                  // {0.955447506641064, -0.410270899466658,  0.225474004890679},
                  // {0.565459993438754, -0.872869311156879,  0.320896396788441}
                  //
                  // Second of 2 degree-6, 10 point rules given by Wissmann.
                  // Either of these will work, I just chose the one with points
                  // slightly further into the element interior.
                  {0.0000000000000000e+00,  8.6983337525005900e-01,  3.9275059096434794e-01},
                  {0.0000000000000000e+00, -4.7940635161211124e-01,  7.5476288124261053e-01},
                  {8.6374282634615388e-01,  8.0283751620765670e-01,  2.0616605058827902e-01},
                  {5.1869052139258234e-01,  2.6214366550805818e-01,  6.8999213848986375e-01},
                  {9.3397254497284950e-01, -3.6309658314806653e-01,  2.6051748873231697e-01},
                  {6.0897753601635630e-01, -8.9660863276245265e-01,  2.6956758608606100e-01}
                };

              _points.resize(10);
              _weights.resize(10);

              wissmann_rule(data, 6);

              return;
            } // end case SIXTH




          case SEVENTH:
            {
              // A degree 7, 12-point rule due to Tyler, can be found in Stroud's book
              //
              // A.H. Stroud, Approximate calculation of multiple integrals,
              // Prentice-Hall, Englewood Cliffs, N.J., 1971.
              //
              // This rule is fully-symmetric and provably minimal in the number of points.
              // A tensor-product rule accurate for "bi-septic" polynomials would have 16 points.
              const Real
                r  = std::sqrt(6.L/7.L),
                s  = std::sqrt( (114.L - 3.L*std::sqrt(583.L)) / 287.L ),
                t  = std::sqrt( (114.L + 3.L*std::sqrt(583.L)) / 287.L ),
                B1 = 196.L / 810.L,
                B2 = 4.L * (178981.L + 2769.L*std::sqrt(583.L)) / 1888920.L,
                B3 = 4.L * (178981.L - 2769.L*std::sqrt(583.L)) / 1888920.L;

              const Real data[3][3] =
                {
                  {r, 0.0, B1}, // 4
                  {s, 0.0, B2}, // 4
                  {t, 0.0, B3}  // 4
                };

              const unsigned int symmetry[3] = {
                3, // Full Symmetry, (x,0)
                2, // Full Symmetry, (x,x)
                2  // Full Symmetry, (x,x)
              };

              _points.resize (12);
              _weights.resize(12);

              stroud_rule(data, symmetry, 3);

              return;
            } // end case SEVENTH




          case EIGHTH:
            {
              // A pair of degree=8 rules for the QUAD "C2" due to
              // Wissmann and Becker. These rules both have 16 points.
              // A tensor product degree-6 Gauss would have 25 points.
              //
              // J. W. Wissmann and T. Becker, Partially symmetric cubature
              // formulas for even degrees of exactness, SIAM J. Numer. Anal.  23
              // (1986), 676--685.
              const Real data[10][3] =
                {
                  // First of 2 degree-8, 16 point rules given by Wissmann
                  // {0.000000000000000,  0.000000000000000,  0.055364705621440},
                  // {0.000000000000000,  0.757629177660505,  0.404389368726076},
                  // {0.000000000000000, -0.236871842255702,  0.533546604952635},
                  // {0.000000000000000, -0.989717929044527,  0.117054188786739},
                  // {0.639091304900370,  0.950520955645667,  0.125614417613747},
                  // {0.937069076924990,  0.663882736885633,  0.136544584733588},
                  // {0.537083530541494,  0.304210681724104,  0.483408479211257},
                  // {0.887188506449625, -0.236496718536120,  0.252528506429544},
                  // {0.494698820670197, -0.698953476086564,  0.361262323882172},
                  // {0.897495818279768, -0.900390774211580,  0.085464254086247}
                  //
                  // Second of 2 degree-8, 16 point rules given by Wissmann.
                  // Either of these will work, I just chose the one with points
                  // further into the element interior.
                  {0.0000000000000000e+00,  6.5956013196034176e-01,  4.5027677630559029e-01},
                  {0.0000000000000000e+00, -9.4914292304312538e-01,  1.6657042677781274e-01},
                  {9.5250946607156228e-01,  7.6505181955768362e-01,  9.8869459933431422e-02},
                  {5.3232745407420624e-01,  9.3697598108841598e-01,  1.5369674714081197e-01},
                  {6.8473629795173504e-01,  3.3365671773574759e-01,  3.9668697607290278e-01},
                  {2.3314324080140552e-01, -7.9583272377396852e-02,  3.5201436794569501e-01},
                  {9.2768331930611748e-01, -2.7224008061253425e-01,  1.8958905457779799e-01},
                  {4.5312068740374942e-01, -6.1373535339802760e-01,  3.7510100114758727e-01},
                  {8.3750364042281223e-01, -8.8847765053597136e-01,  1.2561879164007201e-01}
                };

              _points.resize(16);
              _weights.resize(16);

              wissmann_rule(data, /*10*/ 9);

              return;
            } // end case EIGHTH




          case NINTH:
            {
              // A degree 9, 17-point rule due to Moller.
              //
              // H.M. Moller,  Kubaturformeln mit minimaler Knotenzahl,
              // Numer. Math.  25 (1976), 185--200.
              //
              // This rule is provably minimal in the number of points.
              // A tensor-product rule accurate for "bi-ninth" degree polynomials would have 25 points.
              const Real data[5][3] =
                {
                  {0.0000000000000000e+00, 0.0000000000000000e+00, 5.2674897119341563e-01}, // 1
                  {6.3068011973166885e-01, 9.6884996636197772e-01, 8.8879378170198706e-02}, // 4
                  {9.2796164595956966e-01, 7.5027709997890053e-01, 1.1209960212959648e-01}, // 4
                  {4.5333982113564719e-01, 5.2373582021442933e-01, 3.9828243926207009e-01}, // 4
                  {8.5261572933366230e-01, 7.6208328192617173e-02, 2.6905133763978080e-01}  // 4
                };

              const unsigned int symmetry[5] = {
                0, // Single point
                4, // Rotational Invariant
                4, // Rotational Invariant
                4, // Rotational Invariant
                4  // Rotational Invariant
              };

              _points.resize (17);
              _weights.resize(17);

              stroud_rule(data, symmetry, 5);

              return;
            } // end case NINTH




          case TENTH:
          case ELEVENTH:
            {
              // A degree 11, 24-point rule due to Cools and Haegemans.
              //
              // R. Cools and A. Haegemans, Another step forward in searching for
              // cubature formulae with a minimal number of knots for the square,
              // Computing 40 (1988), 139--146.
              //
              // P. Verlinden and R. Cools, The algebraic construction of a minimal
              // cubature formula of degree 11 for the square, Cubature Formulas
              // and their Applications (Russian) (Krasnoyarsk) (M.V. Noskov, ed.),
              // 1994, pp. 13--23.
              //
              // This rule is provably minimal in the number of points.
              // A tensor-product rule accurate for "bi-tenth" or "bi-eleventh" degree polynomials would have 36 points.
              const Real data[6][3] =
                {
                  {6.9807610454956756e-01, 9.8263922354085547e-01, 4.8020763350723814e-02}, // 4
                  {9.3948638281673690e-01, 8.2577583590296393e-01, 6.6071329164550595e-02}, // 4
                  {9.5353952820153201e-01, 1.8858613871864195e-01, 9.7386777358668164e-02}, // 4
                  {3.1562343291525419e-01, 8.1252054830481310e-01, 2.1173634999894860e-01}, // 4
                  {7.1200191307533630e-01, 5.2532025036454776e-01, 2.2562606172886338e-01}, // 4
                  {4.2484724884866925e-01, 4.1658071912022368e-02, 3.5115871839824543e-01}  // 4
                };

              const unsigned int symmetry[6] = {
                4, // Rotational Invariant
                4, // Rotational Invariant
                4, // Rotational Invariant
                4, // Rotational Invariant
                4, // Rotational Invariant
                4  // Rotational Invariant
              };

              _points.resize (24);
              _weights.resize(24);

              stroud_rule(data, symmetry, 6);

              return;
            } // end case TENTH,ELEVENTH




          case TWELFTH:
          case THIRTEENTH:
            {
              // A degree 13, 33-point rule due to Cools and Haegemans.
              //
              // R. Cools and A. Haegemans, Another step forward in searching for
              // cubature formulae with a minimal number of knots for the square,
              // Computing 40 (1988), 139--146.
              //
              // A tensor-product rule accurate for "bi-12" or "bi-13" degree polynomials would have 49 points.
              const Real data[9][3] =
                {
                  {0.0000000000000000e+00, 0.0000000000000000e+00, 3.0038211543122536e-01}, // 1
                  {9.8348668243987226e-01, 7.7880971155441942e-01, 2.9991838864499131e-02}, // 4
                  {8.5955600564163892e-01, 9.5729769978630736e-01, 3.8174421317083669e-02}, // 4
                  {9.5892517028753485e-01, 1.3818345986246535e-01, 6.0424923817749980e-02}, // 4
                  {3.9073621612946100e-01, 9.4132722587292523e-01, 7.7492738533105339e-02}, // 4
                  {8.5007667369974857e-01, 4.7580862521827590e-01, 1.1884466730059560e-01}, // 4
                  {6.4782163718701073e-01, 7.5580535657208143e-01, 1.2976355037000271e-01}, // 4
                  {7.0741508996444936e-02, 6.9625007849174941e-01, 2.1334158145718938e-01}, // 4
                  {4.0930456169403884e-01, 3.4271655604040678e-01, 2.5687074948196783e-01}  // 4
                };

              const unsigned int symmetry[9] = {
                0, // Single point
                4, // Rotational Invariant
                4, // Rotational Invariant
                4, // Rotational Invariant
                4, // Rotational Invariant
                4, // Rotational Invariant
                4, // Rotational Invariant
                4, // Rotational Invariant
                4  // Rotational Invariant
              };

              _points.resize (33);
              _weights.resize(33);

              stroud_rule(data, symmetry, 9);

              return;
            } // end case TWELFTH,THIRTEENTH




          case FOURTEENTH:
          case FIFTEENTH:
            {
              // A degree-15, 48 point rule originally due to Rabinowitz and Richter,
              // can be found in Cools' 1971 book.
              //
              // A.H. Stroud, Approximate calculation of multiple integrals,
              // Prentice-Hall, Englewood Cliffs, N.J., 1971.
              //
              // The product Gauss rule for this order has 8^2=64 points.
              const Real data[9][3] =
                {
                  {9.915377816777667e-01L, 0.0000000000000000e+00,  3.01245207981210e-02L}, // 4
                  {8.020163879230440e-01L, 0.0000000000000000e+00,  8.71146840209092e-02L}, // 4
                  {5.648674875232742e-01L, 0.0000000000000000e+00, 1.250080294351494e-01L}, // 4
                  {9.354392392539896e-01L, 0.0000000000000000e+00,  2.67651407861666e-02L}, // 4
                  {7.624563338825799e-01L, 0.0000000000000000e+00,  9.59651863624437e-02L}, // 4
                  {2.156164241427213e-01L, 0.0000000000000000e+00, 1.750832998343375e-01L}, // 4
                  {9.769662659711761e-01L, 6.684480048977932e-01L,  2.83136372033274e-02L}, // 4
                  {8.937128379503403e-01L, 3.735205277617582e-01L,  8.66414716025093e-02L}, // 4
                  {6.122485619312083e-01L, 4.078983303613935e-01L, 1.150144605755996e-01L}  // 4
                };

              const unsigned int symmetry[9] = {
                3, // Full Symmetry, (x,0)
                3, // Full Symmetry, (x,0)
                3, // Full Symmetry, (x,0)
                2, // Full Symmetry, (x,x)
                2, // Full Symmetry, (x,x)
                2, // Full Symmetry, (x,x)
                1, // Full Symmetry, (x,y)
                1, // Full Symmetry, (x,y)
                1, // Full Symmetry, (x,y)
              };

              _points.resize (48);
              _weights.resize(48);

              stroud_rule(data, symmetry, 9);

              return;
            } //   case FOURTEENTH, FIFTEENTH:




          case SIXTEENTH:
          case SEVENTEENTH:
            {
              // A degree 17, 60-point rule due to Cools and Haegemans.
              //
              // R. Cools and A. Haegemans, Another step forward in searching for
              // cubature formulae with a minimal number of knots for the square,
              // Computing 40 (1988), 139--146.
              //
              // A tensor-product rule accurate for "bi-14" or "bi-15" degree polynomials would have 64 points.
              // A tensor-product rule accurate for "bi-16" or "bi-17" degree polynomials would have 81 points.
              const Real data[10][3] =
                {
                  {9.8935307451260049e-01, 0.0000000000000000e+00, 2.0614915919990959e-02}, // 4
                  {3.7628520715797329e-01, 0.0000000000000000e+00, 1.2802571617990983e-01}, // 4
                  {9.7884827926223311e-01, 0.0000000000000000e+00, 5.5117395340318905e-03}, // 4
                  {8.8579472916411612e-01, 0.0000000000000000e+00, 3.9207712457141880e-02}, // 4
                  {1.7175612383834817e-01, 0.0000000000000000e+00, 7.6396945079863302e-02}, // 4
                  {5.9049927380600241e-01, 3.1950503663457394e-01, 1.4151372994997245e-01}, // 8
                  {7.9907913191686325e-01, 5.9797245192945738e-01, 8.3903279363797602e-02}, // 8
                  {8.0374396295874471e-01, 5.8344481776550529e-02, 6.0394163649684546e-02}, // 8
                  {9.3650627612749478e-01, 3.4738631616620267e-01, 5.7387752969212695e-02}, // 8
                  {9.8132117980545229e-01, 7.0600028779864611e-01, 2.1922559481863763e-02}, // 8
                };

              const unsigned int symmetry[10] = {
                3, // Fully symmetric (x,0)
                3, // Fully symmetric (x,0)
                2, // Fully symmetric (x,x)
                2, // Fully symmetric (x,x)
                2, // Fully symmetric (x,x)
                1, // Fully symmetric (x,y)
                1, // Fully symmetric (x,y)
                1, // Fully symmetric (x,y)
                1, // Fully symmetric (x,y)
                1  // Fully symmetric (x,y)
              };

              _points.resize (60);
              _weights.resize(60);

              stroud_rule(data, symmetry, 10);

              return;
            } // end case FOURTEENTH through SEVENTEENTH



            // By default: construct and use a Gauss quadrature rule
          default:
            {
              // Break out and fall down into the default: case for the
              // outer switch statement.
              break;
            }

          } // end switch(_order + 2*p)
      } // end case QUAD4/8/9


      // By default: construct and use a Gauss quadrature rule
    default:
      {
        QGauss gauss_rule(2, _order);
        gauss_rule.init(type_in, p);

        // Swap points and weights with the about-to-be destroyed rule.
        _points.swap (gauss_rule.get_points() );
        _weights.swap(gauss_rule.get_weights());

        return;
      }
    } // end switch (type_in)
}

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

More efficient rules for HEXes

Reimplemented from libMesh::QBase.

Definition at line 28 of file quadrature_monomial_3D.C.

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

{

  switch (type_in)
    {
      //---------------------------------------------
      // Hex quadrature rules
    case HEX8:
    case HEX20:
    case HEX27:
      {
        switch(_order + 2*p)
          {

            // The CONSTANT/FIRST rule is the 1-point Gauss "product" rule...we fall
            // through to the default case for this rule.

          case SECOND:
          case THIRD:
            {
              // A degree 3, 6-point, "rotationally-symmetric" rule by
              // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
              //
              // Warning: this rule contains points on the boundary of the reference
              // element, and therefore may be unsuitable for some problems.  The alternative
              // would be a 2x2x2 Gauss product rule.
              const Real data[1][4] =
                {
                  {1.0L, 0.0L, 0.0L, static_cast<Real>(4.0L/3.0L)}
                };

              const unsigned int rule_id[1] = {
                1 // (x,0,0) -> 6 permutations
              };

              _points.resize(6);
              _weights.resize(6);

              kim_rule(data, rule_id, 1);
              return;
            } // end case SECOND,THIRD

          case FOURTH:
          case FIFTH:
            {
              // A degree 5, 13-point rule by Stroud,
              // AH Stroud, "Some Fifth Degree Integration Formulas for Symmetric Regions II.",
              // Numerische Mathematik 9, pp. 460-468 (1967).
              //
              // This rule is provably minimal in the number of points.  The equations given for
              // the n-cube on pg. 466 of the paper for mu/gamma and gamma are wrong, at least for
              // the n=3 case.  The analytical values given here were computed by me [JWP] in Maple.

              // Convenient intermediate values.
              const Real sqrt19 = std::sqrt(19.L);
              const Real tp     = std::sqrt(71440.L + 6802.L*sqrt19);

              // Point data for permutations.
              const Real eta    =  0.00000000000000000000000000000000e+00L;

              const Real lambda =  std::sqrt(1919.L/3285.L - 148.L*sqrt19/3285.L + 4.L*tp/3285.L);
              // 8.8030440669930978047737818209860e-01L;

              const Real xi     = -std::sqrt(1121.L/3285.L +  74.L*sqrt19/3285.L - 2.L*tp/3285.L);
              // -4.9584817142571115281421242364290e-01L;

              const Real mu     =  std::sqrt(1121.L/3285.L +  74.L*sqrt19/3285.L + 2.L*tp/3285.L);
              // 7.9562142216409541542982482567580e-01L;

              const Real gamma  =  std::sqrt(1919.L/3285.L - 148.L*sqrt19/3285.L - 4.L*tp/3285.L);
              // 2.5293711744842581347389255929324e-02L;

              // Weights: the centroid weight is given analytically.  Weight B (resp C) goes
              // with the {lambda,xi} (resp {gamma,mu}) permutation.  The single-precision
              // results reported by Stroud are given for reference.

              const Real A      = 32.0L / 19.0L;
              // Stroud: 0.21052632  * 8.0 = 1.684210560;

              const Real B      = 1.L / ( 260072.L/133225.L  - 1520*sqrt19/133225.L + (133.L - 37.L*sqrt19)*tp/133225.L );
              // 5.4498735127757671684690782180890e-01L; // Stroud: 0.068123420 * 8.0 = 0.544987360;

              const Real C      = 1.L / ( 260072.L/133225.L  - 1520*sqrt19/133225.L - (133.L - 37.L*sqrt19)*tp/133225.L );
              // 5.0764422766979170420572375713840e-01L; // Stroud: 0.063455527 * 8.0 = 0.507644216;

              _points.resize(13);
              _weights.resize(13);

              unsigned int c=0;

              // Point with weight A (origin)
              _points[c] = Point(eta, eta, eta);
              _weights[c++] = A;

              // Points with weight B
              _points[c] = Point(lambda, xi, xi);
              _weights[c++] = B;
              _points[c] = -_points[c-1];
              _weights[c++] = B;

              _points[c] = Point(xi, lambda, xi);
              _weights[c++] = B;
              _points[c] = -_points[c-1];
              _weights[c++] = B;

              _points[c] = Point(xi, xi, lambda);
              _weights[c++] = B;
              _points[c] = -_points[c-1];
              _weights[c++] = B;

              // Points with weight C
              _points[c] = Point(mu, mu, gamma);
              _weights[c++] = C;
              _points[c] = -_points[c-1];
              _weights[c++] = C;

              _points[c] = Point(mu, gamma, mu);
              _weights[c++] = C;
              _points[c] = -_points[c-1];
              _weights[c++] = C;

              _points[c] = Point(gamma, mu, mu);
              _weights[c++] = C;
              _points[c] = -_points[c-1];
              _weights[c++] = C;

              return;


              //       // A degree 5, 14-point, "rotationally-symmetric" rule by
              //       // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
              //       // Was also reported in Stroud's 1971 book.
              //       const Real data[2][4] =
              // {
              //   {7.95822425754221463264548820476135e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 8.86426592797783933518005540166204e-01L},
              //   {7.58786910639328146269034278112267e-01L, 7.58786910639328146269034278112267e-01L, 7.58786910639328146269034278112267e-01L, 3.35180055401662049861495844875346e-01L}
              // };

              //       const unsigned int rule_id[2] = {
              // 1, // (x,0,0) -> 6 permutations
              // 4  // (x,x,x) -> 8 permutations
              //       };

              //       _points.resize(14);
              //       _weights.resize(14);

              //       kim_rule(data, rule_id, 2);
              //       return;
            } // end case FOURTH,FIFTH

          case SIXTH:
          case SEVENTH:
            {
              if (allow_rules_with_negative_weights)
                {
                  // A degree 7, 31-point, "rotationally-symmetric" rule by
                  // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
                  // This rule contains a negative weight, so only use it if such type of
                  // rules are allowed.
                  const Real data[3][4] =
                    {
                      {0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, -1.27536231884057971014492753623188e+00L},
                      {5.85540043769119907612630781744060e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L,  8.71111111111111111111111111111111e-01L},
                      {6.94470135991704766602025803883310e-01L, 9.37161638568208038511047377665396e-01L, 4.15659267604065126239606672567031e-01L,  1.68695652173913043478260869565217e-01L}
                    };

                  const unsigned int rule_id[3] = {
                    0, // (0,0,0) -> 1 permutation
                    1, // (x,0,0) -> 6 permutations
                    6  // (x,y,z) -> 24 permutations
                  };

                  _points.resize(31);
                  _weights.resize(31);

                  kim_rule(data, rule_id, 3);
                  return;
                } // end if (allow_rules_with_negative_weights)


              // A degree 7, 34-point, "fully-symmetric" rule, first published in
              // P.C. Hammer and A.H. Stroud, "Numerical Evaluation of Multiple Integrals II",
              // Mathmatical Tables and Other Aids to Computation, vol 12., no 64, 1958, pp. 272-280
              //
              // This rule happens to fall under the same general
              // construction as the Kim rules, so we've re-used
              // that code here.  Stroud gives 16 digits for his rule,
              // and this is the most accurate version I've found.
              //
              // For comparison, a SEVENTH-order Gauss product rule
              // (which integrates tri-7th order polynomials) would
              // have 4^3=64 points.
              const Real
                r  = std::sqrt(6.L/7.L),
                s  = std::sqrt((960.L - 3.L*std::sqrt(28798.L)) / 2726.L),
                t  = std::sqrt((960.L + 3.L*std::sqrt(28798.L)) / 2726.L),
                B1 = 8624.L / 29160.L,
                B2 = 2744.L / 29160.L,
                B3 = 8.L*(774.L*t*t - 230.L)/(9720.L*(t*t-s*s)),
                B4 = 8.L*(230.L - 774.L*s*s)/(9720.L*(t*t-s*s));

              const Real data[4][4] =
                {
                  {r, 0.L, 0.L, B1},
                  {r,   r, 0.L, B2},
                  {s,   s,   s, B3},
                  {t,   t,   t, B4}
                };

              const unsigned int rule_id[4] = {
                1, // (x,0,0) -> 6 permutations
                2, // (x,x,0) -> 12 permutations
                4, // (x,x,x) -> 8 permutations
                4  // (x,x,x) -> 8 permutations
              };

              _points.resize(34);
              _weights.resize(34);

              kim_rule(data, rule_id, 4);
              return;


              //      // A degree 7, 38-point, "rotationally-symmetric" rule by
              //      // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
              //      //
              //      // This rule is obviously inferior to the 34-point rule above...
              //      const Real data[3][4] =
              //{
              //  {9.01687807821291289082811566285950e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 2.95189738262622903181631100062774e-01L},
              //  {4.08372221499474674069588900002128e-01L, 4.08372221499474674069588900002128e-01L, 4.08372221499474674069588900002128e-01L, 4.04055417266200582425904380777126e-01L},
              //  {8.59523090201054193116477875786220e-01L, 8.59523090201054193116477875786220e-01L, 4.14735913727987720499709244748633e-01L, 1.24850759678944080062624098058597e-01L}
              //};
              //
              //      const unsigned int rule_id[3] = {
              //1, // (x,0,0) -> 6 permutations
              //4, // (x,x,x) -> 8 permutations
              //5  // (x,x,z) -> 24 permutations
              //      };
              //
              //      _points.resize(38);
              //      _weights.resize(38);
              //
              //      kim_rule(data, rule_id, 3);
              //      return;
            } // end case SIXTH,SEVENTH

          case EIGHTH:
            {
              // A degree 8, 47-point, "rotationally-symmetric" rule by
              // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
              //
              // A EIGHTH-order Gauss product rule (which integrates tri-8th order polynomials)
              // would have 5^3=125 points.
              const Real data[5][4] =
                {
                  {0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 4.51903714875199690490763818699555e-01L},
                  {7.82460796435951590652813975429717e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 2.99379177352338919703385618576171e-01L},
                  {4.88094669706366480526729301468686e-01L, 4.88094669706366480526729301468686e-01L, 4.88094669706366480526729301468686e-01L, 3.00876159371240019939698689791164e-01L},
                  {8.62218927661481188856422891110042e-01L, 8.62218927661481188856422891110042e-01L, 8.62218927661481188856422891110042e-01L, 4.94843255877038125738173175714853e-02L},
                  {2.81113909408341856058098281846420e-01L, 9.44196578292008195318687494773744e-01L, 6.97574833707236996779391729948984e-01L, 1.22872389222467338799199767122592e-01L}
                };

              const unsigned int rule_id[5] = {
                0, // (0,0,0) -> 1 permutation
                1, // (x,0,0) -> 6 permutations
                4, // (x,x,x) -> 8 permutations
                4, // (x,x,x) -> 8 permutations
                6  // (x,y,z) -> 24 permutations
              };

              _points.resize(47);
              _weights.resize(47);

              kim_rule(data, rule_id, 5);
              return;
            } // end case EIGHTH


            // By default: construct and use a Gauss quadrature rule
          default:
            {
              // Break out and fall down into the default: case for the
              // outer switch statement.
              break;
            }

          } // end switch(_order + 2*p)
      } // end case HEX8/20/27


      // By default: construct and use a Gauss quadrature rule
    default:
      {
        QGauss gauss_rule(3, _order);
        gauss_rule.init(type_in, p);

        // Swap points and weights with the about-to-be destroyed rule.
        _points.swap (gauss_rule.get_points() );
        _weights.swap(gauss_rule.get_weights());

        return;
      }
    } // end switch (type_in)
}

void libMesh::QMonomial::kim_rule	(	const Real	rule_data[][4],
		const unsigned int *	rule_id,
		const unsigned int	n_pts
	)		`[private]`

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

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

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

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

Definition at line 212 of file quadrature_monomial.C.

References libMesh::QBase::_points, libMesh::QBase::_weights, libMesh::QBase::libmesh_error_msg(), libMesh::Real, and libMesh::x.

Referenced by init_3D().

{
  for (unsigned int i=0, c=0; i<n_pts; ++i)
    {
      const Real
        x=rule_data[i][0],
        y=rule_data[i][1],
        z=rule_data[i][2],
        wt=rule_data[i][3];

      switch(rule_id[i])
        {
        case 0: // (0,0,0) 1 permutation
          {
            _points[c]  = Point( x, y, z);    _weights[c++] = wt;

            break;
          }
        case 1: //  (x,0,0) 6 permutations
          {
            _points[c] = Point( x, 0., 0.);    _weights[c++] = wt;
            _points[c] = Point(0., -x, 0.);    _weights[c++] = wt;
            _points[c] = Point(-x, 0., 0.);    _weights[c++] = wt;
            _points[c] = Point(0.,  x, 0.);    _weights[c++] = wt;
            _points[c] = Point(0., 0., -x);    _weights[c++] = wt;
            _points[c] = Point(0., 0.,  x);    _weights[c++] = wt;

            break;
          }
        case 2: // (x,x,0) 12 permutations
          {
            _points[c] = Point( x,  x, 0.);    _weights[c++] = wt;
            _points[c] = Point( x, -x, 0.);    _weights[c++] = wt;
            _points[c] = Point(-x, -x, 0.);    _weights[c++] = wt;
            _points[c] = Point(-x,  x, 0.);    _weights[c++] = wt;
            _points[c] = Point( x, 0., -x);    _weights[c++] = wt;
            _points[c] = Point( x, 0.,  x);    _weights[c++] = wt;
            _points[c] = Point(0.,  x, -x);    _weights[c++] = wt;
            _points[c] = Point(0.,  x,  x);    _weights[c++] = wt;
            _points[c] = Point(0., -x, -x);    _weights[c++] = wt;
            _points[c] = Point(-x, 0., -x);    _weights[c++] = wt;
            _points[c] = Point(0., -x,  x);    _weights[c++] = wt;
            _points[c] = Point(-x, 0.,  x);    _weights[c++] = wt;

            break;
          }
        case 3: // (x,y,0) 24 permutations
          {
            _points[c] = Point( x,  y, 0.);    _weights[c++] = wt;
            _points[c] = Point( y, -x, 0.);    _weights[c++] = wt;
            _points[c] = Point(-x, -y, 0.);    _weights[c++] = wt;
            _points[c] = Point(-y,  x, 0.);    _weights[c++] = wt;
            _points[c] = Point( x, 0., -y);    _weights[c++] = wt;
            _points[c] = Point( x, -y, 0.);    _weights[c++] = wt;
            _points[c] = Point( x, 0.,  y);    _weights[c++] = wt;
            _points[c] = Point(0.,  y, -x);    _weights[c++] = wt;
            _points[c] = Point(-x,  y, 0.);    _weights[c++] = wt;
            _points[c] = Point(0.,  y,  x);    _weights[c++] = wt;
            _points[c] = Point( y, 0., -x);    _weights[c++] = wt;
            _points[c] = Point(0., -y, -x);    _weights[c++] = wt;
            _points[c] = Point(-y, 0., -x);    _weights[c++] = wt;
            _points[c] = Point( y,  x, 0.);    _weights[c++] = wt;
            _points[c] = Point(-y, -x, 0.);    _weights[c++] = wt;
            _points[c] = Point( y, 0.,  x);    _weights[c++] = wt;
            _points[c] = Point(0., -y,  x);    _weights[c++] = wt;
            _points[c] = Point(-y, 0.,  x);    _weights[c++] = wt;
            _points[c] = Point(-x, 0.,  y);    _weights[c++] = wt;
            _points[c] = Point(0., -x, -y);    _weights[c++] = wt;
            _points[c] = Point(0., -x,  y);    _weights[c++] = wt;
            _points[c] = Point(-x, 0., -y);    _weights[c++] = wt;
            _points[c] = Point(0.,  x,  y);    _weights[c++] = wt;
            _points[c] = Point(0.,  x, -y);    _weights[c++] = wt;

            break;
          }
        case 4: // (x,x,x) 8 permutations
          {
            _points[c] = Point( x,  x,  x);    _weights[c++] = wt;
            _points[c] = Point( x, -x,  x);    _weights[c++] = wt;
            _points[c] = Point(-x, -x,  x);    _weights[c++] = wt;
            _points[c] = Point(-x,  x,  x);    _weights[c++] = wt;
            _points[c] = Point( x,  x, -x);    _weights[c++] = wt;
            _points[c] = Point( x, -x, -x);    _weights[c++] = wt;
            _points[c] = Point(-x,  x, -x);    _weights[c++] = wt;
            _points[c] = Point(-x, -x, -x);    _weights[c++] = wt;

            break;
          }
        case 5: // (x,x,z) 24 permutations
          {
            _points[c] = Point( x,  x,  z);    _weights[c++] = wt;
            _points[c] = Point( x, -x,  z);    _weights[c++] = wt;
            _points[c] = Point(-x, -x,  z);    _weights[c++] = wt;
            _points[c] = Point(-x,  x,  z);    _weights[c++] = wt;
            _points[c] = Point( x,  z, -x);    _weights[c++] = wt;
            _points[c] = Point( x, -x, -z);    _weights[c++] = wt;
            _points[c] = Point( x, -z,  x);    _weights[c++] = wt;
            _points[c] = Point( z,  x, -x);    _weights[c++] = wt;
            _points[c] = Point(-x,  x, -z);    _weights[c++] = wt;
            _points[c] = Point(-z,  x,  x);    _weights[c++] = wt;
            _points[c] = Point( x, -z, -x);    _weights[c++] = wt;
            _points[c] = Point(-z, -x, -x);    _weights[c++] = wt;
            _points[c] = Point(-x,  z, -x);    _weights[c++] = wt;
            _points[c] = Point( x,  x, -z);    _weights[c++] = wt;
            _points[c] = Point(-x, -x, -z);    _weights[c++] = wt;
            _points[c] = Point( x,  z,  x);    _weights[c++] = wt;
            _points[c] = Point( z, -x,  x);    _weights[c++] = wt;
            _points[c] = Point(-x, -z,  x);    _weights[c++] = wt;
            _points[c] = Point(-x,  z,  x);    _weights[c++] = wt;
            _points[c] = Point( z, -x, -x);    _weights[c++] = wt;
            _points[c] = Point(-z, -x,  x);    _weights[c++] = wt;
            _points[c] = Point(-x, -z, -x);    _weights[c++] = wt;
            _points[c] = Point( z,  x,  x);    _weights[c++] = wt;
            _points[c] = Point(-z,  x, -x);    _weights[c++] = wt;

            break;
          }
        case 6: // (x,y,z) 24 permutations
          {
            _points[c] = Point( x,  y,  z);    _weights[c++] = wt;
            _points[c] = Point( y, -x,  z);    _weights[c++] = wt;
            _points[c] = Point(-x, -y,  z);    _weights[c++] = wt;
            _points[c] = Point(-y,  x,  z);    _weights[c++] = wt;
            _points[c] = Point( x,  z, -y);    _weights[c++] = wt;
            _points[c] = Point( x, -y, -z);    _weights[c++] = wt;
            _points[c] = Point( x, -z,  y);    _weights[c++] = wt;
            _points[c] = Point( z,  y, -x);    _weights[c++] = wt;
            _points[c] = Point(-x,  y, -z);    _weights[c++] = wt;
            _points[c] = Point(-z,  y,  x);    _weights[c++] = wt;
            _points[c] = Point( y, -z, -x);    _weights[c++] = wt;
            _points[c] = Point(-z, -y, -x);    _weights[c++] = wt;
            _points[c] = Point(-y,  z, -x);    _weights[c++] = wt;
            _points[c] = Point( y,  x, -z);    _weights[c++] = wt;
            _points[c] = Point(-y, -x, -z);    _weights[c++] = wt;
            _points[c] = Point( y,  z,  x);    _weights[c++] = wt;
            _points[c] = Point( z, -y,  x);    _weights[c++] = wt;
            _points[c] = Point(-y, -z,  x);    _weights[c++] = wt;
            _points[c] = Point(-x,  z,  y);    _weights[c++] = wt;
            _points[c] = Point( z, -x, -y);    _weights[c++] = wt;
            _points[c] = Point(-z, -x,  y);    _weights[c++] = wt;
            _points[c] = Point(-x, -z, -y);    _weights[c++] = wt;
            _points[c] = Point( z,  x,  y);    _weights[c++] = wt;
            _points[c] = Point(-z,  x, -y);    _weights[c++] = wt;

            break;
          }
        default:
          libmesh_error_msg("Unknown rule ID: " << rule_id[i] << "!");
        } // end switch(rule_id[i])
    }
}

libMesh::QBase::libmesh_error_msg ( "ERROR: Seems as if this quadrature rule \nis not implemented for 2D." ) [protected, inherited]

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

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

Definition at line 79 of file reference_counter.h.

References libMesh::ReferenceCounter::_n_objects.

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

  { return _n_objects; }

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

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

Definition at line 118 of file quadrature.h.

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

Referenced by libMesh::ExactSolution::_compute_error(), libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), libMesh::ProjectFEMSolution::operator()(), libMesh::QBase::print_info(), libMesh::QBase::tensor_product_hex(), and libMesh::QBase::tensor_product_prism().

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

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

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

Definition at line 88 of file reference_counter.C.

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

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

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

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

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

Definition at line 372 of file quadrature.h.

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

Referenced by libMesh::operator<<().

{
  libmesh_assert(!_points.empty());
  libmesh_assert(!_weights.empty());

  Real summed_weights=0;
  os << "N_Q_Points=" << this->n_points() << std::endl << std::endl;
  for (unsigned int qpoint=0; qpoint<this->n_points(); qpoint++)
    {
      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;
}

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

Returns:: the $i^{th}$ quadrature point on the reference object.

Definition at line 152 of file quadrature.h.

References libMesh::QBase::_points.

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

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

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

Maps the points of a 1D interval quadrature rule (typically [-1,1]) to any other 1D interval (typically [0,1]) and scales the weights accordingly. The quadrature rule will be mapped from the entries of old_range to the entries of new_range.

Definition at line 93 of file quadrature.C.

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

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

{
  // 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 (unsigned int i=0; i<_points.size(); i++)
    {
      _points[i](0) = new_range.first +
        (_points[i](0) - old_range.first) * scfact;

      // Scale the weights
      _weights[i] *= scfact;
    }
}

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

Returns true if the shape functions need to be recalculated.

This can happen if the number of points or their positions change.

By default this will return false.

Definition at line 212 of file quadrature.h.

{ return false; }

void libMesh::QMonomial::stroud_rule	(	const Real	rule_data[][3],
		const unsigned int *	rule_symmetry,
		const unsigned int	n_pts
	)		`[private]`

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

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

Definition at line 64 of file quadrature_monomial.C.

References libMesh::QBase::_points, libMesh::QBase::_weights, libMesh::QBase::libmesh_error_msg(), libMesh::Real, and libMesh::x.

Referenced by init_2D().

{
  for (unsigned int i=0, c=0; i<n_pts; ++i)
    {
      const Real
        x=rule_data[i][0],
        y=rule_data[i][1],
        wt=rule_data[i][2];

      switch(rule_symmetry[i])
        {
        case 0: // Single point (no symmetry)
          {
            _points[c]  = Point( x, y);
            _weights[c++] = wt;

            break;
          }
        case 1: // Fully-symmetric (x,y)
          {
            _points[c]    = Point( x, y);
            _weights[c++] = wt;

            _points[c]    = Point(-x, y);
            _weights[c++] = wt;

            _points[c]    = Point( x,-y);
            _weights[c++] = wt;

            _points[c]    = Point(-x,-y);
            _weights[c++] = wt;

            _points[c]    = Point( y, x);
            _weights[c++] = wt;

            _points[c]    = Point(-y, x);
            _weights[c++] = wt;

            _points[c]    = Point( y,-x);
            _weights[c++] = wt;

            _points[c]    = Point(-y,-x);
            _weights[c++] = wt;

            break;
          }
        case 2: // Fully-symmetric (x,x)
          {
            _points[c]    = Point( x, x);
            _weights[c++] = wt;

            _points[c]    = Point(-x, x);
            _weights[c++] = wt;

            _points[c]    = Point( x,-x);
            _weights[c++] = wt;

            _points[c]    = Point(-x,-x);
            _weights[c++] = wt;

            break;
          }
        case 3: // Fully-symmetric (x,0)
          {
            libmesh_assert_equal_to (y, 0.0);

            _points[c]    = Point( x,0.);
            _weights[c++] = wt;

            _points[c]    = Point(-x,0.);
            _weights[c++] = wt;

            _points[c]    = Point(0., x);
            _weights[c++] = wt;

            _points[c]    = Point(0.,-x);
            _weights[c++] = wt;

            break;
          }
        case 4: // Rotational invariant
          {
            _points[c]    = Point( x, y);
            _weights[c++] = wt;

            _points[c]    = Point(-x,-y);
            _weights[c++] = wt;

            _points[c]    = Point(-y, x);
            _weights[c++] = wt;

            _points[c]    = Point( y,-x);
            _weights[c++] = wt;

            break;
          }
        case 5: // Partial symmetry (Wissman's rules)
          {
            libmesh_assert_not_equal_to (x, 0.0);

            _points[c]    = Point( x, y);
            _weights[c++] = wt;

            _points[c]    = Point(-x, y);
            _weights[c++] = wt;

            break;
          }
        case 6: // Rectangular symmetry
          {
            _points[c]    = Point( x, y);
            _weights[c++] = wt;

            _points[c]    = Point(-x, y);
            _weights[c++] = wt;

            _points[c]    = Point(-x,-y);
            _weights[c++] = wt;

            _points[c]    = Point( x,-y);
            _weights[c++] = wt;

            break;
          }
        case 7: // Central symmetry
          {
            libmesh_assert_equal_to (x, 0.0);
            libmesh_assert_not_equal_to (y, 0.0);

            _points[c]    = Point(0., y);
            _weights[c++] = wt;

            _points[c]    = Point(0.,-y);
            _weights[c++] = wt;

            break;
          }
        default:
          libmesh_error_msg("Unknown symmetry!");
        } // end switch(rule_symmetry[i])
    }
}

void libMesh::QBase::tensor_product_hex ( const QBase & q1D ) [protected, inherited]

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 154 of file quadrature.C.

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

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

{
  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++;
        }
}

void libMesh::QBase::tensor_product_prism	(	const QBase &	q1D,
		const QBase &	q2D
	)		`[protected, inherited]`

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 181 of file quadrature.C.

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

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

{
  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++;
      }

}

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

Returns:: QMONOMIAL

Implements libMesh::QBase.

Definition at line 83 of file quadrature_monomial.h.

References libMesh::QMONOMIAL.

{ return QMONOMIAL; }

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

Returns:: the $i^{th}$ quadrature weight.

Definition at line 158 of file quadrature.h.

References libMesh::QBase::_weights.

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

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

void libMesh::QMonomial::wissmann_rule	(	const Real	rule_data[][3],
		const unsigned int	n_pts
	)		`[private]`

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

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

Definition at line 44 of file quadrature_monomial.C.

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

Referenced by init_2D().

{
  for (unsigned int i=0, c=0; i<n_pts; ++i)
    {
      _points[c]  = Point( rule_data[i][0], rule_data[i][1] );
      _weights[c++] = rule_data[i][2];

      // This may be an (x1,x2) -> (-x1,x2) point, in which case
      // we will also generate the mirror point using the same weight.
      if (rule_data[i][0] != static_cast<Real>(0.0))
        {
          _points[c]  = Point( -rule_data[i][0], rule_data[i][1] );
          _weights[c++] = rule_data[i][2];
        }
    }
}

Friends And Related Function Documentation

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

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

Definition at line 208 of file quadrature.C.

{
  q.print_info(os);
  return os;
}

Member Data Documentation

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

Actually holds the data.

Definition at line 118 of file reference_counter.h.

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

const unsigned int libMesh::QBase::_dim [protected, inherited]

The dimension

Definition at line 319 of file quadrature.h.

Referenced by libMesh::QBase::get_dim(), libMesh::QBase::init(), libMesh::QGauss::QGauss(), libMesh::QGaussLobatto::QGaussLobatto(), libMesh::QJacobi::QJacobi(), libMesh::QSimpson::QSimpson(), libMesh::QTrap::QTrap(), and libMesh::QBase::scale().

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

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

Definition at line 137 of file reference_counter.h.

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

Threads::spin_mutex libMesh::ReferenceCounter::_mutex [static, protected, inherited]

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 131 of file reference_counter.h.

Threads::atomic< unsigned int > libMesh::ReferenceCounter::_n_objects [static, protected, inherited]

The number of objects. Print the reference count information when the number returns to 0.

Definition at line 126 of file reference_counter.h.

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

const Order libMesh::QBase::_order [protected, inherited]

unsigned int libMesh::QBase::_p_level [protected, inherited]

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

Definition at line 336 of file quadrature.h.

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

std::vector<Point> libMesh::QBase::_points [protected, inherited]

The reference element locations of the quadrature points.

Definition at line 342 of file quadrature.h.

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

ElemType libMesh::QBase::_type [protected, inherited]

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

Definition at line 330 of file quadrature.h.

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

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

The value of the quadrature weights.

Definition at line 347 of file quadrature.h.

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

bool libMesh::QBase::allow_rules_with_negative_weights [inherited]

Flag (default true) controlling the use of quadrature rules with negative weights. Set this to false to ONLY use (potentially) safer but more expensive rules with all positive weights.

Negative weights typically appear in Gaussian quadrature rules over three-dimensional elements. Rules with negative weights can be unsuitable for some problems. For example, it is possible for a rule with negative weights to obtain a negative result when integrating a positive function.

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

Definition at line 229 of file quadrature.h.

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

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

Public Member Functions

Static Public Member Functions

Public Attributes

Protected Types

Protected Member Functions

Protected Attributes

Static Protected Attributes

Private Member Functions

Friends

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Friends And Related Function Documentation

Member Data Documentation