$extrastylesheet
libMesh::QGrundmann_Moller Class Reference

#include <quadrature_gm.h>

Inheritance diagram for libMesh::QGrundmann_Moller:

List of all members.

Public Member Functions

 QGrundmann_Moller (const unsigned int _dim, const Order _order=INVALID_ORDER)
 ~QGrundmann_Moller ()
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< QBasebuild (const std::string &name, const unsigned int _dim, const Order _order=INVALID_ORDER)
static UniquePtr< QBasebuild (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_3D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0)
void init_2D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0)
void gm_rule (unsigned int s, unsigned int dim)
void compose_all (unsigned int s, unsigned int p, std::vector< std::vector< unsigned int > > &result)

Friends

std::ostream & operator<< (std::ostream &os, const QBase &q)

Detailed Description

This class implements the Grundmann-Moller quadrature rules for tetrahedra. The GM rules are well-defined for simplices of arbitrary dimension and to any order, but the rules by Dunavant for two-dimensional simplices are in general superior. This is primarily due to the fact that the GM rules contain a significant proportion of negative weights, making them susceptible to round-off error at high-order.

The GM rules are interesting in 3D because they overlap with the conical product rules at higher order while having significantly fewer evaluation points, making them potentially much more efficient. The table below gives a comparison between the number of points in a conical product (CP) rule and the GM rule of equivalent order. The GM rules are defined to be exact for polynomials of degree d=2*s+1, s=0,1,2,3,... The table also gives the percentage of each GM rule's weights which are negative. The percentage of negative weights appears to approach 50, and the amplification factor (a measure of the effect of round-off) defined as

amp. factor = $ \frac{1}{V} \sum \|w_i\|, $

where V is the volume of the reference element, grows like exp(C*s). (A rule with all positive weights has an amplification factor of 1.0 by definition.)

 * s   degree  n_pts(conical)  n_pts(GM)   % neg wts  amp. factor
 * ------------------------------------------------------------------------
 * 0   1       1               1            0.00      1.00e+00
 * 1   3       8               5           20.00      2.60e+00
 * 2   5       27              15          26.67      5.63e+00
 * 3   7       64              35          31.43      1.19e+01
 * 4   9       125             70          34.29      2.54e+01
 * 5   11      216             126         36.51      5.41e+01
 * 6   13      343             210         38.10      1.16e+02
 * 7   15      512             330         39.39      2.51e+02
 * 8   17      729             495         40.40      5.45e+02
 * 9   19      1000            715         41.26      1.19e+03
 * 10  21      1331            1001        41.96      2.59e+03
 * 11  23      1728            1365        42.56      5.68e+03
 * 12  25      2197            1820        43.08      1.25e+04
 * 13  27      2744            2380        43.53      2.75e+04
 * 14  29      3375            3060        43.92      6.07e+04
 * 15  31      4096            3876        44.27      1.34e+05
 * 16  33      4913            4845        44.58      2.97e+05
 * 17  35      5832            5985        44.86      6.59e+05 <= Conical rule has fewer points for degree >= 34
 * 18  37      6859            7315        45.11      1.46e+06
 * 19  39      8000            8855        45.34      3.25e+06
 * 20  41      9261            10626       45.55      7.23e+06
 * 21  43      10648           12650       45.74      1.61e+07
 * 

Reference: Axel Grundmann and Michael M"{o}ller, "Invariant Integration Formulas for the N-Simplex by Combinatorial Methods," SIAM Journal on Numerical Analysis, Volume 15, Number 2, April 1978, pages 282-290.

Reference LGPL Fortran90 code by John Burkardt can be found here: http://people.scs.fsu.edu/~burkardt/f_src/gm_rules/gm_rules.html

Author:
John W. Peterson, 2008

Definition at line 97 of file quadrature_gm.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::QGrundmann_Moller::QGrundmann_Moller ( const unsigned int  _dim,
const Order  _order = INVALID_ORDER 
)

Constructor. Declares the order of the quadrature rule.

Definition at line 32 of file quadrature_gm.C.

                                                    : QBase(d,o)
{
}

Destructor.

Definition at line 39 of file quadrature_gm.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::QGrundmann_Moller::compose_all ( unsigned int  s,
unsigned int  p,
std::vector< std::vector< unsigned int > > &  result 
) [private]

