Poster of Linux kernelThe best gift for a Linux geek
QGauss

QGauss

Section: C Library Functions (3) Updated: Thu Apr 7 2011
Local index Up
 

NAME

QGauss -  

SYNOPSIS


#include <quadrature_gauss.h>

Inherits QBase.  

Public Member Functions


QGauss (const unsigned int _dim, const Order _order=INVALID_ORDER)

~QGauss ()

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

void init (const ElemType _type=INVALID_ELEM, unsigned int p_level=0)

Order get_order () const

void print_info (std::ostream &os=std::cout) const

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

Static Public Member Functions


static AutoPtr< QBase > build (const std::string &name, const unsigned int _dim, const Order _order=INVALID_ORDER)

static AutoPtr< QBase > build (const QuadratureType _qt, const unsigned int _dim, const Order _order=INVALID_ORDER)

static void print_info ()

static std::string get_info ()

static unsigned int n_objects ()
 

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)

void increment_constructor_count (const std::string &name)

void increment_destructor_count (const std::string &name)
 

Protected Attributes


std::cerr<< 'ERROR: Seems as if this quadrature rule'<< std::endl<< ' is not implemented for 2D.'<< std::endl;libmesh_error();}#endif virtual void init_3D(const ElemType, unsigned int=0)#ifndef DEBUG{}#else{std::cerr<< 'ERROR: Seems as if this quadrature rule'<< std::endl<< ' is not implemented for 3D.'<< std::endl;libmesh_error();}#endif void tensor_product_quad(const QBase &q1D);void tensor_product_hex(const QBase &q1D);void tensor_product_prism(const QBase &q1D, const QBase &q2D);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
 

Private Member Functions


void init_1D (const ElemType _type=INVALID_ELEM, unsigned int p_level=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 dunavant_rule (const Real rule_data[][4], const unsigned int n_pts)

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

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

Friends


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

Detailed Description

This class implemenets specific orders of Gauss quadrature. Gauss quadrature rules of Order p have the property of integrating polynomials of degree p exactly.

Definition at line 40 of file quadrature_gauss.h.  

Member Typedef Documentation

 

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > ReferenceCounter::Counts [protected, inherited]Data structure to log the information. The log is identified by the class name.

Definition at line 105 of file reference_counter.h.  

Constructor & Destructor Documentation

 

QGauss::QGauss (const unsigned int_dim, const Order_order = INVALID_ORDER) [inline]Constructor. Declares the order of the quadrature rule.

Definition at line 102 of file quadrature_gauss.h.

References libMeshEnums::EDGE2, and QBase::init().

                              : QBase(d,o)
{
  // explicitly call the init function in 1D since the
  // other tensor-product rules require this one.
  // note that EDGE will not be used internally, however
  // if we called the function with INVALID_ELEM it would try to
  // be smart and return, thinking it had already done the work.
  if (_dim == 1)
    init(EDGE2);
}
 

QGauss::~QGauss () [inline]Destructor.

Definition at line 118 of file quadrature_gauss.h.

{
}
 

Member Function Documentation

 

AutoPtr< QBase > 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 AutoPtr<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 36 of file quadrature_build.C.

References Utility::string_to_enum< QuadratureType >().

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

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

AutoPtr< QBase > 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 AutoPtr<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 47 of file quadrature_build.C.

References libMeshEnums::FIRST, libMeshEnums::FORTYTHIRD, libMeshEnums::QCLOUGH, libMeshEnums::QGAUSS, libMeshEnums::QJACOBI_1_0, libMeshEnums::QJACOBI_2_0, libMeshEnums::QSIMPSON, libMeshEnums::QTRAP, libMeshEnums::THIRD, and libMeshEnums::TWENTYTHIRD.

{
  switch (_qt)
    {
      
    case QCLOUGH:
      {
#ifdef DEBUG
        if (_order > TWENTYTHIRD)
          {
            std::cout << 'WARNING: Clough quadrature implemented' << std::endl
                      << ' up to TWENTYTHIRD order.' << std::endl;
          }
#endif

        AutoPtr<QBase> ap(new QClough(_dim, _order));
        return ap;
      }

    case QGAUSS:
      {

#ifdef DEBUG
        if (_order > FORTYTHIRD)
          {
            std::cout << 'WARNING: Gauss quadrature implemented' << std::endl
                      << ' up to FORTYTHIRD order.' << std::endl;
          }
#endif

        AutoPtr<QBase> ap(new QGauss(_dim, _order));
        return ap;
      }

    case QJACOBI_1_0:
      {

#ifdef DEBUG
        if (_order > TWENTYTHIRD)
          {
            std::cout << 'WARNING: Jacobi(1,0) quadrature implemented' << std::endl
                      << ' up to TWENTYTHIRD order.' << std::endl;
          }

        if (_dim > 1)
          {
            std::cout << 'WARNING: Jacobi(1,0) quadrature implemented' << std::endl
                      << ' in 1D only.' << std::endl;
          }
#endif

        AutoPtr<QBase> ap(new QJacobi(_dim, _order, 1, 0));
        return ap;
      }

    case QJACOBI_2_0:
      {

#ifdef DEBUG
        if (_order > TWENTYTHIRD)
          {
            std::cout << 'WARNING: Jacobi(2,0) quadrature implemented' << std::endl
                      << ' up to TWENTYTHIRD order.' << std::endl;
          }

        if (_dim > 1)
          {
            std::cout << 'WARNING: Jacobi(2,0) quadrature implemented' << std::endl
                      << ' in 1D only.' << std::endl;
          }
#endif

        AutoPtr<QBase> ap(new QJacobi(_dim, _order, 2, 0));
        return ap;
      }

    case QSIMPSON:
      {

#ifdef DEBUG
        if (_order > THIRD)
          {
            std::cout << 'WARNING: Simpson rule provides only' << std::endl
                      << ' THIRD order!' << std::endl;
          }
#endif

        AutoPtr<QBase> ap(new QSimpson(_dim));
        return ap;
      }

    case QTRAP:
      {

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

        AutoPtr<QBase> ap(new QTrap(_dim));
        return ap;
      }


    default:
      { 
        std::cerr << 'ERROR: Bad qt=' << _qt << std::endl;
        libmesh_error();
      }
    }


  libmesh_error();
  AutoPtr<QBase> ap(NULL);
  return ap;
}
 

void QGauss::dunavant_rule (const Realrule_data[][4], const unsigned intn_pts) [private]The Dunavant rule is for triangles. It takes permutation points and weights in a specific format as input and fills the pre-sized _points and _weights vectors. This function is only used internally by the TRI geometric elements.

Definition at line 265 of file quadrature_gauss.C.

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

Referenced by init_2D().

{
  // The input data array has 4 columns.  The first 3 are the permutation points.
  // The last column is the weights for a given set of permutation points.  A zero
  // in two of the first 3 columns implies the point is a 1-permutation (centroid).
  // A zero in one of the first 3 columns implies the point is a 3-permutation.
  // Otherwise each point is assumed to be a 6-permutation.

  // Always insert into the points & weights vector relative to the offset 
  unsigned int offset=0;

  
  for (unsigned int p=0; p<n_pts; ++p)
    {

      // There must always be a non-zero entry to start the row
      libmesh_assert( rule_data[p][0] != static_cast<Real>(0.0) );
      
      // A zero weight may imply you did not set up the raw data correctly
      libmesh_assert( rule_data[p][3] != static_cast<Real>(0.0) );

      // What kind of point is this?
      // One non-zero entry in first 3 cols   ? 1-perm (centroid) point = 1
      // Two non-zero entries in first 3 cols ? 3-perm point            = 3
      // Three non-zero entries               ? 6-perm point            = 6
      unsigned int pointtype=1;
      
      if (rule_data[p][1] != static_cast<Real>(0.0))      
        {
          if (rule_data[p][2] != static_cast<Real>(0.0))
            pointtype = 6;
          else
            pointtype = 3;
        }
      
      switch (pointtype)
        {
        case 1:
          {
            // Be sure we have enough space to insert this point
            libmesh_assert(offset + 0 < _points.size());
            
            // The point has only a single permutation (the centroid!)
            _points[offset  + 0] = Point(rule_data[p][0], rule_data[p][0]);

            // The weight is always the last entry in the row.
            _weights[offset + 0] = rule_data[p][3];

            offset += 1;
            break;
          }

        case 3:
          {
            // Be sure we have enough space to insert these points
            libmesh_assert(offset + 2 < _points.size());
            
            // Here it's understood the second entry is to be used twice, and
            // thus there are three possible permutations.
            _points[offset + 0] = Point(rule_data[p][0], rule_data[p][1]);
            _points[offset + 1] = Point(rule_data[p][1], rule_data[p][0]);
            _points[offset + 2] = Point(rule_data[p][1], rule_data[p][1]);

            for (unsigned int j=0; j<3; ++j)
              _weights[offset + j] = rule_data[p][3];
            
            offset += 3;
            break;
          }

        case 6:
          {
            // Be sure we have enough space to insert these points
            libmesh_assert(offset + 5 < _points.size());
            
            // Three individual entries with six permutations.
            _points[offset + 0] = Point(rule_data[p][0], rule_data[p][1]);
            _points[offset + 1] = Point(rule_data[p][0], rule_data[p][2]);
            _points[offset + 2] = Point(rule_data[p][1], rule_data[p][0]);
            _points[offset + 3] = Point(rule_data[p][1], rule_data[p][2]);
            _points[offset + 4] = Point(rule_data[p][2], rule_data[p][0]);
            _points[offset + 5] = Point(rule_data[p][2], rule_data[p][1]);

            for (unsigned int j=0; j<6; ++j)
              _weights[offset + j] = rule_data[p][3];
            
            offset += 6;
            break;
          }

        default:
          {
            std::cerr << 'Don't know what to do with this many permutation points!' << std::endl;
            libmesh_error();
          }
        }
    }
}
 

void QGauss::dunavant_rule2 (const Real *wts, const Real *a, const Real *b, const unsigned int *permutation_ids, const unsigned intn_wts) [private]

Definition at line 187 of file quadrature_gauss.C.

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

Referenced by init_2D().

{
  // Figure out how many total points by summing up the entries
  // in the permutation_ids array, and resize the _points and _weights
  // vectors appropriately.
  unsigned int total_pts = 0;
  for (unsigned int p=0; p<n_wts; ++p)
    total_pts += permutation_ids[p];

  // Resize point and weight vectors appropriately.
  _points.resize(total_pts);
  _weights.resize(total_pts);
    
  // Always insert into the points & weights vector relative to the offset 
  unsigned int offset=0;
  
  for (unsigned int p=0; p<n_wts; ++p)
    {
      switch (permutation_ids[p])
        {
        case 1:
          {
            // The point has only a single permutation (the centroid!)
            // So we don't even need to look in the a or b arrays.
            _points [offset  + 0] = Point(1.0L/3.0L, 1.0L/3.0L);
            _weights[offset + 0] = wts[p];
            
            offset += 1;
            break;
          }


        case 3:
          {
            // For this type of rule, don't need to look in the b array.
            _points[offset + 0] = Point(a[p],         a[p]);         // (a,a)
            _points[offset + 1] = Point(a[p],         1.L-2.L*a[p]); // (a,1-2a)
            _points[offset + 2] = Point(1.L-2.L*a[p], a[p]);         // (1-2a,a)

            for (unsigned int j=0; j<3; ++j)
              _weights[offset + j] = wts[p];
            
            offset += 3;
            break;
          }

        case 6:
          {
            // This type of point uses all 3 arrays...
            _points[offset + 0] = Point(a[p], b[p]);
            _points[offset + 1] = Point(b[p], a[p]);
            _points[offset + 2] = Point(a[p], 1.L-a[p]-b[p]);
            _points[offset + 3] = Point(1.L-a[p]-b[p], a[p]);
            _points[offset + 4] = Point(b[p], 1.L-a[p]-b[p]);
            _points[offset + 5] = Point(1.L-a[p]-b[p], b[p]);

            for (unsigned int j=0; j<6; ++j)
              _weights[offset + j] = wts[p];
            
            offset += 6;
            break;
          }
          
        default:
          {
            std::cerr << 'Unknown permutation id: ' << permutation_ids[p] << '!' << std::endl;
            libmesh_error();
          }
        }
    }
  
}
 

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

the dimension of the quadrature rule.

Definition at line 121 of file quadrature.h.

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

{ return _dim;  }
 

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

the current element type we're set up for

Definition at line 103 of file quadrature.h.

    { return _type; }
 

std::string ReferenceCounter::get_info () [static, inherited]Gets a string containing the reference information.

Definition at line 45 of file reference_counter.C.

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

Referenced by ReferenceCounter::print_info().

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

  std::ostringstream out;
  
  out << '
      << ' ---------------------------------------------------------------------------- 
      << '| Reference count information                                                |
      << ' ---------------------------------------------------------------------------- ;
  
  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;

      out << '| ' << name << ' reference count information:
          << '|  Creations:    ' << creations    << '
          << '|  Destructions: ' << destructions << ';
    }
  
  out << ' ---------------------------------------------------------------------------- ;

  return out.str();

#else

  return '';
  
#endif
}
 

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

the order of the quadrature rule.

Definition at line 167 of file quadrature.h.

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

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

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

the current p refinement level we're initialized with

Definition at line 109 of file quadrature.h.

    { return _p_level; }
 

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

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

Definition at line 127 of file quadrature.h.

References QBase::_points.

Referenced by FE< Dim, T >::edge_reinit(), QClough::init_1D(), QMonomial::init_2D(), init_2D(), QClough::init_2D(), QMonomial::init_3D(), init_3D(), InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), InfFE< Dim, T_radial, T_map >::reinit(), and REINIT_ERROR().

{ return _points;  }
 

std::vector<Point>& 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 133 of file quadrature.h.

References QBase::_points.

{ return _points;  }
 

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

a std::vector containing the quadrature weights.

Definition at line 138 of file quadrature.h.

References QBase::_weights.

Referenced by FE< Dim, T >::edge_reinit(), QClough::init_1D(), QMonomial::init_2D(), init_2D(), QClough::init_2D(), QMonomial::init_3D(), init_3D(), InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), and REINIT_ERROR().

{ return _weights; }
 

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

a std::vector containing the quadrature weights.

Definition at line 143 of file quadrature.h.

References QBase::_weights.

{ return _weights; }
 

void 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 149 of file reference_counter.h.

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

Referenced by ReferenceCountedObject< Value >::ReferenceCountedObject().

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

  p.first++;
}
 

void 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 167 of file reference_counter.h.

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

Referenced by ReferenceCountedObject< Value >::~ReferenceCountedObject().

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

  p.second++;
}
 

void QBase::init (const ElemType_type = INVALID_ELEM, unsigned intp_level = 0) [inherited]Initializes the data structures to contain a quadrature rule for an object of type type.

Definition at line 26 of file quadrature.C.

References QBase::init_0D(), QBase::init_1D(), and QBase::init_2D().

Referenced by FE< Dim, T >::edge_reinit(), QClough::init_1D(), QTrap::init_2D(), QSimpson::init_2D(), QMonomial::init_2D(), QGrid::init_2D(), init_2D(), QClough::init_2D(), QTrap::init_3D(), QSimpson::init_3D(), QMonomial::init_3D(), QGrid::init_3D(), init_3D(), InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), QGauss(), QJacobi::QJacobi(), QSimpson::QSimpson(), QTrap::QTrap(), InfFE< Dim, T_radial, T_map >::reinit(), and REINIT_ERROR().

{
  // 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();
    }
}
 

