libMesh::QGrundmann_Moller Class Reference

#include <quadrature_gm.h>

Inheritance diagram for libMesh::QGrundmann_Moller:

List of all members.

Public Member Functions

 QGrundmann_Moller (const unsigned int _dim, const Order _order=INVALID_ORDER)
 ~QGrundmann_Moller ()
QuadratureType type () const
ElemType get_elem_type () const
unsigned int get_p_level () const
unsigned int n_points () const
unsigned int get_dim () const
const std::vector< Point > & get_points () const
std::vector< Point > & get_points ()
const std::vector< Real > & get_weights () const
std::vector< Real > & get_weights ()
Point qp (const unsigned int i) const
Real w (const unsigned int i) const
void init (const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Order get_order () const
void print_info (std::ostream &os=libMesh::out) const
void scale (std::pair< Real, Real > old_range, std::pair< Real, Real > new_range)
virtual bool shapes_need_reinit ()

Static Public Member Functions

static AutoPtr< QBasebuild (const std::string &name, const unsigned int _dim, const Order _order=INVALID_ORDER)
static AutoPtr< QBasebuild (const QuadratureType _qt, const unsigned int _dim, const Order _order=INVALID_ORDER)
static void print_info (std::ostream &out=libMesh::out)
static std::string get_info ()
static unsigned int n_objects ()
static void enable_print_counter_info ()
static void disable_print_counter_info ()

Public Attributes

bool allow_rules_with_negative_weights

Protected Types

typedef std::map< std::string,
std::pair< unsigned int,
unsigned int > > 
Counts

Protected Member Functions

virtual void init_0D (const ElemType type=INVALID_ELEM, unsigned int p_level=0)
virtual void init_2D (const ElemType, unsigned int=0)
void increment_constructor_count (const std::string &name)
void increment_destructor_count (const std::string &name)

Protected Attributes

libMesh::err<< "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{libMesh::err<< "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
static bool _enable_print_counter = true

Private Member Functions

void init_1D (const ElemType, unsigned int=0)
void init_3D (const ElemType _type=INVALID_ELEM, unsigned int p_level=0)
void gm_rule (unsigned int s)
void compose_all (unsigned int s, unsigned int p, std::vector< std::vector< unsigned int > > &result)

Friends

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

Detailed Description

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

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

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

where V is the volume of the reference element, does grow quickly. (A rule with all positive has has an amplification factor of 1.0 by definition.)

  s  | d     | N. CP        | N. GM   | % neg wts | amp. factor
-----------------------------------------------------------------
  0  | 1     |              | 1       |           |
  1  | 2-3   |              | 5       |           |
  2  | 4-5   |              | 15      |           |
  3  | 6-7   |              | 35      | 31.43     |   11.94
  4  | 8-9   |  5^3=125     | 70      | 34.29     |   25.35
  5  | 10-11 |  6^3=216     | 126     | 36.51     |   54.14
  6  | 12-13 |  7^3=343     | 210     | 38.10     |  116.30
  7  | 14-15 |  8^3=512     | 330     | 39.39     |  251.10
  8  | 16-17 |  9^3=729     | 495     | 40.40     |  544.68
  9  | 18-19 | 10^3=1,000   | 715     | 41.26     | 1186.16
 10  | 20-21 | 11^3=1,331   | 1,001   | 41.96     | 2591.97
 11  | 22-23 | 12^3=1,728   | 1,365   | 42.56     | 5680.75
 ...
 16  | 32-33 | 17^3=4,913   | 4,845   |
 17  | 34-35 | 18^3=5,832   | 5,985   | <= Cross-over point, CP has fewer points for d >= 34
 18  | 36-37 | 19^3=6,859   | 7,315   |
 ...
 21  | 42-43 | 22^3=10,648  | 12,650  |

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

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

Author:
John W. Peterson, 2008

Definition at line 100 of file quadrature_gm.h.


Member Typedef Documentation

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

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

Definition at line 113 of file reference_counter.h.


Constructor & Destructor Documentation

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

Constructor. Declares the order of the quadrature rule.

Definition at line 32 of file quadrature_gm.C.

00033                                                     : QBase(d,o)
00034 {
00035 }

libMesh::QGrundmann_Moller::~QGrundmann_Moller (  ) 

Destructor.

Definition at line 39 of file quadrature_gm.C.

00040 {
00041 }


Member Function Documentation

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

Builds a specific quadrature rule, identified through the QuadratureType. An 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 48 of file quadrature_build.C.

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

00051 {
00052   switch (_qt)
00053     {
00054 
00055     case QCLOUGH:
00056       {
00057 #ifdef DEBUG
00058         if (_order > TWENTYTHIRD)
00059           {
00060             libMesh::out << "WARNING: Clough quadrature implemented" << std::endl
00061                           << " up to TWENTYTHIRD order." << std::endl;
00062           }
00063 #endif
00064 
00065         AutoPtr<QBase> ap(new QClough(_dim, _order));
00066         return ap;
00067       }
00068 
00069     case QGAUSS:
00070       {
00071 
00072 #ifdef DEBUG
00073         if (_order > FORTYTHIRD)
00074           {
00075             libMesh::out << "WARNING: Gauss quadrature implemented" << std::endl
00076                           << " up to FORTYTHIRD order." << std::endl;
00077           }
00078 #endif
00079 
00080         AutoPtr<QBase> ap(new QGauss(_dim, _order));
00081         return ap;
00082       }
00083 
00084     case QJACOBI_1_0:
00085       {
00086 
00087 #ifdef DEBUG
00088         if (_order > TWENTYTHIRD)
00089           {
00090             libMesh::out << "WARNING: Jacobi(1,0) quadrature implemented" << std::endl
00091                           << " up to TWENTYTHIRD order." << std::endl;
00092           }
00093 
00094         if (_dim > 1)
00095           {
00096             libMesh::out << "WARNING: Jacobi(1,0) quadrature implemented" << std::endl
00097                           << " in 1D only." << std::endl;
00098           }
00099 #endif
00100 
00101         AutoPtr<QBase> ap(new QJacobi(_dim, _order, 1, 0));
00102         return ap;
00103       }
00104 
00105     case QJACOBI_2_0:
00106       {
00107 
00108 #ifdef DEBUG
00109         if (_order > TWENTYTHIRD)
00110           {
00111             libMesh::out << "WARNING: Jacobi(2,0) quadrature implemented" << std::endl
00112                           << " up to TWENTYTHIRD order." << std::endl;
00113           }
00114 
00115         if (_dim > 1)
00116           {
00117             libMesh::out << "WARNING: Jacobi(2,0) quadrature implemented" << std::endl
00118                           << " in 1D only." << std::endl;
00119           }
00120 #endif
00121 
00122         AutoPtr<QBase> ap(new QJacobi(_dim, _order, 2, 0));
00123         return ap;
00124       }
00125 
00126     case QSIMPSON:
00127       {
00128 
00129 #ifdef DEBUG
00130         if (_order > THIRD)
00131           {
00132             libMesh::out << "WARNING: Simpson rule provides only" << std::endl
00133                           << " THIRD order!" << std::endl;
00134           }
00135 #endif
00136 
00137         AutoPtr<QBase> ap(new QSimpson(_dim));
00138         return ap;
00139       }
00140 
00141     case QTRAP:
00142       {
00143 
00144 #ifdef DEBUG
00145         if (_order > FIRST)
00146           {
00147             libMesh::out << "WARNING: Trapezoidal rule provides only" << std::endl
00148                           << " FIRST order!" << std::endl;
00149           }
00150 #endif
00151 
00152         AutoPtr<QBase> ap(new QTrap(_dim));
00153         return ap;
00154       }
00155 
00156 
00157     default:
00158       {
00159         libMesh::err << "ERROR: Bad qt=" << _qt << std::endl;
00160         libmesh_error();
00161       }
00162     }
00163 
00164 
00165   libmesh_error();
00166   AutoPtr<QBase> ap(NULL);
00167   return ap;
00168 }

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

Builds a specific quadrature rule, identified through the name string. An 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 37 of file quadrature_build.C.

References libMesh::Utility::string_to_enum< QuadratureType >().

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

00040 {
00041   return QBase::build (Utility::string_to_enum<QuadratureType> (type),
00042                        _dim,
00043                        _order);
00044 }

void libMesh::QGrundmann_Moller::compose_all ( unsigned int  s,
unsigned int  p,
std::vector< std::vector< unsigned int > > &  result 
) [private]

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

Definition at line 144 of file quadrature_gm.C.

Referenced by gm_rule().

00147 {
00148   // Clear out results remaining from previous calls
00149   result.clear();
00150 
00151   // Allocate storage for a workspace.  The workspace will periodically
00152   // be copied into the result container.
00153   std::vector<unsigned int> workspace(p);
00154 
00155   // The first result is always (s,0,...,0)
00156   workspace[0] = s;
00157   result.push_back(workspace);
00158 
00159   // the value of the first non-zero entry
00160   unsigned int head_value=s;
00161 
00162   // When head_index=-1, it refers to "off the front" of the array.  Therefore,
00163   // this needs to be a regular int rather than unsigned.  I initially tried to
00164   // do this with head_index unsigned and an else statement below, but then there
00165   // is the special case: (1,0,...,0) which does not work correctly.
00166   int head_index = -1;
00167 
00168   // At the end, all the entries will be in the final slot of workspace
00169   while (workspace.back() != s)
00170     {
00171       // Uncomment for debugging
00172       //libMesh::out << "previous head_value=" << head_value << " -> ";
00173 
00174       // If the previous head value is still larger than 1, reset the index
00175       // to "off the front" of the array
00176       if (head_value > 1)
00177         head_index = -1;
00178 
00179       // Either move the index onto the front of the array or on to
00180       // the next value.
00181       head_index++;
00182 
00183       // Get current value of the head entry
00184       head_value = workspace[head_index];
00185 
00186       // Uncomment for debugging
00187       //std::copy(workspace.begin(), workspace.end(), std::ostream_iterator<int>(libMesh::out, " "));
00188       //libMesh::out << ", head_index=" << head_index;
00189       //libMesh::out << ", head_value=" << head_value << " -> ";
00190 
00191       // Put a zero into the head_index of the array.  If head_index==0,
00192       // this will be overwritten in the next line with head_value-1.
00193       workspace[head_index] = 0;
00194 
00195       // The initial entry gets the current head value, minus 1.
00196       // If head_value > 1, the next loop iteration will start back
00197       // at workspace[0] again.
00198       libmesh_assert_greater (head_value, 0);
00199       workspace[0] = head_value - 1;
00200 
00201       // Increment the head+1 value
00202       workspace[head_index+1] += 1;
00203 
00204       // Save this composition in the results
00205       result.push_back(workspace);
00206 
00207       // Uncomment for debugging
00208       //std::copy(workspace.begin(), workspace.end(), std::ostream_iterator<int>(libMesh::out, " "));
00209       //libMesh::out<<"\n";
00210     }
00211 }

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

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

00107 {
00108   _enable_print_counter = false;
00109   return;
00110 }

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

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

Definition at line 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

00101 {
00102   _enable_print_counter = true;
00103   return;
00104 }

unsigned int libMesh::QBase::get_dim (  )  const [inline, inherited]
ElemType libMesh::QBase::get_elem_type (  )  const [inline, inherited]
Returns:
the current element type we're set up for

Definition at line 104 of file quadrature.h.

00105     { return _type; }

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

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

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

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

00048 {
00049 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
00050 
00051   std::ostringstream oss;
00052 
00053   oss << '\n'
00054       << " ---------------------------------------------------------------------------- \n"
00055       << "| Reference count information                                                |\n"
00056       << " ---------------------------------------------------------------------------- \n";
00057 
00058   for (Counts::iterator it = _counts.begin();
00059        it != _counts.end(); ++it)
00060     {
00061       const std::string name(it->first);
00062       const unsigned int creations    = it->second.first;
00063       const unsigned int destructions = it->second.second;
00064 
00065       oss << "| " << name << " reference count information:\n"
00066           << "|  Creations:    " << creations    << '\n'
00067           << "|  Destructions: " << destructions << '\n';
00068     }
00069 
00070   oss << " ---------------------------------------------------------------------------- \n";
00071 
00072   return oss.str();
00073 
00074 #else
00075 
00076   return "";
00077 
00078 #endif
00079 }

Order libMesh::QBase::get_order (  )  const [inline, inherited]
Returns:
the order of the quadrature rule.

Definition at line 169 of file quadrature.h.

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

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

unsigned int libMesh::QBase::get_p_level (  )  const [inline, inherited]
Returns:
the current p refinement level we're initialized with

Definition at line 110 of file quadrature.h.

00111     { return _p_level; }

std::vector<Point>& libMesh::QBase::get_points (  )  [inline, inherited]
Returns:
a std::vector containing the quadrature point locations on a reference object as a writeable reference.

Definition at line 135 of file quadrature.h.

References libMesh::QBase::_points.

00135 { return _points;  }

std::vector<Real>& libMesh::QBase::get_weights (  )  [inline, inherited]
Returns:
a std::vector containing the quadrature weights.

Definition at line 145 of file quadrature.h.

References libMesh::QBase::_weights.

00145 { return _weights; }

void libMesh::QGrundmann_Moller::gm_rule ( unsigned int  s  )  [private]

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

Definition at line 47 of file quadrature_gm.C.

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

Referenced by init_3D().

00048 {
00049   // A GM rule of index s can integrate polynomials of degree 2*s+1 exactly
00050   const unsigned int degree = 2*s+1;
00051 
00052   // Here we are considering only tetrahedra rules, so dim==3
00053   const unsigned int dim = 3;
00054 
00055   // The number of points for rule of index s is
00056   // (dim+1+s)! / (dim+1)! / s!
00057   // In 3D, this is = 1/24 * P_{i=1}^4 (s+i)
00058   const unsigned int n_pts = (s+4)*(s+3)*(s+2)*(s+1) / 24;
00059   //libMesh::out << "n_pts=" << n_pts << std::endl;
00060 
00061   // Allocate space for points and weights
00062   _points.resize(n_pts);
00063   _weights.resize(n_pts);
00064 
00065   // (-1)^i -> This one flips sign at each iteration of the i-loop below.
00066   int one_pm=1;
00067 
00068   // Where we store all the integer point compositions/permutations
00069   std::vector<std::vector<unsigned int> > permutations;
00070 
00071   // Index into the vector where we should start adding the next round of points/weights
00072   std::size_t offset=0;
00073 
00074   // Implement the GM formula 4.1 on page 286 of the paper
00075   for (unsigned int i=0; i<=s; ++i)
00076     {
00077       // Get all the ordered compositions (and their permutations)
00078       // of |beta| = s-i into dim+1=4 parts
00079       compose_all(s-i, dim+1, permutations);
00080       //libMesh::out << "n. permutations=" << permutations.size() << std::endl;
00081 
00082       for (unsigned int p=0; p<permutations.size(); ++p)
00083         {
00084           // We use the first dim=3 entries of each permutation to
00085           // construct an integration point.
00086           for (unsigned int j=0; j<3; ++j)
00087             _points[offset+p](j) =
00088               static_cast<Real>(2.*permutations[p][j] + 1.) /
00089               static_cast<Real>(  degree + dim - 2.*i     );
00090         }
00091 
00092       // Compute the weight for this i, being careful to avoid overflow.
00093       // This technique is borrowed from Burkardt's code as well.
00094       // Use once for each of the points obtained from the permutations array.
00095       Real weight = one_pm;
00096 
00097       // This for loop needs to run for dim, degree, or dim+degree-i iterations,
00098       // whichever is largest.
00099       const unsigned int weight_loop_index =
00100         std::max(dim, std::max(degree, degree+dim-i));
00101 
00102       for (unsigned int j=1; j<=weight_loop_index; ++j)
00103         {
00104           if (j <= degree) // Accumulate (d+n-2i)^d term
00105             weight *= static_cast<Real>(degree+dim-2*i);
00106 
00107           if (j <= 2*s) // Accumulate 2^{-2s}
00108             weight *= 0.5;
00109 
00110           if (j <= i) // Accumulate (i!)^{-1}
00111             weight /= static_cast<Real>(j);
00112 
00113           if (j <= degree+dim-i) // Accumulate ( (d+n-i)! )^{-1}
00114             weight /= static_cast<Real>(j);
00115         }
00116 
00117       // This is the weight for each of the points computed previously
00118       for (unsigned int j=0; j<permutations.size(); ++j)
00119         _weights[offset+j] = weight;
00120 
00121       // Change sign for next iteration
00122       one_pm = -one_pm;
00123 
00124       // Update offset for the next set of points
00125       offset += permutations.size();
00126     }
00127 }

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

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

Definition at line 163 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, and libMesh::Threads::spin_mtx.

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

00164 {
00165   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
00166   std::pair<unsigned int, unsigned int>& p = _counts[name];
00167 
00168   p.first++;
00169 }

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

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

Definition at line 176 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, and libMesh::Threads::spin_mtx.

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

00177 {
00178   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
00179   std::pair<unsigned int, unsigned int>& p = _counts[name];
00180 
00181   p.second++;
00182 }

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

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

Definition at line 27 of file quadrature.C.

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

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

00029 {
00030   // check to see if we have already
00031   // done the work for this quadrature rule
00032   if (t == _type && p == _p_level)
00033     return;
00034   else
00035     {
00036       _type = t;
00037       _p_level = p;
00038     }
00039 
00040 
00041 
00042   switch(_dim)
00043     {
00044     case 0:
00045       this->init_0D(_type,_p_level);
00046 
00047       return;
00048 
00049     case 1:
00050       this->init_1D(_type,_p_level);
00051 
00052       return;
00053 
00054     case 2:
00055       this->init_2D(_type,_p_level);
00056 
00057       return;
00058 
00059     case 3:
00060       this->init_3D(_type,_p_level);
00061 
00062       return;
00063 
00064     default:
00065       libmesh_error();
00066     }
00067 }

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

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

Definition at line 71 of file quadrature.C.

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

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

00073 {
00074   _points.resize(1);
00075   _weights.resize(1);
00076   _points[0] = Point(0.);
00077   _weights[0] = 1.0;
00078 }

void libMesh::QGrundmann_Moller::init_1D ( const   type,
unsigned  p_level = 0 
) [inline, private, virtual]

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

Implements libMesh::QBase.

Definition at line 123 of file quadrature_gm.h.

00125   {
00126     // See about making this non-pure virtual in the base class
00127     libmesh_error();
00128   }

virtual void libMesh::QBase::init_2D ( const   ElemType,
unsigned int  = 0 
) [inline, protected, virtual, inherited]

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 in libMesh::QClough, libMesh::QConical, libMesh::QGauss, libMesh::QGrid, libMesh::QMonomial, libMesh::QSimpson, and libMesh::QTrap.

Definition at line 246 of file quadrature.h.

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

00249   {}

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

The GM rules are only defined for 3D since better 2D rules for simplexes are available.

Definition at line 28 of file quadrature_gm_3D.C.

References libMesh::QBase::allow_rules_with_negative_weights, libMesh::err, gm_rule(), libMeshEnums::TET10, and libMeshEnums::TET4.

00030 {
00031   // Nearly all GM rules contain negative weights, so if you are not
00032   // allowing rules with negative weights, we cannot continue!
00033   if (!allow_rules_with_negative_weights)
00034     {
00035       libMesh::err << "You requested a Grundmann-Moller rule but\n"
00036                     << "are not allowing rules with negative weights!\n"
00037                     << "Either select a different quadrature class or\n"
00038                     << "set allow_rules_with_negative_weights==true."
00039                     << std::endl;
00040 
00041       libmesh_error();
00042     }
00043 
00044   switch (type_in)
00045     {
00046     case TET4:
00047     case TET10:
00048       {
00049         // Untested above _order=23 but should work...
00050         gm_rule( (_order + 2*p)/2 );
00051         return;
00052 
00053       } // end case TET4, TET10
00054 
00055 
00056 
00057       //---------------------------------------------
00058       // Unsupported element type
00059     default:
00060       {
00061         libMesh::err << "ERROR: Unsupported element type: " << type_in << std::endl;
00062         libmesh_error();
00063       }
00064     } // end switch (type_in)
00065 
00066   // We must have returned or errored-out by this point.  If not,
00067   // throw an error now.
00068   libmesh_error();
00069   return;
00070 }

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

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

Definition at line 79 of file reference_counter.h.

References libMesh::ReferenceCounter::_n_objects.

00080   { return _n_objects; }

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

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

Definition at line 88 of file reference_counter.C.

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

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

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

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

Definition at line 362 of file quadrature.h.

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

Referenced by libMesh::operator<<().

00363 {
00364   libmesh_assert(!_points.empty());
00365   libmesh_assert(!_weights.empty());
00366 
00367   os << "N_Q_Points=" << this->n_points() << std::endl << std::endl;
00368   for (unsigned int qpoint=0; qpoint<this->n_points(); qpoint++)
00369     {
00370       os << " Point " << qpoint << ":\n"
00371          << "  "
00372          << _points[qpoint]
00373          << " Weight:\n "
00374          << "  w=" << _weights[qpoint] << "\n" << std::endl;
00375     }
00376 }

Point libMesh::QBase::qp ( const unsigned int  i  )  const [inline, inherited]
Returns:
the $ i^{th} $ quadrature point on the reference object.

Definition at line 150 of file quadrature.h.

References libMesh::QBase::_points.

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

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

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

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

Definition at line 82 of file quadrature.C.

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

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

00084 {
00085   // Make sure we are in 1D
00086   libmesh_assert_equal_to (_dim, 1);
00087 
00088   // Make sure that we have sane ranges
00089   libmesh_assert_greater (new_range.second, new_range.first);
00090   libmesh_assert_greater (old_range.second, old_range.first);
00091 
00092   // Make sure there are some points
00093   libmesh_assert_greater (_points.size(), 0);
00094 
00095   // We're mapping from old_range -> new_range
00096   for (unsigned int i=0; i<_points.size(); i++)
00097     {
00098       _points[i](0) =
00099         (_points[i](0) - old_range.first) *
00100         (new_range.second - new_range.first) /
00101         (old_range.second - old_range.first) +
00102         new_range.first;
00103     }
00104 
00105   // Compute the scale factor and scale the weights
00106   const Real scfact = (new_range.second - new_range.first) /
00107                       (old_range.second - old_range.first);
00108 
00109   for (unsigned int i=0; i<_points.size(); i++)
00110     _weights[i] *= scfact;
00111 }

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

Returns true if the shape functions need to be recalculated.

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

By default this will return false.

Definition at line 198 of file quadrature.h.

Referenced by libMesh::FE< Dim, T >::edge_reinit(), and libMesh::FE< Dim, T >::reinit().

00198 { return false; }

QuadratureType libMesh::QGrundmann_Moller::type (  )  const [inline, virtual]
Returns:
QGRUNDMANN_MOLLER

Implements libMesh::QBase.

Definition at line 118 of file quadrature_gm.h.

References libMeshEnums::QGRUNDMANN_MOLLER.

00118 { return QGRUNDMANN_MOLLER; }

Real libMesh::QBase::w ( const unsigned int  i  )  const [inline, inherited]
Returns:
the $ i^{th} $ quadrature weight.

Definition at line 156 of file quadrature.h.

References libMesh::QBase::_weights.

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

00157     { libmesh_assert_less (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.


Member Data Documentation

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

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

Definition at line 137 of file reference_counter.h.

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

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 131 of file reference_counter.h.

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

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

Definition at line 126 of file reference_counter.h.

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

libMesh::err<< "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 { libMesh::err << "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> libMesh::QBase::_points [protected, 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 215 of file quadrature.h.

Referenced by libMesh::QMonomial::init_3D(), init_3D(), and libMesh::QGauss::init_3D().


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

Site Created By: libMesh Developers
Last modified: February 05 2013 19:55:35 UTC

Hosted By:
SourceForge.net Logo