Routine which generates p-compositions of a given order, s, as well as permutations thereof. This routine is called internally by the gm_rule() routine, you should not call this yourself!

Definition at line 145 of file quadrature_gm.C.

Referenced by gm_rule().

{
  // Clear out results remaining from previous calls
  result.clear();

  // Allocate storage for a workspace.  The workspace will periodically
  // be copied into the result container.
  std::vector<unsigned int> workspace(p);

  // The first result is always (s,0,...,0)
  workspace[0] = s;
  result.push_back(workspace);

  // the value of the first non-zero entry
  unsigned int head_value=s;

  // When head_index=-1, it refers to "off the front" of the array.  Therefore,
  // this needs to be a regular int rather than unsigned.  I initially tried to
  // do this with head_index unsigned and an else statement below, but then there
  // is the special case: (1,0,...,0) which does not work correctly.
  int head_index = -1;

  // At the end, all the entries will be in the final slot of workspace
  while (workspace.back() != s)
    {
      // Uncomment for debugging
      //libMesh::out << "previous head_value=" << head_value << " -> ";

      // If the previous head value is still larger than 1, reset the index
      // to "off the front" of the array
      if (head_value > 1)
        head_index = -1;

      // Either move the index onto the front of the array or on to
      // the next value.
      head_index++;

      // Get current value of the head entry
      head_value = workspace[head_index];

      // Uncomment for debugging
      //std::copy(workspace.begin(), workspace.end(), std::ostream_iterator<int>(libMesh::out, " "));
      //libMesh::out << ", head_index=" << head_index;
      //libMesh::out << ", head_value=" << head_value << " -> ";

      // Put a zero into the head_index of the array.  If head_index==0,
      // this will be overwritten in the next line with head_value-1.
      workspace[head_index] = 0;

      // The initial entry gets the current head value, minus 1.
      // If head_value > 1, the next loop iteration will start back
      // at workspace[0] again.
      libmesh_assert_greater (head_value, 0);
      workspace[0] = head_value - 1;

      // Increment the head+1 value
      workspace[head_index+1] += 1;

      // Save this composition in the results
      result.push_back(workspace);

      // Uncomment for debugging
      //std::copy(workspace.begin(), workspace.end(), std::ostream_iterator<int>(libMesh::out, " "));
      //libMesh::out<<"\n";
    }
}

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]
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]
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;  }
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::QGrundmann_Moller::gm_rule ( unsigned int  s,
unsigned int  dim 
) [private]

This routine is called from init_2D() and init_3D(). It actually fills the _points and _weights vectors for a given rule index, s and dimension, dim.

Definition at line 47 of file quadrature_gm.C.

References libMesh::QBase::_points, libMesh::QBase::_weights, compose_all(), libMesh::libmesh_assert(), std::max(), libMesh::Real, and libMesh::MeshTools::weight().

Referenced by init_2D(), and init_3D().