void QBase::init_0D (const ElemType_type = INVALID_ELEM, unsigned intp_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 70 of file quadrature.C.

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

Referenced by QBase::init().

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

void QGauss::init_1D (const ElemType_type = INVALID_ELEM, unsigned intp_level = 0) [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 QBase.

Definition at line 29 of file quadrature_gauss_1D.C.

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

{
  //----------------------------------------------------------------------
  // 1D quadrature rules
  switch(_order + 2*p)
    {
    case CONSTANT:
    case FIRST:
      {
        _points.resize (1);
        _weights.resize(1);
        
        _points[0](0)  = 0.;
        
        _weights[0]    = 2.;
        
        return;
      }
    case SECOND:
    case THIRD:
      {
        _points.resize (2);
        _weights.resize(2);
              
        _points[0](0) = -5.7735026918962576450914878050196e-01L; // -sqrt(3)/3
        _points[1]    = -_points[0];

        _weights[0]   = 1.;
        _weights[1]   = _weights[0];

        return;
      }
    case FOURTH:
    case FIFTH:
      {
        _points.resize (3);
        _weights.resize(3);
              
        _points[ 0](0) = -7.7459666924148337703585307995648e-01L;
        _points[ 1](0) = 0.;
        _points[ 2]    = -_points[0];

        _weights[ 0]   = 5.5555555555555555555555555555556e-01L;
        _weights[ 1]   = 8.8888888888888888888888888888889e-01L;
        _weights[ 2]   = _weights[0];

        return;
      }
    case SIXTH:
    case SEVENTH:
      {
        _points.resize (4);
        _weights.resize(4);
              
        _points[ 0](0) = -8.6113631159405257522394648889281e-01L;
        _points[ 1](0) = -3.3998104358485626480266575910324e-01L;
        _points[ 2]    = -_points[1];
        _points[ 3]    = -_points[0];

        _weights[ 0]   = 3.4785484513745385737306394922200e-01L;
        _weights[ 1]   = 6.5214515486254614262693605077800e-01L;
        _weights[ 2]   = _weights[1];
        _weights[ 3]   = _weights[0];

        return;
      }
    case EIGHTH:
    case NINTH:
      {
        _points.resize (5);
        _weights.resize(5);
              
        _points[ 0](0) = -9.0617984593866399279762687829939e-01L;
        _points[ 1](0) = -5.3846931010568309103631442070021e-01L;
        _points[ 2](0) = 0.;
        _points[ 3]    = -_points[1];
        _points[ 4]    = -_points[0];

        _weights[ 0]   = 2.3692688505618908751426404071992e-01L;
        _weights[ 1]   = 4.7862867049936646804129151483564e-01L;
        _weights[ 2]   = 5.6888888888888888888888888888889e-01L;
        _weights[ 3]   = _weights[1];
        _weights[ 4]   = _weights[0];

        return;
      }
    case TENTH:
    case ELEVENTH:
      {
        _points.resize (6);
        _weights.resize(6);
              
        _points[ 0](0) = -9.3246951420315202781230155449399e-01L;
        _points[ 1](0) = -6.6120938646626451366139959501991e-01L;
        _points[ 2](0) = -2.3861918608319690863050172168071e-01L;
        _points[ 3]    = -_points[2];
        _points[ 4]    = -_points[1];
        _points[ 5]    = -_points[0];

        _weights[ 0]   = 1.7132449237917034504029614217273e-01L;
        _weights[ 1]   = 3.6076157304813860756983351383772e-01L;
        _weights[ 2]   = 4.6791393457269104738987034398955e-01L;
        _weights[ 3]   = _weights[2];
        _weights[ 4]   = _weights[1];
        _weights[ 5]   = _weights[0];

        return;
      }
    case TWELFTH:
    case THIRTEENTH:
      {
        _points.resize (7);
        _weights.resize(7);
              
        _points[ 0](0) = -9.4910791234275852452618968404785e-01L;
        _points[ 1](0) = -7.4153118559939443986386477328079e-01L;
        _points[ 2](0) = -4.0584515137739716690660641207696e-01L;
        _points[ 3](0) = 0.;
        _points[ 4]    = -_points[2];
        _points[ 5]    = -_points[1];
        _points[ 6]    = -_points[0];

        _weights[ 0]   = 1.2948496616886969327061143267908e-01L;
        _weights[ 1]   = 2.7970539148927666790146777142378e-01L;
        _weights[ 2]   = 3.8183005050511894495036977548898e-01L;
        _weights[ 3]   = 4.1795918367346938775510204081633e-01L;
        _weights[ 4]   = _weights[2];
        _weights[ 5]   = _weights[1];
        _weights[ 6]   = _weights[0];

        return;
      }
    case FOURTEENTH:
    case FIFTEENTH:
      {
        _points.resize (8);
        _weights.resize(8);
              
        _points[ 0](0) = -9.6028985649753623168356086856947e-01L;
        _points[ 1](0) = -7.9666647741362673959155393647583e-01L;
        _points[ 2](0) = -5.2553240991632898581773904918925e-01L;
        _points[ 3](0) = -1.8343464249564980493947614236018e-01L;
        _points[ 4]    = -_points[3];
        _points[ 5]    = -_points[2];
        _points[ 6]    = -_points[1];
        _points[ 7]    = -_points[0];

        _weights[ 0]   = 1.0122853629037625915253135430996e-01L;
        _weights[ 1]   = 2.2238103445337447054435599442624e-01L;
        _weights[ 2]   = 3.1370664587788728733796220198660e-01L;
        _weights[ 3]   = 3.6268378337836198296515044927720e-01L;
        _weights[ 4]   = _weights[3];
        _weights[ 5]   = _weights[2];
        _weights[ 6]   = _weights[1];
        _weights[ 7]   = _weights[0];

        return;
      }
    case SIXTEENTH:
    case SEVENTEENTH:
      {
        _points.resize (9);
        _weights.resize(9);

        _points[ 0](0) = -9.6816023950762608983557620290367e-01L;
        _points[ 1](0) = -8.3603110732663579429942978806973e-01L;
        _points[ 2](0) = -6.1337143270059039730870203934147e-01L;
        _points[ 3](0) = -3.2425342340380892903853801464334e-01L;
        _points[ 4](0) = 0.;
        _points[ 5]    = -_points[3];
        _points[ 6]    = -_points[2];
        _points[ 7]    = -_points[1];
        _points[ 8]    = -_points[0];

        _weights[ 0]   = 8.1274388361574411971892158110524e-02L;
        _weights[ 1]   = 1.8064816069485740405847203124291e-01L;
        _weights[ 2]   = 2.6061069640293546231874286941863e-01L;
        _weights[ 3]   = 3.1234707704000284006863040658444e-01L;
        _weights[ 4]   = 3.3023935500125976316452506928697e-01L;
        _weights[ 5]   = _weights[3];
        _weights[ 6]   = _weights[2];
        _weights[ 7]   = _weights[1];
        _weights[ 8]   = _weights[0];

        return;
      }
    case EIGHTTEENTH:
    case NINTEENTH:
      {
        _points.resize (10);
        _weights.resize(10);

        _points[ 0](0) = -9.7390652851717172007796401208445e-01L;
        _points[ 1](0) = -8.6506336668898451073209668842349e-01L;
        _points[ 2](0) = -6.7940956829902440623432736511487e-01L;
        _points[ 3](0) = -4.3339539412924719079926594316578e-01L;
        _points[ 4](0) = -1.4887433898163121088482600112972e-01L;
        _points[ 5]    = -_points[4];
        _points[ 6]    = -_points[3];
        _points[ 7]    = -_points[2];
        _points[ 8]    = -_points[1];
        _points[ 9]    = -_points[0];

        _weights[ 0]   = 6.6671344308688137593568809893332e-02L;
        _weights[ 1]   = 1.4945134915058059314577633965770e-01L;
        _weights[ 2]   = 2.1908636251598204399553493422816e-01L;
        _weights[ 3]   = 2.6926671930999635509122692156947e-01L;
        _weights[ 4]   = 2.9552422471475287017389299465134e-01L;
        _weights[ 5]   = _weights[4];
        _weights[ 6]   = _weights[3];
        _weights[ 7]   = _weights[2];
        _weights[ 8]   = _weights[1];
        _weights[ 9]   = _weights[0];

        return;
      }      

    case TWENTIETH:
    case TWENTYFIRST:
      {
        _points.resize (11);
        _weights.resize(11);

        _points[ 0](0) = -9.7822865814605699280393800112286e-01L;
        _points[ 1](0) = -8.8706259976809529907515776930393e-01L;
        _points[ 2](0) = -7.3015200557404932409341625203115e-01L;
        _points[ 3](0) = -5.1909612920681181592572566945861e-01L;
        _points[ 4](0) = -2.6954315595234497233153198540086e-01L;
        _points[ 5](0) = 0.;
        _points[ 6]    = -_points[4];
        _points[ 7]    = -_points[3];
        _points[ 8]    = -_points[2];
        _points[ 9]    = -_points[1];
        _points[10]    = -_points[0];

        _weights[ 0]   = 5.5668567116173666482753720442549e-02L;
        _weights[ 1]   = 1.2558036946490462463469429922394e-01L;
        _weights[ 2]   = 1.8629021092773425142609764143166e-01L;
        _weights[ 3]   = 2.3319376459199047991852370484318e-01L;
        _weights[ 4]   = 2.6280454451024666218068886989051e-01L;
        _weights[ 5]   = 2.7292508677790063071448352833634e-01L;
        _weights[ 6]   = _weights[4];
        _weights[ 7]   = _weights[3];
        _weights[ 8]   = _weights[2];
        _weights[ 9]   = _weights[1];
        _weights[10]   = _weights[0];

        return;
      }

    case TWENTYSECOND:
    case TWENTYTHIRD:
      {
        _points.resize (12);
        _weights.resize(12);

        _points[ 0](0) = -9.8156063424671925069054909014928e-01L;
        _points[ 1](0) = -9.0411725637047485667846586611910e-01L;
        _points[ 2](0) = -7.6990267419430468703689383321282e-01L;
        _points[ 3](0) = -5.8731795428661744729670241894053e-01L;
        _points[ 4](0) = -3.6783149899818019375269153664372e-01L;
        _points[ 5](0) = -1.2523340851146891547244136946385e-01L;
        _points[ 6]    = -_points[5];
        _points[ 7]    = -_points[4];
        _points[ 8]    = -_points[3];
        _points[ 9]    = -_points[2];
        _points[10]    = -_points[1];
        _points[11]    = -_points[0];

        _weights[ 0]   = 4.7175336386511827194615961485017e-02L;
        _weights[ 1]   = 1.0693932599531843096025471819400e-01L;
        _weights[ 2]   = 1.6007832854334622633465252954336e-01L;
        _weights[ 3]   = 2.0316742672306592174906445580980e-01L;
        _weights[ 4]   = 2.3349253653835480876084989892483e-01L;
        _weights[ 5]   = 2.4914704581340278500056243604295e-01L;
        _weights[ 6]   = _weights[5];
        _weights[ 7]   = _weights[4];
        _weights[ 8]   = _weights[3];
        _weights[ 9]   = _weights[2];
        _weights[10]   = _weights[1];
        _weights[11]   = _weights[0];

        return;
      }      

    case TWENTYFOURTH:
    case TWENTYFIFTH:
      {
        _points.resize (13);
        _weights.resize(13);

        _points[ 0](0) = -9.8418305471858814947282944880711e-01L;
        _points[ 1](0) = -9.1759839922297796520654783650072e-01L;
        _points[ 2](0) = -8.0157809073330991279420648958286e-01L;
        _points[ 3](0) = -6.4234933944034022064398460699552e-01L;
        _points[ 4](0) = -4.4849275103644685287791285212764e-01L;
        _points[ 5](0) = -2.3045831595513479406552812109799e-01L;
        _points[ 6](0) = 0.;
        _points[ 7]    = -_points[5];
        _points[ 8]    = -_points[4];
        _points[ 9]    = -_points[3];
        _points[10]    = -_points[2];
        _points[11]    = -_points[1];
        _points[12]    = -_points[0];

        _weights[ 0]   = 4.0484004765315879520021592200986e-02L;
        _weights[ 1]   = 9.2121499837728447914421775953797e-02L;
        _weights[ 2]   = 1.3887351021978723846360177686887e-01L;
        _weights[ 3]   = 1.7814598076194573828004669199610e-01L;
        _weights[ 4]   = 2.0781604753688850231252321930605e-01L;
        _weights[ 5]   = 2.2628318026289723841209018603978e-01L;
        _weights[ 6]   = 2.3255155323087391019458951526884e-01L;
        _weights[ 7]   = _weights[5];
        _weights[ 8]   = _weights[4];
        _weights[ 9]   = _weights[3];
        _weights[10]   = _weights[2];
        _weights[11]   = _weights[1];
        _weights[12]   = _weights[0];

        return;
      }

    case TWENTYSIXTH:
    case TWENTYSEVENTH:
      {
        _points.resize (14);
        _weights.resize(14);

        _points[ 0](0) = -9.8628380869681233884159726670405e-01L;
        _points[ 1](0) = -9.2843488366357351733639113937787e-01L;
        _points[ 2](0) = -8.2720131506976499318979474265039e-01L;
        _points[ 3](0) = -6.8729290481168547014801980301933e-01L;
        _points[ 4](0) = -5.1524863635815409196529071855119e-01L;
        _points[ 5](0) = -3.1911236892788976043567182416848e-01L;
        _points[ 6](0) = -1.0805494870734366206624465021983e-01L;
        _points[ 7]    = -_points[6];
        _points[ 8]    = -_points[5];
        _points[ 9]    = -_points[4];
        _points[10]    = -_points[3];
        _points[11]    = -_points[2];
        _points[12]    = -_points[1];
        _points[13]    = -_points[0];

        _weights[ 0]   = 3.5119460331751863031832876138192e-02L;
        _weights[ 1]   = 8.0158087159760209805633277062854e-02L;
        _weights[ 2]   = 1.2151857068790318468941480907248e-01L;
        _weights[ 3]   = 1.5720316715819353456960193862384e-01L;
        _weights[ 4]   = 1.8553839747793781374171659012516e-01L;
        _weights[ 5]   = 2.0519846372129560396592406566122e-01L;
        _weights[ 6]   = 2.1526385346315779019587644331626e-01L;
        _weights[ 7]   = _weights[6];
        _weights[ 8]   = _weights[5];
        _weights[ 9]   = _weights[4];
        _weights[10]   = _weights[3];
        _weights[11]   = _weights[2];
        _weights[12]   = _weights[1];
        _weights[13]   = _weights[0];

        return;
      }

    case TWENTYEIGHTH:
    case TWENTYNINTH:
      {
        _points.resize (15);
        _weights.resize(15);

        _points[ 0](0) = -9.8799251802048542848956571858661e-01L;
        _points[ 1](0) = -9.3727339240070590430775894771021e-01L;
        _points[ 2](0) = -8.4820658341042721620064832077422e-01L;
        _points[ 3](0) = -7.2441773136017004741618605461394e-01L;
        _points[ 4](0) = -5.7097217260853884753722673725391e-01L;
        _points[ 5](0) = -3.9415134707756336989720737098105e-01L;
        _points[ 6](0) = -2.0119409399743452230062830339460e-01L;
        _points[ 7](0) = 0.;
        _points[ 8]    = -_points[6];
        _points[ 9]    = -_points[5];
        _points[10]    = -_points[4];
        _points[11]    = -_points[3];
        _points[12]    = -_points[2];
        _points[13]    = -_points[1];
        _points[14]    = -_points[0];

        _weights[ 0]   = 3.0753241996117268354628393577204e-02L;
        _weights[ 1]   = 7.0366047488108124709267416450667e-02L;
        _weights[ 2]   = 1.0715922046717193501186954668587e-01L;
        _weights[ 3]   = 1.3957067792615431444780479451103e-01L;
        _weights[ 4]   = 1.6626920581699393355320086048121e-01L;
        _weights[ 5]   = 1.8616100001556221102680056186642e-01L;
        _weights[ 6]   = 1.9843148532711157645611832644384e-01L;
        _weights[ 7]   = 2.0257824192556127288062019996752e-01L;
        _weights[ 8]   = _weights[6];
        _weights[ 9]   = _weights[5];
        _weights[10]   = _weights[4];
        _weights[11]   = _weights[3];
        _weights[12]   = _weights[2];
        _weights[13]   = _weights[1];
        _weights[14]   = _weights[0];

        return;
      }

    case THIRTIETH:
    case THIRTYFIRST:
      {
        _points.resize (16);
        _weights.resize(16);

        _points[ 0](0) = -9.8940093499164993259615417345033e-01L;
        _points[ 1](0) = -9.4457502307323257607798841553461e-01L;
        _points[ 2](0) = -8.6563120238783174388046789771239e-01L;
        _points[ 3](0) = -7.5540440835500303389510119484744e-01L;
        _points[ 4](0) = -6.1787624440264374844667176404879e-01L;
        _points[ 5](0) = -4.5801677765722738634241944298358e-01L;
        _points[ 6](0) = -2.8160355077925891323046050146050e-01L;
        _points[ 7](0) = -9.5012509837637440185319335424958e-02L;
        _points[ 8]    = -_points[7];
        _points[ 9]    = -_points[6];
        _points[10]    = -_points[5];
        _points[11]    = -_points[4];
        _points[12]    = -_points[3];
        _points[13]    = -_points[2];
        _points[14]    = -_points[1];
        _points[15]    = -_points[0];

        _weights[ 0]   = 2.7152459411754094851780572456018e-02L;
        _weights[ 1]   = 6.2253523938647892862843836994378e-02L;
        _weights[ 2]   = 9.5158511682492784809925107602246e-02L;
        _weights[ 3]   = 1.2462897125553387205247628219202e-01L;
        _weights[ 4]   = 1.4959598881657673208150173054748e-01L;
        _weights[ 5]   = 1.6915651939500253818931207903033e-01L;
        _weights[ 6]   = 1.8260341504492358886676366796922e-01L;
        _weights[ 7]   = 1.8945061045506849628539672320828e-01L;
        _weights[ 8]   = _weights[7];
        _weights[ 9]   = _weights[6];
        _weights[10]   = _weights[5];
        _weights[11]   = _weights[4];
        _weights[12]   = _weights[3];
        _weights[13]   = _weights[2];
        _weights[14]   = _weights[1];
        _weights[15]   = _weights[0];

        return;
      }

    case THIRTYSECOND:
    case THIRTYTHIRD:
      {
        _points.resize (17);
        _weights.resize(17);

        _points[ 0](0) = -9.9057547531441733567543401994067e-01L;
        _points[ 1](0) = -9.5067552176876776122271695789580e-01L;
        _points[ 2](0) = -8.8023915372698590212295569448816e-01L;
        _points[ 3](0) = -7.8151400389680140692523005552048e-01L;
        _points[ 4](0) = -6.5767115921669076585030221664300e-01L;
        _points[ 5](0) = -5.1269053708647696788624656862955e-01L;
        _points[ 6](0) = -3.5123176345387631529718551709535e-01L;
        _points[ 7](0) = -1.7848418149584785585067749365407e-01L;
        _points[ 8](0) = 0.;
        _points[ 9]    = -_points[7];
        _points[10]    = -_points[6];
        _points[11]    = -_points[5];
        _points[12]    = -_points[4];
        _points[13]    = -_points[3];
        _points[14]    = -_points[2];
        _points[15]    = -_points[1];
        _points[16]    = -_points[0];

        _weights[ 0]   = 2.4148302868547931960110026287565e-02L;
        _weights[ 1]   = 5.5459529373987201129440165358245e-02L;
        _weights[ 2]   = 8.5036148317179180883535370191062e-02L;
        _weights[ 3]   = 1.1188384719340397109478838562636e-01L;
        _weights[ 4]   = 1.3513636846852547328631998170235e-01L;
        _weights[ 5]   = 1.5404576107681028808143159480196e-01L;
        _weights[ 6]   = 1.6800410215645004450997066378832e-01L;
        _weights[ 7]   = 1.7656270536699264632527099011320e-01L;
        _weights[ 8]   = 1.7944647035620652545826564426189e-01L;
        _weights[ 9]   = _weights[7];
        _weights[10]   = _weights[6];
        _weights[11]   = _weights[5];
        _weights[12]   = _weights[4];
        _weights[13]   = _weights[3];
        _weights[14]   = _weights[2];
        _weights[15]   = _weights[1];
        _weights[16]   = _weights[0];

        return;
      }

    case THIRTYFOURTH:
    case THIRTYFIFTH:
      {
        _points.resize (18);
        _weights.resize(18);

        _points[ 0](0) = -9.9156516842093094673001600470615e-01L;
        _points[ 1](0) = -9.5582394957139775518119589292978e-01L;
        _points[ 2](0) = -8.9260246649755573920606059112715e-01L;
        _points[ 3](0) = -8.0370495897252311568241745501459e-01L;
        _points[ 4](0) = -6.9168704306035320787489108128885e-01L;
        _points[ 5](0) = -5.5977083107394753460787154852533e-01L;
        _points[ 6](0) = -4.1175116146284264603593179383305e-01L;
        _points[ 7](0) = -2.5188622569150550958897285487791e-01L;
        _points[ 8](0) = -8.4775013041735301242261852935784e-02L;
        _points[ 9]    = -_points[8];
        _points[10]    = -_points[7];
        _points[11]    = -_points[6];
        _points[12]    = -_points[5];
        _points[13]    = -_points[4];
        _points[14]    = -_points[3];
        _points[15]    = -_points[2];
        _points[16]    = -_points[1];
        _points[17]    = -_points[0];

        _weights[ 0]   = 2.1616013526483310313342710266452e-02L;
        _weights[ 1]   = 4.9714548894969796453334946202639e-02L;
        _weights[ 2]   = 7.6425730254889056529129677616637e-02L;
        _weights[ 3]   = 1.0094204410628716556281398492483e-01L;
        _weights[ 4]   = 1.2255520671147846018451912680020e-01L;
        _weights[ 5]   = 1.4064291467065065120473130375195e-01L;
        _weights[ 6]   = 1.5468467512626524492541800383637e-01L;
        _weights[ 7]   = 1.6427648374583272298605377646593e-01L;
        _weights[ 8]   = 1.6914238296314359184065647013499e-01L;
        _weights[ 9]   = _weights[8];
        _weights[10]   = _weights[7];
        _weights[11]   = _weights[6];
        _weights[12]   = _weights[5];
        _weights[13]   = _weights[4];
        _weights[14]   = _weights[3];
        _weights[15]   = _weights[2];
        _weights[16]   = _weights[1];
        _weights[17]   = _weights[0];

        return;
      }

    case THIRTYSIXTH:
    case THIRTYSEVENTH:
      {
        _points.resize (19);
        _weights.resize(19);

        _points[ 0](0) = -9.9240684384358440318901767025326e-01L;
        _points[ 1](0) = -9.6020815213483003085277884068765e-01L;
        _points[ 2](0) = -9.0315590361481790164266092853231e-01L;
        _points[ 3](0) = -8.2271465653714282497892248671271e-01L;
        _points[ 4](0) = -7.2096617733522937861709586082378e-01L;
        _points[ 5](0) = -6.0054530466168102346963816494624e-01L;
        _points[ 6](0) = -4.6457074137596094571726714810410e-01L;
        _points[ 7](0) = -3.1656409996362983199011732884984e-01L;
        _points[ 8](0) = -1.6035864564022537586809611574074e-01L;
        _points[ 9](0) = 0.;
        _points[10]    = -_points[8];
        _points[11]    = -_points[7];
        _points[12]    = -_points[6];
        _points[13]    = -_points[5];
        _points[14]    = -_points[4];
        _points[15]    = -_points[3];
        _points[16]    = -_points[2];
        _points[17]    = -_points[1];
        _points[18]    = -_points[0];

        _weights[ 0]   = 1.9461788229726477036312041464438e-02L;
        _weights[ 1]   = 4.4814226765699600332838157401994e-02L;
        _weights[ 2]   = 6.9044542737641226580708258006013e-02L;
        _weights[ 3]   = 9.1490021622449999464462094123840e-02L;
        _weights[ 4]   = 1.1156664554733399471602390168177e-01L;
        _weights[ 5]   = 1.2875396253933622767551578485688e-01L;
        _weights[ 6]   = 1.4260670217360661177574610944190e-01L;
        _weights[ 7]   = 1.5276604206585966677885540089766e-01L;
        _weights[ 8]   = 1.5896884339395434764995643946505e-01L;
        _weights[ 9]   = 1.6105444984878369597916362532092e-01L;
        _weights[10]   = _weights[8];
        _weights[11]   = _weights[7];
        _weights[12]   = _weights[6];
        _weights[13]   = _weights[5];
        _weights[14]   = _weights[4];
        _weights[15]   = _weights[3];
        _weights[16]   = _weights[2];
        _weights[17]   = _weights[1];
        _weights[18]   = _weights[0];

        return;
      }

    case THIRTYEIGHTH:
    case THIRTYNINTH:
      {
        _points.resize (20);
        _weights.resize(20);

        _points[ 0](0) = -9.9312859918509492478612238847132e-01L;
        _points[ 1](0) = -9.6397192727791379126766613119728e-01L;
        _points[ 2](0) = -9.1223442825132590586775244120330e-01L;
        _points[ 3](0) = -8.3911697182221882339452906170152e-01L;
        _points[ 4](0) = -7.4633190646015079261430507035564e-01L;
        _points[ 5](0) = -6.3605368072651502545283669622629e-01L;
        _points[ 6](0) = -5.1086700195082709800436405095525e-01L;
        _points[ 7](0) = -3.7370608871541956067254817702493e-01L;
        _points[ 8](0) = -2.2778585114164507808049619536857e-01L;
        _points[ 9](0) = -7.6526521133497333754640409398838e-02L;
        _points[10]    = -_points[9];
        _points[11]    = -_points[8];
        _points[12]    = -_points[7];
        _points[13]    = -_points[6];
        _points[14]    = -_points[5];
        _points[15]    = -_points[4];
        _points[16]    = -_points[3];
        _points[17]    = -_points[2];
        _points[18]    = -_points[1];
        _points[19]    = -_points[0];

        _weights[ 0]   = 1.7614007139152118311861962351853e-02L;
        _weights[ 1]   = 4.0601429800386941331039952274932e-02L;
        _weights[ 2]   = 6.2672048334109063569506535187042e-02L;
        _weights[ 3]   = 8.3276741576704748724758143222046e-02L;
        _weights[ 4]   = 1.0193011981724043503675013548035e-01L;
        _weights[ 5]   = 1.1819453196151841731237737771138e-01L;
        _weights[ 6]   = 1.3168863844917662689849449974816e-01L;
        _weights[ 7]   = 1.4209610931838205132929832506716e-01L;
        _weights[ 8]   = 1.4917298647260374678782873700197e-01L;
        _weights[ 9]   = 1.5275338713072585069808433195510e-01L;
        _weights[10]   = _weights[9];
        _weights[11]   = _weights[8];
        _weights[12]   = _weights[7];
        _weights[13]   = _weights[6];
        _weights[14]   = _weights[5];
        _weights[15]   = _weights[4];
        _weights[16]   = _weights[3];
        _weights[17]   = _weights[2];
        _weights[18]   = _weights[1];
        _weights[19]   = _weights[0];

        return;
      }

    case FORTIETH:
    case FORTYFIRST:
      {
        _points.resize (21);
        _weights.resize(21);

        _points[ 0](0) = -9.9375217062038950026024203593794e-01L;
        _points[ 1](0) = -9.6722683856630629431662221490770e-01L;
        _points[ 2](0) = -9.2009933415040082879018713371497e-01L;
        _points[ 3](0) = -8.5336336458331728364725063858757e-01L;
        _points[ 4](0) = -7.6843996347567790861587785130623e-01L;
        _points[ 5](0) = -6.6713880419741231930596666999034e-01L;
        _points[ 6](0) = -5.5161883588721980705901879672431e-01L;
        _points[ 7](0) = -4.2434212020743878357366888854379e-01L;
        _points[ 8](0) = -2.8802131680240109660079251606460e-01L;
        _points[ 9](0) = -1.4556185416089509093703098233869e-01L;
        _points[10](0) = 0.;
        _points[11]    = -_points[9];
        _points[12]    = -_points[8];
        _points[13]    = -_points[7];
        _points[14]    = -_points[6];
        _points[15]    = -_points[5];
        _points[16]    = -_points[4];
        _points[17]    = -_points[3];
        _points[18]    = -_points[2];
        _points[19]    = -_points[1];
        _points[20]    = -_points[0];

        _weights[ 0]   = 1.6017228257774333324224616858471e-02L;
        _weights[ 1]   = 3.6953789770852493799950668299330e-02L;
        _weights[ 2]   = 5.7134425426857208283635826472448e-02L;
        _weights[ 3]   = 7.6100113628379302017051653300183e-02L;
        _weights[ 4]   = 9.3444423456033861553289741113932e-02L;
        _weights[ 5]   = 1.0879729916714837766347457807011e-01L;
        _weights[ 6]   = 1.2183141605372853419536717712572e-01L;
        _weights[ 7]   = 1.3226893863333746178105257449678e-01L;
        _weights[ 8]   = 1.3988739479107315472213342386758e-01L;
        _weights[ 9]   = 1.4452440398997005906382716655375e-01L;
        _weights[10]   = 1.4608113364969042719198514768337e-01L;
        _weights[11]   = _weights[9];
        _weights[12]   = _weights[8];
        _weights[13]   = _weights[7];
        _weights[14]   = _weights[6];
        _weights[15]   = _weights[5];
        _weights[16]   = _weights[4];
        _weights[17]   = _weights[3];
        _weights[18]   = _weights[2];
        _weights[19]   = _weights[1];
        _weights[20]   = _weights[0];

        return;
      }

    case FORTYSECOND:
    case FORTYTHIRD:
      {
        _points.resize (22);
        _weights.resize(22);

        _points[ 0](0) = -9.9429458548239929207303142116130e-01L;
        _points[ 1](0) = -9.7006049783542872712395098676527e-01L;
        _points[ 2](0) = -9.2695677218717400052069293925905e-01L;
        _points[ 3](0) = -8.6581257772030013653642563701938e-01L;
        _points[ 4](0) = -7.8781680597920816200427795540835e-01L;
        _points[ 5](0) = -6.9448726318668278005068983576226e-01L;
        _points[ 6](0) = -5.8764040350691159295887692763865e-01L;
        _points[ 7](0) = -4.6935583798675702640633071096641e-01L;
        _points[ 8](0) = -3.4193582089208422515814742042738e-01L;
        _points[ 9](0) = -2.0786042668822128547884653391955e-01L;
        _points[10](0) = -6.9739273319722221213841796118628e-02L;
        _points[11]    = -_points[10];
        _points[12]    = -_points[9];
        _points[13]    = -_points[8];
        _points[14]    = -_points[7];
        _points[15]    = -_points[6];
        _points[16]    = -_points[5];
        _points[17]    = -_points[4];
        _points[18]    = -_points[3];
        _points[19]    = -_points[2];
        _points[20]    = -_points[1];
        _points[21]    = -_points[0];

        _weights[ 0]   = 1.4627995298272200684991098047185e-02L;
        _weights[ 1]   = 3.3774901584814154793302246865913e-02L;
        _weights[ 2]   = 5.2293335152683285940312051273211e-02L;
        _weights[ 3]   = 6.9796468424520488094961418930218e-02L;
        _weights[ 4]   = 8.5941606217067727414443681372703e-02L;
        _weights[ 5]   = 1.0041414444288096493207883783054e-01L;
        _weights[ 6]   = 1.1293229608053921839340060742175e-01L;
        _weights[ 7]   = 1.2325237681051242428556098615481e-01L;
        _weights[ 8]   = 1.3117350478706237073296499253031e-01L;
        _weights[ 9]   = 1.3654149834601517135257383123152e-01L;
        _weights[10]   = 1.3925187285563199337541024834181e-01L;
        _weights[11]   = _weights[10];
        _weights[12]   = _weights[9];
        _weights[13]   = _weights[8];
        _weights[14]   = _weights[7];
        _weights[15]   = _weights[6];
        _weights[16]   = _weights[5];
        _weights[17]   = _weights[4];
        _weights[18]   = _weights[3];
        _weights[19]   = _weights[2];
        _weights[20]   = _weights[1];
        _weights[21]   = _weights[0];

        return;
      }


    default:
      {
        std::cerr << 'Quadrature rule ' << _order
                  << ' not supported!' << std::endl;
              
        libmesh_error();
      }
    }



  return;
}
 

void QGauss::init_2D (const ElemType = INVALID_ELEM, unsigned int = 0) [private, virtual]Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate values. The order of the rule will be defined by the implementing class. Should not be pure virtual since a derived quadrature rule may only be defined in 1D. If not redefined, gives an error (when DEBUG defined) when called.

Reimplemented from QBase.

Definition at line 27 of file quadrature_gauss_2D.C.

References QBase::_points, QBase::_weights, libMeshEnums::CONSTANT, dunavant_rule(), dunavant_rule2(), libMeshEnums::EDGE2, libMeshEnums::EIGHTH, libMeshEnums::EIGHTTEENTH, libMeshEnums::ELEVENTH, libMeshEnums::FIFTEENTH, libMeshEnums::FIFTH, libMeshEnums::FIRST, libMeshEnums::FOURTEENTH, libMeshEnums::FOURTH, QBase::get_points(), QBase::get_weights(), QBase::init(), libMeshEnums::NINTEENTH, libMeshEnums::NINTH, libMeshEnums::QUAD4, libMeshEnums::QUAD8, libMeshEnums::QUAD9, libMeshEnums::SECOND, libMeshEnums::SEVENTEENTH, libMeshEnums::SEVENTH, libMeshEnums::SIXTEENTH, libMeshEnums::SIXTH, libMeshEnums::TENTH, libMeshEnums::THIRD, libMeshEnums::THIRTEENTH, libMeshEnums::THIRTIETH, libMeshEnums::TRI3, libMeshEnums::TRI6, libMeshEnums::TWELFTH, libMeshEnums::TWENTIETH, libMeshEnums::TWENTYEIGHTH, libMeshEnums::TWENTYFIFTH, libMeshEnums::TWENTYFOURTH, libMeshEnums::TWENTYNINTH, libMeshEnums::TWENTYSECOND, libMeshEnums::TWENTYSEVENTH, libMeshEnums::TWENTYSIXTH, and libMeshEnums::TWENTYTHIRD.

{
#if LIBMESH_DIM > 1
  
  //-----------------------------------------------------------------------
  // 2D quadrature rules
  switch (_type)
    {


      //---------------------------------------------
      // Quadrilateral quadrature rules
    case QUAD4:
    case QUAD8:
    case QUAD9:
      {
        // We compute the 2D quadrature rule as a tensor
        // product of the 1D quadrature rule.
        //
        // For QUADs, a quadrature rule of order 'p' must be able to integrate
        // bilinear (p=1), biquadratic (p=2), bicubic (p=3), etc. polynomials of the form
        //
        // (x^p + x^{p-1} + ... + 1) * (y^p + y^{p-1} + ... + 1)
        //
        // These polynomials have terms *up to* degree 2p but they are *not* complete
        // polynomials of degree 2p. For example, when p=2 we have
        //        1
        //     x      y
        // x^2    xy     y^2
        //    yx^2   xy^2
        //       x^2y^2
        QGauss q1D(1,_order);
        q1D.init(EDGE2,p);
        tensor_product_quad( q1D );
        return;
      }

            
      //---------------------------------------------
      // Triangle quadrature rules
    case TRI3:
    case TRI6:
      {
        switch(_order + 2*p)
          {
          case CONSTANT:
          case FIRST:
            {
              // Exact for linears
              _points.resize(1);
              _weights.resize(1);
                  
              _points[0](0) = 1.0L/3.0L;
              _points[0](1) = 1.0L/3.0L;

              _weights[0] = 0.5;

              return;
            }
          case SECOND:
            {
              // Exact for quadratics
              _points.resize(3);
              _weights.resize(3);

              // Alternate rule with points on ref. elt. boundaries.
              // Not ideal for problems with material coefficient discontinuities
              // aligned along element boundaries.
              // _points[0](0) = .5;
              // _points[0](1) = .5;
              // _points[1](0) = 0.;
              // _points[1](1) = .5;
              // _points[2](0) = .5;
              // _points[2](1) = .0;
              
              _points[0](0) = 2.0L/3.0L;
              _points[0](1) = 1.0L/6.0L;

              _points[1](0) = 1.0L/6.0L;
              _points[1](1) = 2.0L/3.0L;

              _points[2](0) = 1.0L/6.0L;
              _points[2](1) = 1.0L/6.0L;


              _weights[0] = 1.0L/6.0L;
              _weights[1] = 1.0L/6.0L;
              _weights[2] = 1.0L/6.0L;

              return;
            }
          case THIRD:
            {
              // Exact for cubics
              _points.resize(4);
              _weights.resize(4);

              // This rule is formed from a tensor product of
              // appropriately-scaled Gauss and Jacobi rules.  (See
              // also the QConical quadrature class, this is a
              // hard-coded version of one of those rules.)  For high
              // orders these rules generally have too many points,
              // but at extremely low order they are competitive and
              // have the additional benefit of having all positive
              // weights.
              _points[0](0) = 1.5505102572168219018027159252941e-01L;
              _points[0](1) = 1.7855872826361642311703513337422e-01L;
              _points[1](0) = 6.4494897427831780981972840747059e-01L;
              _points[1](1) = 7.5031110222608118177475598324603e-02L;
              _points[2](0) = 1.5505102572168219018027159252941e-01L;
              _points[2](1) = 6.6639024601470138670269327409637e-01L;
              _points[3](0) = 6.4494897427831780981972840747059e-01L;
              _points[3](1) = 2.8001991549907407200279599420481e-01L;

              _weights[0] = 1.5902069087198858469718450103758e-01L;
              _weights[1] = 9.0979309128011415302815498962418e-02L;
              _weights[2] = 1.5902069087198858469718450103758e-01L;
              _weights[3] = 9.0979309128011415302815498962418e-02L;

              return;

              
              // The following third-order rule is quite commonly cited
              // in the literature and most likely works fine.  However,
              // we generally prefer a rule with all positive weights
              // and an equal number of points, when available.
              //
              //  (allow_rules_with_negative_weights)
              // {
              //   // Exact for cubics
              //   _points.resize(4);
              //   _weights.resize(4);
              //   
              //   _points[0](0) = .33333333333333333333333333333333;
              //   _points[0](1) = .33333333333333333333333333333333;
              // 
              //   _points[1](0) = .2;
              //   _points[1](1) = .6;
              // 
              //   _points[2](0) = .2;
              //   _points[2](1) = .2;
              // 
              //   _points[3](0) = .6;
              //   _points[3](1) = .2;
              // 
              // 
              //   _weights[0] = -27./96.;
              //   _weights[1] =  25./96.;
              //   _weights[2] =  25./96.;
              //   _weights[3] =  25./96.;
              // 
              //   return;
              // } // end if (allow_rules_with_negative_weights)
              // Note: if !allow_rules_with_negative_weights, fall through to next case.
            }


            
            // A degree 4 rule with six points.  This rule can be found in many places
            // including:
            //
            // J.N. Lyness and D. Jespersen, Moderate degree symmetric
            // quadrature rules for the triangle, J. Inst. Math. Appl.  15 (1975),
            // 19--32.
            //
            // We used the code in: 
            // L. Zhang, T. Cui, and H. Liu. 'A set of symmetric quadrature rules
            // on triangles and tetrahedra'  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
            // to generate additional precision.
          case FOURTH:
            {
              const unsigned int n_wts = 2;
              const Real wts[n_wts] =
                {
                  1.1169079483900573284750350421656140e-01L,
                  5.4975871827660933819163162450105264e-02L
                };

              const Real a[n_wts] =
                {
                  4.4594849091596488631832925388305199e-01L,
                  9.1576213509770743459571463402201508e-02L
                };

              const Real b[n_wts] = {0., 0.}; // not used
              const unsigned int permutation_ids[n_wts] = {3, 3};

              dunavant_rule2(wts, a, b, permutation_ids, n_wts); // 6 total points
              
              return;
            }


            
            // Exact for quintics
            // Can be found in 'Quadrature on Simplices of Arbitrary
            // Dimension' by Walkington.
          case FIFTH:
            {
              const unsigned int n_wts = 3;
              const Real wts[n_wts] =
                {
                  9.0L/80.0L,
                  31.0L/480.0L + std::sqrt(15.0L)/2400.0L,
                  31.0L/480.0L - std::sqrt(15.0L)/2400.0L
                };

              const Real a[n_wts] =
                {
                  0., // 'a' parameter not used for origin permutation
                  2.0L/7.0L + std::sqrt(15.0L)/21.0L,
                  2.0L/7.0L - std::sqrt(15.0L)/21.0L
                };

              const Real b[n_wts] = {0., 0., 0.}; // not used
              const unsigned int permutation_ids[n_wts] = {1, 3, 3};

              dunavant_rule2(wts, a, b, permutation_ids, n_wts); // 6 total points

              return;
            }


            
            // A degree 6 rule with 12 points.  This rule can be found in many places
            // including:
            //
            // J.N. Lyness and D. Jespersen, Moderate degree symmetric
            // quadrature rules for the triangle, J. Inst. Math. Appl.  15 (1975),
            // 19--32.
            //
            // We used the code in: 
            // L. Zhang, T. Cui, and H. Liu. 'A set of symmetric quadrature rules
            // on triangles and tetrahedra'  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
            // to generate additional precision.
            //
            // Note that the following 7th-order Ro3-invariant rule also has only 12 points,
            // which technically makes it the superior rule.  This one is here for completeness.
          case SIXTH:
            {
              const unsigned int n_wts = 3;
              const Real wts[n_wts] =
                {
                  5.8393137863189683012644805692789721e-02L,
                  2.5422453185103408460468404553434492e-02L,
                  4.1425537809186787596776728210221227e-02L
                };

              const Real a[n_wts] =
                {
                  2.4928674517091042129163855310701908e-01L,
                  6.3089014491502228340331602870819157e-02L,
                  3.1035245103378440541660773395655215e-01L
                };

              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  6.3650249912139864723014259441204970e-01L               
                };
              
              const unsigned int permutation_ids[n_wts] = {3, 3, 6}; // 12 total points

              dunavant_rule2(wts, a, b, permutation_ids, n_wts);

              return;
            }
            
            
            // A degree 7 rule with 12 points.  This rule can be found in:
            //
            // K. Gatermann, The construction of symmetric cubature
            // formulas for the square and the triangle, Computing 40
            // (1988), 229--240.
            //
            // This rule, which is provably minimal in the number of
            // integration points, is said to be 'Ro3 invariant' which
            // means that a given set of barycentric coordinates
            // (z1,z2,z3) implies the quadrature points (z1,z2),
            // (z3,z1), (z2,z3) which are formed by taking the first
            // two entries in cyclic permutations of the barycentric
            // point.  Barycentric coordinates are related in the
            // sense that: z3 = 1 - z1 - z2.
            //
            // The 12-point sixth-order rule for triangles given in
            // Flaherty's (http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps)
            // lecture notes has been removed in favor of this rule
            // which is higher-order (for the same number of
            // quadrature points) and has a few more digits of
            // precision in the points and weights.  Some 10-point
            // degree 6 rules exist for the triangle but they have
            // quadrature points outside the region of integration.
          case SEVENTH:
            {
              _points.resize (12);
              _weights.resize(12);
              
              const unsigned int nrows=4;
              
              // In each of the rows below, the first two entries are (z1, z2) which imply
              // z3.  The third entry is the weight for each of the points in the cyclic permutation.
              const Real p[nrows][3] = {
                {6.2382265094402118e-02, 6.7517867073916085e-02, 2.6517028157436251e-02}, // group A
                {5.5225456656926611e-02, 3.2150249385198182e-01, 4.3881408714446055e-02}, // group B
                {3.4324302945097146e-02, 6.6094919618673565e-01, 2.8775042784981585e-02}, // group C
                {5.1584233435359177e-01, 2.7771616697639178e-01, 6.7493187009802774e-02}  // group D
              };

              for (unsigned int i=0, offset=0; i<nrows; ++i)
                {
                  _points[offset + 0] = Point(p[i][0],            p[i][1]); // (z1,z2)
                  _points[offset + 1] = Point(1.-p[i][0]-p[i][1], p[i][0]); // (z3,z1)
                  _points[offset + 2] = Point(p[i][1], 1.-p[i][0]-p[i][1]); // (z2,z3)

                  // All these points get the same weight
                  _weights[offset + 0] = p[i][2];
                  _weights[offset + 1] = p[i][2];
                  _weights[offset + 2] = p[i][2];

                  // Increment offset
                  offset += 3;
                }

              return;

              
//            // The following is an inferior 7th-order Lyness-style rule with 15 points.
//            // It's here only for completeness and the Ro3-invariant rule above should
//            // be used instead!
//            const unsigned int n_wts = 3;
//            const Real wts[n_wts] =
//              {
//                2.6538900895116205835977487499847719e-02L,
//                3.5426541846066783659206291623201826e-02L,
//                3.4637341039708446756138297960207647e-02L
//              };
// 
//            const Real a[n_wts] =
//              {
//                6.4930513159164863078379776030396538e-02L,
//                2.8457558424917033519741605734978046e-01L,
//                3.1355918438493150795585190219862865e-01L
//              };
// 
//            const Real b[n_wts] =
//              {
//                0.,
//                1.9838447668150671917987659863332941e-01L,
//                4.3863471792372471511798695971295936e-02L               
//              };
//            
//            const unsigned int permutation_ids[n_wts] = {3, 6, 6}; // 15 total points
// 
//            dunavant_rule2(wts, a, b, permutation_ids, n_wts);
// 
//            return;
            }



            
            // Another Dunavant rule.  This one has all positive weights.  This rule has
            // 16 points while a comparable conical product rule would have 5*5=25.
            //
            // It was copied 23rd June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
            //
            // Additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. 'A set of symmetric quadrature rules
            // on triangles and tetrahedra'  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
          case EIGHTH:
            {
              const unsigned int n_wts = 5;
              const Real wts[n_wts] =
                {
                  7.2157803838893584125545555244532310e-02L,
                  4.7545817133642312396948052194292159e-02L,
                  5.1608685267359125140895775146064515e-02L,
                  1.6229248811599040155462964170890299e-02L,
                  1.3615157087217497132422345036954462e-02L
                };

              const Real a[n_wts] =
                {
                  0.0, // 'a' parameter not used for origin permutation
                  4.5929258829272315602881551449416932e-01L,
                  1.7056930775176020662229350149146450e-01L,
                  5.0547228317030975458423550596598947e-02L,
                  2.6311282963463811342178578628464359e-01L,
                };

              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  0.,
                  0.,
                  7.2849239295540428124100037917606196e-01L
                };
              
              const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 6}; // 16 total points

              dunavant_rule2(wts, a, b, permutation_ids, n_wts);

              return;
            }


            
            // Another Dunavant rule.  This one has all positive weights.  This rule has 19
            // points. The comparable conical product rule would have 25.
            // It was copied 23rd June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
            //
            // Additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. 'A set of symmetric quadrature rules
            // on triangles and tetrahedra'  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
          case NINTH:
            {
              const unsigned int n_wts = 6;
              const Real wts[n_wts] =
                {
                  4.8567898141399416909620991253644315e-02L,
                  1.5667350113569535268427415643604658e-02L, 
                  1.2788837829349015630839399279499912e-02L, 
                  3.8913770502387139658369678149701978e-02L, 
                  3.9823869463605126516445887132022637e-02L, 
                  2.1641769688644688644688644688644689e-02L
                };

              const Real a[n_wts] =
                {
                  0.0, // 'a' parameter not used for origin permutation
                  4.8968251919873762778370692483619280e-01L,
                  4.4729513394452709865106589966276365e-02L,
                  4.3708959149293663726993036443535497e-01L,
                  1.8820353561903273024096128046733557e-01L,
                  2.2196298916076569567510252769319107e-01L
                };

              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  0.,
                  0.,
                  0.,
                  7.4119859878449802069007987352342383e-01L
                };
              
              const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 3, 6}; // 19 total points

              dunavant_rule2(wts, a, b, permutation_ids, n_wts);

              return;
            }

            
            // Another Dunavant rule with all positive weights.  This rule has 25
            // points. The comparable conical product rule would have 36.
            // It was copied 23rd June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
            //
            // Additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. 'A set of symmetric quadrature rules
            // on triangles and tetrahedra'  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
          case TENTH:
            {
              const unsigned int n_wts = 6;
              const Real wts[n_wts] =
                {
                  4.5408995191376790047643297550014267e-02L,
                  1.8362978878233352358503035945683300e-02L,
                  2.2660529717763967391302822369298659e-02L,
                  3.6378958422710054302157588309680344e-02L,
                  1.4163621265528742418368530791049552e-02L,
                  4.7108334818664117299637354834434138e-03L
                };

              const Real a[n_wts] =
                {
                  0.0, // 'a' parameter not used for origin permutation
                  4.8557763338365737736750753220812615e-01L,
                  1.0948157548503705479545863134052284e-01L,
                  3.0793983876412095016515502293063162e-01L,
                  2.4667256063990269391727646541117681e-01L,
                  6.6803251012200265773540212762024737e-02L
                };

              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  0.,
                  5.5035294182099909507816172659300821e-01L,
                  7.2832390459741092000873505358107866e-01L,
                  9.2365593358750027664630697761508843e-01L
                };
              
              const unsigned int permutation_ids[n_wts] = {1, 3, 3, 6, 6, 6}; // 25 total points

              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
              
              return;
            }


            // Dunavant's 11th-order rule contains points outside the region of
            // integration, and is thus unacceptable for our FEM calculations.
            //
            // This 30-point, 11th-order rule was obtained by me [JWP] using the code in
            //
            // Additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. 'A set of symmetric quadrature rules
            // on triangles and tetrahedra'  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
            //
            // Note: the 28-point 11th-order rule obtained by Zhang in the paper above
            // does not appear to be unique.  It is a solution in the sense that it
            // minimizes the error in the least-squares minimization problem, but
            // it involves too many unknowns and the Jacobian is therefore singular
            // when attempting to improve the solution via Newton's method. 
          case ELEVENTH:
            {
              const unsigned int n_wts = 6;
              const Real wts[n_wts] =
                {
                  3.6089021198604635216985338480426484e-02L,
                  2.1607717807680420303346736867931050e-02L,
                  3.1144524293927978774861144478241807e-03L,
                  2.9086855161081509446654185084988077e-02L,
                  8.4879241614917017182977532679947624e-03L,
                  1.3795732078224796530729242858347546e-02L
                };

              const Real a[n_wts] =
                {
                  3.9355079629947969884346551941969960e-01L,
                  4.7979065808897448654107733982929214e-01L,
                  5.1003445645828061436081405648347852e-03L,
                  2.6597620190330158952732822450744488e-01L,
                  2.8536418538696461608233522814483715e-01L,
                  1.3723536747817085036455583801851025e-01L
                };

              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  5.6817155788572446538150614865768991e-02L,
                  1.2539956353662088473247489775203396e-01L,
                  1.2409970153698532116262152247041742e-02L,
                  5.2792057988217708934207928630851643e-02L               
                };
              
              const unsigned int permutation_ids[n_wts] = {3, 3, 6, 6, 6, 6}; // 30 total points

              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
              
              return;
            }



            
            // Another Dunavant rule with all positive weights.  This rule has 33
            // points. The comparable conical product rule would have 36 (ELEVENTH) or 49 (TWELFTH).
            //
            // It was copied 23rd June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
            //
            // Additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. 'A set of symmetric quadrature rules
            // on triangles and tetrahedra'  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
          case TWELFTH:
            {
              const unsigned int n_wts = 8;
              const Real wts[n_wts] =
                {
                  3.0831305257795086169332418926151771e-03L,
                  3.1429112108942550177135256546441273e-02L,
                  1.7398056465354471494664198647499687e-02L,
                  2.1846272269019201067728631278737487e-02L,
                  1.2865533220227667708895461535782215e-02L,
                  1.1178386601151722855919538351159995e-02L,
                  8.6581155543294461858210504055170332e-03L,
                  2.0185778883190464758914349626118386e-02L
                };

              const Real a[n_wts] =
                {
                  2.1317350453210370246856975515728246e-02L,
                  2.7121038501211592234595134039689474e-01L,
                  1.2757614554158592467389632515428357e-01L,
                  4.3972439229446027297973662348436108e-01L,
                  4.8821738977380488256466206525881104e-01L,
                  2.8132558098993954824813069297455275e-01L,
                  1.1625191590759714124135414784260182e-01L,
                  2.7571326968551419397479634607976398e-01L
                };

              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  0.,
                  0.,
                  0.,
                  6.9583608678780342214163552323607254e-01L,
                  8.5801403354407263059053661662617818e-01L,
                  6.0894323577978780685619243776371007e-01L
                };
              
              const unsigned int permutation_ids[n_wts] = {3, 3, 3, 3, 3, 6, 6, 6}; // 33 total points

              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
              
              return;
            }

            
            // Another Dunavant rule with all positive weights.  This rule has 37
            // points. The comparable conical product rule would have 49 points.
            //
            // It was copied 23rd June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
            //
            // A second rule with additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. 'A set of symmetric quadrature rules
            // on triangles and tetrahedra'  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
          case THIRTEENTH:
            {
              const unsigned int n_wts = 9;
              const Real wts[n_wts] =
                {
                  3.3980018293415822140887212340442440e-02L,
                  2.7800983765226664353628733005230734e-02L,
                  2.9139242559599990702383541756669905e-02L,
                  3.0261685517695859208964000161454122e-03L,
                  1.1997200964447365386855399725479827e-02L,
                  1.7320638070424185232993414255459110e-02L,
                  7.4827005525828336316229285664517190e-03L,
                  1.2089519905796909568722872786530380e-02L,
                  4.7953405017716313612975450830554457e-03L
                };

              const Real a[n_wts] =
                {
                  0., // 'a' parameter not used for origin permutation
                  4.2694141425980040602081253503137421e-01L,
                  2.2137228629183290065481255470507908e-01L,
                  2.1509681108843183869291313534052083e-02L,
                  4.8907694645253934990068971909020439e-01L,
                  3.0844176089211777465847185254124531e-01L,
                  1.1092204280346339541286954522167452e-01L,
                  1.6359740106785048023388790171095725e-01L,
                  2.7251581777342966618005046435408685e-01L
                };

              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  0.,
                  0.,
                  0.,
                  6.2354599555367557081585435318623659e-01L,
                  8.6470777029544277530254595089569318e-01L,
                  7.4850711589995219517301859578870965e-01L,
                  7.2235779312418796526062013230478405e-01L 
                };
              
              const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 3, 6, 6, 6, 6}; // 37 total points

              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
              
              return;
            }

            
            // Another Dunavant rule.  This rule has 42 points, while
            // a comparable conical product rule would have 64.
            //
            // It was copied 23rd June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
            //
            // Additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. 'A set of symmetric quadrature rules
            // on triangles and tetrahedra'  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
          case FOURTEENTH:
            {
              const unsigned int n_wts = 10;
              const Real wts[n_wts] =
                {
                  1.0941790684714445320422472981662986e-02L,
                  1.6394176772062675320655489369312672e-02L,
                  2.5887052253645793157392455083198201e-02L,
                  2.1081294368496508769115218662093065e-02L,
                  7.2168498348883338008549607403266583e-03L,
                  2.4617018012000408409130117545210774e-03L,
                  1.2332876606281836981437622591818114e-02L,
                  1.9285755393530341614244513905205430e-02L,
                  7.2181540567669202480443459995079017e-03L,
                  2.5051144192503358849300465412445582e-03L
                };

              const Real a[n_wts] =
                {
                  4.8896391036217863867737602045239024e-01L,
                  4.1764471934045392250944082218564344e-01L,
                  2.7347752830883865975494428326269856e-01L,
                  1.7720553241254343695661069046505908e-01L,
                  6.1799883090872601267478828436935788e-02L,
                  1.9390961248701048178250095054529511e-02L,
                  1.7226668782135557837528960161365733e-01L,
                  3.3686145979634500174405519708892539e-01L,
                  2.9837288213625775297083151805961273e-01L,
                  1.1897449769695684539818196192990548e-01L
                };

              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  0.,
                  0.,
                  0.,
                  0.,             
                  7.7060855477499648258903327416742796e-01L,
                  5.7022229084668317349769621336235426e-01L,
                  6.8698016780808783735862715402031306e-01L,
                  8.7975717137017112951457163697460183e-01L 
                };
              
              const unsigned int permutation_ids[n_wts]
                = {3, 3, 3, 3, 3, 3, 6, 6, 6, 6}; // 42 total points

              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
              
              return;
            }

            
            // This 49-point rule was found by me [JWP] using the code in:
            //
            // L. Zhang, T. Cui, and H. Liu. 'A set of symmetric quadrature rules
            // on triangles and tetrahedra'  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
            //
            // A 54-point, 15th-order rule is reported by
            //
            // Stephen Wandzura, Hong Xiao,
            // Symmetric Quadrature Rules on a Triangle,
            // Computers and Mathematics with Applications,
            // Volume 45, Number 12, June 2003, pages 1829-1840.
            //
            // can be found here:
            // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
            //
            // but this 49-point rule is superior.
          case FIFTEENTH:
            {
              const unsigned int n_wts = 11;
              const Real wts[n_wts] =
                {
                  2.4777380743035579804788826970198951e-02L,
                  9.2433943023307730591540642828347660e-03L,
                  2.2485768962175402793245929133296627e-03L,
                  6.7052581900064143760518398833360903e-03L,
                  1.9011381726930579256700190357527956e-02L,
                  1.4605445387471889398286155981802858e-02L,
                  1.5087322572773133722829435011138258e-02L,
                  1.5630213780078803020711746273129099e-02L,
                  6.1808086085778203192616856133701233e-03L,
                  3.2209366452594664857296985751120513e-03L,
                  5.8747373242569702667677969985668817e-03L
                };

              const Real a[n_wts] =
                {
                  0.0, // 'a' parameter not used for origin
                  7.9031013655541635005816956762252155e-02L,
                  1.8789501810770077611247984432284226e-02L,
                  4.9250168823249670532514526605352905e-01L,
                  4.0886316907744105975059040108092775e-01L,
                  5.3877851064220142445952549348423733e-01L,
                  2.0250549804829997692885033941362673e-01L,
                  5.5349674918711643207148086558288110e-01L,
                  7.8345022567320812359258882143250181e-01L,
                  8.9514624528794883409864566727625002e-01L,
                  3.2515745241110782862789881780746490e-01L
                };

              const Real b[n_wts] =
                {
                  0.,
                  0.,
                  0.,
                  0.,
                  0.,
                  1.9412620368774630292701241080996842e-01L,
                  9.8765911355712115933807754318089099e-02L,
                  7.7663767064308164090246588765178087e-02L,
                  2.1594628433980258573654682690950798e-02L,
                  1.2563596287784997705599005477153617e-02L,
                  1.5082654870922784345283124845552190e-02L               
                };
              
              const unsigned int permutation_ids[n_wts]
                = {1, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6}; // 49 total points

              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
              
              return;
            }



            
            // Dunavant's 16th-order rule contains points outside the region of
            // integration, and is thus unacceptable for our FEM calculations.
            //
            // This 55-point, 16th-order rule was obtained by me [JWP] using the code in
            //
            // Additional precision obtained from the code in:
            // L. Zhang, T. Cui, and H. Liu. 'A set of symmetric quadrature rules
            // on triangles and tetrahedra'  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
            //
            // Note: the 55-point 16th-order rule obtained by Zhang in the paper above
            // does not appear to be unique.  It is a solution in the sense that it
            // minimizes the error in the least-squares minimization problem, but
            // it involves too many unknowns and the Jacobian is therefore singular
            // when attempting to improve the solution via Newton's method. 
          case SIXTEENTH:    
            {
              const unsigned int n_wts = 12;
              const Real wts[n_wts] =
                {
                  2.2668082505910087151996321171534230e-02L,
                  8.4043060714818596159798961899306135e-03L,
                  1.0850949634049747713966288634484161e-03L,
                  7.2252773375423638869298219383808751e-03L,
                  1.2997715227338366024036316182572871e-02L,
                  2.0054466616677715883228810959112227e-02L,
                  9.7299841600417010281624372720122710e-03L,
                  1.1651974438298104227427176444311766e-02L,
                  9.1291185550484450744725847363097389e-03L,
                  3.5568614040947150231712567900113671e-03L,
                  5.8355861686234326181790822005304303e-03L,
                  4.7411314396804228041879331486234396e-03L
                };

              const Real a[n_wts] =
                {
                  0.0, // 'a' parameter not used for centroid weight
                  8.5402539407933203673769900926355911e-02L,
                  1.2425572001444092841183633409631260e-02L,
                  4.9174838341891594024701017768490960e-01L,
                  4.5669426695387464162068900231444462e-01L,
                  4.8506759880447437974189793537259677e-01L,
                  2.0622099278664205707909858461264083e-01L,
                  3.2374950270039093446805340265853956e-01L,
                  7.3834330556606586255186213302750029e-01L,
                  9.1210673061680792565673823935174611e-01L,
                  6.6129919222598721544966837350891531e-01L,
                  1.7807138906021476039088828811346122e-01L
                };

              const Real b[n_wts] =
                {
                  0.0,
                  0.0,
                  0.0,
                  0.0,
                  0.0,
                  3.2315912848634384647700266402091638e-01L,
                  1.5341553679414688425981898952416987e-01L,
                  7.4295478991330687632977899141707872e-02L,
                  7.1278762832147862035977841733532020e-02L,
                  1.6623223223705792825395256602140459e-02L,
                  1.4160772533794791868984026749196156e-02L,
                  1.4539694958941854654807449467759690e-02L               
                };
              
              const unsigned int permutation_ids[n_wts]
                = {1, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6, 6}; // 55 total points

              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
              
              return;
            }
            
            
            // Dunavant's 17th-order rule has 61 points, while a
            // comparable conical product rule would have 81 (16th and 17th orders).
            //
            // It can be found here:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
            //
            // Zhang reports an identical rule in:
            // L. Zhang, T. Cui, and H. Liu. 'A set of symmetric quadrature rules
            // on triangles and tetrahedra'  Journal of Computational Mathematics,
            // v. 27, no. 1, 2009, pp. 89-96.
            //
            // Note: the 61-point 17th-order rule obtained by Dunavant and Zhang
            // does not appear to be unique.  It is a solution in the sense that it
            // minimizes the error in the least-squares minimization problem, but
            // it involves too many unknowns and the Jacobian is therefore singular
            // when attempting to improve the solution via Newton's method.
            //
            // Therefore, we prefer the following 63-point rule which
            // I [JWP] found.  It appears to be more accurate than the
            // rule reported by Dunavant and Zhang, even though it has
            // a few more points.
          case SEVENTEENTH:
            {
              const unsigned int n_wts = 12;
              const Real wts[n_wts] =
                {
                  1.7464603792572004485690588092246146e-02L,
                  5.9429003555801725246549713984660076e-03L,
                  1.2490753345169579649319736639588729e-02L,
                  1.5386987188875607593083456905596468e-02L,
                  1.1185807311917706362674684312990270e-02L,
                  1.0301845740670206831327304917180007e-02L,
                  1.1767783072977049696840016810370464e-02L,
                  3.8045312849431209558329128678945240e-03L,
                  4.5139302178876351271037137230354382e-03L,
                  2.2178812517580586419412547665472893e-03L,
                  5.2216271537483672304731416553063103e-03L,
                  9.8381136389470256422419930926212114e-04L
                };

              const Real a[n_wts] =
                {
                  2.8796825754667362165337965123570514e-01L,
                  4.9216175986208465345536805750663939e-01L,
                  4.6252866763171173685916780827044612e-01L,
                  1.6730292951631792248498303276090273e-01L,
                  1.5816335500814652972296428532213019e-01L,
                  1.6352252138387564873002458959679529e-01L,
                  6.2447680488959768233910286168417367e-01L,
                  8.7317249935244454285263604347964179e-01L,
                  3.4428164322282694677972239461699271e-01L,
                  9.1584484467813674010523309855340209e-02L,
                  2.0172088013378989086826623852040632e-01L,
                  9.6538762758254643474731509845084691e-01L
                };

              const Real b[n_wts] =
                {
                  0.0,
                  0.0,
                  0.0,
                  3.4429160695501713926320695771253348e-01L,
                  2.2541623431550639817203145525444726e-01L,
                  8.0670083153531811694942222940484991e-02L,
                  6.5967451375050925655738829747288190e-02L,
                  4.5677879890996762665044366994439565e-02L,
                  1.1528411723154215812386518751976084e-02L,
                  9.3057714323900610398389176844165892e-03L,
                  1.5916814107619812717966560404970160e-02L,
                  1.0734733163764032541125434215228937e-02L
                };
              
              const unsigned int permutation_ids[n_wts]
                = {3, 3, 3, 6, 6, 6, 6, 6, 6, 6, 6, 6}; // 63 total points

              dunavant_rule2(wts, a, b, permutation_ids, n_wts);
              
              return;
   
//            _points.resize (61);
//            _weights.resize(61);

//            // The raw data for the quadrature rule.
//            const Real p[15][4] = {
//              {                1./3.,                    0.,                    0., 0.033437199290803e+00 / 2.0}, // 1-perm
//              {0.005658918886452e+00, 0.497170540556774e+00,                    0., 0.005093415440507e+00 / 2.0}, // 3-perm
//              {0.035647354750751e+00, 0.482176322624625e+00,                    0., 0.014670864527638e+00 / 2.0}, // 3-perm
//              {0.099520061958437e+00, 0.450239969020782e+00,                    0., 0.024350878353672e+00 / 2.0}, // 3-perm
//              {0.199467521245206e+00, 0.400266239377397e+00,                    0., 0.031107550868969e+00 / 2.0}, // 3-perm
//              {0.495717464058095e+00, 0.252141267970953e+00,                    0., 0.031257111218620e+00 / 2.0}, // 3-perm
//              {0.675905990683077e+00, 0.162047004658461e+00,                    0., 0.024815654339665e+00 / 2.0}, // 3-perm
//              {0.848248235478508e+00, 0.075875882260746e+00,                    0., 0.014056073070557e+00 / 2.0}, // 3-perm
//              {0.968690546064356e+00, 0.015654726967822e+00,                    0., 0.003194676173779e+00 / 2.0}, // 3-perm
//              {0.010186928826919e+00, 0.334319867363658e+00, 0.655493203809423e+00, 0.008119655318993e+00 / 2.0}, // 6-perm
//              {0.135440871671036e+00, 0.292221537796944e+00, 0.572337590532020e+00, 0.026805742283163e+00 / 2.0}, // 6-perm
//              {0.054423924290583e+00, 0.319574885423190e+00, 0.626001190286228e+00, 0.018459993210822e+00 / 2.0}, // 6-perm 
//              {0.012868560833637e+00, 0.190704224192292e+00, 0.796427214974071e+00, 0.008476868534328e+00 / 2.0}, // 6-perm
//              {0.067165782413524e+00, 0.180483211648746e+00, 0.752351005937729e+00, 0.018292796770025e+00 / 2.0}, // 6-perm 
//              {0.014663182224828e+00, 0.080711313679564e+00, 0.904625504095608e+00, 0.006665632004165e+00 / 2.0}  // 6-perm
//            };

              
//            // Now call the dunavant routine to generate _points and _weights
//            dunavant_rule(p, 15);

//            return;
            }


            
            // Dunavant's 18th-order rule contains points outside the region and is therefore unsuitable
            // for our FEM calculations.  His 19th-order rule has 73 points, compared with 100 points for
            // a comparable-order conical product rule.
            //
            // It was copied 23rd June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
          case EIGHTTEENTH:  
          case NINTEENTH:
            {
              _points.resize (73);
              _weights.resize(73);

              // The raw data for the quadrature rule.
              const Real p[17][4] = {
                {                1./3.,                    0.,                    0., 0.032906331388919e+00 / 2.0}, // 1-perm
                {0.020780025853987e+00, 0.489609987073006e+00,                    0., 0.010330731891272e+00 / 2.0}, // 3-perm
                {0.090926214604215e+00, 0.454536892697893e+00,                    0., 0.022387247263016e+00 / 2.0}, // 3-perm
                {0.197166638701138e+00, 0.401416680649431e+00,                    0., 0.030266125869468e+00 / 2.0}, // 3-perm
                {0.488896691193805e+00, 0.255551654403098e+00,                    0., 0.030490967802198e+00 / 2.0}, // 3-perm
                {0.645844115695741e+00, 0.177077942152130e+00,                    0., 0.024159212741641e+00 / 2.0}, // 3-perm
                {0.779877893544096e+00, 0.110061053227952e+00,                    0., 0.016050803586801e+00 / 2.0}, // 3-perm
                {0.888942751496321e+00, 0.055528624251840e+00,                    0., 0.008084580261784e+00 / 2.0}, // 3-perm
                {0.974756272445543e+00, 0.012621863777229e+00,                    0., 0.002079362027485e+00 / 2.0}, // 3-perm
                {0.003611417848412e+00, 0.395754787356943e+00, 0.600633794794645e+00, 0.003884876904981e+00 / 2.0}, // 6-perm
                {0.134466754530780e+00, 0.307929983880436e+00, 0.557603261588784e+00, 0.025574160612022e+00 / 2.0}, // 6-perm
                {0.014446025776115e+00, 0.264566948406520e+00, 0.720987025817365e+00, 0.008880903573338e+00 / 2.0}, // 6-perm 
                {0.046933578838178e+00, 0.358539352205951e+00, 0.594527068955871e+00, 0.016124546761731e+00 / 2.0}, // 6-perm
                {0.002861120350567e+00, 0.157807405968595e+00, 0.839331473680839e+00, 0.002491941817491e+00 / 2.0}, // 6-perm 
                {0.223861424097916e+00, 0.075050596975911e+00, 0.701087978926173e+00, 0.018242840118951e+00 / 2.0}, // 6-perm
                {0.034647074816760e+00, 0.142421601113383e+00, 0.822931324069857e+00, 0.010258563736199e+00 / 2.0}, // 6-perm 
                {0.010161119296278e+00, 0.065494628082938e+00, 0.924344252620784e+00, 0.003799928855302e+00 / 2.0}  // 6-perm
              };

              
              // Now call the dunavant routine to generate _points and _weights
              dunavant_rule(p, 17);

              return;
            }

            
            // 20th-order rule by Wandzura.  
            //
            // Stephen Wandzura, Hong Xiao,
            // Symmetric Quadrature Rules on a Triangle,
            // Computers and Mathematics with Applications,
            // Volume 45, Number 12, June 2003, pages 1829-1840.
            //
            // Wandzura's work extends the work of Dunavant by providing degree
            // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
            //
            // Copied on 3rd July 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
          case TWENTIETH:
            {
              // The equivalent concial product rule would have 121 points
              _points.resize (85);
              _weights.resize(85);

              // The raw data for the quadrature rule.
              const Real p[19][4] = {
                {0.33333333333333e+00,                  0.0,                  0.0, 0.2761042699769952e-01 / 2.0}, // 1-perm
                {0.00150064932443e+00, 0.49924967533779e+00,                  0.0, 0.1779029547326740e-02 / 2.0}, // 3-perm 
                {0.09413975193895e+00, 0.45293012403052e+00,                  0.0, 0.2011239811396117e-01 / 2.0}, // 3-perm
                {0.20447212408953e+00, 0.39776393795524e+00,                  0.0, 0.2681784725933157e-01 / 2.0}, // 3-perm
                {0.47099959493443e+00, 0.26450020253279e+00,                  0.0, 0.2452313380150201e-01 / 2.0}, // 3-perm
                {0.57796207181585e+00, 0.21101896409208e+00,                  0.0, 0.1639457841069539e-01 / 2.0}, // 3-perm
                {0.78452878565746e+00, 0.10773560717127e+00,                  0.0, 0.1479590739864960e-01 / 2.0}, // 3-perm
                {0.92186182432439e+00, 0.03906908783780e+00,                  0.0, 0.4579282277704251e-02 / 2.0}, // 3-perm
                {0.97765124054134e+00, 0.01117437972933e+00,                  0.0, 0.1651826515576217e-02 / 2.0}, // 3-perm
                {0.00534961818734e+00, 0.06354966590835e+00, 0.93110071590431e+00, 0.2349170908575584e-02 / 2.0}, // 6-perm
                {0.00795481706620e+00, 0.15710691894071e+00, 0.83493826399309e+00, 0.4465925754181793e-02 / 2.0}, // 6-perm
                {0.01042239828126e+00, 0.39564211436437e+00, 0.59393548735436e+00, 0.6099566807907972e-02 / 2.0}, // 6-perm
                {0.01096441479612e+00, 0.27316757071291e+00, 0.71586801449097e+00, 0.6891081327188203e-02 / 2.0}, // 6-perm
                {0.03856671208546e+00, 0.10178538248502e+00, 0.85964790542952e+00, 0.7997475072478163e-02 / 2.0}, // 6-perm
                {0.03558050781722e+00, 0.44665854917641e+00, 0.51776094300637e+00, 0.7386134285336024e-02 / 2.0}, // 6-perm
                {0.04967081636276e+00, 0.19901079414950e+00, 0.75131838948773e+00, 0.1279933187864826e-01 / 2.0}, // 6-perm
                {0.05851972508433e+00, 0.32426118369228e+00, 0.61721909122339e+00, 0.1725807117569655e-01 / 2.0}, // 6-perm
                {0.12149778700439e+00, 0.20853136321013e+00, 0.66997084978547e+00, 0.1867294590293547e-01 / 2.0}, // 6-perm
                {0.14071084494394e+00, 0.32317056653626e+00, 0.53611858851980e+00, 0.2281822405839526e-01 / 2.0}  // 6-perm
              };

              
              // Now call the dunavant routine to generate _points and _weights
              dunavant_rule(p, 19);

              return;
            }



            // 25th-order rule by Wandzura.  
            //
            // Stephen Wandzura, Hong Xiao,
            // Symmetric Quadrature Rules on a Triangle,
            // Computers and Mathematics with Applications,
            // Volume 45, Number 12, June 2003, pages 1829-1840.
            //
            // Wandzura's work extends the work of Dunavant by providing degree
            // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
            //
            // Copied on 3rd July 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
            // case TWENTYFIRST: // fall through to 121 point conical product rule below
          case TWENTYSECOND:  
          case TWENTYTHIRD:   
          case TWENTYFOURTH:        
          case TWENTYFIFTH:
            {
              // The equivalent concial product rule would have 169 points
              _points.resize (126);
              _weights.resize(126);

              // The raw data for the quadrature rule.
              const Real p[26][4] = {
                {0.02794648307317e+00, 0.48602675846341e+00,                  0.0, 0.8005581880020417e-02 / 2.0},  // 3-perm
                {0.13117860132765e+00, 0.43441069933617e+00,                  0.0, 0.1594707683239050e-01 / 2.0},  // 3-perm
                {0.22022172951207e+00, 0.38988913524396e+00,                  0.0, 0.1310914123079553e-01 / 2.0},  // 3-perm
                {0.40311353196039e+00, 0.29844323401980e+00,                  0.0, 0.1958300096563562e-01 / 2.0},  // 3-perm
                {0.53191165532526e+00, 0.23404417233737e+00,                  0.0, 0.1647088544153727e-01 / 2.0},  // 3-perm
                {0.69706333078196e+00, 0.15146833460902e+00,                  0.0, 0.8547279074092100e-02 / 2.0},  // 3-perm
                {0.77453221290801e+00, 0.11273389354599e+00,                  0.0, 0.8161885857226492e-02 / 2.0},  // 3-perm
                {0.84456861581695e+00, 0.07771569209153e+00,                  0.0, 0.6121146539983779e-02 / 2.0},  // 3-perm
                {0.93021381277141e+00, 0.03489309361430e+00,                  0.0, 0.2908498264936665e-02 / 2.0},  // 3-perm
                {0.98548363075813e+00, 0.00725818462093e+00,                  0.0, 0.6922752456619963e-03 / 2.0},  // 3-perm
                {0.00129235270444e+00, 0.22721445215336e+00, 0.77149319514219e+00, 0.1248289199277397e-02 / 2.0},  // 6-perm
                {0.00539970127212e+00, 0.43501055485357e+00, 0.55958974387431e+00, 0.3404752908803022e-02 / 2.0},  // 6-perm
                {0.00638400303398e+00, 0.32030959927220e+00, 0.67330639769382e+00, 0.3359654326064051e-02 / 2.0},  // 6-perm
                {0.00502821150199e+00, 0.09175032228001e+00, 0.90322146621800e+00, 0.1716156539496754e-02 / 2.0},  // 6-perm
                {0.00682675862178e+00, 0.03801083585872e+00, 0.95516240551949e+00, 0.1480856316715606e-02 / 2.0},  // 6-perm
                {0.01001619963993e+00, 0.15742521848531e+00, 0.83255858187476e+00, 0.3511312610728685e-02 / 2.0},  // 6-perm
                {0.02575781317339e+00, 0.23988965977853e+00, 0.73435252704808e+00, 0.7393550149706484e-02 / 2.0},  // 6-perm
                {0.03022789811992e+00, 0.36194311812606e+00, 0.60782898375402e+00, 0.7983087477376558e-02 / 2.0},  // 6-perm
                {0.03050499010716e+00, 0.08355196095483e+00, 0.88594304893801e+00, 0.4355962613158041e-02 / 2.0},  // 6-perm
                {0.04595654736257e+00, 0.14844322073242e+00, 0.80560023190501e+00, 0.7365056701417832e-02 / 2.0},  // 6-perm
                {0.06744280054028e+00, 0.28373970872753e+00, 0.64881749073219e+00, 0.1096357284641955e-01 / 2.0},  // 6-perm
                {0.07004509141591e+00, 0.40689937511879e+00, 0.52305553346530e+00, 0.1174996174354112e-01 / 2.0},  // 6-perm
                {0.08391152464012e+00, 0.19411398702489e+00, 0.72197448833499e+00, 0.1001560071379857e-01 / 2.0},  // 6-perm
                {0.12037553567715e+00, 0.32413434700070e+00, 0.55549011732214e+00, 0.1330964078762868e-01 / 2.0},  // 6-perm
                {0.14806689915737e+00, 0.22927748355598e+00, 0.62265561728665e+00, 0.1415444650522614e-01 / 2.0},  // 6-perm
                {0.19177186586733e+00, 0.32561812259598e+00, 0.48261001153669e+00, 0.1488137956116801e-01 / 2.0}   // 6-perm
              };

              
              // Now call the dunavant routine to generate _points and _weights
              dunavant_rule(p, 26);

              return;
            }



            // 30th-order rule by Wandzura.  
            //
            // Stephen Wandzura, Hong Xiao,
            // Symmetric Quadrature Rules on a Triangle,
            // Computers and Mathematics with Applications,
            // Volume 45, Number 12, June 2003, pages 1829-1840.
            //
            // Wandzura's work extends the work of Dunavant by providing degree
            // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
            //
            // Copied on 3rd July 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
          case TWENTYSIXTH:   
          case TWENTYSEVENTH: 
          case TWENTYEIGHTH:  
          case TWENTYNINTH:         
          case THIRTIETH:
            {
              // The equivalent concial product rule would have 256 points
              _points.resize (175);
              _weights.resize(175);

              // The raw data for the quadrature rule.
              const Real p[36][4] = {
                {0.33333333333333e+00,                  0.0,                  0.0, 0.1557996020289920e-01 / 2.0}, // 1-perm 
                {0.00733011643277e+00, 0.49633494178362e+00,                  0.0, 0.3177233700534134e-02 / 2.0}, // 3-perm  
                {0.08299567580296e+00, 0.45850216209852e+00,                  0.0, 0.1048342663573077e-01 / 2.0}, // 3-perm  
                {0.15098095612541e+00, 0.42450952193729e+00,                  0.0, 0.1320945957774363e-01 / 2.0}, // 3-perm  
                {0.23590585989217e+00, 0.38204707005392e+00,                  0.0, 0.1497500696627150e-01 / 2.0}, // 3-perm  
                {0.43802430840785e+00, 0.28098784579608e+00,                  0.0, 0.1498790444338419e-01 / 2.0}, // 3-perm  
                {0.54530204829193e+00, 0.22734897585403e+00,                  0.0, 0.1333886474102166e-01 / 2.0}, // 3-perm  
                {0.65088177698254e+00, 0.17455911150873e+00,                  0.0, 0.1088917111390201e-01 / 2.0}, // 3-perm  
                {0.75348314559713e+00, 0.12325842720144e+00,                  0.0, 0.8189440660893461e-02 / 2.0}, // 3-perm  
                {0.83983154221561e+00, 0.08008422889220e+00,                  0.0, 0.5575387588607785e-02 / 2.0}, // 3-perm  
                {0.90445106518420e+00, 0.04777446740790e+00,                  0.0, 0.3191216473411976e-02 / 2.0}, // 3-perm  
                {0.95655897063972e+00, 0.02172051468014e+00,                  0.0, 0.1296715144327045e-02 / 2.0}, // 3-perm  
                {0.99047064476913e+00, 0.00476467761544e+00,                  0.0, 0.2982628261349172e-03 / 2.0}, // 3-perm  
                {0.00092537119335e+00, 0.41529527091331e+00, 0.58377935789334e+00, 0.9989056850788964e-03 / 2.0}, // 6-perm  
                {0.00138592585556e+00, 0.06118990978535e+00, 0.93742416435909e+00, 0.4628508491732533e-03 / 2.0}, // 6-perm  
                {0.00368241545591e+00, 0.16490869013691e+00, 0.83140889440718e+00, 0.1234451336382413e-02 / 2.0}, // 6-perm  
                {0.00390322342416e+00, 0.02503506223200e+00, 0.97106171434384e+00, 0.5707198522432062e-03 / 2.0}, // 6-perm  
                {0.00323324815501e+00, 0.30606446515110e+00, 0.69070228669389e+00, 0.1126946125877624e-02 / 2.0}, // 6-perm  
                {0.00646743211224e+00, 0.10707328373022e+00, 0.88645928415754e+00, 0.1747866949407337e-02 / 2.0}, // 6-perm  
                {0.00324747549133e+00, 0.22995754934558e+00, 0.76679497516308e+00, 0.1182818815031657e-02 / 2.0}, // 6-perm  
                {0.00867509080675e+00, 0.33703663330578e+00, 0.65428827588746e+00, 0.1990839294675034e-02 / 2.0}, // 6-perm  
                {0.01559702646731e+00, 0.05625657618206e+00, 0.92814639735063e+00, 0.1900412795035980e-02 / 2.0}, // 6-perm  
                {0.01797672125369e+00, 0.40245137521240e+00, 0.57957190353391e+00, 0.4498365808817451e-02 / 2.0}, // 6-perm  
                {0.01712424535389e+00, 0.24365470201083e+00, 0.73922105263528e+00, 0.3478719460274719e-02 / 2.0}, // 6-perm  
                {0.02288340534658e+00, 0.16538958561453e+00, 0.81172700903888e+00, 0.4102399036723953e-02 / 2.0}, // 6-perm  
                {0.03273759728777e+00, 0.09930187449585e+00, 0.86796052821639e+00, 0.4021761549744162e-02 / 2.0}, // 6-perm  
                {0.03382101234234e+00, 0.30847833306905e+00, 0.65770065458860e+00, 0.6033164660795066e-02 / 2.0}, // 6-perm  
                {0.03554761446002e+00, 0.46066831859211e+00, 0.50378406694787e+00, 0.3946290302129598e-02 / 2.0}, // 6-perm  
                {0.05053979030687e+00, 0.21881529945393e+00, 0.73064491023920e+00, 0.6644044537680268e-02 / 2.0}, // 6-perm  
                {0.05701471491573e+00, 0.37920955156027e+00, 0.56377573352399e+00, 0.8254305856078458e-02 / 2.0}, // 6-perm  
                {0.06415280642120e+00, 0.14296081941819e+00, 0.79288637416061e+00, 0.6496056633406411e-02 / 2.0}, // 6-perm  
                {0.08050114828763e+00, 0.28373128210592e+00, 0.63576756960645e+00, 0.9252778144146602e-02 / 2.0}, // 6-perm  
                {0.10436706813453e+00, 0.19673744100444e+00, 0.69889549086103e+00, 0.9164920726294280e-02 / 2.0}, // 6-perm  
                {0.11384489442875e+00, 0.35588914121166e+00, 0.53026596435959e+00, 0.1156952462809767e-01 / 2.0}, // 6-perm  
                {0.14536348771552e+00, 0.25981868535191e+00, 0.59481782693256e+00, 0.1176111646760917e-01 / 2.0}, // 6-perm  
                {0.18994565282198e+00, 0.32192318123130e+00, 0.48813116594672e+00, 0.1382470218216540e-01 / 2.0}  // 6-perm   
              };

              
              // Now call the dunavant routine to generate _points and _weights
              dunavant_rule(p, 36);

              return;
            }
            
            
            // By default, we fall back on the conical product rules.  If the user
            // requests an order higher than what is currently available in the 1D
            // rules, an error will be thrown from the respective 1D code.
          default:
            {
              // The following quadrature rules are generated as
              // conical products.  These tend to be non-optimal
              // (use too many points, cluster points in certain
              // regions of the domain) but they are quite easy to
              // automatically generate using a 1D Gauss rule on
              // [0,1] and two 1D Jacobi-Gauss rules on [0,1].
              QConical conical_rule(2, _order);
              conical_rule.init(_type, p);

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

            
      //---------------------------------------------
      // Unsupported type
    default:
      {
        std::cerr << 'Element type not supported!:' << _type << std::endl;
        libmesh_error();
      }
    }

  libmesh_error();

  return;