{
  // Only dim==2 or dim==3 is allowed
  libmesh_assert(dim==2 || dim==3);

  // A GM rule of index s can integrate polynomials of degree 2*s+1 exactly
  const unsigned int degree = 2*s+1;

  // The number of points for rule of index s is
  // (dim+1+s)! / (dim+1)! / s!
  // In 3D, this is = 1/24 * Pi_{i=1}^4 (s+i)
  // In 2D, this is =  1/6 * Pi_{i=1}^3 (s+i)
  const unsigned int n_pts = dim==2 ? (s+3)*(s+2)*(s+1) / 6 : (s+4)*(s+3)*(s+2)*(s+1) / 24;
  //libMesh::out << "n_pts=" << n_pts << std::endl;

  // Allocate space for points and weights
  _points.resize(n_pts);
  _weights.resize(n_pts);

  // (-1)^i -> This one flips sign at each iteration of the i-loop below.
  int one_pm=1;

  // Where we store all the integer point compositions/permutations
  std::vector<std::vector<unsigned int> > permutations;

  // Index into the vector where we should start adding the next round of points/weights
  std::size_t offset=0;

  // Implement the GM formula 4.1 on page 286 of the paper
  for (unsigned int i=0; i<=s; ++i)
    {
      // Get all the ordered compositions (and their permutations)
      // of |beta| = s-i into dim+1 parts
      compose_all(s-i, dim+1, permutations);
      //libMesh::out << "n. permutations=" << permutations.size() << std::endl;

      for (unsigned int p=0; p<permutations.size(); ++p)
        {
          // We use the first dim entries of each permutation to
          // construct an integration point.
          for (unsigned int j=0; j<dim; ++j)
            _points[offset+p](j) =
              static_cast<Real>(2.*permutations[p][j] + 1.) /
              static_cast<Real>(  degree + dim - 2.*i     );
        }

      // Compute the weight for this i, being careful to avoid overflow.
      // This technique is borrowed from Burkardt's code as well.
      // Use once for each of the points obtained from the permutations array.
      Real weight = one_pm;

      // This for loop needs to run for dim, degree, or dim+degree-i iterations,
      // whichever is largest.
      const unsigned int weight_loop_index =
        std::max(dim, std::max(degree, degree+dim-i))+1;

      for (unsigned int j=1; j<weight_loop_index; ++j)
        {
          if (j <= degree) // Accumulate (d+n-2i)^d term
            weight *= static_cast<Real>(degree+dim-2*i);

          if (j <= 2*s) // Accumulate 2^{-2s}
            weight *= 0.5;

          if (j <= i) // Accumulate (i!)^{-1}
            weight /= static_cast<Real>(j);

          if (j <= degree+dim-i) // Accumulate ( (d+n-i)! )^{-1}
            weight /= static_cast<Real>(j);
        }

      // This is the weight for each of the points computed previously
      for (unsigned int j=0; j<permutations.size(); ++j)
        _weights[offset+j] = weight;

      // Change sign for next iteration
      one_pm = -one_pm;

      // Update offset for the next set of points
      offset += permutations.size();
    }
}
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]
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::QGrundmann_Moller::init_1D ( const ElemType  type,
unsigned  p_level = 0 
) [inline, private, virtual]

Initializes the 1D quadrature rule by filling the points and weights vectors with the appropriate values. The order of the rule will be defined by the implementing class. It is assumed that derived quadrature rules will at least define the init_1D function, therefore it is pure virtual.

Implements libMesh::QBase.

Definition at line 120 of file quadrature_gm.h.

  {
    // See about making this non-pure virtual in the base class
    libmesh_not_implemented();
  }
void libMesh::QGrundmann_Moller::init_2D ( const ElemType  _type = INVALID_ELEM,
unsigned int  p_level = 0 
) [private, virtual]

Initialize a 2D GM rule. Only makes sense for Tris.

Reimplemented from libMesh::QBase.

Definition at line 28 of file quadrature_gm_2D.C.

References libMesh::QBase::_order, libMesh::QBase::allow_rules_with_negative_weights, gm_rule(), libMesh::QBase::libmesh_error_msg(), libMesh::TRI3, and libMesh::TRI6.

{
  // Nearly all GM rules contain negative weights, so if you are not
  // allowing rules with negative weights, we cannot continue!
  if (!allow_rules_with_negative_weights)
    libmesh_error_msg("You requested a Grundmann-Moller rule but\n"        \
                      << "are not allowing rules with negative weights!\n" \
                      << "Either select a different quadrature class or\n" \
                      << "set allow_rules_with_negative_weights==true.");

  switch (type_in)
    {
    case TRI3:
    case TRI6:
      {
        switch(_order + 2*p)
          {

          default:
            {
              // Untested above _order=23 but should work...
              gm_rule((_order + 2*p)/2, /*dim=*/2);
              return;
            }
          } // end switch (order)
      } // end case TRI3, TRI6

    default:
      libmesh_error_msg("ERROR: Unsupported element type: " << type_in);
    } // end switch (type_in)
}
void libMesh::QGrundmann_Moller::init_3D ( const ElemType  _type = INVALID_ELEM,
unsigned int  p_level = 0 
) [private, virtual]

Initialize a 3D GM rule. Only makes sense for Tets.

Reimplemented from libMesh::QBase.

Definition at line 28 of file quadrature_gm_3D.C.

References libMesh::QBase::_order, libMesh::QBase::_points, libMesh::QBase::_weights, libMesh::QBase::allow_rules_with_negative_weights, libMesh::CONSTANT, libMesh::FIFTH, libMesh::FIRST, libMesh::FOURTH, gm_rule(), libMesh::QBase::libmesh_error_msg(), libMesh::SECOND, libMesh::SEVENTH, libMesh::SIXTH, libMesh::TET10, libMesh::TET4, and libMesh::THIRD.