#endif
}
 

void QGauss::init_3D (const ElemType_type = INVALID_ELEM, unsigned intp_level = 0) [private]

Definition at line 27 of file quadrature_gauss_3D.C.

References QBase::_points, QBase::_weights, QBase::allow_rules_with_negative_weights, libMeshEnums::CONSTANT, libMeshEnums::EDGE2, libMeshEnums::EIGHTH, libMeshEnums::FIFTH, libMeshEnums::FIRST, libMeshEnums::FOURTH, QBase::get_points(), QBase::get_weights(), libMeshEnums::HEX20, libMeshEnums::HEX27, libMeshEnums::HEX8, QBase::init(), keast_rule(), libMeshEnums::PRISM15, libMeshEnums::PRISM18, libMeshEnums::PRISM6, libMeshEnums::PYRAMID5, libMeshEnums::SECOND, libMeshEnums::SEVENTH, libMeshEnums::SIXTH, libMeshEnums::TET10, libMeshEnums::TET4, libMeshEnums::THIRD, libMeshEnums::TRI3, and QBase::w().

{
#if LIBMESH_DIM == 3
  
  //-----------------------------------------------------------------------
  // 3D quadrature rules
  switch (_type)
    {
      //---------------------------------------------
      // Hex quadrature rules
    case HEX8:
    case HEX20:
    case HEX27:
      {
        // We compute the 3D quadrature rule as a tensor
        // product of the 1D quadrature rule.
        QGauss q1D(1,_order);
        q1D.init(EDGE2,p);

        tensor_product_hex( q1D );
        
        return;
      }


      
      //---------------------------------------------
      // Tetrahedral quadrature rules
    case TET4:
    case TET10:
      {
        switch(_order + 2*p)
          {
            // Taken from pg. 222 of 'The finite element method,' vol. 1
            // ed. 5 by Zienkiewicz & Taylor
          case CONSTANT:
          case FIRST:
            {
              // Exact for linears
              _points.resize(1);
              _weights.resize(1);
                    
                    
              _points[0](0) = .25;
              _points[0](1) = .25;
              _points[0](2) = .25;
                    
              _weights[0] = .1666666666666666666666666666666666666666666667L;
                    
              return;
            }
          case SECOND:
            {
              // Exact for quadratics
              _points.resize(4);
              _weights.resize(4);
                    
                    
              const Real a = .585410196624969;
              const Real b = .138196601125011;
                    
              _points[0](0) = a;
              _points[0](1) = b;
              _points[0](2) = b;
                    
              _points[1](0) = b;
              _points[1](1) = a;
              _points[1](2) = b;
                    
              _points[2](0) = b;
              _points[2](1) = b;
              _points[2](2) = a;
                    
              _points[3](0) = b;
              _points[3](1) = b;
              _points[3](2) = b;
                    
                    
                    
              _weights[0] = .0416666666666666666666666666666666666666666667;
              _weights[1] = _weights[0];
              _weights[2] = _weights[0];
              _weights[3] = _weights[0];
                    
              return;
            }


            
            // Can be found in the class notes
            // http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps
            // by Flaherty.
            //
            // Caution: this rule has a negative weight and may be
            // unsuitable for some problems.
            // Exact for cubics.
            //
            // Note: Keast (see ref. elsewhere in this file) also gives
            // a third-order rule with positive weights, but it contains points
            // on the ref. elt. boundary, making it less suitable for FEM calculations.
          case THIRD:
            {
              if (allow_rules_with_negative_weights)
                {
                  _points.resize(5);
                  _weights.resize(5);
                    
                    
                  _points[0](0) = .25;
                  _points[0](1) = .25;
                  _points[0](2) = .25;
                    
                  _points[1](0) = .5;
                  _points[1](1) = .16666666666666666666666666666666666666666667;
                  _points[1](2) = .16666666666666666666666666666666666666666667;
                    
                  _points[2](0) = .16666666666666666666666666666666666666666667;
                  _points[2](1) = .5;
                  _points[2](2) = .16666666666666666666666666666666666666666667;
                    
                  _points[3](0) = .16666666666666666666666666666666666666666667;
                  _points[3](1) = .16666666666666666666666666666666666666666667;
                  _points[3](2) = .5;
                    
                  _points[4](0) = .16666666666666666666666666666666666666666667;
                  _points[4](1) = .16666666666666666666666666666666666666666667;
                  _points[4](2) = .16666666666666666666666666666666666666666667;
                    
                    
                  _weights[0] = -.133333333333333333333333333333333333333333333;
                  _weights[1] = .075;
                  _weights[2] = _weights[1];
                  _weights[3] = _weights[1];
                  _weights[4] = _weights[1];
                    
                  return;
                } // end if (allow_rules_with_negative_weights)
              else
                {
                  // If a rule with postive weights is required, the 2x2x2 conical
                  // product rule is third-order accurate and has less points than
                  // the next-available positive-weight rule at FIFTH order. 
                  QConical conical_rule(3, _order);
                  conical_rule.init(_type, p);

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

                  return;
                }
              // Note: if !allow_rules_with_negative_weights, fall through to next case.
            }

            

            // Originally a Keast rule,
            //    Patrick Keast,
            //    Moderate Degree Tetrahedral Quadrature Formulas,
            //    Computer Methods in Applied Mechanics and Engineering,
            //    Volume 55, Number 3, May 1986, pages 339-348.
            //
            // Can also be found the class notes
            // http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps
            // by Flaherty.
            //
            // Caution: this rule has a negative weight and may be
            // unsuitable for some problems.
          case FOURTH:
            {
              if (allow_rules_with_negative_weights)
                {
                  _points.resize(11);
                  _weights.resize(11);
                    
                  // The raw data for the quadrature rule.
                  const Real p[3][4] = {
                    {0.250000000000000000e+00,                         0.,                            0.,  -0.131555555555555556e-01},  // 1
                    {0.785714285714285714e+00,   0.714285714285714285e-01,                            0.,   0.762222222222222222e-02},  // 4
                    {0.399403576166799219e+00,                         0.,      0.100596423833200785e+00,   0.248888888888888889e-01}   // 6
                  };

              
                  // Now call the keast routine to generate _points and _weights
                  keast_rule(p, 3);

                  return;
                } // end if (allow_rules_with_negative_weights)
              // Note: if !allow_rules_with_negative_weights, fall through to next case.
            }



            
            // Walkington's fifth-order 14-point rule from
            // 'Quadrature on Simplices of Arbitrary Dimension'
            //
            // We originally had a Keast rule here, but this rule had
            // more points than an equivalent rule by Walkington and
            // also contained points on the boundary of the ref. elt,
            // making it less suitable for FEM calculations.
          case FIFTH:
            {
              _points.resize(14);
              _weights.resize(14);

              // permutations of these points and suitably-modified versions of
              // these points are the quadrature point locations
              const Real a[3] = {0.31088591926330060980,    // a1 from the paper
                                 0.092735250310891226402,   // a2 from the paper
                                 0.045503704125649649492};  // a3 from the paper

              // weights.  a[] and w[] are the only floating-point inputs required
              // for this rule.
              const Real w[3] = {0.018781320953002641800,    // w1 from the paper
                                 0.012248840519393658257,    // w2 from the paper
                                 0.0070910034628469110730};  // w3 from the paper

              // The first two sets of 4 points are formed in a similar manner
              for (unsigned int i=0; i<2; ++i)
                {
                  // Where we will insert values into _points and _weights
                  const unsigned int offset=4*i;

                  // Stuff points and weights values into their arrays
                  const Real b = 1. - 3.*a[i];

                  // Here are the permutations.  Order of these is not important,
                  // all have the same weight
                  _points[offset + 0] = Point(a[i], a[i], a[i]);
                  _points[offset + 1] = Point(a[i],    b, a[i]);
                  _points[offset + 2] = Point(   b, a[i], a[i]);
                  _points[offset + 3] = Point(a[i], a[i],    b);
                                            
                  // These 4 points all have the same weights 
                  for (unsigned int j=0; j<4; ++j)
                    _weights[offset + j] = w[i];
                } // end for


              {
                // The third set contains 6 points and is formed a little differently
                const unsigned int offset = 8;
                const Real b = 0.5*(1. - 2.*a[2]);

                // Here are the permutations.  Order of these is not important,
                // all have the same weight
                _points[offset + 0] = Point(b   ,    b, a[2]);
                _points[offset + 1] = Point(b   , a[2], a[2]);
                _points[offset + 2] = Point(a[2], a[2],    b);
                _points[offset + 3] = Point(a[2],    b, a[2]);
                _points[offset + 4] = Point(   b, a[2],    b);
                _points[offset + 5] = Point(a[2],    b,    b);
                  
                // These 6 points all have the same weights 
                for (unsigned int j=0; j<6; ++j)
                  _weights[offset + j] = w[2];
              }

              
              // Original rule by Keast, unsuitable because it has points on the
              // reference element boundary!
              //              _points.resize(15);
              //              _weights.resize(15);
                    
              //              _points[0](0) = 0.25;
              //              _points[0](1) = 0.25;
              //              _points[0](2) = 0.25;

              //              {
              //                const Real a = 0.;
              //                const Real b = 0.333333333333333333333333333333333333333;
                
              //                _points[1](0) = a;
              //                _points[1](1) = b;
              //                _points[1](2) = b;
                
              //                _points[2](0) = b;
              //                _points[2](1) = a;
              //                _points[2](2) = b;
                
              //                _points[3](0) = b;
              //                _points[3](1) = b;
              //                _points[3](2) = a;
                
              //                _points[4](0) = b;
              //                _points[4](1) = b;
              //                _points[4](2) = b;
              //              }
              //              {
              //                const Real a = 0.7272727272727272727272727272727272727272727272727272727;
              //                const Real b = 0.0909090909090909090909090909090909090909090909090909091;
                
              //                _points[5](0) = a;
              //                _points[5](1) = b;
              //                _points[5](2) = b;
                
              //                _points[6](0) = b;
              //                _points[6](1) = a;
              //                _points[6](2) = b;
                
              //                _points[7](0) = b;
              //                _points[7](1) = b;
              //                _points[7](2) = a;
                
              //                _points[8](0) = b;
              //                _points[8](1) = b;
              //                _points[8](2) = b;
              //              }
              //              {
              //                const Real a = 0.066550153573664;
              //                const Real b = 0.433449846426336;
                
              //                _points[9](0) = b;
              //                _points[9](1) = a;
              //                _points[9](2) = a;
                
              //                _points[10](0) = a;
              //                _points[10](1) = a;
              //                _points[10](2) = b;
                
              //                _points[11](0) = a;
              //                _points[11](1) = b;
              //                _points[11](2) = b;
                
              //                _points[12](0) = b;
              //                _points[12](1) = a;
              //                _points[12](2) = b;
                
              //                _points[13](0) = b;
              //                _points[13](1) = b;
              //                _points[13](2) = a;
                
              //                _points[14](0) = a;
              //                _points[14](1) = b;
              //                _points[14](2) = a;             
              //              }
              
              //              _weights[0]  = 0.030283678097089;
              //              _weights[1]  = 0.006026785714286;
              //              _weights[2]  = _weights[1];
              //              _weights[3]  = _weights[1];
              //              _weights[4]  = _weights[1];
              //              _weights[5]  = 0.011645249086029;
              //              _weights[6]  = _weights[5];
              //              _weights[7]  = _weights[5];
              //              _weights[8]  = _weights[5];
              //              _weights[9]  = 0.010949141561386;
              //              _weights[10] = _weights[9];
              //              _weights[11] = _weights[9];
              //              _weights[12] = _weights[9];
              //              _weights[13] = _weights[9];
              //              _weights[14] = _weights[9];
              
              return;
            }



            
            // This rule is originally from Keast:
            //    Patrick Keast,
            //    Moderate Degree Tetrahedral Quadrature Formulas,
            //    Computer Methods in Applied Mechanics and Engineering,
            //    Volume 55, Number 3, May 1986, pages 339-348.
            //
            // It is accurate on 6th-degree polynomials and has 24 points
            // vs. 64 for the comparable conical product rule.
            //
            // Values copied 24th June 2008 from:
            // http://people.scs.fsu.edu/~burkardt/f_src/keast/keast.f90
          case SIXTH:
            {
              _points.resize (24);
              _weights.resize(24);

              // The raw data for the quadrature rule.
              const Real p[4][4] = {
                {0.356191386222544953e+00 , 0.214602871259151684e+00 ,                       0., 0.00665379170969464506e+00}, // 4
                {0.877978124396165982e+00 , 0.0406739585346113397e+00,                       0., 0.00167953517588677620e+00}, // 4
                {0.0329863295731730594e+00, 0.322337890142275646e+00 ,                       0., 0.00922619692394239843e+00}, // 4
                {0.0636610018750175299e+00, 0.269672331458315867e+00 , 0.603005664791649076e+00, 0.00803571428571428248e+00}  // 12    
              };

              
              // Now call the keast routine to generate _points and _weights
              keast_rule(p, 4);

              return;
            }


            // Keast's 31 point, 7th-order rule contains points on the reference
            // element boundary, so we've decided not to include it here.
            //
            // Keast's 8th-order rule has 45 points.  and a negative
            // weight, so if you've explicitly disallowed such rules
            // you will fall through to the conical product rules
            // below.
          case SEVENTH:
          case EIGHTH:
            {
              if (allow_rules_with_negative_weights)
                {
                  _points.resize (45);
                  _weights.resize(45);

                  // The raw data for the quadrature rule.
                  const Real p[7][4] = {
                    {0.250000000000000000e+00,                        0.,                        0.,   -0.393270066412926145e-01},  // 1
                    {0.617587190300082967e+00,  0.127470936566639015e+00,                        0.,    0.408131605934270525e-02},  // 4
                    {0.903763508822103123e+00,  0.320788303926322960e-01,                        0.,    0.658086773304341943e-03},  // 4
                    {0.497770956432810185e-01,                        0.,  0.450222904356718978e+00,    0.438425882512284693e-02},  // 6
                    {0.183730447398549945e+00,                        0.,  0.316269552601450060e+00,    0.138300638425098166e-01},  // 6
                    {0.231901089397150906e+00,  0.229177878448171174e-01,  0.513280033360881072e+00,    0.424043742468372453e-02},  // 12
                    {0.379700484718286102e-01,  0.730313427807538396e+00,  0.193746475248804382e+00,    0.223873973961420164e-02}   // 12
                  };

              
                  // Now call the keast routine to generate _points and _weights
                  keast_rule(p, 7);

                  return;
                } // end if (allow_rules_with_negative_weights)
              // Note: if !allow_rules_with_negative_weights, fall through to next case.
            }

            // Fall back on Grundmann-Moller or Conical Product rules at high orders.
          default:
            {
              if ((allow_rules_with_negative_weights) && (_order + 2*p < 34))
                {
                  // The Grundmann-Moller rules are defined to arbitrary order and
                  // can have significantly fewer evaluation points than concial product
                  // rules.  If you allow rules with negative weights, the GM rules
                  // will be more efficient for degree up to 33 (for degree 34 and
                  // higher, CP is more efficient!) but may be more susceptible
                  // to round-off error.  Safest is to disallow rules with negative
                  // weights, but this decision should be made on a case-by-case basis.
                  QGrundmann_Moller gm_rule(3, _order);
                  gm_rule.init(_type, p);

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

                  return;
                }

              else
                {
                  // The following quadrature rules are generated as
                  // conical products.  These tend to be non-optimal
                  // (use too many points, cluster points in certain
                  // regions of the domain) but they are quite easy to
                  // automatically generate using a 1D Gauss rule on
                  // [0,1] and two 1D Jacobi-Gauss rules on [0,1].
                  QConical conical_rule(3, _order);
                  conical_rule.init(_type, p);

                  // Swap points and weights with the about-to-be destroyed rule.
                  _points.swap (conical_rule.get_points() );
                  _weights.swap(conical_rule.get_weights());
                  
                  return;
                }
            }
          }
      } // end case TET4,TET10

      
            
      //---------------------------------------------
      // Prism quadrature rules
    case PRISM6:
    case PRISM15:
    case PRISM18:
      {
        // We compute the 3D quadrature rule as a tensor
        // product of the 1D quadrature rule and a 2D
        // triangle quadrature rule
            
        QGauss q1D(1,_order);
        QGauss q2D(2,_order);

        // Initialize 
        q1D.init(EDGE2,p);
        q2D.init(TRI3,p);

        tensor_product_prism(q1D, q2D);
        
        return;
      }
      

      
      //---------------------------------------------
      // Pyramid
    case PYRAMID5:
      {
        // We compute the Pyramid rule as a conical product of a
        // Jacobi rule with alpha==2 on the interval [0,1] two 1D
        // Gauss rules with suitably modified points.  The idea comes
        // from: Stroud, A.H. 'Approximate Calculation of Multiple
        // Integrals.'
        QConical conical_rule(3, _order);
        conical_rule.init(_type, p);

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

      }


      
      //---------------------------------------------
      // Unsupported type
    default:
      {
        std::cerr << 'ERROR: Unsupported type: ' << _type << std::endl;
        libmesh_error();
      }
    }

  libmesh_error();

  return;
  
#endif
}
 

void QGauss::keast_rule (const Realrule_data[][4], const unsigned intn_pts) [private]The Keast rule is for tets. It takes permutation points and weights in a specific format as input and fills the pre-sized _points and _weights vectors. This function is only used internally by the TET geometric elements.

Definition at line 32 of file quadrature_gauss.C.

References QBase::_points, QBase::_weights, and QBase::w().

Referenced by init_3D().

{
  // Like the Dunavant rule, the input data should have 4 columns.  These columns
  // have the following format and implied permutations (w=weight).
  // {a, 0, 0, w} = 1-permutation  (a,a,a)
  // {a, b, 0, w} = 4-permutation  (a,b,b), (b,a,b), (b,b,a), (b,b,b)
  // {a, 0, b, w} = 6-permutation  (a,a,b), (a,b,b), (b,b,a), (b,a,b), (b,a,a), (a,b,a)
  // {a, b, c, w} = 12-permutation (a,a,b), (a,a,c), (b,a,a), (c,a,a), (a,b,a), (a,c,a)
  //                               (a,b,c), (a,c,b), (b,a,c), (b,c,a), (c,a,b), (c,b,a)

  // Always insert into the points & weights vector relative to the offset 
  unsigned int offset=0;
  
  
  for (unsigned int p=0; p<n_pts; ++p)
    {

      // There must always be a non-zero entry to start the row
      libmesh_assert(rule_data[p][0] != static_cast<Real>(0.0));
      
      // A zero weight may imply you did not set up the raw data correctly
      libmesh_assert(rule_data[p][3] != static_cast<Real>(0.0));

      // What kind of point is this?
      // One non-zero entry in first 3 cols   ? 1-perm (centroid) point = 1
      // Two non-zero entries in first 3 cols ? 3-perm point            = 3
      // Three non-zero entries               ? 6-perm point            = 6
      unsigned int pointtype=1;

      if (rule_data[p][1] != static_cast<Real>(0.0))      
        {
          if (rule_data[p][2] != static_cast<Real>(0.0))
            pointtype = 12;
          else
            pointtype = 4;
        }
      else
        {
          // The second entry is zero.  What about the third?
          if (rule_data[p][2] != static_cast<Real>(0.0))
            pointtype = 6;
        }

      
      switch (pointtype)
        {
        case 1:
          {
            // Be sure we have enough space to insert this point
            libmesh_assert(offset + 0 < _points.size());

            const Real a = rule_data[p][0];
            
            // The point has only a single permutation (the centroid!)
            _points[offset  + 0] = Point(a,a,a);

            // The weight is always the last entry in the row.
            _weights[offset + 0] = rule_data[p][3];

            offset += pointtype;
            break;
          }

        case 4:
          {
            // Be sure we have enough space to insert these points
            libmesh_assert(offset + 3 < _points.size());

            const Real a = rule_data[p][0];
            const Real b = rule_data[p][1];
            const Real w = rule_data[p][3];
            
            // Here it's understood the second entry is to be used twice, and
            // thus there are three possible permutations.
            _points[offset + 0] = Point(a,b,b);
            _points[offset + 1] = Point(b,a,b);
            _points[offset + 2] = Point(b,b,a);
            _points[offset + 3] = Point(b,b,b);

            for (unsigned int j=0; j<pointtype; ++j)
              _weights[offset + j] = w;
            
            offset += pointtype;
            break;
          }

        case 6:
          {
            // Be sure we have enough space to insert these points
            libmesh_assert(offset + 5 < _points.size());

            const Real a = rule_data[p][0];
            const Real b = rule_data[p][2];
            const Real w = rule_data[p][3];
            
            // Three individual entries with six permutations.
            _points[offset + 0] = Point(a,a,b);
            _points[offset + 1] = Point(a,b,b);
            _points[offset + 2] = Point(b,b,a);
            _points[offset + 3] = Point(b,a,b);
            _points[offset + 4] = Point(b,a,a);
            _points[offset + 5] = Point(a,b,a);

            for (unsigned int j=0; j<pointtype; ++j)
              _weights[offset + j] = w;
            
            offset += pointtype;
            break;
          }

          
        case 12:
          {
            // Be sure we have enough space to insert these points
            libmesh_assert(offset + 11 < _points.size());

            const Real a = rule_data[p][0];
            const Real b = rule_data[p][1];
            const Real c = rule_data[p][2];
            const Real w = rule_data[p][3];
            
            // Three individual entries with six permutations.
            _points[offset + 0] = Point(a,a,b);  _points[offset + 6]  = Point(a,b,c);
            _points[offset + 1] = Point(a,a,c);  _points[offset + 7]  = Point(a,c,b);
            _points[offset + 2] = Point(b,a,a);  _points[offset + 8]  = Point(b,a,c);
            _points[offset + 3] = Point(c,a,a);  _points[offset + 9]  = Point(b,c,a);
            _points[offset + 4] = Point(a,b,a);  _points[offset + 10] = Point(c,a,b);
            _points[offset + 5] = Point(a,c,a);  _points[offset + 11] = Point(c,b,a);

            for (unsigned int j=0; j<pointtype; ++j)
              _weights[offset + j] = w;
            
            offset += pointtype;
            break;
          }

        default:
          {
            std::cerr << 'Don't know what to do with this many permutation points!' << std::endl;
            libmesh_error();
          }
        }
      
    }
  
}
 