{
  // Nearly all GM rules contain negative weights, so if you are not
  // allowing rules with negative weights, we cannot continue!
  if (!allow_rules_with_negative_weights)
    libmesh_error_msg("You requested a Grundmann-Moller rule but\n"        \
                      << "are not allowing rules with negative weights!\n" \
                      << "Either select a different quadrature class or\n" \
                      << "set allow_rules_with_negative_weights==true.");

  switch (type_in)
    {
    case TET4:
    case TET10:
      {
        switch(_order + 2*p)
          {
            // We hard-code the first few orders based on output from
            // the mp-quadrature library:
            //
            // https://code.google.com/p/mp-quadrature
            //
            // The points are ordered in such a way that the amount of
            // round-off error in the quadrature calculations is
            // (hopefully) minimized.  These orderings were obtained
            // via a simple permutation optimization strategy designed
            // to produce the smallest overall average error while
            // integrating all polynomials of degree <= d.
          case CONSTANT:
          case FIRST:
            {
              _points.resize(1);
              _weights.resize(1);

              _points[0](0) = 1.L/4.L;
              _points[0](1) = 1.L/4.L;
              _points[0](2) = 1.L/4.L;

              _weights[0] = 1.L/6.L;
              return;
            }

          case SECOND:
          case THIRD:
            {
              _points.resize(5);
              _weights.resize(5);

              _points[0](0) = 1.L/6.L;  _points[0](1) = 1.L/6.L;  _points[0](2) = 1.L/2.L;  _weights[0] =  3.L/40.L;
              _points[1](0) = 1.L/6.L;  _points[1](1) = 1.L/6.L;  _points[1](2) = 1.L/6.L;  _weights[1] =  3.L/40.L;
              _points[2](0) = 1.L/4.L;  _points[2](1) = 1.L/4.L;  _points[2](2) = 1.L/4.L;  _weights[2] = -2.L/15.L;
              _points[3](0) = 1.L/2.L;  _points[3](1) = 1.L/6.L;  _points[3](2) = 1.L/6.L;  _weights[3] =  3.L/40.L;
              _points[4](0) = 1.L/6.L;  _points[4](1) = 1.L/2.L;  _points[4](2) = 1.L/6.L;  _weights[4] =  3.L/40.L;
              return;
            }

          case FOURTH:
          case FIFTH:
            {
              _points.resize(15);
              _weights.resize(15);

              _points[ 0](0) = 1.L/8.L;  _points[ 0](1) = 3.L/8.L;  _points[ 0](2) = 1.L/8.L;  _weights[ 0] =  16.L/315.L;
              _points[ 1](0) = 1.L/8.L;  _points[ 1](1) = 1.L/8.L;  _points[ 1](2) = 5.L/8.L;  _weights[ 1] =  16.L/315.L;
              _points[ 2](0) = 3.L/8.L;  _points[ 2](1) = 1.L/8.L;  _points[ 2](2) = 1.L/8.L;  _weights[ 2] =  16.L/315.L;
              _points[ 3](0) = 1.L/6.L;  _points[ 3](1) = 1.L/2.L;  _points[ 3](2) = 1.L/6.L;  _weights[ 3] = -27.L/280.L;
              _points[ 4](0) = 3.L/8.L;  _points[ 4](1) = 1.L/8.L;  _points[ 4](2) = 3.L/8.L;  _weights[ 4] =  16.L/315.L;
              _points[ 5](0) = 1.L/8.L;  _points[ 5](1) = 3.L/8.L;  _points[ 5](2) = 3.L/8.L;  _weights[ 5] =  16.L/315.L;
              _points[ 6](0) = 1.L/6.L;  _points[ 6](1) = 1.L/6.L;  _points[ 6](2) = 1.L/6.L;  _weights[ 6] = -27.L/280.L;
              _points[ 7](0) = 1.L/6.L;  _points[ 7](1) = 1.L/6.L;  _points[ 7](2) = 1.L/2.L;  _weights[ 7] = -27.L/280.L;
              _points[ 8](0) = 1.L/8.L;  _points[ 8](1) = 1.L/8.L;  _points[ 8](2) = 1.L/8.L;  _weights[ 8] =  16.L/315.L;
              _points[ 9](0) = 1.L/4.L;  _points[ 9](1) = 1.L/4.L;  _points[ 9](2) = 1.L/4.L;  _weights[ 9] =    2.L/45.L;
              _points[10](0) = 1.L/8.L;  _points[10](1) = 5.L/8.L;  _points[10](2) = 1.L/8.L;  _weights[10] =  16.L/315.L;
              _points[11](0) = 1.L/2.L;  _points[11](1) = 1.L/6.L;  _points[11](2) = 1.L/6.L;  _weights[11] = -27.L/280.L;
              _points[12](0) = 1.L/8.L;  _points[12](1) = 1.L/8.L;  _points[12](2) = 3.L/8.L;  _weights[12] =  16.L/315.L;
              _points[13](0) = 5.L/8.L;  _points[13](1) = 1.L/8.L;  _points[13](2) = 1.L/8.L;  _weights[13] =  16.L/315.L;
              _points[14](0) = 3.L/8.L;  _points[14](1) = 3.L/8.L;  _points[14](2) = 1.L/8.L;  _weights[14] =  16.L/315.L;
              return;
            }

          case SIXTH:
          case SEVENTH:
            {
              _points.resize(35);
              _weights.resize(35);

              _points[ 0](0) = 3.L/10.L;  _points[ 0](1) = 1.L/10.L;  _points[ 0](2) = 3.L/10.L;  _weights[ 0] = 3125.L/72576.L;
              _points[ 1](0) =  1.L/6.L;  _points[ 1](1) =  1.L/2.L;  _points[ 1](2) =  1.L/6.L;  _weights[ 1] =   243.L/4480.L;
              _points[ 2](0) =  1.L/6.L;  _points[ 2](1) =  1.L/6.L;  _points[ 2](2) =  1.L/2.L;  _weights[ 2] =   243.L/4480.L;
              _points[ 3](0) =  1.L/2.L;  _points[ 3](1) = 1.L/10.L;  _points[ 3](2) = 1.L/10.L;  _weights[ 3] = 3125.L/72576.L;
              _points[ 4](0) = 3.L/10.L;  _points[ 4](1) = 1.L/10.L;  _points[ 4](2) = 1.L/10.L;  _weights[ 4] = 3125.L/72576.L;
              _points[ 5](0) = 3.L/10.L;  _points[ 5](1) = 3.L/10.L;  _points[ 5](2) = 1.L/10.L;  _weights[ 5] = 3125.L/72576.L;
              _points[ 6](0) =  1.L/2.L;  _points[ 6](1) =  1.L/6.L;  _points[ 6](2) =  1.L/6.L;  _weights[ 6] =   243.L/4480.L;
              _points[ 7](0) = 3.L/10.L;  _points[ 7](1) = 1.L/10.L;  _points[ 7](2) =  1.L/2.L;  _weights[ 7] = 3125.L/72576.L;
              _points[ 8](0) = 7.L/10.L;  _points[ 8](1) = 1.L/10.L;  _points[ 8](2) = 1.L/10.L;  _weights[ 8] = 3125.L/72576.L;
              _points[ 9](0) =  3.L/8.L;  _points[ 9](1) =  1.L/8.L;  _points[ 9](2) =  1.L/8.L;  _weights[ 9] =  -256.L/2835.L;
              _points[10](0) = 3.L/10.L;  _points[10](1) = 3.L/10.L;  _points[10](2) = 3.L/10.L;  _weights[10] = 3125.L/72576.L;
              _points[11](0) = 1.L/10.L;  _points[11](1) =  1.L/2.L;  _points[11](2) = 3.L/10.L;  _weights[11] = 3125.L/72576.L;
              _points[12](0) = 1.L/10.L;  _points[12](1) = 1.L/10.L;  _points[12](2) = 7.L/10.L;  _weights[12] = 3125.L/72576.L;
              _points[13](0) =  1.L/2.L;  _points[13](1) = 1.L/10.L;  _points[13](2) = 3.L/10.L;  _weights[13] = 3125.L/72576.L;
              _points[14](0) = 1.L/10.L;  _points[14](1) = 7.L/10.L;  _points[14](2) = 1.L/10.L;  _weights[14] = 3125.L/72576.L;
              _points[15](0) = 1.L/10.L;  _points[15](1) =  1.L/2.L;  _points[15](2) = 1.L/10.L;  _weights[15] = 3125.L/72576.L;
              _points[16](0) =  1.L/6.L;  _points[16](1) =  1.L/6.L;  _points[16](2) =  1.L/6.L;  _weights[16] =   243.L/4480.L;
              _points[17](0) =  3.L/8.L;  _points[17](1) =  1.L/8.L;  _points[17](2) =  3.L/8.L;  _weights[17] =  -256.L/2835.L;
              _points[18](0) =  1.L/8.L;  _points[18](1) =  1.L/8.L;  _points[18](2) =  5.L/8.L;  _weights[18] =  -256.L/2835.L;
              _points[19](0) = 1.L/10.L;  _points[19](1) = 1.L/10.L;  _points[19](2) = 3.L/10.L;  _weights[19] = 3125.L/72576.L;
              _points[20](0) =  1.L/8.L;  _points[20](1) =  3.L/8.L;  _points[20](2) =  3.L/8.L;  _weights[20] =  -256.L/2835.L;
              _points[21](0) =  5.L/8.L;  _points[21](1) =  1.L/8.L;  _points[21](2) =  1.L/8.L;  _weights[21] =  -256.L/2835.L;
              _points[22](0) =  1.L/8.L;  _points[22](1) =  5.L/8.L;  _points[22](2) =  1.L/8.L;  _weights[22] =  -256.L/2835.L;
              _points[23](0) = 1.L/10.L;  _points[23](1) = 3.L/10.L;  _points[23](2) = 1.L/10.L;  _weights[23] = 3125.L/72576.L;
              _points[24](0) =  1.L/4.L;  _points[24](1) =  1.L/4.L;  _points[24](2) =  1.L/4.L;  _weights[24] =     -8.L/945.L;
              _points[25](0) =  1.L/8.L;  _points[25](1) =  1.L/8.L;  _points[25](2) =  3.L/8.L;  _weights[25] =  -256.L/2835.L;
              _points[26](0) =  3.L/8.L;  _points[26](1) =  3.L/8.L;  _points[26](2) =  1.L/8.L;  _weights[26] =  -256.L/2835.L;
              _points[27](0) =  1.L/8.L;  _points[27](1) =  3.L/8.L;  _points[27](2) =  1.L/8.L;  _weights[27] =  -256.L/2835.L;
              _points[28](0) = 1.L/10.L;  _points[28](1) = 3.L/10.L;  _points[28](2) =  1.L/2.L;  _weights[28] = 3125.L/72576.L;
              _points[29](0) = 3.L/10.L;  _points[29](1) =  1.L/2.L;  _points[29](2) = 1.L/10.L;  _weights[29] = 3125.L/72576.L;
              _points[30](0) = 1.L/10.L;  _points[30](1) = 1.L/10.L;  _points[30](2) =  1.L/2.L;  _weights[30] = 3125.L/72576.L;
              _points[31](0) =  1.L/2.L;  _points[31](1) = 3.L/10.L;  _points[31](2) = 1.L/10.L;  _weights[31] = 3125.L/72576.L;
              _points[32](0) =  1.L/8.L;  _points[32](1) =  1.L/8.L;  _points[32](2) =  1.L/8.L;  _weights[32] =  -256.L/2835.L;
              _points[33](0) = 1.L/10.L;  _points[33](1) = 3.L/10.L;  _points[33](2) = 3.L/10.L;  _weights[33] = 3125.L/72576.L;
              _points[34](0) = 1.L/10.L;  _points[34](1) = 1.L/10.L;  _points[34](2) = 1.L/10.L;  _weights[34] = 3125.L/72576.L;
              return;
            }

          default:
            {
              // Untested above _order=23 but should work...
              gm_rule((_order + 2*p)/2, /*dim=*/3);
              return;
            }
          } // end switch (order)
      } // end case TET4, TET10

    default:
      libmesh_error_msg("ERROR: Unsupported element type: " << type_in);
    } // end switch (type_in)
}
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; }
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().

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::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::QGrundmann_Moller::type ( ) const [inline, virtual]
Returns:
QGRUNDMANN_MOLLER

Implements libMesh::QBase.

Definition at line 115 of file quadrature_gm.h.

References libMesh::QGRUNDMANN_MOLLER.

{ return QGRUNDMANN_MOLLER; }
Real libMesh::QBase::w ( const unsigned int  i) const [inline, inherited]

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

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

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

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

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

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


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