static unsigned int ReferenceCounter::n_objects () [inline, static, inherited]Prints the number of outstanding (created, but not yet destroyed) objects.

Definition at line 76 of file reference_counter.h.

References ReferenceCounter::_n_objects.

Referenced by System::read_serialized_blocked_dof_objects(), and System::write_serialized_blocked_dof_objects().

  { return _n_objects; }
 

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

the number of points associated with the quadrature rule.

Definition at line 115 of file quadrature.h.

References QBase::_points.

Referenced by FEBase::coarsened_dof_values(), QConical::conical_product_pyramid(), QConical::conical_product_tet(), QConical::conical_product_tri(), FEMSystem::eulerian_residual(), InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), FEMSystem::mass_residual(), and QBase::print_info().

    { libmesh_assert (!_points.empty()); return _points.size(); }
 

void QBase::print_info (std::ostream &os = std::cout) const [inline, inherited]Prints information relevant to the quadrature rule.

Definition at line 350 of file quadrature.h.

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

Referenced by operator<<().

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

  os << 'N_Q_Points=' << this->n_points() << std::endl << std::endl;
  for (unsigned int qp=0; qp<this->n_points(); qp++)
    {
      os << ' Point ' << qp << ':
         << '  '
         << _points[qp]
         << ' Weight:'
         << '  w=' << _weights[qp] << ' << std::endl;
    }
}
 

void ReferenceCounter::print_info () [static, inherited]Prints the reference information to std::cout.

Definition at line 83 of file reference_counter.C.

References ReferenceCounter::get_info().

{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
  
  std::cout << ReferenceCounter::get_info();
  
#endif
}
 

Point QBase::qp (const unsigned inti) const [inline, inherited]Returns:

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

Definition at line 148 of file quadrature.h.

References QBase::_points.

Referenced by QConical::conical_product_pyramid(), QConical::conical_product_tet(), QConical::conical_product_tri(), and QBase::print_info().

    { libmesh_assert (i < _points.size()); return _points[i]; }
 

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

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

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

{
  // Make sure we are in 1D
  libmesh_assert(_dim == 1);
  
  // Make sure that we have sane ranges 
  libmesh_assert(new_range.second > new_range.first);
  libmesh_assert(old_range.second > old_range.first);

  // Make sure there are some points
  libmesh_assert(_points.size() > 0);

  // We're mapping from old_range -> new_range 
  for (unsigned int i=0; i<_points.size(); i++)
    {
      _points[i](0) =
        (_points[i](0) - old_range.first) *
        (new_range.second - new_range.first) /
        (old_range.second - old_range.first) +
        new_range.first;
    }

  // Compute the scale factor and scale the weights
  const Real scfact = (new_range.second - new_range.first) /
                      (old_range.second - old_range.first);

  for (unsigned int i=0; i<_points.size(); i++)
    _weights[i] *= scfact;
}
 

QuadratureType QGauss::type () const [inline, virtual]Returns:

QGAUSS

Implements QBase.

Definition at line 58 of file quadrature_gauss.h.

References libMeshEnums::QGAUSS.

{ return QGAUSS; }
 

Real QBase::w (const unsigned inti) const [inline, inherited]Returns:

the $ i^{th} $ quadrature weight.

Definition at line 154 of file quadrature.h.

References QBase::_weights.

Referenced by QConical::conical_product_pyramid(), QConical::conical_product_tet(), QConical::conical_product_tri(), init_3D(), keast_rule(), QMonomial::kim_rule(), and QMonomial::stroud_rule().

    { libmesh_assert (i < _weights.size()); return _weights[i]; }
 

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

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

Member Data Documentation

 

ReferenceCounter::Counts ReferenceCounter::_counts [static, protected, inherited]Actually holds the data.

Definition at line 110 of file reference_counter.h.

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

Threads::spin_mutex ReferenceCounter::_mutex [static, protected, inherited]Mutual exclusion object to enable thread-safe reference counting.

Definition at line 123 of file reference_counter.h.  

Threads::atomic< unsigned int > ReferenceCounter::_n_objects [static, protected, inherited]The number of objects. Print the reference count information when the number returns to 0.

Definition at line 118 of file reference_counter.h.

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

std::cerr<< 'ERROR: Seems as if this quadrature rule' << std::endl << ' is not implemented for 2D.' << std::endl; libmesh_error(); }#endif virtual void init_3D (const ElemType, unsigned int =0)#ifndef DEBUG {}#else { std::cerr << 'ERROR: Seems as if this quadrature rule' << std::endl << ' is not implemented for 3D.' << std::endl; libmesh_error(); }#endif void tensor_product_quad (const QBase& q1D); void tensor_product_hex (const QBase& q1D); void tensor_product_prism (const QBase& q1D, const QBase& q2D); const unsigned int _dim; const Order _order; ElemType _type; unsigned int _p_level; std::vector<Point> QBase::_points [protected, inherited]

Definition at line 320 of file quadrature.h.

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

std::vector<Real> QBase::_weights [protected, inherited]The value of the quadrature weights.

Definition at line 325 of file quadrature.h.

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

bool 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 203 of file quadrature.h.

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

 

Author

Generated automatically by Doxygen for libMesh from the source code.


 

Index

NAME
SYNOPSIS
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
typedef std::map<std::string, std::pair<unsigned int, unsigned int> > ReferenceCounter::Counts [protected, inherited]Data structure to log the information. The log is identified by the class name.
Constructor & Destructor Documentation
QGauss::QGauss (const unsigned int_dim, const Order_order = INVALID_ORDER) [inline]Constructor. Declares the order of the quadrature rule.
QGauss::~QGauss () [inline]Destructor.
Member Function Documentation
AutoPtr< QBase > 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 AutoPtr<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.
AutoPtr< QBase > 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 AutoPtr<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.
void QGauss::dunavant_rule (const Realrule_data[][4], const unsigned intn_pts) [private]The Dunavant rule is for triangles. It takes permutation points and weights in a specific format as input and fills the pre-sized _points and _weights vectors. This function is only used internally by the TRI geometric elements.
void QGauss::dunavant_rule2 (const Real *wts, const Real *a, const Real *b, const unsigned int *permutation_ids, const unsigned intn_wts) [private]
unsigned int QBase::get_dim () const [inline, inherited]Returns:
ElemType QBase::get_elem_type () const [inline, inherited]Returns:
std::string ReferenceCounter::get_info () [static, inherited]Gets a string containing the reference information.
Order QBase::get_order () const [inline, inherited]Returns:
unsigned int QBase::get_p_level () const [inline, inherited]Returns:
const std::vector<Point>& QBase::get_points () const [inline, inherited]Returns:
std::vector<Point>& QBase::get_points () [inline, inherited]Returns:
const std::vector<Real>& QBase::get_weights () const [inline, inherited]Returns:
std::vector<Real>& QBase::get_weights () [inline, inherited]Returns:
void 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.
void 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.
void QBase::init (const ElemType_type = INVALID_ELEM, unsigned intp_level = 0) [inherited]Initializes the data structures to contain a quadrature rule for an object of type type.
void QBase::init_0D (const ElemType_type = INVALID_ELEM, unsigned intp_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.
void QGauss::init_1D (const ElemType_type = INVALID_ELEM, unsigned intp_level = 0) [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.
void QGauss::init_2D (const ElemType = INVALID_ELEM, unsigned int = 0) [private, virtual]Initializes the 2D quadrature rule by filling the points and weights vectors with the appropriate values. The order of the rule will be defined by the implementing class. Should not be pure virtual since a derived quadrature rule may only be defined in 1D. If not redefined, gives an error (when DEBUG defined) when called.
void QGauss::init_3D (const ElemType_type = INVALID_ELEM, unsigned intp_level = 0) [private]
void QGauss::keast_rule (const Realrule_data[][4], const unsigned intn_pts) [private]The Keast rule is for tets. It takes permutation points and weights in a specific format as input and fills the pre-sized _points and _weights vectors. This function is only used internally by the TET geometric elements.
static unsigned int ReferenceCounter::n_objects () [inline, static, inherited]Prints the number of outstanding (created, but not yet destroyed) objects.
unsigned int QBase::n_points () const [inline, inherited]Returns:
void QBase::print_info (std::ostream &os = std::cout) const [inline, inherited]Prints information relevant to the quadrature rule.
void ReferenceCounter::print_info () [static, inherited]Prints the reference information to std::cout.
Point QBase::qp (const unsigned inti) const [inline, inherited]Returns:
void 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.
QuadratureType QGauss::type () const [inline, virtual]Returns:
Real QBase::w (const unsigned inti) const [inline, inherited]Returns:
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.
Member Data Documentation
ReferenceCounter::Counts ReferenceCounter::_counts [static, protected, inherited]Actually holds the data.
Threads::spin_mutex ReferenceCounter::_mutex [static, protected, inherited]Mutual exclusion object to enable thread-safe reference counting.
Threads::atomic< unsigned int > ReferenceCounter::_n_objects [static, protected, inherited]The number of objects. Print the reference count information when the number returns to 0.
std::cerr<< 'ERROR: Seems as if this quadrature rule' << std::endl << ' is not implemented for 2D.' << std::endl; libmesh_error(); }#endif virtual void init_3D (const ElemType, unsigned int =0)#ifndef DEBUG {}#else { std::cerr << 'ERROR: Seems as if this quadrature rule' << std::endl << ' is not implemented for 3D.' << std::endl; libmesh_error(); }#endif void tensor_product_quad (const QBase& q1D); void tensor_product_hex (const QBase& q1D); void tensor_product_prism (const QBase& q1D, const QBase& q2D); const unsigned int _dim; const Order _order; ElemType _type; unsigned int _p_level; std::vector<Point> QBase::_points [protected, inherited]
std::vector<Real> QBase::_weights [protected, inherited]The value of the quadrature weights.
bool 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.
Author

This document was created by man2html, using the manual pages.
Time: 21:53:32 GMT, April 16, 2011