libMesh::QGauss Class Reference
#include <quadrature_gauss.h>
Inheritance diagram for libMesh::QGauss:
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=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< 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 (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) |
| 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 _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 43 of file quadrature_gauss.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::QGauss::QGauss | ( | const unsigned int | _dim, | |
| const Order | _order = INVALID_ORDER | |||
| ) | [inline] |
Constructor. Declares the order of the quadrature rule.
Definition at line 105 of file quadrature_gauss.h.
References libMeshEnums::EDGE2, and libMesh::QBase::init().
00106 : QBase(d,o) 00107 { 00108 // explicitly call the init function in 1D since the 00109 // other tensor-product rules require this one. 00110 // note that EDGE will not be used internally, however 00111 // if we called the function with INVALID_ELEM it would try to 00112 // be smart and return, thinking it had already done the work. 00113 if (_dim == 1) 00114 init(EDGE2); 00115 }
| libMesh::QGauss::~QGauss | ( | ) | [inline] |
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::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::QGauss::dunavant_rule | ( | const Real | rule_data[][4], | |
| const unsigned int | n_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 266 of file quadrature_gauss.C.
References libMesh::QBase::_points, libMesh::QBase::_weights, and libMesh::err.
Referenced by init_2D().
00268 { 00269 // The input data array has 4 columns. The first 3 are the permutation points. 00270 // The last column is the weights for a given set of permutation points. A zero 00271 // in two of the first 3 columns implies the point is a 1-permutation (centroid). 00272 // A zero in one of the first 3 columns implies the point is a 3-permutation. 00273 // Otherwise each point is assumed to be a 6-permutation. 00274 00275 // Always insert into the points & weights vector relative to the offset 00276 unsigned int offset=0; 00277 00278 00279 for (unsigned int p=0; p<n_pts; ++p) 00280 { 00281 00282 // There must always be a non-zero entry to start the row 00283 libmesh_assert_not_equal_to ( rule_data[p][0], static_cast<Real>(0.0) ); 00284 00285 // A zero weight may imply you did not set up the raw data correctly 00286 libmesh_assert_not_equal_to ( rule_data[p][3], static_cast<Real>(0.0) ); 00287 00288 // What kind of point is this? 00289 // One non-zero entry in first 3 cols ? 1-perm (centroid) point = 1 00290 // Two non-zero entries in first 3 cols ? 3-perm point = 3 00291 // Three non-zero entries ? 6-perm point = 6 00292 unsigned int pointtype=1; 00293 00294 if (rule_data[p][1] != static_cast<Real>(0.0)) 00295 { 00296 if (rule_data[p][2] != static_cast<Real>(0.0)) 00297 pointtype = 6; 00298 else 00299 pointtype = 3; 00300 } 00301 00302 switch (pointtype) 00303 { 00304 case 1: 00305 { 00306 // Be sure we have enough space to insert this point 00307 libmesh_assert_less (offset + 0, _points.size()); 00308 00309 // The point has only a single permutation (the centroid!) 00310 _points[offset + 0] = Point(rule_data[p][0], rule_data[p][0]); 00311 00312 // The weight is always the last entry in the row. 00313 _weights[offset + 0] = rule_data[p][3]; 00314 00315 offset += 1; 00316 break; 00317 } 00318 00319 case 3: 00320 { 00321 // Be sure we have enough space to insert these points 00322 libmesh_assert_less (offset + 2, _points.size()); 00323 00324 // Here it's understood the second entry is to be used twice, and 00325 // thus there are three possible permutations. 00326 _points[offset + 0] = Point(rule_data[p][0], rule_data[p][1]); 00327 _points[offset + 1] = Point(rule_data[p][1], rule_data[p][0]); 00328 _points[offset + 2] = Point(rule_data[p][1], rule_data[p][1]); 00329 00330 for (unsigned int j=0; j<3; ++j) 00331 _weights[offset + j] = rule_data[p][3]; 00332 00333 offset += 3; 00334 break; 00335 } 00336 00337 case 6: 00338 { 00339 // Be sure we have enough space to insert these points 00340 libmesh_assert_less (offset + 5, _points.size()); 00341 00342 // Three individual entries with six permutations. 00343 _points[offset + 0] = Point(rule_data[p][0], rule_data[p][1]); 00344 _points[offset + 1] = Point(rule_data[p][0], rule_data[p][2]); 00345 _points[offset + 2] = Point(rule_data[p][1], rule_data[p][0]); 00346 _points[offset + 3] = Point(rule_data[p][1], rule_data[p][2]); 00347 _points[offset + 4] = Point(rule_data[p][2], rule_data[p][0]); 00348 _points[offset + 5] = Point(rule_data[p][2], rule_data[p][1]); 00349 00350 for (unsigned int j=0; j<6; ++j) 00351 _weights[offset + j] = rule_data[p][3]; 00352 00353 offset += 6; 00354 break; 00355 } 00356 00357 default: 00358 { 00359 libMesh::err << "Don't know what to do with this many permutation points!" << std::endl; 00360 libmesh_error(); 00361 } 00362 } 00363 } 00364 }
| void libMesh::QGauss::dunavant_rule2 | ( | const Real * | wts, | |
| const Real * | a, | |||
| const Real * | b, | |||
| const unsigned int * | permutation_ids, | |||
| const unsigned int | n_wts | |||
| ) | [private] |
Definition at line 188 of file quadrature_gauss.C.
References libMesh::QBase::_points, libMesh::QBase::_weights, and libMesh::err.
Referenced by init_2D().
00193 { 00194 // Figure out how many total points by summing up the entries 00195 // in the permutation_ids array, and resize the _points and _weights 00196 // vectors appropriately. 00197 unsigned int total_pts = 0; 00198 for (unsigned int p=0; p<n_wts; ++p) 00199 total_pts += permutation_ids[p]; 00200 00201 // Resize point and weight vectors appropriately. 00202 _points.resize(total_pts); 00203 _weights.resize(total_pts); 00204 00205 // Always insert into the points & weights vector relative to the offset 00206 unsigned int offset=0; 00207 00208 for (unsigned int p=0; p<n_wts; ++p) 00209 { 00210 switch (permutation_ids[p]) 00211 { 00212 case 1: 00213 { 00214 // The point has only a single permutation (the centroid!) 00215 // So we don't even need to look in the a or b arrays. 00216 _points [offset + 0] = Point(1.0L/3.0L, 1.0L/3.0L); 00217 _weights[offset + 0] = wts[p]; 00218 00219 offset += 1; 00220 break; 00221 } 00222 00223 00224 case 3: 00225 { 00226 // For this type of rule, don't need to look in the b array. 00227 _points[offset + 0] = Point(a[p], a[p]); // (a,a) 00228 _points[offset + 1] = Point(a[p], 1.L-2.L*a[p]); // (a,1-2a) 00229 _points[offset + 2] = Point(1.L-2.L*a[p], a[p]); // (1-2a,a) 00230 00231 for (unsigned int j=0; j<3; ++j) 00232 _weights[offset + j] = wts[p]; 00233 00234 offset += 3; 00235 break; 00236 } 00237 00238 case 6: 00239 { 00240 // This type of point uses all 3 arrays... 00241 _points[offset + 0] = Point(a[p], b[p]); 00242 _points[offset + 1] = Point(b[p], a[p]); 00243 _points[offset + 2] = Point(a[p], 1.L-a[p]-b[p]); 00244 _points[offset + 3] = Point(1.L-a[p]-b[p], a[p]); 00245 _points[offset + 4] = Point(b[p], 1.L-a[p]-b[p]); 00246 _points[offset + 5] = Point(1.L-a[p]-b[p], b[p]); 00247 00248 for (unsigned int j=0; j<6; ++j) 00249 _weights[offset + j] = wts[p]; 00250 00251 offset += 6; 00252 break; 00253 } 00254 00255 default: 00256 { 00257 libMesh::err << "Unknown permutation id: " << permutation_ids[p] << "!" << std::endl; 00258 libmesh_error(); 00259 } 00260 } 00261 } 00262 00263 }
| 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] |
- Returns:
- the dimension of the quadrature rule.
Definition at line 123 of file quadrature.h.
Referenced by libMesh::InfFE< Dim, T_radial, T_map >::attach_quadrature_rule(), libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), and libMesh::QConical::conical_product_tri().
| 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.
| 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.
| std::vector<Point>& libMesh::QBase::get_points | ( | ) | [inline, inherited] |
- Returns:
- a
std::vectorcontaining 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; }
| const std::vector<Point>& libMesh::QBase::get_points | ( | ) | const [inline, inherited] |
- Returns:
- a
std::vectorcontaining the quadrature point locations on a reference object.
Definition at line 129 of file quadrature.h.
References libMesh::QBase::_points.
Referenced by libMesh::FE< Dim, T >::edge_reinit(), libMesh::QClough::init_1D(), libMesh::QMonomial::init_2D(), init_2D(), libMesh::QClough::init_2D(), libMesh::QMonomial::init_3D(), init_3D(), libMesh::InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), libMesh::InfFE< Dim, T_radial, T_map >::init_shape_functions(), libMesh::InfFE< Dim, T_radial, T_map >::reinit(), libMesh::FEXYZ< Dim >::reinit(), and libMesh::FE< Dim, T >::reinit().
00129 { return _points; }
| std::vector<Real>& libMesh::QBase::get_weights | ( | ) | [inline, inherited] |
- Returns:
- a
std::vectorcontaining the quadrature weights.
Definition at line 145 of file quadrature.h.
References libMesh::QBase::_weights.
00145 { return _weights; }
| const std::vector<Real>& libMesh::QBase::get_weights | ( | ) | const [inline, inherited] |
- Returns:
- a
std::vectorcontaining the quadrature weights.
Definition at line 140 of file quadrature.h.
References libMesh::QBase::_weights.
Referenced by libMesh::FE< Dim, T >::edge_reinit(), libMesh::QClough::init_1D(), libMesh::QMonomial::init_2D(), init_2D(), libMesh::QClough::init_2D(), libMesh::QMonomial::init_3D(), init_3D(), libMesh::InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), libMesh::InfFE< Dim, T_radial, T_map >::init_shape_functions(), libMesh::FEXYZ< Dim >::reinit(), and libMesh::FE< Dim, T >::reinit().
00140 { return _weights; }
| void libMesh::ReferenceCounter::increment_constructor_count | ( | const std::string & | name | ) | [inline, protected, inherited] |
Increments the construction counter. Should be called in the constructor of any derived class that will be reference counted.
Definition at line 163 of file reference_counter.h.
References libMesh::ReferenceCounter::_counts, 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(), init_2D(), libMesh::QClough::init_2D(), libMesh::QTrap::init_3D(), libMesh::QSimpson::init_3D(), libMesh::QMonomial::init_3D(), libMesh::QGrid::init_3D(), init_3D(), libMesh::InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), 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().
| void libMesh::QGauss::init_1D | ( | const ElemType | type = INVALID_ELEM, |
|
| unsigned int | p_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 libMesh::QBase.
Definition at line 30 of file quadrature_gauss_1D.C.
References libMesh::QBase::_points, libMesh::QBase::_weights, libMeshEnums::CONSTANT, libMeshEnums::EIGHTH, libMeshEnums::EIGHTTEENTH, libMeshEnums::ELEVENTH, libMesh::err, 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.
00032 { 00033 //---------------------------------------------------------------------- 00034 // 1D quadrature rules 00035 switch(_order + 2*p) 00036 { 00037 case CONSTANT: 00038 case FIRST: 00039 { 00040 _points.resize (1); 00041 _weights.resize(1); 00042 00043 _points[0](0) = 0.; 00044 00045 _weights[0] = 2.; 00046 00047 return; 00048 } 00049 case SECOND: 00050 case THIRD: 00051 { 00052 _points.resize (2); 00053 _weights.resize(2); 00054 00055 _points[0](0) = -5.7735026918962576450914878050196e-01L; // -sqrt(3)/3 00056 _points[1] = -_points[0]; 00057 00058 _weights[0] = 1.; 00059 _weights[1] = _weights[0]; 00060 00061 return; 00062 } 00063 case FOURTH: 00064 case FIFTH: 00065 { 00066 _points.resize (3); 00067 _weights.resize(3); 00068 00069 _points[ 0](0) = -7.7459666924148337703585307995648e-01L; 00070 _points[ 1](0) = 0.; 00071 _points[ 2] = -_points[0]; 00072 00073 _weights[ 0] = 5.5555555555555555555555555555556e-01L; 00074 _weights[ 1] = 8.8888888888888888888888888888889e-01L; 00075 _weights[ 2] = _weights[0]; 00076 00077 return; 00078 } 00079 case SIXTH: 00080 case SEVENTH: 00081 { 00082 _points.resize (4); 00083 _weights.resize(4); 00084 00085 _points[ 0](0) = -8.6113631159405257522394648889281e-01L; 00086 _points[ 1](0) = -3.3998104358485626480266575910324e-01L; 00087 _points[ 2] = -_points[1]; 00088 _points[ 3] = -_points[0]; 00089 00090 _weights[ 0] = 3.4785484513745385737306394922200e-01L; 00091 _weights[ 1] = 6.5214515486254614262693605077800e-01L; 00092 _weights[ 2] = _weights[1]; 00093 _weights[ 3] = _weights[0]; 00094 00095 return; 00096 } 00097 case EIGHTH: 00098 case NINTH: 00099 { 00100 _points.resize (5); 00101 _weights.resize(5); 00102 00103 _points[ 0](0) = -9.0617984593866399279762687829939e-01L; 00104 _points[ 1](0) = -5.3846931010568309103631442070021e-01L; 00105 _points[ 2](0) = 0.; 00106 _points[ 3] = -_points[1]; 00107 _points[ 4] = -_points[0]; 00108 00109 _weights[ 0] = 2.3692688505618908751426404071992e-01L; 00110 _weights[ 1] = 4.7862867049936646804129151483564e-01L; 00111 _weights[ 2] = 5.6888888888888888888888888888889e-01L; 00112 _weights[ 3] = _weights[1]; 00113 _weights[ 4] = _weights[0]; 00114 00115 return; 00116 } 00117 case TENTH: 00118 case ELEVENTH: 00119 { 00120 _points.resize (6); 00121 _weights.resize(6); 00122 00123 _points[ 0](0) = -9.3246951420315202781230155449399e-01L; 00124 _points[ 1](0) = -6.6120938646626451366139959501991e-01L; 00125 _points[ 2](0) = -2.3861918608319690863050172168071e-01L; 00126 _points[ 3] = -_points[2]; 00127 _points[ 4] = -_points[1]; 00128 _points[ 5] = -_points[0]; 00129 00130 _weights[ 0] = 1.7132449237917034504029614217273e-01L; 00131 _weights[ 1] = 3.6076157304813860756983351383772e-01L; 00132 _weights[ 2] = 4.6791393457269104738987034398955e-01L; 00133 _weights[ 3] = _weights[2]; 00134 _weights[ 4] = _weights[1]; 00135 _weights[ 5] = _weights[0]; 00136 00137 return; 00138 } 00139 case TWELFTH: 00140 case THIRTEENTH: 00141 { 00142 _points.resize (7); 00143 _weights.resize(7); 00144 00145 _points[ 0](0) = -9.4910791234275852452618968404785e-01L; 00146 _points[ 1](0) = -7.4153118559939443986386477328079e-01L; 00147 _points[ 2](0) = -4.0584515137739716690660641207696e-01L; 00148 _points[ 3](0) = 0.; 00149 _points[ 4] = -_points[2]; 00150 _points[ 5] = -_points[1]; 00151 _points[ 6] = -_points[0]; 00152 00153 _weights[ 0] = 1.2948496616886969327061143267908e-01L; 00154 _weights[ 1] = 2.7970539148927666790146777142378e-01L; 00155 _weights[ 2] = 3.8183005050511894495036977548898e-01L; 00156 _weights[ 3] = 4.1795918367346938775510204081633e-01L; 00157 _weights[ 4] = _weights[2]; 00158 _weights[ 5] = _weights[1]; 00159 _weights[ 6] = _weights[0]; 00160 00161 return; 00162 } 00163 case FOURTEENTH: 00164 case FIFTEENTH: 00165 { 00166 _points.resize (8); 00167 _weights.resize(8); 00168 00169 _points[ 0](0) = -9.6028985649753623168356086856947e-01L; 00170 _points[ 1](0) = -7.9666647741362673959155393647583e-01L; 00171 _points[ 2](0) = -5.2553240991632898581773904918925e-01L; 00172 _points[ 3](0) = -1.8343464249564980493947614236018e-01L; 00173 _points[ 4] = -_points[3]; 00174 _points[ 5] = -_points[2]; 00175 _points[ 6] = -_points[1]; 00176 _points[ 7] = -_points[0]; 00177 00178 _weights[ 0] = 1.0122853629037625915253135430996e-01L; 00179 _weights[ 1] = 2.2238103445337447054435599442624e-01L; 00180 _weights[ 2] = 3.1370664587788728733796220198660e-01L; 00181 _weights[ 3] = 3.6268378337836198296515044927720e-01L; 00182 _weights[ 4] = _weights[3]; 00183 _weights[ 5] = _weights[2]; 00184 _weights[ 6] = _weights[1]; 00185 _weights[ 7] = _weights[0]; 00186 00187 return; 00188 } 00189 case SIXTEENTH: 00190 case SEVENTEENTH: 00191 { 00192 _points.resize (9); 00193 _weights.resize(9); 00194 00195 _points[ 0](0) = -9.6816023950762608983557620290367e-01L; 00196 _points[ 1](0) = -8.3603110732663579429942978806973e-01L; 00197 _points[ 2](0) = -6.1337143270059039730870203934147e-01L; 00198 _points[ 3](0) = -3.2425342340380892903853801464334e-01L; 00199 _points[ 4](0) = 0.; 00200 _points[ 5] = -_points[3]; 00201 _points[ 6] = -_points[2]; 00202 _points[ 7] = -_points[1]; 00203 _points[ 8] = -_points[0]; 00204 00205 _weights[ 0] = 8.1274388361574411971892158110524e-02L; 00206 _weights[ 1] = 1.8064816069485740405847203124291e-01L; 00207 _weights[ 2] = 2.6061069640293546231874286941863e-01L; 00208 _weights[ 3] = 3.1234707704000284006863040658444e-01L; 00209 _weights[ 4] = 3.3023935500125976316452506928697e-01L; 00210 _weights[ 5] = _weights[3]; 00211 _weights[ 6] = _weights[2]; 00212 _weights[ 7] = _weights[1]; 00213 _weights[ 8] = _weights[0]; 00214 00215 return; 00216 } 00217 case EIGHTTEENTH: 00218 case NINTEENTH: 00219 { 00220 _points.resize (10); 00221 _weights.resize(10); 00222 00223 _points[ 0](0) = -9.7390652851717172007796401208445e-01L; 00224 _points[ 1](0) = -8.6506336668898451073209668842349e-01L; 00225 _points[ 2](0) = -6.7940956829902440623432736511487e-01L; 00226 _points[ 3](0) = -4.3339539412924719079926594316578e-01L; 00227 _points[ 4](0) = -1.4887433898163121088482600112972e-01L; 00228 _points[ 5] = -_points[4]; 00229 _points[ 6] = -_points[3]; 00230 _points[ 7] = -_points[2]; 00231 _points[ 8] = -_points[1]; 00232 _points[ 9] = -_points[0]; 00233 00234 _weights[ 0] = 6.6671344308688137593568809893332e-02L; 00235 _weights[ 1] = 1.4945134915058059314577633965770e-01L; 00236 _weights[ 2] = 2.1908636251598204399553493422816e-01L; 00237 _weights[ 3] = 2.6926671930999635509122692156947e-01L; 00238 _weights[ 4] = 2.9552422471475287017389299465134e-01L; 00239 _weights[ 5] = _weights[4]; 00240 _weights[ 6] = _weights[3]; 00241 _weights[ 7] = _weights[2]; 00242 _weights[ 8] = _weights[1]; 00243 _weights[ 9] = _weights[0]; 00244 00245 return; 00246 } 00247 00248 case TWENTIETH: 00249 case TWENTYFIRST: 00250 { 00251 _points.resize (11); 00252 _weights.resize(11); 00253 00254 _points[ 0](0) = -9.7822865814605699280393800112286e-01L; 00255 _points[ 1](0) = -8.8706259976809529907515776930393e-01L; 00256 _points[ 2](0) = -7.3015200557404932409341625203115e-01L; 00257 _points[ 3](0) = -5.1909612920681181592572566945861e-01L; 00258 _points[ 4](0) = -2.6954315595234497233153198540086e-01L; 00259 _points[ 5](0) = 0.; 00260 _points[ 6] = -_points[4]; 00261 _points[ 7] = -_points[3]; 00262 _points[ 8] = -_points[2]; 00263 _points[ 9] = -_points[1]; 00264 _points[10] = -_points[0]; 00265 00266 _weights[ 0] = 5.5668567116173666482753720442549e-02L; 00267 _weights[ 1] = 1.2558036946490462463469429922394e-01L; 00268 _weights[ 2] = 1.8629021092773425142609764143166e-01L; 00269 _weights[ 3] = 2.3319376459199047991852370484318e-01L; 00270 _weights[ 4] = 2.6280454451024666218068886989051e-01L; 00271 _weights[ 5] = 2.7292508677790063071448352833634e-01L; 00272 _weights[ 6] = _weights[4]; 00273 _weights[ 7] = _weights[3]; 00274 _weights[ 8] = _weights[2]; 00275 _weights[ 9] = _weights[1]; 00276 _weights[10] = _weights[0]; 00277 00278 return; 00279 } 00280 00281 case TWENTYSECOND: 00282 case TWENTYTHIRD: 00283 { 00284 _points.resize (12); 00285 _weights.resize(12); 00286 00287 _points[ 0](0) = -9.8156063424671925069054909014928e-01L; 00288 _points[ 1](0) = -9.0411725637047485667846586611910e-01L; 00289 _points[ 2](0) = -7.6990267419430468703689383321282e-01L; 00290 _points[ 3](0) = -5.8731795428661744729670241894053e-01L; 00291 _points[ 4](0) = -3.6783149899818019375269153664372e-01L; 00292 _points[ 5](0) = -1.2523340851146891547244136946385e-01L; 00293 _points[ 6] = -_points[5]; 00294 _points[ 7] = -_points[4]; 00295 _points[ 8] = -_points[3]; 00296 _points[ 9] = -_points[2]; 00297 _points[10] = -_points[1]; 00298 _points[11] = -_points[0]; 00299 00300 _weights[ 0] = 4.7175336386511827194615961485017e-02L; 00301 _weights[ 1] = 1.0693932599531843096025471819400e-01L; 00302 _weights[ 2] = 1.6007832854334622633465252954336e-01L; 00303 _weights[ 3] = 2.0316742672306592174906445580980e-01L; 00304 _weights[ 4] = 2.3349253653835480876084989892483e-01L; 00305 _weights[ 5] = 2.4914704581340278500056243604295e-01L; 00306 _weights[ 6] = _weights[5]; 00307 _weights[ 7] = _weights[4]; 00308 _weights[ 8] = _weights[3]; 00309 _weights[ 9] = _weights[2]; 00310 _weights[10] = _weights[1]; 00311 _weights[11] = _weights[0]; 00312 00313 return; 00314 } 00315 00316 case TWENTYFOURTH: 00317 case TWENTYFIFTH: 00318 { 00319 _points.resize (13); 00320 _weights.resize(13); 00321 00322 _points[ 0](0) = -9.8418305471858814947282944880711e-01L; 00323 _points[ 1](0) = -9.1759839922297796520654783650072e-01L; 00324 _points[ 2](0) = -8.0157809073330991279420648958286e-01L; 00325 _points[ 3](0) = -6.4234933944034022064398460699552e-01L; 00326 _points[ 4](0) = -4.4849275103644685287791285212764e-01L; 00327 _points[ 5](0) = -2.3045831595513479406552812109799e-01L; 00328 _points[ 6](0) = 0.; 00329 _points[ 7] = -_points[5]; 00330 _points[ 8] = -_points[4]; 00331 _points[ 9] = -_points[3]; 00332 _points[10] = -_points[2]; 00333 _points[11] = -_points[1]; 00334 _points[12] = -_points[0]; 00335 00336 _weights[ 0] = 4.0484004765315879520021592200986e-02L; 00337 _weights[ 1] = 9.2121499837728447914421775953797e-02L; 00338 _weights[ 2] = 1.3887351021978723846360177686887e-01L; 00339 _weights[ 3] = 1.7814598076194573828004669199610e-01L; 00340 _weights[ 4] = 2.0781604753688850231252321930605e-01L; 00341 _weights[ 5] = 2.2628318026289723841209018603978e-01L; 00342 _weights[ 6] = 2.3255155323087391019458951526884e-01L; 00343 _weights[ 7] = _weights[5]; 00344 _weights[ 8] = _weights[4]; 00345 _weights[ 9] = _weights[3]; 00346 _weights[10] = _weights[2]; 00347 _weights[11] = _weights[1]; 00348 _weights[12] = _weights[0]; 00349 00350 return; 00351 } 00352 00353 case TWENTYSIXTH: 00354 case TWENTYSEVENTH: 00355 { 00356 _points.resize (14); 00357 _weights.resize(14); 00358 00359 _points[ 0](0) = -9.8628380869681233884159726670405e-01L; 00360 _points[ 1](0) = -9.2843488366357351733639113937787e-01L; 00361 _points[ 2](0) = -8.2720131506976499318979474265039e-01L; 00362 _points[ 3](0) = -6.8729290481168547014801980301933e-01L; 00363 _points[ 4](0) = -5.1524863635815409196529071855119e-01L; 00364 _points[ 5](0) = -3.1911236892788976043567182416848e-01L; 00365 _points[ 6](0) = -1.0805494870734366206624465021983e-01L; 00366 _points[ 7] = -_points[6]; 00367 _points[ 8] = -_points[5]; 00368 _points[ 9] = -_points[4]; 00369 _points[10] = -_points[3]; 00370 _points[11] = -_points[2]; 00371 _points[12] = -_points[1]; 00372 _points[13] = -_points[0]; 00373 00374 _weights[ 0] = 3.5119460331751863031832876138192e-02L; 00375 _weights[ 1] = 8.0158087159760209805633277062854e-02L; 00376 _weights[ 2] = 1.2151857068790318468941480907248e-01L; 00377 _weights[ 3] = 1.5720316715819353456960193862384e-01L; 00378 _weights[ 4] = 1.8553839747793781374171659012516e-01L; 00379 _weights[ 5] = 2.0519846372129560396592406566122e-01L; 00380 _weights[ 6] = 2.1526385346315779019587644331626e-01L; 00381 _weights[ 7] = _weights[6]; 00382 _weights[ 8] = _weights[5]; 00383 _weights[ 9] = _weights[4]; 00384 _weights[10] = _weights[3]; 00385 _weights[11] = _weights[2]; 00386 _weights[12] = _weights[1]; 00387 _weights[13] = _weights[0]; 00388 00389 return; 00390 } 00391 00392 case TWENTYEIGHTH: 00393 case TWENTYNINTH: 00394 { 00395 _points.resize (15); 00396 _weights.resize(15); 00397 00398 _points[ 0](0) = -9.8799251802048542848956571858661e-01L; 00399 _points[ 1](0) = -9.3727339240070590430775894771021e-01L; 00400 _points[ 2](0) = -8.4820658341042721620064832077422e-01L; 00401 _points[ 3](0) = -7.2441773136017004741618605461394e-01L; 00402 _points[ 4](0) = -5.7097217260853884753722673725391e-01L; 00403 _points[ 5](0) = -3.9415134707756336989720737098105e-01L; 00404 _points[ 6](0) = -2.0119409399743452230062830339460e-01L; 00405 _points[ 7](0) = 0.; 00406 _points[ 8] = -_points[6]; 00407 _points[ 9] = -_points[5]; 00408 _points[10] = -_points[4]; 00409 _points[11] = -_points[3]; 00410 _points[12] = -_points[2]; 00411 _points[13] = -_points[1]; 00412 _points[14] = -_points[0]; 00413 00414 _weights[ 0] = 3.0753241996117268354628393577204e-02L; 00415 _weights[ 1] = 7.0366047488108124709267416450667e-02L; 00416 _weights[ 2] = 1.0715922046717193501186954668587e-01L; 00417 _weights[ 3] = 1.3957067792615431444780479451103e-01L; 00418 _weights[ 4] = 1.6626920581699393355320086048121e-01L; 00419 _weights[ 5] = 1.8616100001556221102680056186642e-01L; 00420 _weights[ 6] = 1.9843148532711157645611832644384e-01L; 00421 _weights[ 7] = 2.0257824192556127288062019996752e-01L; 00422 _weights[ 8] = _weights[6]; 00423 _weights[ 9] = _weights[5]; 00424 _weights[10] = _weights[4]; 00425 _weights[11] = _weights[3]; 00426 _weights[12] = _weights[2]; 00427 _weights[13] = _weights[1]; 00428 _weights[14] = _weights[0]; 00429 00430 return; 00431 } 00432 00433 case THIRTIETH: 00434 case THIRTYFIRST: 00435 { 00436 _points.resize (16); 00437 _weights.resize(16); 00438 00439 _points[ 0](0) = -9.8940093499164993259615417345033e-01L; 00440 _points[ 1](0) = -9.4457502307323257607798841553461e-01L; 00441 _points[ 2](0) = -8.6563120238783174388046789771239e-01L; 00442 _points[ 3](0) = -7.5540440835500303389510119484744e-01L; 00443 _points[ 4](0) = -6.1787624440264374844667176404879e-01L; 00444 _points[ 5](0) = -4.5801677765722738634241944298358e-01L; 00445 _points[ 6](0) = -2.8160355077925891323046050146050e-01L; 00446 _points[ 7](0) = -9.5012509837637440185319335424958e-02L; 00447 _points[ 8] = -_points[7]; 00448 _points[ 9] = -_points[6]; 00449 _points[10] = -_points[5]; 00450 _points[11] = -_points[4]; 00451 _points[12] = -_points[3]; 00452 _points[13] = -_points[2]; 00453 _points[14] = -_points[1]; 00454 _points[15] = -_points[0]; 00455 00456 _weights[ 0] = 2.7152459411754094851780572456018e-02L; 00457 _weights[ 1] = 6.2253523938647892862843836994378e-02L; 00458 _weights[ 2] = 9.5158511682492784809925107602246e-02L; 00459 _weights[ 3] = 1.2462897125553387205247628219202e-01L; 00460 _weights[ 4] = 1.4959598881657673208150173054748e-01L; 00461 _weights[ 5] = 1.6915651939500253818931207903033e-01L; 00462 _weights[ 6] = 1.8260341504492358886676366796922e-01L; 00463 _weights[ 7] = 1.8945061045506849628539672320828e-01L; 00464 _weights[ 8] = _weights[7]; 00465 _weights[ 9] = _weights[6]; 00466 _weights[10] = _weights[5]; 00467 _weights[11] = _weights[4]; 00468 _weights[12] = _weights[3]; 00469 _weights[13] = _weights[2]; 00470 _weights[14] = _weights[1]; 00471 _weights[15] = _weights[0]; 00472 00473 return; 00474 } 00475 00476 case THIRTYSECOND: 00477 case THIRTYTHIRD: 00478 { 00479 _points.resize (17); 00480 _weights.resize(17); 00481 00482 _points[ 0](0) = -9.9057547531441733567543401994067e-01L; 00483 _points[ 1](0) = -9.5067552176876776122271695789580e-01L; 00484 _points[ 2](0) = -8.8023915372698590212295569448816e-01L; 00485 _points[ 3](0) = -7.8151400389680140692523005552048e-01L; 00486 _points[ 4](0) = -6.5767115921669076585030221664300e-01L; 00487 _points[ 5](0) = -5.1269053708647696788624656862955e-01L; 00488 _points[ 6](0) = -3.5123176345387631529718551709535e-01L; 00489 _points[ 7](0) = -1.7848418149584785585067749365407e-01L; 00490 _points[ 8](0) = 0.; 00491 _points[ 9] = -_points[7]; 00492 _points[10] = -_points[6]; 00493 _points[11] = -_points[5]; 00494 _points[12] = -_points[4]; 00495 _points[13] = -_points[3]; 00496 _points[14] = -_points[2]; 00497 _points[15] = -_points[1]; 00498 _points[16] = -_points[0]; 00499 00500 _weights[ 0] = 2.4148302868547931960110026287565e-02L; 00501 _weights[ 1] = 5.5459529373987201129440165358245e-02L; 00502 _weights[ 2] = 8.5036148317179180883535370191062e-02L; 00503 _weights[ 3] = 1.1188384719340397109478838562636e-01L; 00504 _weights[ 4] = 1.3513636846852547328631998170235e-01L; 00505 _weights[ 5] = 1.5404576107681028808143159480196e-01L; 00506 _weights[ 6] = 1.6800410215645004450997066378832e-01L; 00507 _weights[ 7] = 1.7656270536699264632527099011320e-01L; 00508 _weights[ 8] = 1.7944647035620652545826564426189e-01L; 00509 _weights[ 9] = _weights[7]; 00510 _weights[10] = _weights[6]; 00511 _weights[11] = _weights[5]; 00512 _weights[12] = _weights[4]; 00513 _weights[13] = _weights[3]; 00514 _weights[14] = _weights[2]; 00515 _weights[15] = _weights[1]; 00516 _weights[16] = _weights[0]; 00517 00518 return; 00519 } 00520 00521 case THIRTYFOURTH: 00522 case THIRTYFIFTH: 00523 { 00524 _points.resize (18); 00525 _weights.resize(18); 00526 00527 _points[ 0](0) = -9.9156516842093094673001600470615e-01L; 00528 _points[ 1](0) = -9.5582394957139775518119589292978e-01L; 00529 _points[ 2](0) = -8.9260246649755573920606059112715e-01L; 00530 _points[ 3](0) = -8.0370495897252311568241745501459e-01L; 00531 _points[ 4](0) = -6.9168704306035320787489108128885e-01L; 00532 _points[ 5](0) = -5.5977083107394753460787154852533e-01L; 00533 _points[ 6](0) = -4.1175116146284264603593179383305e-01L; 00534 _points[ 7](0) = -2.5188622569150550958897285487791e-01L; 00535 _points[ 8](0) = -8.4775013041735301242261852935784e-02L; 00536 _points[ 9] = -_points[8]; 00537 _points[10] = -_points[7]; 00538 _points[11] = -_points[6]; 00539 _points[12] = -_points[5]; 00540 _points[13] = -_points[4]; 00541 _points[14] = -_points[3]; 00542 _points[15] = -_points[2]; 00543 _points[16] = -_points[1]; 00544 _points[17] = -_points[0]; 00545 00546 _weights[ 0] = 2.1616013526483310313342710266452e-02L; 00547 _weights[ 1] = 4.9714548894969796453334946202639e-02L; 00548 _weights[ 2] = 7.6425730254889056529129677616637e-02L; 00549 _weights[ 3] = 1.0094204410628716556281398492483e-01L; 00550 _weights[ 4] = 1.2255520671147846018451912680020e-01L; 00551 _weights[ 5] = 1.4064291467065065120473130375195e-01L; 00552 _weights[ 6] = 1.5468467512626524492541800383637e-01L; 00553 _weights[ 7] = 1.6427648374583272298605377646593e-01L; 00554 _weights[ 8] = 1.6914238296314359184065647013499e-01L; 00555 _weights[ 9] = _weights[8]; 00556 _weights[10] = _weights[7]; 00557 _weights[11] = _weights[6]; 00558 _weights[12] = _weights[5]; 00559 _weights[13] = _weights[4]; 00560 _weights[14] = _weights[3]; 00561 _weights[15] = _weights[2]; 00562 _weights[16] = _weights[1]; 00563 _weights[17] = _weights[0]; 00564 00565 return; 00566 } 00567 00568 case THIRTYSIXTH: 00569 case THIRTYSEVENTH: 00570 { 00571 _points.resize (19); 00572 _weights.resize(19); 00573 00574 _points[ 0](0) = -9.9240684384358440318901767025326e-01L; 00575 _points[ 1](0) = -9.6020815213483003085277884068765e-01L; 00576 _points[ 2](0) = -9.0315590361481790164266092853231e-01L; 00577 _points[ 3](0) = -8.2271465653714282497892248671271e-01L; 00578 _points[ 4](0) = -7.2096617733522937861709586082378e-01L; 00579 _points[ 5](0) = -6.0054530466168102346963816494624e-01L; 00580 _points[ 6](0) = -4.6457074137596094571726714810410e-01L; 00581 _points[ 7](0) = -3.1656409996362983199011732884984e-01L; 00582 _points[ 8](0) = -1.6035864564022537586809611574074e-01L; 00583 _points[ 9](0) = 0.; 00584 _points[10] = -_points[8]; 00585 _points[11] = -_points[7]; 00586 _points[12] = -_points[6]; 00587 _points[13] = -_points[5]; 00588 _points[14] = -_points[4]; 00589 _points[15] = -_points[3]; 00590 _points[16] = -_points[2]; 00591 _points[17] = -_points[1]; 00592 _points[18] = -_points[0]; 00593 00594 _weights[ 0] = 1.9461788229726477036312041464438e-02L; 00595 _weights[ 1] = 4.4814226765699600332838157401994e-02L; 00596 _weights[ 2] = 6.9044542737641226580708258006013e-02L; 00597 _weights[ 3] = 9.1490021622449999464462094123840e-02L; 00598 _weights[ 4] = 1.1156664554733399471602390168177e-01L; 00599 _weights[ 5] = 1.2875396253933622767551578485688e-01L; 00600 _weights[ 6] = 1.4260670217360661177574610944190e-01L; 00601 _weights[ 7] = 1.5276604206585966677885540089766e-01L; 00602 _weights[ 8] = 1.5896884339395434764995643946505e-01L; 00603 _weights[ 9] = 1.6105444984878369597916362532092e-01L; 00604 _weights[10] = _weights[8]; 00605 _weights[11] = _weights[7]; 00606 _weights[12] = _weights[6]; 00607 _weights[13] = _weights[5]; 00608 _weights[14] = _weights[4]; 00609 _weights[15] = _weights[3]; 00610 _weights[16] = _weights[2]; 00611 _weights[17] = _weights[1]; 00612 _weights[18] = _weights[0]; 00613 00614 return; 00615 } 00616 00617 case THIRTYEIGHTH: 00618 case THIRTYNINTH: 00619 { 00620 _points.resize (20); 00621 _weights.resize(20); 00622 00623 _points[ 0](0) = -9.9312859918509492478612238847132e-01L; 00624 _points[ 1](0) = -9.6397192727791379126766613119728e-01L; 00625 _points[ 2](0) = -9.1223442825132590586775244120330e-01L; 00626 _points[ 3](0) = -8.3911697182221882339452906170152e-01L; 00627 _points[ 4](0) = -7.4633190646015079261430507035564e-01L; 00628 _points[ 5](0) = -6.3605368072651502545283669622629e-01L; 00629 _points[ 6](0) = -5.1086700195082709800436405095525e-01L; 00630 _points[ 7](0) = -3.7370608871541956067254817702493e-01L; 00631 _points[ 8](0) = -2.2778585114164507808049619536857e-01L; 00632 _points[ 9](0) = -7.6526521133497333754640409398838e-02L; 00633 _points[10] = -_points[9]; 00634 _points[11] = -_points[8]; 00635 _points[12] = -_points[7]; 00636 _points[13] = -_points[6]; 00637 _points[14] = -_points[5]; 00638 _points[15] = -_points[4]; 00639 _points[16] = -_points[3]; 00640 _points[17] = -_points[2]; 00641 _points[18] = -_points[1]; 00642 _points[19] = -_points[0]; 00643 00644 _weights[ 0] = 1.7614007139152118311861962351853e-02L; 00645 _weights[ 1] = 4.0601429800386941331039952274932e-02L; 00646 _weights[ 2] = 6.2672048334109063569506535187042e-02L; 00647 _weights[ 3] = 8.3276741576704748724758143222046e-02L; 00648 _weights[ 4] = 1.0193011981724043503675013548035e-01L; 00649 _weights[ 5] = 1.1819453196151841731237737771138e-01L; 00650 _weights[ 6] = 1.3168863844917662689849449974816e-01L; 00651 _weights[ 7] = 1.4209610931838205132929832506716e-01L; 00652 _weights[ 8] = 1.4917298647260374678782873700197e-01L; 00653 _weights[ 9] = 1.5275338713072585069808433195510e-01L; 00654 _weights[10] = _weights[9]; 00655 _weights[11] = _weights[8]; 00656 _weights[12] = _weights[7]; 00657 _weights[13] = _weights[6]; 00658 _weights[14] = _weights[5]; 00659 _weights[15] = _weights[4]; 00660 _weights[16] = _weights[3]; 00661 _weights[17] = _weights[2]; 00662 _weights[18] = _weights[1]; 00663 _weights[19] = _weights[0]; 00664 00665 return; 00666 } 00667 00668 case FORTIETH: 00669 case FORTYFIRST: 00670 { 00671 _points.resize (21); 00672 _weights.resize(21); 00673 00674 _points[ 0](0) = -9.9375217062038950026024203593794e-01L; 00675 _points[ 1](0) = -9.6722683856630629431662221490770e-01L; 00676 _points[ 2](0) = -9.2009933415040082879018713371497e-01L; 00677 _points[ 3](0) = -8.5336336458331728364725063858757e-01L; 00678 _points[ 4](0) = -7.6843996347567790861587785130623e-01L; 00679 _points[ 5](0) = -6.6713880419741231930596666999034e-01L; 00680 _points[ 6](0) = -5.5161883588721980705901879672431e-01L; 00681 _points[ 7](0) = -4.2434212020743878357366888854379e-01L; 00682 _points[ 8](0) = -2.8802131680240109660079251606460e-01L; 00683 _points[ 9](0) = -1.4556185416089509093703098233869e-01L; 00684 _points[10](0) = 0.; 00685 _points[11] = -_points[9]; 00686 _points[12] = -_points[8]; 00687 _points[13] = -_points[7]; 00688 _points[14] = -_points[6]; 00689 _points[15] = -_points[5]; 00690 _points[16] = -_points[4]; 00691 _points[17] = -_points[3]; 00692 _points[18] = -_points[2]; 00693 _points[19] = -_points[1]; 00694 _points[20] = -_points[0]; 00695 00696 _weights[ 0] = 1.6017228257774333324224616858471e-02L; 00697 _weights[ 1] = 3.6953789770852493799950668299330e-02L; 00698 _weights[ 2] = 5.7134425426857208283635826472448e-02L; 00699 _weights[ 3] = 7.6100113628379302017051653300183e-02L; 00700 _weights[ 4] = 9.3444423456033861553289741113932e-02L; 00701 _weights[ 5] = 1.0879729916714837766347457807011e-01L; 00702 _weights[ 6] = 1.2183141605372853419536717712572e-01L; 00703 _weights[ 7] = 1.3226893863333746178105257449678e-01L; 00704 _weights[ 8] = 1.3988739479107315472213342386758e-01L; 00705 _weights[ 9] = 1.4452440398997005906382716655375e-01L; 00706 _weights[10] = 1.4608113364969042719198514768337e-01L; 00707 _weights[11] = _weights[9]; 00708 _weights[12] = _weights[8]; 00709 _weights[13] = _weights[7]; 00710 _weights[14] = _weights[6]; 00711 _weights[15] = _weights[5]; 00712 _weights[16] = _weights[4]; 00713 _weights[17] = _weights[3]; 00714 _weights[18] = _weights[2]; 00715 _weights[19] = _weights[1]; 00716 _weights[20] = _weights[0]; 00717 00718 return; 00719 } 00720 00721 case FORTYSECOND: 00722 case FORTYTHIRD: 00723 { 00724 _points.resize (22); 00725 _weights.resize(22); 00726 00727 _points[ 0](0) = -9.9429458548239929207303142116130e-01L; 00728 _points[ 1](0) = -9.7006049783542872712395098676527e-01L; 00729 _points[ 2](0) = -9.2695677218717400052069293925905e-01L; 00730 _points[ 3](0) = -8.6581257772030013653642563701938e-01L; 00731 _points[ 4](0) = -7.8781680597920816200427795540835e-01L; 00732 _points[ 5](0) = -6.9448726318668278005068983576226e-01L; 00733 _points[ 6](0) = -5.8764040350691159295887692763865e-01L; 00734 _points[ 7](0) = -4.6935583798675702640633071096641e-01L; 00735 _points[ 8](0) = -3.4193582089208422515814742042738e-01L; 00736 _points[ 9](0) = -2.0786042668822128547884653391955e-01L; 00737 _points[10](0) = -6.9739273319722221213841796118628e-02L; 00738 _points[11] = -_points[10]; 00739 _points[12] = -_points[9]; 00740 _points[13] = -_points[8]; 00741 _points[14] = -_points[7]; 00742 _points[15] = -_points[6]; 00743 _points[16] = -_points[5]; 00744 _points[17] = -_points[4]; 00745 _points[18] = -_points[3]; 00746 _points[19] = -_points[2]; 00747 _points[20] = -_points[1]; 00748 _points[21] = -_points[0]; 00749 00750 _weights[ 0] = 1.4627995298272200684991098047185e-02L; 00751 _weights[ 1] = 3.3774901584814154793302246865913e-02L; 00752 _weights[ 2] = 5.2293335152683285940312051273211e-02L; 00753 _weights[ 3] = 6.9796468424520488094961418930218e-02L; 00754 _weights[ 4] = 8.5941606217067727414443681372703e-02L; 00755 _weights[ 5] = 1.0041414444288096493207883783054e-01L; 00756 _weights[ 6] = 1.1293229608053921839340060742175e-01L; 00757 _weights[ 7] = 1.2325237681051242428556098615481e-01L; 00758 _weights[ 8] = 1.3117350478706237073296499253031e-01L; 00759 _weights[ 9] = 1.3654149834601517135257383123152e-01L; 00760 _weights[10] = 1.3925187285563199337541024834181e-01L; 00761 _weights[11] = _weights[10]; 00762 _weights[12] = _weights[9]; 00763 _weights[13] = _weights[8]; 00764 _weights[14] = _weights[7]; 00765 _weights[15] = _weights[6]; 00766 _weights[16] = _weights[5]; 00767 _weights[17] = _weights[4]; 00768 _weights[18] = _weights[3]; 00769 _weights[19] = _weights[2]; 00770 _weights[20] = _weights[1]; 00771 _weights[21] = _weights[0]; 00772 00773 return; 00774 } 00775 00776 00777 default: 00778 { 00779 libMesh::err << "Quadrature rule " << _order 00780 << " not supported!" << std::endl; 00781 00782 libmesh_error(); 00783 } 00784 } 00785 00786 00787 00788 return; 00789 }
| void libMesh::QGauss::init_2D | ( | const ElemType | 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 libMesh::QBase.
Definition at line 28 of file quadrature_gauss_2D.C.
References libMesh::QBase::_points, libMesh::QBase::_weights, libMeshEnums::CONSTANT, dunavant_rule(), dunavant_rule2(), libMeshEnums::EDGE2, libMeshEnums::EIGHTH, libMeshEnums::EIGHTTEENTH, libMeshEnums::ELEVENTH, libMesh::err, libMeshEnums::FIFTEENTH, libMeshEnums::FIFTH, libMeshEnums::FIRST, libMeshEnums::FOURTEENTH, libMeshEnums::FOURTH, libMesh::QBase::get_points(), libMesh::QBase::get_weights(), libMesh::QBase::init(), libMeshEnums::NINTEENTH, libMeshEnums::NINTH, libMeshEnums::QUAD4, libMeshEnums::QUAD8, libMeshEnums::QUAD9, libMesh::Real, 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.
00030 { 00031 #if LIBMESH_DIM > 1 00032 00033 //----------------------------------------------------------------------- 00034 // 2D quadrature rules 00035 switch (type_in) 00036 { 00037 00038 00039 //--------------------------------------------- 00040 // Quadrilateral quadrature rules 00041 case QUAD4: 00042 case QUAD8: 00043 case QUAD9: 00044 { 00045 // We compute the 2D quadrature rule as a tensor 00046 // product of the 1D quadrature rule. 00047 // 00048 // For QUADs, a quadrature rule of order 'p' must be able to integrate 00049 // bilinear (p=1), biquadratic (p=2), bicubic (p=3), etc. polynomials of the form 00050 // 00051 // (x^p + x^{p-1} + ... + 1) * (y^p + y^{p-1} + ... + 1) 00052 // 00053 // These polynomials have terms *up to* degree 2p but they are *not* complete 00054 // polynomials of degree 2p. For example, when p=2 we have 00055 // 1 00056 // x y 00057 // x^2 xy y^2 00058 // yx^2 xy^2 00059 // x^2y^2 00060 QGauss q1D(1,_order); 00061 q1D.init(EDGE2,p); 00062 tensor_product_quad( q1D ); 00063 return; 00064 } 00065 00066 00067 //--------------------------------------------- 00068 // Triangle quadrature rules 00069 case TRI3: 00070 case TRI6: 00071 { 00072 switch(_order + 2*p) 00073 { 00074 case CONSTANT: 00075 case FIRST: 00076 { 00077 // Exact for linears 00078 _points.resize(1); 00079 _weights.resize(1); 00080 00081 _points[0](0) = 1.0L/3.0L; 00082 _points[0](1) = 1.0L/3.0L; 00083 00084 _weights[0] = 0.5; 00085 00086 return; 00087 } 00088 case SECOND: 00089 { 00090 // Exact for quadratics 00091 _points.resize(3); 00092 _weights.resize(3); 00093 00094 // Alternate rule with points on ref. elt. boundaries. 00095 // Not ideal for problems with material coefficient discontinuities 00096 // aligned along element boundaries. 00097 // _points[0](0) = .5; 00098 // _points[0](1) = .5; 00099 // _points[1](0) = 0.; 00100 // _points[1](1) = .5; 00101 // _points[2](0) = .5; 00102 // _points[2](1) = .0; 00103 00104 _points[0](0) = 2.0L/3.0L; 00105 _points[0](1) = 1.0L/6.0L; 00106 00107 _points[1](0) = 1.0L/6.0L; 00108 _points[1](1) = 2.0L/3.0L; 00109 00110 _points[2](0) = 1.0L/6.0L; 00111 _points[2](1) = 1.0L/6.0L; 00112 00113 00114 _weights[0] = 1.0L/6.0L; 00115 _weights[1] = 1.0L/6.0L; 00116 _weights[2] = 1.0L/6.0L; 00117 00118 return; 00119 } 00120 case THIRD: 00121 { 00122 // Exact for cubics 00123 _points.resize(4); 00124 _weights.resize(4); 00125 00126 // This rule is formed from a tensor product of 00127 // appropriately-scaled Gauss and Jacobi rules. (See 00128 // also the QConical quadrature class, this is a 00129 // hard-coded version of one of those rules.) For high 00130 // orders these rules generally have too many points, 00131 // but at extremely low order they are competitive and 00132 // have the additional benefit of having all positive 00133 // weights. 00134 _points[0](0) = 1.5505102572168219018027159252941e-01L; 00135 _points[0](1) = 1.7855872826361642311703513337422e-01L; 00136 _points[1](0) = 6.4494897427831780981972840747059e-01L; 00137 _points[1](1) = 7.5031110222608118177475598324603e-02L; 00138 _points[2](0) = 1.5505102572168219018027159252941e-01L; 00139 _points[2](1) = 6.6639024601470138670269327409637e-01L; 00140 _points[3](0) = 6.4494897427831780981972840747059e-01L; 00141 _points[3](1) = 2.8001991549907407200279599420481e-01L; 00142 00143 _weights[0] = 1.5902069087198858469718450103758e-01L; 00144 _weights[1] = 9.0979309128011415302815498962418e-02L; 00145 _weights[2] = 1.5902069087198858469718450103758e-01L; 00146 _weights[3] = 9.0979309128011415302815498962418e-02L; 00147 00148 return; 00149 00150 00151 // The following third-order rule is quite commonly cited 00152 // in the literature and most likely works fine. However, 00153 // we generally prefer a rule with all positive weights 00154 // and an equal number of points, when available. 00155 // 00156 // (allow_rules_with_negative_weights) 00157 // { 00158 // // Exact for cubics 00159 // _points.resize(4); 00160 // _weights.resize(4); 00161 // 00162 // _points[0](0) = .33333333333333333333333333333333; 00163 // _points[0](1) = .33333333333333333333333333333333; 00164 // 00165 // _points[1](0) = .2; 00166 // _points[1](1) = .6; 00167 // 00168 // _points[2](0) = .2; 00169 // _points[2](1) = .2; 00170 // 00171 // _points[3](0) = .6; 00172 // _points[3](1) = .2; 00173 // 00174 // 00175 // _weights[0] = -27./96.; 00176 // _weights[1] = 25./96.; 00177 // _weights[2] = 25./96.; 00178 // _weights[3] = 25./96.; 00179 // 00180 // return; 00181 // } // end if (allow_rules_with_negative_weights) 00182 // Note: if !allow_rules_with_negative_weights, fall through to next case. 00183 } 00184 00185 00186 00187 // A degree 4 rule with six points. This rule can be found in many places 00188 // including: 00189 // 00190 // J.N. Lyness and D. Jespersen, Moderate degree symmetric 00191 // quadrature rules for the triangle, J. Inst. Math. Appl. 15 (1975), 00192 // 19--32. 00193 // 00194 // We used the code in: 00195 // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules 00196 // on triangles and tetrahedra" Journal of Computational Mathematics, 00197 // v. 27, no. 1, 2009, pp. 89-96. 00198 // to generate additional precision. 00199 case FOURTH: 00200 { 00201 const unsigned int n_wts = 2; 00202 const Real wts[n_wts] = 00203 { 00204 1.1169079483900573284750350421656140e-01L, 00205 5.4975871827660933819163162450105264e-02L 00206 }; 00207 00208 const Real a[n_wts] = 00209 { 00210 4.4594849091596488631832925388305199e-01L, 00211 9.1576213509770743459571463402201508e-02L 00212 }; 00213 00214 const Real b[n_wts] = {0., 0.}; // not used 00215 const unsigned int permutation_ids[n_wts] = {3, 3}; 00216 00217 dunavant_rule2(wts, a, b, permutation_ids, n_wts); // 6 total points 00218 00219 return; 00220 } 00221 00222 00223 00224 // Exact for quintics 00225 // Can be found in "Quadrature on Simplices of Arbitrary 00226 // Dimension" by Walkington. 00227 case FIFTH: 00228 { 00229 const unsigned int n_wts = 3; 00230 const Real wts[n_wts] = 00231 { 00232 static_cast<Real>(9.0L/80.0L), 00233 static_cast<Real>(31.0L/480.0L + std::sqrt(15.0L)/2400.0L), 00234 static_cast<Real>(31.0L/480.0L - std::sqrt(15.0L)/2400.0L) 00235 }; 00236 00237 const Real a[n_wts] = 00238 { 00239 0., // 'a' parameter not used for origin permutation 00240 static_cast<Real>(2.0L/7.0L + std::sqrt(15.0L)/21.0L), 00241 static_cast<Real>(2.0L/7.0L - std::sqrt(15.0L)/21.0L) 00242 }; 00243 00244 const Real b[n_wts] = {0., 0., 0.}; // not used 00245 const unsigned int permutation_ids[n_wts] = {1, 3, 3}; 00246 00247 dunavant_rule2(wts, a, b, permutation_ids, n_wts); // 7 total points 00248 00249 return; 00250 } 00251 00252 00253 00254 // A degree 6 rule with 12 points. This rule can be found in many places 00255 // including: 00256 // 00257 // J.N. Lyness and D. Jespersen, Moderate degree symmetric 00258 // quadrature rules for the triangle, J. Inst. Math. Appl. 15 (1975), 00259 // 19--32. 00260 // 00261 // We used the code in: 00262 // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules 00263 // on triangles and tetrahedra" Journal of Computational Mathematics, 00264 // v. 27, no. 1, 2009, pp. 89-96. 00265 // to generate additional precision. 00266 // 00267 // Note that the following 7th-order Ro3-invariant rule also has only 12 points, 00268 // which technically makes it the superior rule. This one is here for completeness. 00269 case SIXTH: 00270 { 00271 const unsigned int n_wts = 3; 00272 const Real wts[n_wts] = 00273 { 00274 5.8393137863189683012644805692789721e-02L, 00275 2.5422453185103408460468404553434492e-02L, 00276 4.1425537809186787596776728210221227e-02L 00277 }; 00278 00279 const Real a[n_wts] = 00280 { 00281 2.4928674517091042129163855310701908e-01L, 00282 6.3089014491502228340331602870819157e-02L, 00283 3.1035245103378440541660773395655215e-01L 00284 }; 00285 00286 const Real b[n_wts] = 00287 { 00288 0., 00289 0., 00290 6.3650249912139864723014259441204970e-01L 00291 }; 00292 00293 const unsigned int permutation_ids[n_wts] = {3, 3, 6}; // 12 total points 00294 00295 dunavant_rule2(wts, a, b, permutation_ids, n_wts); 00296 00297 return; 00298 } 00299 00300 00301 // A degree 7 rule with 12 points. This rule can be found in: 00302 // 00303 // K. Gatermann, The construction of symmetric cubature 00304 // formulas for the square and the triangle, Computing 40 00305 // (1988), 229--240. 00306 // 00307 // This rule, which is provably minimal in the number of 00308 // integration points, is said to be 'Ro3 invariant' which 00309 // means that a given set of barycentric coordinates 00310 // (z1,z2,z3) implies the quadrature points (z1,z2), 00311 // (z3,z1), (z2,z3) which are formed by taking the first 00312 // two entries in cyclic permutations of the barycentric 00313 // point. Barycentric coordinates are related in the 00314 // sense that: z3 = 1 - z1 - z2. 00315 // 00316 // The 12-point sixth-order rule for triangles given in 00317 // Flaherty's (http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps) 00318 // lecture notes has been removed in favor of this rule 00319 // which is higher-order (for the same number of 00320 // quadrature points) and has a few more digits of 00321 // precision in the points and weights. Some 10-point 00322 // degree 6 rules exist for the triangle but they have 00323 // quadrature points outside the region of integration. 00324 case SEVENTH: 00325 { 00326 _points.resize (12); 00327 _weights.resize(12); 00328 00329 const unsigned int nrows=4; 00330 00331 // In each of the rows below, the first two entries are (z1, z2) which imply 00332 // z3. The third entry is the weight for each of the points in the cyclic permutation. 00333 const Real rule_data[nrows][3] = { 00334 {6.2382265094402118e-02, 6.7517867073916085e-02, 2.6517028157436251e-02}, // group A 00335 {5.5225456656926611e-02, 3.2150249385198182e-01, 4.3881408714446055e-02}, // group B 00336 {3.4324302945097146e-02, 6.6094919618673565e-01, 2.8775042784981585e-02}, // group C 00337 {5.1584233435359177e-01, 2.7771616697639178e-01, 6.7493187009802774e-02} // group D 00338 }; 00339 00340 for (unsigned int i=0, offset=0; i<nrows; ++i) 00341 { 00342 _points[offset + 0] = Point(rule_data[i][0], rule_data[i][1]); // (z1,z2) 00343 _points[offset + 1] = Point(1.-rule_data[i][0]-rule_data[i][1], rule_data[i][0]); // (z3,z1) 00344 _points[offset + 2] = Point(rule_data[i][1], 1.-rule_data[i][0]-rule_data[i][1]); // (z2,z3) 00345 00346 // All these points get the same weight 00347 _weights[offset + 0] = rule_data[i][2]; 00348 _weights[offset + 1] = rule_data[i][2]; 00349 _weights[offset + 2] = rule_data[i][2]; 00350 00351 // Increment offset 00352 offset += 3; 00353 } 00354 00355 return; 00356 00357 00358 // // The following is an inferior 7th-order Lyness-style rule with 15 points. 00359 // // It's here only for completeness and the Ro3-invariant rule above should 00360 // // be used instead! 00361 // const unsigned int n_wts = 3; 00362 // const Real wts[n_wts] = 00363 // { 00364 // 2.6538900895116205835977487499847719e-02L, 00365 // 3.5426541846066783659206291623201826e-02L, 00366 // 3.4637341039708446756138297960207647e-02L 00367 // }; 00368 // 00369 // const Real a[n_wts] = 00370 // { 00371 // 6.4930513159164863078379776030396538e-02L, 00372 // 2.8457558424917033519741605734978046e-01L, 00373 // 3.1355918438493150795585190219862865e-01L 00374 // }; 00375 // 00376 // const Real b[n_wts] = 00377 // { 00378 // 0., 00379 // 1.9838447668150671917987659863332941e-01L, 00380 // 4.3863471792372471511798695971295936e-02L 00381 // }; 00382 // 00383 // const unsigned int permutation_ids[n_wts] = {3, 6, 6}; // 15 total points 00384 // 00385 // dunavant_rule2(wts, a, b, permutation_ids, n_wts); 00386 // 00387 // return; 00388 } 00389 00390 00391 00392 00393 // Another Dunavant rule. This one has all positive weights. This rule has 00394 // 16 points while a comparable conical product rule would have 5*5=25. 00395 // 00396 // It was copied 23rd June 2008 from: 00397 // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90 00398 // 00399 // Additional precision obtained from the code in: 00400 // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules 00401 // on triangles and tetrahedra" Journal of Computational Mathematics, 00402 // v. 27, no. 1, 2009, pp. 89-96. 00403 case EIGHTH: 00404 { 00405 const unsigned int n_wts = 5; 00406 const Real wts[n_wts] = 00407 { 00408 7.2157803838893584125545555244532310e-02L, 00409 4.7545817133642312396948052194292159e-02L, 00410 5.1608685267359125140895775146064515e-02L, 00411 1.6229248811599040155462964170890299e-02L, 00412 1.3615157087217497132422345036954462e-02L 00413 }; 00414 00415 const Real a[n_wts] = 00416 { 00417 0.0, // 'a' parameter not used for origin permutation 00418 4.5929258829272315602881551449416932e-01L, 00419 1.7056930775176020662229350149146450e-01L, 00420 5.0547228317030975458423550596598947e-02L, 00421 2.6311282963463811342178578628464359e-01L, 00422 }; 00423 00424 const Real b[n_wts] = 00425 { 00426 0., 00427 0., 00428 0., 00429 0., 00430 7.2849239295540428124100037917606196e-01L 00431 }; 00432 00433 const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 6}; // 16 total points 00434 00435 dunavant_rule2(wts, a, b, permutation_ids, n_wts); 00436 00437 return; 00438 } 00439 00440 00441 00442 // Another Dunavant rule. This one has all positive weights. This rule has 19 00443 // points. The comparable conical product rule would have 25. 00444 // It was copied 23rd June 2008 from: 00445 // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90 00446 // 00447 // Additional precision obtained from the code in: 00448 // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules 00449 // on triangles and tetrahedra" Journal of Computational Mathematics, 00450 // v. 27, no. 1, 2009, pp. 89-96. 00451 case NINTH: 00452 { 00453 const unsigned int n_wts = 6; 00454 const Real wts[n_wts] = 00455 { 00456 4.8567898141399416909620991253644315e-02L, 00457 1.5667350113569535268427415643604658e-02L, 00458 1.2788837829349015630839399279499912e-02L, 00459 3.8913770502387139658369678149701978e-02L, 00460 3.9823869463605126516445887132022637e-02L, 00461 2.1641769688644688644688644688644689e-02L 00462 }; 00463 00464 const Real a[n_wts] = 00465 { 00466 0.0, // 'a' parameter not used for origin permutation 00467 4.8968251919873762778370692483619280e-01L, 00468 4.4729513394452709865106589966276365e-02L, 00469 4.3708959149293663726993036443535497e-01L, 00470 1.8820353561903273024096128046733557e-01L, 00471 2.2196298916076569567510252769319107e-01L 00472 }; 00473 00474 const Real b[n_wts] = 00475 { 00476 0., 00477 0., 00478 0., 00479 0., 00480 0., 00481 7.4119859878449802069007987352342383e-01L 00482 }; 00483 00484 const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 3, 6}; // 19 total points 00485 00486 dunavant_rule2(wts, a, b, permutation_ids, n_wts); 00487 00488 return; 00489 } 00490 00491 00492 // Another Dunavant rule with all positive weights. This rule has 25 00493 // points. The comparable conical product rule would have 36. 00494 // It was copied 23rd June 2008 from: 00495 // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90 00496 // 00497 // Additional precision obtained from the code in: 00498 // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules 00499 // on triangles and tetrahedra" Journal of Computational Mathematics, 00500 // v. 27, no. 1, 2009, pp. 89-96. 00501 case TENTH: 00502 { 00503 const unsigned int n_wts = 6; 00504 const Real wts[n_wts] = 00505 { 00506 4.5408995191376790047643297550014267e-02L, 00507 1.8362978878233352358503035945683300e-02L, 00508 2.2660529717763967391302822369298659e-02L, 00509 3.6378958422710054302157588309680344e-02L, 00510 1.4163621265528742418368530791049552e-02L, 00511 4.7108334818664117299637354834434138e-03L 00512 }; 00513 00514 const Real a[n_wts] = 00515 { 00516 0.0, // 'a' parameter not used for origin permutation 00517 4.8557763338365737736750753220812615e-01L, 00518 1.0948157548503705479545863134052284e-01L, 00519 3.0793983876412095016515502293063162e-01L, 00520 2.4667256063990269391727646541117681e-01L, 00521 6.6803251012200265773540212762024737e-02L 00522 }; 00523 00524 const Real b[n_wts] = 00525 { 00526 0., 00527 0., 00528 0., 00529 5.5035294182099909507816172659300821e-01L, 00530 7.2832390459741092000873505358107866e-01L, 00531 9.2365593358750027664630697761508843e-01L 00532 }; 00533 00534 const unsigned int permutation_ids[n_wts] = {1, 3, 3, 6, 6, 6}; // 25 total points 00535 00536 dunavant_rule2(wts, a, b, permutation_ids, n_wts); 00537 00538 return; 00539 } 00540 00541 00542 // Dunavant's 11th-order rule contains points outside the region of 00543 // integration, and is thus unacceptable for our FEM calculations. 00544 // 00545 // This 30-point, 11th-order rule was obtained by me [JWP] using the code in 00546 // 00547 // Additional precision obtained from the code in: 00548 // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules 00549 // on triangles and tetrahedra" Journal of Computational Mathematics, 00550 // v. 27, no. 1, 2009, pp. 89-96. 00551 // 00552 // Note: the 28-point 11th-order rule obtained by Zhang in the paper above 00553 // does not appear to be unique. It is a solution in the sense that it 00554 // minimizes the error in the least-squares minimization problem, but 00555 // it involves too many unknowns and the Jacobian is therefore singular 00556 // when attempting to improve the solution via Newton's method. 00557 case ELEVENTH: 00558 { 00559 const unsigned int n_wts = 6; 00560 const Real wts[n_wts] = 00561 { 00562 3.6089021198604635216985338480426484e-02L, 00563 2.1607717807680420303346736867931050e-02L, 00564 3.1144524293927978774861144478241807e-03L, 00565 2.9086855161081509446654185084988077e-02L, 00566 8.4879241614917017182977532679947624e-03L, 00567 1.3795732078224796530729242858347546e-02L 00568 }; 00569 00570 const Real a[n_wts] = 00571 { 00572 3.9355079629947969884346551941969960e-01L, 00573 4.7979065808897448654107733982929214e-01L, 00574 5.1003445645828061436081405648347852e-03L, 00575 2.6597620190330158952732822450744488e-01L, 00576 2.8536418538696461608233522814483715e-01L, 00577 1.3723536747817085036455583801851025e-01L 00578 }; 00579 00580 const Real b[n_wts] = 00581 { 00582 0., 00583 0., 00584 5.6817155788572446538150614865768991e-02L, 00585 1.2539956353662088473247489775203396e-01L, 00586 1.2409970153698532116262152247041742e-02L, 00587 5.2792057988217708934207928630851643e-02L 00588 }; 00589 00590 const unsigned int permutation_ids[n_wts] = {3, 3, 6, 6, 6, 6}; // 30 total points 00591 00592 dunavant_rule2(wts, a, b, permutation_ids, n_wts); 00593 00594 return; 00595 } 00596 00597 00598 00599 00600 // Another Dunavant rule with all positive weights. This rule has 33 00601 // points. The comparable conical product rule would have 36 (ELEVENTH) or 49 (TWELFTH). 00602 // 00603 // It was copied 23rd June 2008 from: 00604 // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90 00605 // 00606 // Additional precision obtained from the code in: 00607 // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules 00608 // on triangles and tetrahedra" Journal of Computational Mathematics, 00609 // v. 27, no. 1, 2009, pp. 89-96. 00610 case TWELFTH: 00611 { 00612 const unsigned int n_wts = 8; 00613 const Real wts[n_wts] = 00614 { 00615 3.0831305257795086169332418926151771e-03L, 00616 3.1429112108942550177135256546441273e-02L, 00617 1.7398056465354471494664198647499687e-02L, 00618 2.1846272269019201067728631278737487e-02L, 00619 1.2865533220227667708895461535782215e-02L, 00620 1.1178386601151722855919538351159995e-02L, 00621 8.6581155543294461858210504055170332e-03L, 00622 2.0185778883190464758914349626118386e-02L 00623 }; 00624 00625 const Real a[n_wts] = 00626 { 00627 2.1317350453210370246856975515728246e-02L, 00628 2.7121038501211592234595134039689474e-01L, 00629 1.2757614554158592467389632515428357e-01L, 00630 4.3972439229446027297973662348436108e-01L, 00631 4.8821738977380488256466206525881104e-01L, 00632 2.8132558098993954824813069297455275e-01L, 00633 1.1625191590759714124135414784260182e-01L, 00634 2.7571326968551419397479634607976398e-01L 00635 }; 00636 00637 const Real b[n_wts] = 00638 { 00639 0., 00640 0., 00641 0., 00642 0., 00643 0., 00644 6.9583608678780342214163552323607254e-01L, 00645 8.5801403354407263059053661662617818e-01L, 00646 6.0894323577978780685619243776371007e-01L 00647 }; 00648 00649 const unsigned int permutation_ids[n_wts] = {3, 3, 3, 3, 3, 6, 6, 6}; // 33 total points 00650 00651 dunavant_rule2(wts, a, b, permutation_ids, n_wts); 00652 00653 return; 00654 } 00655 00656 00657 // Another Dunavant rule with all positive weights. This rule has 37 00658 // points. The comparable conical product rule would have 49 points. 00659 // 00660 // It was copied 23rd June 2008 from: 00661 // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90 00662 // 00663 // A second rule with additional precision obtained from the code in: 00664 // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules 00665 // on triangles and tetrahedra" Journal of Computational Mathematics, 00666 // v. 27, no. 1, 2009, pp. 89-96. 00667 case THIRTEENTH: 00668 { 00669 const unsigned int n_wts = 9; 00670 const Real wts[n_wts] = 00671 { 00672 3.3980018293415822140887212340442440e-02L, 00673 2.7800983765226664353628733005230734e-02L, 00674 2.9139242559599990702383541756669905e-02L, 00675 3.0261685517695859208964000161454122e-03L, 00676 1.1997200964447365386855399725479827e-02L, 00677 1.7320638070424185232993414255459110e-02L, 00678 7.4827005525828336316229285664517190e-03L, 00679 1.2089519905796909568722872786530380e-02L, 00680 4.7953405017716313612975450830554457e-03L 00681 }; 00682 00683 const Real a[n_wts] = 00684 { 00685 0., // 'a' parameter not used for origin permutation 00686 4.2694141425980040602081253503137421e-01L, 00687 2.2137228629183290065481255470507908e-01L, 00688 2.1509681108843183869291313534052083e-02L, 00689 4.8907694645253934990068971909020439e-01L, 00690 3.0844176089211777465847185254124531e-01L, 00691 1.1092204280346339541286954522167452e-01L, 00692 1.6359740106785048023388790171095725e-01L, 00693 2.7251581777342966618005046435408685e-01L 00694 }; 00695 00696 const Real b[n_wts] = 00697 { 00698 0., 00699 0., 00700 0., 00701 0., 00702 0., 00703 6.2354599555367557081585435318623659e-01L, 00704 8.6470777029544277530254595089569318e-01L, 00705 7.4850711589995219517301859578870965e-01L, 00706 7.2235779312418796526062013230478405e-01L 00707 }; 00708 00709 const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 3, 6, 6, 6, 6}; // 37 total points 00710 00711 dunavant_rule2(wts, a, b, permutation_ids, n_wts); 00712 00713 return; 00714 } 00715 00716 00717 // Another Dunavant rule. This rule has 42 points, while 00718 // a comparable conical product rule would have 64. 00719 // 00720 // It was copied 23rd June 2008 from: 00721 // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90 00722 // 00723 // Additional precision obtained from the code in: 00724 // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules 00725 // on triangles and tetrahedra" Journal of Computational Mathematics, 00726 // v. 27, no. 1, 2009, pp. 89-96. 00727 case FOURTEENTH: 00728 { 00729 const unsigned int n_wts = 10; 00730 const Real wts[n_wts] = 00731 { 00732 1.0941790684714445320422472981662986e-02L, 00733 1.6394176772062675320655489369312672e-02L, 00734 2.5887052253645793157392455083198201e-02L, 00735 2.1081294368496508769115218662093065e-02L, 00736 7.2168498348883338008549607403266583e-03L, 00737 2.4617018012000408409130117545210774e-03L, 00738 1.2332876606281836981437622591818114e-02L, 00739 1.9285755393530341614244513905205430e-02L, 00740 7.2181540567669202480443459995079017e-03L, 00741 2.5051144192503358849300465412445582e-03L 00742 }; 00743 00744 const Real a[n_wts] = 00745 { 00746 4.8896391036217863867737602045239024e-01L, 00747 4.1764471934045392250944082218564344e-01L, 00748 2.7347752830883865975494428326269856e-01L, 00749 1.7720553241254343695661069046505908e-01L, 00750 6.1799883090872601267478828436935788e-02L, 00751 1.9390961248701048178250095054529511e-02L, 00752 1.7226668782135557837528960161365733e-01L, 00753 3.3686145979634500174405519708892539e-01L, 00754 2.9837288213625775297083151805961273e-01L, 00755 1.1897449769695684539818196192990548e-01L 00756 }; 00757 00758 const Real b[n_wts] = 00759 { 00760 0., 00761 0., 00762 0., 00763 0., 00764 0., 00765 0., 00766 7.7060855477499648258903327416742796e-01L, 00767 5.7022229084668317349769621336235426e-01L, 00768 6.8698016780808783735862715402031306e-01L, 00769 8.7975717137017112951457163697460183e-01L 00770 }; 00771 00772 const unsigned int permutation_ids[n_wts] 00773 = {3, 3, 3, 3, 3, 3, 6, 6, 6, 6}; // 42 total points 00774 00775 dunavant_rule2(wts, a, b, permutation_ids, n_wts); 00776 00777 return; 00778 } 00779 00780 00781 // This 49-point rule was found by me [JWP] using the code in: 00782 // 00783 // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules 00784 // on triangles and tetrahedra" Journal of Computational Mathematics, 00785 // v. 27, no. 1, 2009, pp. 89-96. 00786 // 00787 // A 54-point, 15th-order rule is reported by 00788 // 00789 // Stephen Wandzura, Hong Xiao, 00790 // Symmetric Quadrature Rules on a Triangle, 00791 // Computers and Mathematics with Applications, 00792 // Volume 45, Number 12, June 2003, pages 1829-1840. 00793 // 00794 // can be found here: 00795 // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90 00796 // 00797 // but this 49-point rule is superior. 00798 case FIFTEENTH: 00799 { 00800 const unsigned int n_wts = 11; 00801 const Real wts[n_wts] = 00802 { 00803 2.4777380743035579804788826970198951e-02L, 00804 9.2433943023307730591540642828347660e-03L, 00805 2.2485768962175402793245929133296627e-03L, 00806 6.7052581900064143760518398833360903e-03L, 00807 1.9011381726930579256700190357527956e-02L, 00808 1.4605445387471889398286155981802858e-02L, 00809 1.5087322572773133722829435011138258e-02L, 00810 1.5630213780078803020711746273129099e-02L, 00811 6.1808086085778203192616856133701233e-03L, 00812 3.2209366452594664857296985751120513e-03L, 00813 5.8747373242569702667677969985668817e-03L 00814 }; 00815 00816 const Real a[n_wts] = 00817 { 00818 0.0, // 'a' parameter not used for origin 00819 7.9031013655541635005816956762252155e-02L, 00820 1.8789501810770077611247984432284226e-02L, 00821 4.9250168823249670532514526605352905e-01L, 00822 4.0886316907744105975059040108092775e-01L, 00823 5.3877851064220142445952549348423733e-01L, 00824 2.0250549804829997692885033941362673e-01L, 00825 5.5349674918711643207148086558288110e-01L, 00826 7.8345022567320812359258882143250181e-01L, 00827 8.9514624528794883409864566727625002e-01L, 00828 3.2515745241110782862789881780746490e-01L 00829 }; 00830 00831 const Real b[n_wts] = 00832 { 00833 0., 00834 0., 00835 0., 00836 0., 00837 0., 00838 1.9412620368774630292701241080996842e-01L, 00839 9.8765911355712115933807754318089099e-02L, 00840 7.7663767064308164090246588765178087e-02L, 00841 2.1594628433980258573654682690950798e-02L, 00842 1.2563596287784997705599005477153617e-02L, 00843 1.5082654870922784345283124845552190e-02L 00844 }; 00845 00846 const unsigned int permutation_ids[n_wts] 00847 = {1, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6}; // 49 total points 00848 00849 dunavant_rule2(wts, a, b, permutation_ids, n_wts); 00850 00851 return; 00852 } 00853 00854 00855 00856 00857 // Dunavant's 16th-order rule contains points outside the region of 00858 // integration, and is thus unacceptable for our FEM calculations. 00859 // 00860 // This 55-point, 16th-order rule was obtained by me [JWP] using the code in 00861 // 00862 // Additional precision obtained from the code in: 00863 // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules 00864 // on triangles and tetrahedra" Journal of Computational Mathematics, 00865 // v. 27, no. 1, 2009, pp. 89-96. 00866 // 00867 // Note: the 55-point 16th-order rule obtained by Zhang in the paper above 00868 // does not appear to be unique. It is a solution in the sense that it 00869 // minimizes the error in the least-squares minimization problem, but 00870 // it involves too many unknowns and the Jacobian is therefore singular 00871 // when attempting to improve the solution via Newton's method. 00872 case SIXTEENTH: 00873 { 00874 const unsigned int n_wts = 12; 00875 const Real wts[n_wts] = 00876 { 00877 2.2668082505910087151996321171534230e-02L, 00878 8.4043060714818596159798961899306135e-03L, 00879 1.0850949634049747713966288634484161e-03L, 00880 7.2252773375423638869298219383808751e-03L, 00881 1.2997715227338366024036316182572871e-02L, 00882 2.0054466616677715883228810959112227e-02L, 00883 9.7299841600417010281624372720122710e-03L, 00884 1.1651974438298104227427176444311766e-02L, 00885 9.1291185550484450744725847363097389e-03L, 00886 3.5568614040947150231712567900113671e-03L, 00887 5.8355861686234326181790822005304303e-03L, 00888 4.7411314396804228041879331486234396e-03L 00889 }; 00890 00891 const Real a[n_wts] = 00892 { 00893 0.0, // 'a' parameter not used for centroid weight 00894 8.5402539407933203673769900926355911e-02L, 00895 1.2425572001444092841183633409631260e-02L, 00896 4.9174838341891594024701017768490960e-01L, 00897 4.5669426695387464162068900231444462e-01L, 00898 4.8506759880447437974189793537259677e-01L, 00899 2.0622099278664205707909858461264083e-01L, 00900 3.2374950270039093446805340265853956e-01L, 00901 7.3834330556606586255186213302750029e-01L, 00902 9.1210673061680792565673823935174611e-01L, 00903 6.6129919222598721544966837350891531e-01L, 00904 1.7807138906021476039088828811346122e-01L 00905 }; 00906 00907 const Real b[n_wts] = 00908 { 00909 0.0, 00910 0.0, 00911 0.0, 00912 0.0, 00913 0.0, 00914 3.2315912848634384647700266402091638e-01L, 00915 1.5341553679414688425981898952416987e-01L, 00916 7.4295478991330687632977899141707872e-02L, 00917 7.1278762832147862035977841733532020e-02L, 00918 1.6623223223705792825395256602140459e-02L, 00919 1.4160772533794791868984026749196156e-02L, 00920 1.4539694958941854654807449467759690e-02L 00921 }; 00922 00923 const unsigned int permutation_ids[n_wts] 00924 = {1, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6, 6}; // 55 total points 00925 00926 dunavant_rule2(wts, a, b, permutation_ids, n_wts); 00927 00928 return; 00929 } 00930 00931 00932 // Dunavant's 17th-order rule has 61 points, while a 00933 // comparable conical product rule would have 81 (16th and 17th orders). 00934 // 00935 // It can be found here: 00936 // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90 00937 // 00938 // Zhang reports an identical rule in: 00939 // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules 00940 // on triangles and tetrahedra" Journal of Computational Mathematics, 00941 // v. 27, no. 1, 2009, pp. 89-96. 00942 // 00943 // Note: the 61-point 17th-order rule obtained by Dunavant and Zhang 00944 // does not appear to be unique. It is a solution in the sense that it 00945 // minimizes the error in the least-squares minimization problem, but 00946 // it involves too many unknowns and the Jacobian is therefore singular 00947 // when attempting to improve the solution via Newton's method. 00948 // 00949 // Therefore, we prefer the following 63-point rule which 00950 // I [JWP] found. It appears to be more accurate than the 00951 // rule reported by Dunavant and Zhang, even though it has 00952 // a few more points. 00953 case SEVENTEENTH: 00954 { 00955 const unsigned int n_wts = 12; 00956 const Real wts[n_wts] = 00957 { 00958 1.7464603792572004485690588092246146e-02L, 00959 5.9429003555801725246549713984660076e-03L, 00960 1.2490753345169579649319736639588729e-02L, 00961 1.5386987188875607593083456905596468e-02L, 00962 1.1185807311917706362674684312990270e-02L, 00963 1.0301845740670206831327304917180007e-02L, 00964 1.1767783072977049696840016810370464e-02L, 00965 3.8045312849431209558329128678945240e-03L, 00966 4.5139302178876351271037137230354382e-03L, 00967 2.2178812517580586419412547665472893e-03L, 00968 5.2216271537483672304731416553063103e-03L, 00969 9.8381136389470256422419930926212114e-04L 00970 }; 00971 00972 const Real a[n_wts] = 00973 { 00974 2.8796825754667362165337965123570514e-01L, 00975 4.9216175986208465345536805750663939e-01L, 00976 4.6252866763171173685916780827044612e-01L, 00977 1.6730292951631792248498303276090273e-01L, 00978 1.5816335500814652972296428532213019e-01L, 00979 1.6352252138387564873002458959679529e-01L, 00980 6.2447680488959768233910286168417367e-01L, 00981 8.7317249935244454285263604347964179e-01L, 00982 3.4428164322282694677972239461699271e-01L, 00983 9.1584484467813674010523309855340209e-02L, 00984 2.0172088013378989086826623852040632e-01L, 00985 9.6538762758254643474731509845084691e-01L 00986 }; 00987 00988 const Real b[n_wts] = 00989 { 00990 0.0, 00991 0.0, 00992 0.0, 00993 3.4429160695501713926320695771253348e-01L, 00994 2.2541623431550639817203145525444726e-01L, 00995 8.0670083153531811694942222940484991e-02L, 00996 6.5967451375050925655738829747288190e-02L, 00997 4.5677879890996762665044366994439565e-02L, 00998 1.1528411723154215812386518751976084e-02L, 00999 9.3057714323900610398389176844165892e-03L, 01000 1.5916814107619812717966560404970160e-02L, 01001 1.0734733163764032541125434215228937e-02L 01002 }; 01003 01004 const unsigned int permutation_ids[n_wts] 01005 = {3, 3, 3, 6, 6, 6, 6, 6, 6, 6, 6, 6}; // 63 total points 01006 01007 dunavant_rule2(wts, a, b, permutation_ids, n_wts); 01008 01009 return; 01010 01011 // _points.resize (61); 01012 // _weights.resize(61); 01013 01014 // // The raw data for the quadrature rule. 01015 // const Real p[15][4] = { 01016 // { 1./3., 0., 0., 0.033437199290803e+00 / 2.0}, // 1-perm 01017 // {0.005658918886452e+00, 0.497170540556774e+00, 0., 0.005093415440507e+00 / 2.0}, // 3-perm 01018 // {0.035647354750751e+00, 0.482176322624625e+00, 0., 0.014670864527638e+00 / 2.0}, // 3-perm 01019 // {0.099520061958437e+00, 0.450239969020782e+00, 0., 0.024350878353672e+00 / 2.0}, // 3-perm 01020 // {0.199467521245206e+00, 0.400266239377397e+00, 0., 0.031107550868969e+00 / 2.0}, // 3-perm 01021 // {0.495717464058095e+00, 0.252141267970953e+00, 0., 0.031257111218620e+00 / 2.0}, // 3-perm 01022 // {0.675905990683077e+00, 0.162047004658461e+00, 0., 0.024815654339665e+00 / 2.0}, // 3-perm 01023 // {0.848248235478508e+00, 0.075875882260746e+00, 0., 0.014056073070557e+00 / 2.0}, // 3-perm 01024 // {0.968690546064356e+00, 0.015654726967822e+00, 0., 0.003194676173779e+00 / 2.0}, // 3-perm 01025 // {0.010186928826919e+00, 0.334319867363658e+00, 0.655493203809423e+00, 0.008119655318993e+00 / 2.0}, // 6-perm 01026 // {0.135440871671036e+00, 0.292221537796944e+00, 0.572337590532020e+00, 0.026805742283163e+00 / 2.0}, // 6-perm 01027 // {0.054423924290583e+00, 0.319574885423190e+00, 0.626001190286228e+00, 0.018459993210822e+00 / 2.0}, // 6-perm 01028 // {0.012868560833637e+00, 0.190704224192292e+00, 0.796427214974071e+00, 0.008476868534328e+00 / 2.0}, // 6-perm 01029 // {0.067165782413524e+00, 0.180483211648746e+00, 0.752351005937729e+00, 0.018292796770025e+00 / 2.0}, // 6-perm 01030 // {0.014663182224828e+00, 0.080711313679564e+00, 0.904625504095608e+00, 0.006665632004165e+00 / 2.0} // 6-perm 01031 // }; 01032 01033 01034 // // Now call the dunavant routine to generate _points and _weights 01035 // dunavant_rule(p, 15); 01036 01037 // return; 01038 } 01039 01040 01041 01042 // Dunavant's 18th-order rule contains points outside the region and is therefore unsuitable 01043 // for our FEM calculations. His 19th-order rule has 73 points, compared with 100 points for 01044 // a comparable-order conical product rule. 01045 // 01046 // It was copied 23rd June 2008 from: 01047 // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90 01048 case EIGHTTEENTH: 01049 case NINTEENTH: 01050 { 01051 _points.resize (73); 01052 _weights.resize(73); 01053 01054 // The raw data for the quadrature rule. 01055 const Real rule_data[17][4] = { 01056 { 1./3., 0., 0., 0.032906331388919e+00 / 2.0}, // 1-perm 01057 {0.020780025853987e+00, 0.489609987073006e+00, 0., 0.010330731891272e+00 / 2.0}, // 3-perm 01058 {0.090926214604215e+00, 0.454536892697893e+00, 0., 0.022387247263016e+00 / 2.0}, // 3-perm 01059 {0.197166638701138e+00, 0.401416680649431e+00, 0., 0.030266125869468e+00 / 2.0}, // 3-perm 01060 {0.488896691193805e+00, 0.255551654403098e+00, 0., 0.030490967802198e+00 / 2.0}, // 3-perm 01061 {0.645844115695741e+00, 0.177077942152130e+00, 0., 0.024159212741641e+00 / 2.0}, // 3-perm 01062 {0.779877893544096e+00, 0.110061053227952e+00, 0., 0.016050803586801e+00 / 2.0}, // 3-perm 01063 {0.888942751496321e+00, 0.055528624251840e+00, 0., 0.008084580261784e+00 / 2.0}, // 3-perm 01064 {0.974756272445543e+00, 0.012621863777229e+00, 0., 0.002079362027485e+00 / 2.0}, // 3-perm 01065 {0.003611417848412e+00, 0.395754787356943e+00, 0.600633794794645e+00, 0.003884876904981e+00 / 2.0}, // 6-perm 01066 {0.134466754530780e+00, 0.307929983880436e+00, 0.557603261588784e+00, 0.025574160612022e+00 / 2.0}, // 6-perm 01067 {0.014446025776115e+00, 0.264566948406520e+00, 0.720987025817365e+00, 0.008880903573338e+00 / 2.0}, // 6-perm 01068 {0.046933578838178e+00, 0.358539352205951e+00, 0.594527068955871e+00, 0.016124546761731e+00 / 2.0}, // 6-perm 01069 {0.002861120350567e+00, 0.157807405968595e+00, 0.839331473680839e+00, 0.002491941817491e+00 / 2.0}, // 6-perm 01070 {0.223861424097916e+00, 0.075050596975911e+00, 0.701087978926173e+00, 0.018242840118951e+00 / 2.0}, // 6-perm 01071 {0.034647074816760e+00, 0.142421601113383e+00, 0.822931324069857e+00, 0.010258563736199e+00 / 2.0}, // 6-perm 01072 {0.010161119296278e+00, 0.065494628082938e+00, 0.924344252620784e+00, 0.003799928855302e+00 / 2.0} // 6-perm 01073 }; 01074 01075 01076 // Now call the dunavant routine to generate _points and _weights 01077 dunavant_rule(rule_data, 17); 01078 01079 return; 01080 } 01081 01082 01083 // 20th-order rule by Wandzura. 01084 // 01085 // Stephen Wandzura, Hong Xiao, 01086 // Symmetric Quadrature Rules on a Triangle, 01087 // Computers and Mathematics with Applications, 01088 // Volume 45, Number 12, June 2003, pages 1829-1840. 01089 // 01090 // Wandzura's work extends the work of Dunavant by providing degree 01091 // 5,10,15,20,25, and 30 rules with positive weights for the triangle. 01092 // 01093 // Copied on 3rd July 2008 from: 01094 // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90 01095 case TWENTIETH: 01096 { 01097 // The equivalent concial product rule would have 121 points 01098 _points.resize (85); 01099 _weights.resize(85); 01100 01101 // The raw data for the quadrature rule. 01102 const Real rule_data[19][4] = { 01103 {0.33333333333333e+00, 0.0, 0.0, 0.2761042699769952e-01 / 2.0}, // 1-perm 01104 {0.00150064932443e+00, 0.49924967533779e+00, 0.0, 0.1779029547326740e-02 / 2.0}, // 3-perm 01105 {0.09413975193895e+00, 0.45293012403052e+00, 0.0, 0.2011239811396117e-01 / 2.0}, // 3-perm 01106 {0.20447212408953e+00, 0.39776393795524e+00, 0.0, 0.2681784725933157e-01 / 2.0}, // 3-perm 01107 {0.47099959493443e+00, 0.26450020253279e+00, 0.0, 0.2452313380150201e-01 / 2.0}, // 3-perm 01108 {0.57796207181585e+00, 0.21101896409208e+00, 0.0, 0.1639457841069539e-01 / 2.0}, // 3-perm 01109 {0.78452878565746e+00, 0.10773560717127e+00, 0.0, 0.1479590739864960e-01 / 2.0}, // 3-perm 01110 {0.92186182432439e+00, 0.03906908783780e+00, 0.0, 0.4579282277704251e-02 / 2.0}, // 3-perm 01111 {0.97765124054134e+00, 0.01117437972933e+00, 0.0, 0.1651826515576217e-02 / 2.0}, // 3-perm 01112 {0.00534961818734e+00, 0.06354966590835e+00, 0.93110071590431e+00, 0.2349170908575584e-02 / 2.0}, // 6-perm 01113 {0.00795481706620e+00, 0.15710691894071e+00, 0.83493826399309e+00, 0.4465925754181793e-02 / 2.0}, // 6-perm 01114 {0.01042239828126e+00, 0.39564211436437e+00, 0.59393548735436e+00, 0.6099566807907972e-02 / 2.0}, // 6-perm 01115 {0.01096441479612e+00, 0.27316757071291e+00, 0.71586801449097e+00, 0.6891081327188203e-02 / 2.0}, // 6-perm 01116 {0.03856671208546e+00, 0.10178538248502e+00, 0.85964790542952e+00, 0.7997475072478163e-02 / 2.0}, // 6-perm 01117 {0.03558050781722e+00, 0.44665854917641e+00, 0.51776094300637e+00, 0.7386134285336024e-02 / 2.0}, // 6-perm 01118 {0.04967081636276e+00, 0.19901079414950e+00, 0.75131838948773e+00, 0.1279933187864826e-01 / 2.0}, // 6-perm 01119 {0.05851972508433e+00, 0.32426118369228e+00, 0.61721909122339e+00, 0.1725807117569655e-01 / 2.0}, // 6-perm 01120 {0.12149778700439e+00, 0.20853136321013e+00, 0.66997084978547e+00, 0.1867294590293547e-01 / 2.0}, // 6-perm 01121 {0.14071084494394e+00, 0.32317056653626e+00, 0.53611858851980e+00, 0.2281822405839526e-01 / 2.0} // 6-perm 01122 }; 01123 01124 01125 // Now call the dunavant routine to generate _points and _weights 01126 dunavant_rule(rule_data, 19); 01127 01128 return; 01129 } 01130 01131 01132 01133 // 25th-order rule by Wandzura. 01134 // 01135 // Stephen Wandzura, Hong Xiao, 01136 // Symmetric Quadrature Rules on a Triangle, 01137 // Computers and Mathematics with Applications, 01138 // Volume 45, Number 12, June 2003, pages 1829-1840. 01139 // 01140 // Wandzura's work extends the work of Dunavant by providing degree 01141 // 5,10,15,20,25, and 30 rules with positive weights for the triangle. 01142 // 01143 // Copied on 3rd July 2008 from: 01144 // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90 01145 // case TWENTYFIRST: // fall through to 121 point conical product rule below 01146 case TWENTYSECOND: 01147 case TWENTYTHIRD: 01148 case TWENTYFOURTH: 01149 case TWENTYFIFTH: 01150 { 01151 // The equivalent concial product rule would have 169 points 01152 _points.resize (126); 01153 _weights.resize(126); 01154 01155 // The raw data for the quadrature rule. 01156 const Real rule_data[26][4] = { 01157 {0.02794648307317e+00, 0.48602675846341e+00, 0.0, 0.8005581880020417e-02 / 2.0}, // 3-perm 01158 {0.13117860132765e+00, 0.43441069933617e+00, 0.0, 0.1594707683239050e-01 / 2.0}, // 3-perm 01159 {0.22022172951207e+00, 0.38988913524396e+00, 0.0, 0.1310914123079553e-01 / 2.0}, // 3-perm 01160 {0.40311353196039e+00, 0.29844323401980e+00, 0.0, 0.1958300096563562e-01 / 2.0}, // 3-perm 01161 {0.53191165532526e+00, 0.23404417233737e+00, 0.0, 0.1647088544153727e-01 / 2.0}, // 3-perm 01162 {0.69706333078196e+00, 0.15146833460902e+00, 0.0, 0.8547279074092100e-02 / 2.0}, // 3-perm 01163 {0.77453221290801e+00, 0.11273389354599e+00, 0.0, 0.8161885857226492e-02 / 2.0}, // 3-perm 01164 {0.84456861581695e+00, 0.07771569209153e+00, 0.0, 0.6121146539983779e-02 / 2.0}, // 3-perm 01165 {0.93021381277141e+00, 0.03489309361430e+00, 0.0, 0.2908498264936665e-02 / 2.0}, // 3-perm 01166 {0.98548363075813e+00, 0.00725818462093e+00, 0.0, 0.6922752456619963e-03 / 2.0}, // 3-perm 01167 {0.00129235270444e+00, 0.22721445215336e+00, 0.77149319514219e+00, 0.1248289199277397e-02 / 2.0}, // 6-perm 01168 {0.00539970127212e+00, 0.43501055485357e+00, 0.55958974387431e+00, 0.3404752908803022e-02 / 2.0}, // 6-perm 01169 {0.00638400303398e+00, 0.32030959927220e+00, 0.67330639769382e+00, 0.3359654326064051e-02 / 2.0}, // 6-perm 01170 {0.00502821150199e+00, 0.09175032228001e+00, 0.90322146621800e+00, 0.1716156539496754e-02 / 2.0}, // 6-perm 01171 {0.00682675862178e+00, 0.03801083585872e+00, 0.95516240551949e+00, 0.1480856316715606e-02 / 2.0}, // 6-perm 01172 {0.01001619963993e+00, 0.15742521848531e+00, 0.83255858187476e+00, 0.3511312610728685e-02 / 2.0}, // 6-perm 01173 {0.02575781317339e+00, 0.23988965977853e+00, 0.73435252704808e+00, 0.7393550149706484e-02 / 2.0}, // 6-perm 01174 {0.03022789811992e+00, 0.36194311812606e+00, 0.60782898375402e+00, 0.7983087477376558e-02 / 2.0}, // 6-perm 01175 {0.03050499010716e+00, 0.08355196095483e+00, 0.88594304893801e+00, 0.4355962613158041e-02 / 2.0}, // 6-perm 01176 {0.04595654736257e+00, 0.14844322073242e+00, 0.80560023190501e+00, 0.7365056701417832e-02 / 2.0}, // 6-perm 01177 {0.06744280054028e+00, 0.28373970872753e+00, 0.64881749073219e+00, 0.1096357284641955e-01 / 2.0}, // 6-perm 01178 {0.07004509141591e+00, 0.40689937511879e+00, 0.52305553346530e+00, 0.1174996174354112e-01 / 2.0}, // 6-perm 01179 {0.08391152464012e+00, 0.19411398702489e+00, 0.72197448833499e+00, 0.1001560071379857e-01 / 2.0}, // 6-perm 01180 {0.12037553567715e+00, 0.32413434700070e+00, 0.55549011732214e+00, 0.1330964078762868e-01 / 2.0}, // 6-perm 01181 {0.14806689915737e+00, 0.22927748355598e+00, 0.62265561728665e+00, 0.1415444650522614e-01 / 2.0}, // 6-perm 01182 {0.19177186586733e+00, 0.32561812259598e+00, 0.48261001153669e+00, 0.1488137956116801e-01 / 2.0} // 6-perm 01183 }; 01184 01185 01186 // Now call the dunavant routine to generate _points and _weights 01187 dunavant_rule(rule_data, 26); 01188 01189 return; 01190 } 01191 01192 01193 01194 // 30th-order rule by Wandzura. 01195 // 01196 // Stephen Wandzura, Hong Xiao, 01197 // Symmetric Quadrature Rules on a Triangle, 01198 // Computers and Mathematics with Applications, 01199 // Volume 45, Number 12, June 2003, pages 1829-1840. 01200 // 01201 // Wandzura's work extends the work of Dunavant by providing degree 01202 // 5,10,15,20,25, and 30 rules with positive weights for the triangle. 01203 // 01204 // Copied on 3rd July 2008 from: 01205 // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90 01206 case TWENTYSIXTH: 01207 case TWENTYSEVENTH: 01208 case TWENTYEIGHTH: 01209 case TWENTYNINTH: 01210 case THIRTIETH: 01211 { 01212 // The equivalent concial product rule would have 256 points 01213 _points.resize (175); 01214 _weights.resize(175); 01215 01216 // The raw data for the quadrature rule. 01217 const Real rule_data[36][4] = { 01218 {0.33333333333333e+00, 0.0, 0.0, 0.1557996020289920e-01 / 2.0}, // 1-perm 01219 {0.00733011643277e+00, 0.49633494178362e+00, 0.0, 0.3177233700534134e-02 / 2.0}, // 3-perm 01220 {0.08299567580296e+00, 0.45850216209852e+00, 0.0, 0.1048342663573077e-01 / 2.0}, // 3-perm 01221 {0.15098095612541e+00, 0.42450952193729e+00, 0.0, 0.1320945957774363e-01 / 2.0}, // 3-perm 01222 {0.23590585989217e+00, 0.38204707005392e+00, 0.0, 0.1497500696627150e-01 / 2.0}, // 3-perm 01223 {0.43802430840785e+00, 0.28098784579608e+00, 0.0, 0.1498790444338419e-01 / 2.0}, // 3-perm 01224 {0.54530204829193e+00, 0.22734897585403e+00, 0.0, 0.1333886474102166e-01 / 2.0}, // 3-perm 01225 {0.65088177698254e+00, 0.17455911150873e+00, 0.0, 0.1088917111390201e-01 / 2.0}, // 3-perm 01226 {0.75348314559713e+00, 0.12325842720144e+00, 0.0, 0.8189440660893461e-02 / 2.0}, // 3-perm 01227 {0.83983154221561e+00, 0.08008422889220e+00, 0.0, 0.5575387588607785e-02 / 2.0}, // 3-perm 01228 {0.90445106518420e+00, 0.04777446740790e+00, 0.0, 0.3191216473411976e-02 / 2.0}, // 3-perm 01229 {0.95655897063972e+00, 0.02172051468014e+00, 0.0, 0.1296715144327045e-02 / 2.0}, // 3-perm 01230 {0.99047064476913e+00, 0.00476467761544e+00, 0.0, 0.2982628261349172e-03 / 2.0}, // 3-perm 01231 {0.00092537119335e+00, 0.41529527091331e+00, 0.58377935789334e+00, 0.9989056850788964e-03 / 2.0}, // 6-perm 01232 {0.00138592585556e+00, 0.06118990978535e+00, 0.93742416435909e+00, 0.4628508491732533e-03 / 2.0}, // 6-perm 01233 {0.00368241545591e+00, 0.16490869013691e+00, 0.83140889440718e+00, 0.1234451336382413e-02 / 2.0}, // 6-perm 01234 {0.00390322342416e+00, 0.02503506223200e+00, 0.97106171434384e+00, 0.5707198522432062e-03 / 2.0}, // 6-perm 01235 {0.00323324815501e+00, 0.30606446515110e+00, 0.69070228669389e+00, 0.1126946125877624e-02 / 2.0}, // 6-perm 01236 {0.00646743211224e+00, 0.10707328373022e+00, 0.88645928415754e+00, 0.1747866949407337e-02 / 2.0}, // 6-perm 01237 {0.00324747549133e+00, 0.22995754934558e+00, 0.76679497516308e+00, 0.1182818815031657e-02 / 2.0}, // 6-perm 01238 {0.00867509080675e+00, 0.33703663330578e+00, 0.65428827588746e+00, 0.1990839294675034e-02 / 2.0}, // 6-perm 01239 {0.01559702646731e+00, 0.05625657618206e+00, 0.92814639735063e+00, 0.1900412795035980e-02 / 2.0}, // 6-perm 01240 {0.01797672125369e+00, 0.40245137521240e+00, 0.57957190353391e+00, 0.4498365808817451e-02 / 2.0}, // 6-perm 01241 {0.01712424535389e+00, 0.24365470201083e+00, 0.73922105263528e+00, 0.3478719460274719e-02 / 2.0}, // 6-perm 01242 {0.02288340534658e+00, 0.16538958561453e+00, 0.81172700903888e+00, 0.4102399036723953e-02 / 2.0}, // 6-perm 01243 {0.03273759728777e+00, 0.09930187449585e+00, 0.86796052821639e+00, 0.4021761549744162e-02 / 2.0}, // 6-perm 01244 {0.03382101234234e+00, 0.30847833306905e+00, 0.65770065458860e+00, 0.6033164660795066e-02 / 2.0}, // 6-perm 01245 {0.03554761446002e+00, 0.46066831859211e+00, 0.50378406694787e+00, 0.3946290302129598e-02 / 2.0}, // 6-perm 01246 {0.05053979030687e+00, 0.21881529945393e+00, 0.73064491023920e+00, 0.6644044537680268e-02 / 2.0}, // 6-perm 01247 {0.05701471491573e+00, 0.37920955156027e+00, 0.56377573352399e+00, 0.8254305856078458e-02 / 2.0}, // 6-perm 01248 {0.06415280642120e+00, 0.14296081941819e+00, 0.79288637416061e+00, 0.6496056633406411e-02 / 2.0}, // 6-perm 01249 {0.08050114828763e+00, 0.28373128210592e+00, 0.63576756960645e+00, 0.9252778144146602e-02 / 2.0}, // 6-perm 01250 {0.10436706813453e+00, 0.19673744100444e+00, 0.69889549086103e+00, 0.9164920726294280e-02 / 2.0}, // 6-perm 01251 {0.11384489442875e+00, 0.35588914121166e+00, 0.53026596435959e+00, 0.1156952462809767e-01 / 2.0}, // 6-perm 01252 {0.14536348771552e+00, 0.25981868535191e+00, 0.59481782693256e+00, 0.1176111646760917e-01 / 2.0}, // 6-perm 01253 {0.18994565282198e+00, 0.32192318123130e+00, 0.48813116594672e+00, 0.1382470218216540e-01 / 2.0} // 6-perm 01254 }; 01255 01256 01257 // Now call the dunavant routine to generate _points and _weights 01258 dunavant_rule(rule_data, 36); 01259 01260 return; 01261 } 01262 01263 01264 // By default, we fall back on the conical product rules. If the user 01265 // requests an order higher than what is currently available in the 1D 01266 // rules, an error will be thrown from the respective 1D code. 01267 default: 01268 { 01269 // The following quadrature rules are generated as 01270 // conical products. These tend to be non-optimal 01271 // (use too many points, cluster points in certain 01272 // regions of the domain) but they are quite easy to 01273 // automatically generate using a 1D Gauss rule on 01274 // [0,1] and two 1D Jacobi-Gauss rules on [0,1]. 01275 QConical conical_rule(2, _order); 01276 conical_rule.init(type_in, p); 01277 01278 // Swap points and weights with the about-to-be destroyed rule. 01279 _points.swap (conical_rule.get_points() ); 01280 _weights.swap(conical_rule.get_weights()); 01281 01282 return; 01283 } 01284 } 01285 } 01286 01287 01288 //--------------------------------------------- 01289 // Unsupported type 01290 default: 01291 { 01292 libMesh::err << "Element type not supported!:" << type_in << std::endl; 01293 libmesh_error(); 01294 } 01295 } 01296 01297 libmesh_error(); 01298 01299 return; 01300 01301 #endif 01302 }
| void libMesh::QGauss::init_3D | ( | const ElemType | _type = INVALID_ELEM, |
|
| unsigned int | p_level = 0 | |||
| ) | [private] |
Definition at line 28 of file quadrature_gauss_3D.C.
References libMesh::QBase::_points, libMesh::QBase::_weights, libMesh::QBase::allow_rules_with_negative_weights, libMeshEnums::CONSTANT, libMeshEnums::EDGE2, libMeshEnums::EIGHTH, libMesh::err, libMeshEnums::FIFTH, libMeshEnums::FIRST, libMeshEnums::FOURTH, libMesh::QBase::get_points(), libMesh::QBase::get_weights(), libMeshEnums::HEX20, libMeshEnums::HEX27, libMeshEnums::HEX8, libMesh::QBase::init(), keast_rule(), libMeshEnums::PRISM15, libMeshEnums::PRISM18, libMeshEnums::PRISM6, libMeshEnums::PYRAMID5, libMesh::Real, libMeshEnums::SECOND, libMeshEnums::SEVENTH, libMeshEnums::SIXTH, libMeshEnums::TET10, libMeshEnums::TET4, libMeshEnums::THIRD, and libMeshEnums::TRI3.
00030 { 00031 #if LIBMESH_DIM == 3 00032 00033 //----------------------------------------------------------------------- 00034 // 3D quadrature rules 00035 switch (type_in) 00036 { 00037 //--------------------------------------------- 00038 // Hex quadrature rules 00039 case HEX8: 00040 case HEX20: 00041 case HEX27: 00042 { 00043 // We compute the 3D quadrature rule as a tensor 00044 // product of the 1D quadrature rule. 00045 QGauss q1D(1,_order); 00046 q1D.init(EDGE2,p); 00047 00048 tensor_product_hex( q1D ); 00049 00050 return; 00051 } 00052 00053 00054 00055 //--------------------------------------------- 00056 // Tetrahedral quadrature rules 00057 case TET4: 00058 case TET10: 00059 { 00060 switch(_order + 2*p) 00061 { 00062 // Taken from pg. 222 of "The finite element method," vol. 1 00063 // ed. 5 by Zienkiewicz & Taylor 00064 case CONSTANT: 00065 case FIRST: 00066 { 00067 // Exact for linears 00068 _points.resize(1); 00069 _weights.resize(1); 00070 00071 00072 _points[0](0) = .25; 00073 _points[0](1) = .25; 00074 _points[0](2) = .25; 00075 00076 _weights[0] = .1666666666666666666666666666666666666666666667L; 00077 00078 return; 00079 } 00080 case SECOND: 00081 { 00082 // Exact for quadratics 00083 _points.resize(4); 00084 _weights.resize(4); 00085 00086 00087 const Real a = .585410196624969; 00088 const Real b = .138196601125011; 00089 00090 _points[0](0) = a; 00091 _points[0](1) = b; 00092 _points[0](2) = b; 00093 00094 _points[1](0) = b; 00095 _points[1](1) = a; 00096 _points[1](2) = b; 00097 00098 _points[2](0) = b; 00099 _points[2](1) = b; 00100 _points[2](2) = a; 00101 00102 _points[3](0) = b; 00103 _points[3](1) = b; 00104 _points[3](2) = b; 00105 00106 00107 00108 _weights[0] = .0416666666666666666666666666666666666666666667; 00109 _weights[1] = _weights[0]; 00110 _weights[2] = _weights[0]; 00111 _weights[3] = _weights[0]; 00112 00113 return; 00114 } 00115 00116 00117 00118 // Can be found in the class notes 00119 // http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps 00120 // by Flaherty. 00121 // 00122 // Caution: this rule has a negative weight and may be 00123 // unsuitable for some problems. 00124 // Exact for cubics. 00125 // 00126 // Note: Keast (see ref. elsewhere in this file) also gives 00127 // a third-order rule with positive weights, but it contains points 00128 // on the ref. elt. boundary, making it less suitable for FEM calculations. 00129 case THIRD: 00130 { 00131 if (allow_rules_with_negative_weights) 00132 { 00133 _points.resize(5); 00134 _weights.resize(5); 00135 00136 00137 _points[0](0) = .25; 00138 _points[0](1) = .25; 00139 _points[0](2) = .25; 00140 00141 _points[1](0) = .5; 00142 _points[1](1) = .16666666666666666666666666666666666666666667; 00143 _points[1](2) = .16666666666666666666666666666666666666666667; 00144 00145 _points[2](0) = .16666666666666666666666666666666666666666667; 00146 _points[2](1) = .5; 00147 _points[2](2) = .16666666666666666666666666666666666666666667; 00148 00149 _points[3](0) = .16666666666666666666666666666666666666666667; 00150 _points[3](1) = .16666666666666666666666666666666666666666667; 00151 _points[3](2) = .5; 00152 00153 _points[4](0) = .16666666666666666666666666666666666666666667; 00154 _points[4](1) = .16666666666666666666666666666666666666666667; 00155 _points[4](2) = .16666666666666666666666666666666666666666667; 00156 00157 00158 _weights[0] = -.133333333333333333333333333333333333333333333; 00159 _weights[1] = .075; 00160 _weights[2] = _weights[1]; 00161 _weights[3] = _weights[1]; 00162 _weights[4] = _weights[1]; 00163 00164 return; 00165 } // end if (allow_rules_with_negative_weights) 00166 else 00167 { 00168 // If a rule with postive weights is required, the 2x2x2 conical 00169 // product rule is third-order accurate and has less points than 00170 // the next-available positive-weight rule at FIFTH order. 00171 QConical conical_rule(3, _order); 00172 conical_rule.init(type_in, p); 00173 00174 // Swap points and weights with the about-to-be destroyed rule. 00175 _points.swap (conical_rule.get_points() ); 00176 _weights.swap(conical_rule.get_weights()); 00177 00178 return; 00179 } 00180 // Note: if !allow_rules_with_negative_weights, fall through to next case. 00181 } 00182 00183 00184 00185 // Originally a Keast rule, 00186 // Patrick Keast, 00187 // Moderate Degree Tetrahedral Quadrature Formulas, 00188 // Computer Methods in Applied Mechanics and Engineering, 00189 // Volume 55, Number 3, May 1986, pages 339-348. 00190 // 00191 // Can also be found the class notes 00192 // http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps 00193 // by Flaherty. 00194 // 00195 // Caution: this rule has a negative weight and may be 00196 // unsuitable for some problems. 00197 case FOURTH: 00198 { 00199 if (allow_rules_with_negative_weights) 00200 { 00201 _points.resize(11); 00202 _weights.resize(11); 00203 00204 // The raw data for the quadrature rule. 00205 const Real rule_data[3][4] = { 00206 {0.250000000000000000e+00, 0., 0., -0.131555555555555556e-01}, // 1 00207 {0.785714285714285714e+00, 0.714285714285714285e-01, 0., 0.762222222222222222e-02}, // 4 00208 {0.399403576166799219e+00, 0., 0.100596423833200785e+00, 0.248888888888888889e-01} // 6 00209 }; 00210 00211 00212 // Now call the keast routine to generate _points and _weights 00213 keast_rule(rule_data, 3); 00214 00215 return; 00216 } // end if (allow_rules_with_negative_weights) 00217 // Note: if !allow_rules_with_negative_weights, fall through to next case. 00218 } 00219 00220 00221 00222 00223 // Walkington's fifth-order 14-point rule from 00224 // "Quadrature on Simplices of Arbitrary Dimension" 00225 // 00226 // We originally had a Keast rule here, but this rule had 00227 // more points than an equivalent rule by Walkington and 00228 // also contained points on the boundary of the ref. elt, 00229 // making it less suitable for FEM calculations. 00230 case FIFTH: 00231 { 00232 _points.resize(14); 00233 _weights.resize(14); 00234 00235 // permutations of these points and suitably-modified versions of 00236 // these points are the quadrature point locations 00237 const Real a[3] = {0.31088591926330060980, // a1 from the paper 00238 0.092735250310891226402, // a2 from the paper 00239 0.045503704125649649492}; // a3 from the paper 00240 00241 // weights. a[] and wt[] are the only floating-point inputs required 00242 // for this rule. 00243 const Real wt[3] = {0.018781320953002641800, // w1 from the paper 00244 0.012248840519393658257, // w2 from the paper 00245 0.0070910034628469110730}; // w3 from the paper 00246 00247 // The first two sets of 4 points are formed in a similar manner 00248 for (unsigned int i=0; i<2; ++i) 00249 { 00250 // Where we will insert values into _points and _weights 00251 const unsigned int offset=4*i; 00252 00253 // Stuff points and weights values into their arrays 00254 const Real b = 1. - 3.*a[i]; 00255 00256 // Here are the permutations. Order of these is not important, 00257 // all have the same weight 00258 _points[offset + 0] = Point(a[i], a[i], a[i]); 00259 _points[offset + 1] = Point(a[i], b, a[i]); 00260 _points[offset + 2] = Point( b, a[i], a[i]); 00261 _points[offset + 3] = Point(a[i], a[i], b); 00262 00263 // These 4 points all have the same weights 00264 for (unsigned int j=0; j<4; ++j) 00265 _weights[offset + j] = wt[i]; 00266 } // end for 00267 00268 00269 { 00270 // The third set contains 6 points and is formed a little differently 00271 const unsigned int offset = 8; 00272 const Real b = 0.5*(1. - 2.*a[2]); 00273 00274 // Here are the permutations. Order of these is not important, 00275 // all have the same weight 00276 _points[offset + 0] = Point(b , b, a[2]); 00277 _points[offset + 1] = Point(b , a[2], a[2]); 00278 _points[offset + 2] = Point(a[2], a[2], b); 00279 _points[offset + 3] = Point(a[2], b, a[2]); 00280 _points[offset + 4] = Point( b, a[2], b); 00281 _points[offset + 5] = Point(a[2], b, b); 00282 00283 // These 6 points all have the same weights 00284 for (unsigned int j=0; j<6; ++j) 00285 _weights[offset + j] = wt[2]; 00286 } 00287 00288 00289 // Original rule by Keast, unsuitable because it has points on the 00290 // reference element boundary! 00291 // _points.resize(15); 00292 // _weights.resize(15); 00293 00294 // _points[0](0) = 0.25; 00295 // _points[0](1) = 0.25; 00296 // _points[0](2) = 0.25; 00297 00298 // { 00299 // const Real a = 0.; 00300 // const Real b = 0.333333333333333333333333333333333333333; 00301 00302 // _points[1](0) = a; 00303 // _points[1](1) = b; 00304 // _points[1](2) = b; 00305 00306 // _points[2](0) = b; 00307 // _points[2](1) = a; 00308 // _points[2](2) = b; 00309 00310 // _points[3](0) = b; 00311 // _points[3](1) = b; 00312 // _points[3](2) = a; 00313 00314 // _points[4](0) = b; 00315 // _points[4](1) = b; 00316 // _points[4](2) = b; 00317 // } 00318 // { 00319 // const Real a = 0.7272727272727272727272727272727272727272727272727272727; 00320 // const Real b = 0.0909090909090909090909090909090909090909090909090909091; 00321 00322 // _points[5](0) = a; 00323 // _points[5](1) = b; 00324 // _points[5](2) = b; 00325 00326 // _points[6](0) = b; 00327 // _points[6](1) = a; 00328 // _points[6](2) = b; 00329 00330 // _points[7](0) = b; 00331 // _points[7](1) = b; 00332 // _points[7](2) = a; 00333 00334 // _points[8](0) = b; 00335 // _points[8](1) = b; 00336 // _points[8](2) = b; 00337 // } 00338 // { 00339 // const Real a = 0.066550153573664; 00340 // const Real b = 0.433449846426336; 00341 00342 // _points[9](0) = b; 00343 // _points[9](1) = a; 00344 // _points[9](2) = a; 00345 00346 // _points[10](0) = a; 00347 // _points[10](1) = a; 00348 // _points[10](2) = b; 00349 00350 // _points[11](0) = a; 00351 // _points[11](1) = b; 00352 // _points[11](2) = b; 00353 00354 // _points[12](0) = b; 00355 // _points[12](1) = a; 00356 // _points[12](2) = b; 00357 00358 // _points[13](0) = b; 00359 // _points[13](1) = b; 00360 // _points[13](2) = a; 00361 00362 // _points[14](0) = a; 00363 // _points[14](1) = b; 00364 // _points[14](2) = a; 00365 // } 00366 00367 // _weights[0] = 0.030283678097089; 00368 // _weights[1] = 0.006026785714286; 00369 // _weights[2] = _weights[1]; 00370 // _weights[3] = _weights[1]; 00371 // _weights[4] = _weights[1]; 00372 // _weights[5] = 0.011645249086029; 00373 // _weights[6] = _weights[5]; 00374 // _weights[7] = _weights[5]; 00375 // _weights[8] = _weights[5]; 00376 // _weights[9] = 0.010949141561386; 00377 // _weights[10] = _weights[9]; 00378 // _weights[11] = _weights[9]; 00379 // _weights[12] = _weights[9]; 00380 // _weights[13] = _weights[9]; 00381 // _weights[14] = _weights[9]; 00382 00383 return; 00384 } 00385 00386 00387 00388 00389 // This rule is originally from Keast: 00390 // Patrick Keast, 00391 // Moderate Degree Tetrahedral Quadrature Formulas, 00392 // Computer Methods in Applied Mechanics and Engineering, 00393 // Volume 55, Number 3, May 1986, pages 339-348. 00394 // 00395 // It is accurate on 6th-degree polynomials and has 24 points 00396 // vs. 64 for the comparable conical product rule. 00397 // 00398 // Values copied 24th June 2008 from: 00399 // http://people.scs.fsu.edu/~burkardt/f_src/keast/keast.f90 00400 case SIXTH: 00401 { 00402 _points.resize (24); 00403 _weights.resize(24); 00404 00405 // The raw data for the quadrature rule. 00406 const Real rule_data[4][4] = { 00407 {0.356191386222544953e+00 , 0.214602871259151684e+00 , 0., 0.00665379170969464506e+00}, // 4 00408 {0.877978124396165982e+00 , 0.0406739585346113397e+00, 0., 0.00167953517588677620e+00}, // 4 00409 {0.0329863295731730594e+00, 0.322337890142275646e+00 , 0., 0.00922619692394239843e+00}, // 4 00410 {0.0636610018750175299e+00, 0.269672331458315867e+00 , 0.603005664791649076e+00, 0.00803571428571428248e+00} // 12 00411 }; 00412 00413 00414 // Now call the keast routine to generate _points and _weights 00415 keast_rule(rule_data, 4); 00416 00417 return; 00418 } 00419 00420 00421 // Keast's 31 point, 7th-order rule contains points on the reference 00422 // element boundary, so we've decided not to include it here. 00423 // 00424 // Keast's 8th-order rule has 45 points. and a negative 00425 // weight, so if you've explicitly disallowed such rules 00426 // you will fall through to the conical product rules 00427 // below. 00428 case SEVENTH: 00429 case EIGHTH: 00430 { 00431 if (allow_rules_with_negative_weights) 00432 { 00433 _points.resize (45); 00434 _weights.resize(45); 00435 00436 // The raw data for the quadrature rule. 00437 const Real rule_data[7][4] = { 00438 {0.250000000000000000e+00, 0., 0., -0.393270066412926145e-01}, // 1 00439 {0.617587190300082967e+00, 0.127470936566639015e+00, 0., 0.408131605934270525e-02}, // 4 00440 {0.903763508822103123e+00, 0.320788303926322960e-01, 0., 0.658086773304341943e-03}, // 4 00441 {0.497770956432810185e-01, 0., 0.450222904356718978e+00, 0.438425882512284693e-02}, // 6 00442 {0.183730447398549945e+00, 0., 0.316269552601450060e+00, 0.138300638425098166e-01}, // 6 00443 {0.231901089397150906e+00, 0.229177878448171174e-01, 0.513280033360881072e+00, 0.424043742468372453e-02}, // 12 00444 {0.379700484718286102e-01, 0.730313427807538396e+00, 0.193746475248804382e+00, 0.223873973961420164e-02} // 12 00445 }; 00446 00447 00448 // Now call the keast routine to generate _points and _weights 00449 keast_rule(rule_data, 7); 00450 00451 return; 00452 } // end if (allow_rules_with_negative_weights) 00453 // Note: if !allow_rules_with_negative_weights, fall through to next case. 00454 } 00455 00456 // Fall back on Grundmann-Moller or Conical Product rules at high orders. 00457 default: 00458 { 00459 if ((allow_rules_with_negative_weights) && (_order + 2*p < 34)) 00460 { 00461 // The Grundmann-Moller rules are defined to arbitrary order and 00462 // can have significantly fewer evaluation points than concial product 00463 // rules. If you allow rules with negative weights, the GM rules 00464 // will be more efficient for degree up to 33 (for degree 34 and 00465 // higher, CP is more efficient!) but may be more susceptible 00466 // to round-off error. Safest is to disallow rules with negative 00467 // weights, but this decision should be made on a case-by-case basis. 00468 QGrundmann_Moller gm_rule(3, _order); 00469 gm_rule.init(type_in, p); 00470 00471 // Swap points and weights with the about-to-be destroyed rule. 00472 _points.swap (gm_rule.get_points() ); 00473 _weights.swap(gm_rule.get_weights()); 00474 00475 return; 00476 } 00477 00478 else 00479 { 00480 // The following quadrature rules are generated as 00481 // conical products. These tend to be non-optimal 00482 // (use too many points, cluster points in certain 00483 // regions of the domain) but they are quite easy to 00484 // automatically generate using a 1D Gauss rule on 00485 // [0,1] and two 1D Jacobi-Gauss rules on [0,1]. 00486 QConical conical_rule(3, _order); 00487 conical_rule.init(type_in, p); 00488 00489 // Swap points and weights with the about-to-be destroyed rule. 00490 _points.swap (conical_rule.get_points() ); 00491 _weights.swap(conical_rule.get_weights()); 00492 00493 return; 00494 } 00495 } 00496 } 00497 } // end case TET4,TET10 00498 00499 00500 00501 //--------------------------------------------- 00502 // Prism quadrature rules 00503 case PRISM6: 00504 case PRISM15: 00505 case PRISM18: 00506 { 00507 // We compute the 3D quadrature rule as a tensor 00508 // product of the 1D quadrature rule and a 2D 00509 // triangle quadrature rule 00510 00511 QGauss q1D(1,_order); 00512 QGauss q2D(2,_order); 00513 00514 // Initialize 00515 q1D.init(EDGE2,p); 00516 q2D.init(TRI3,p); 00517 00518 tensor_product_prism(q1D, q2D); 00519 00520 return; 00521 } 00522 00523 00524 00525 //--------------------------------------------- 00526 // Pyramid 00527 case PYRAMID5: 00528 { 00529 // We compute the Pyramid rule as a conical product of a 00530 // Jacobi rule with alpha==2 on the interval [0,1] two 1D 00531 // Gauss rules with suitably modified points. The idea comes 00532 // from: Stroud, A.H. "Approximate Calculation of Multiple 00533 // Integrals." 00534 QConical conical_rule(3, _order); 00535 conical_rule.init(type_in, p); 00536 00537 // Swap points and weights with the about-to-be destroyed rule. 00538 _points.swap (conical_rule.get_points() ); 00539 _weights.swap(conical_rule.get_weights()); 00540 00541 return; 00542 00543 } 00544 00545 00546 00547 //--------------------------------------------- 00548 // Unsupported type 00549 default: 00550 { 00551 libMesh::err << "ERROR: Unsupported type: " << type_in << std::endl; 00552 libmesh_error(); 00553 } 00554 } 00555 00556 libmesh_error(); 00557 00558 return; 00559 00560 #endif 00561 }
| void libMesh::QGauss::keast_rule | ( | const Real | rule_data[][4], | |
| const unsigned int | n_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 33 of file quadrature_gauss.C.
References libMesh::QBase::_points, libMesh::QBase::_weights, libMesh::err, and libMesh::Real.
Referenced by init_3D().
00035 { 00036 // Like the Dunavant rule, the input data should have 4 columns. These columns 00037 // have the following format and implied permutations (w=weight). 00038 // {a, 0, 0, w} = 1-permutation (a,a,a) 00039 // {a, b, 0, w} = 4-permutation (a,b,b), (b,a,b), (b,b,a), (b,b,b) 00040 // {a, 0, b, w} = 6-permutation (a,a,b), (a,b,b), (b,b,a), (b,a,b), (b,a,a), (a,b,a) 00041 // {a, b, c, w} = 12-permutation (a,a,b), (a,a,c), (b,a,a), (c,a,a), (a,b,a), (a,c,a) 00042 // (a,b,c), (a,c,b), (b,a,c), (b,c,a), (c,a,b), (c,b,a) 00043 00044 // Always insert into the points & weights vector relative to the offset 00045 unsigned int offset=0; 00046 00047 00048 for (unsigned int p=0; p<n_pts; ++p) 00049 { 00050 00051 // There must always be a non-zero entry to start the row 00052 libmesh_assert_not_equal_to (rule_data[p][0], static_cast<Real>(0.0)); 00053 00054 // A zero weight may imply you did not set up the raw data correctly 00055 libmesh_assert_not_equal_to (rule_data[p][3], static_cast<Real>(0.0)); 00056 00057 // What kind of point is this? 00058 // One non-zero entry in first 3 cols ? 1-perm (centroid) point = 1 00059 // Two non-zero entries in first 3 cols ? 3-perm point = 3 00060 // Three non-zero entries ? 6-perm point = 6 00061 unsigned int pointtype=1; 00062 00063 if (rule_data[p][1] != static_cast<Real>(0.0)) 00064 { 00065 if (rule_data[p][2] != static_cast<Real>(0.0)) 00066 pointtype = 12; 00067 else 00068 pointtype = 4; 00069 } 00070 else 00071 { 00072 // The second entry is zero. What about the third? 00073 if (rule_data[p][2] != static_cast<Real>(0.0)) 00074 pointtype = 6; 00075 } 00076 00077 00078 switch (pointtype) 00079 { 00080 case 1: 00081 { 00082 // Be sure we have enough space to insert this point 00083 libmesh_assert_less (offset + 0, _points.size()); 00084 00085 const Real a = rule_data[p][0]; 00086 00087 // The point has only a single permutation (the centroid!) 00088 _points[offset + 0] = Point(a,a,a); 00089 00090 // The weight is always the last entry in the row. 00091 _weights[offset + 0] = rule_data[p][3]; 00092 00093 offset += pointtype; 00094 break; 00095 } 00096 00097 case 4: 00098 { 00099 // Be sure we have enough space to insert these points 00100 libmesh_assert_less (offset + 3, _points.size()); 00101 00102 const Real a = rule_data[p][0]; 00103 const Real b = rule_data[p][1]; 00104 const Real wt = rule_data[p][3]; 00105 00106 // Here it's understood the second entry is to be used twice, and 00107 // thus there are three possible permutations. 00108 _points[offset + 0] = Point(a,b,b); 00109 _points[offset + 1] = Point(b,a,b); 00110 _points[offset + 2] = Point(b,b,a); 00111 _points[offset + 3] = Point(b,b,b); 00112 00113 for (unsigned int j=0; j<pointtype; ++j) 00114 _weights[offset + j] = wt; 00115 00116 offset += pointtype; 00117 break; 00118 } 00119 00120 case 6: 00121 { 00122 // Be sure we have enough space to insert these points 00123 libmesh_assert_less (offset + 5, _points.size()); 00124 00125 const Real a = rule_data[p][0]; 00126 const Real b = rule_data[p][2]; 00127 const Real wt = rule_data[p][3]; 00128 00129 // Three individual entries with six permutations. 00130 _points[offset + 0] = Point(a,a,b); 00131 _points[offset + 1] = Point(a,b,b); 00132 _points[offset + 2] = Point(b,b,a); 00133 _points[offset + 3] = Point(b,a,b); 00134 _points[offset + 4] = Point(b,a,a); 00135 _points[offset + 5] = Point(a,b,a); 00136 00137 for (unsigned int j=0; j<pointtype; ++j) 00138 _weights[offset + j] = wt; 00139 00140 offset += pointtype; 00141 break; 00142 } 00143 00144 00145 case 12: 00146 { 00147 // Be sure we have enough space to insert these points 00148 libmesh_assert_less (offset + 11, _points.size()); 00149 00150 const Real a = rule_data[p][0]; 00151 const Real b = rule_data[p][1]; 00152 const Real c = rule_data[p][2]; 00153 const Real wt = rule_data[p][3]; 00154 00155 // Three individual entries with six permutations. 00156 _points[offset + 0] = Point(a,a,b); _points[offset + 6] = Point(a,b,c); 00157 _points[offset + 1] = Point(a,a,c); _points[offset + 7] = Point(a,c,b); 00158 _points[offset + 2] = Point(b,a,a); _points[offset + 8] = Point(b,a,c); 00159 _points[offset + 3] = Point(c,a,a); _points[offset + 9] = Point(b,c,a); 00160 _points[offset + 4] = Point(a,b,a); _points[offset + 10] = Point(c,a,b); 00161 _points[offset + 5] = Point(a,c,a); _points[offset + 11] = Point(c,b,a); 00162 00163 for (unsigned int j=0; j<pointtype; ++j) 00164 _weights[offset + j] = wt; 00165 00166 offset += pointtype; 00167 break; 00168 } 00169 00170 default: 00171 { 00172 libMesh::err << "Don't know what to do with this many permutation points!" << std::endl; 00173 libmesh_error(); 00174 } 00175 } 00176 00177 } 00178 00179 }
| 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; }
| unsigned int libMesh::QBase::n_points | ( | ) | const [inline, inherited] |
- Returns:
- the number of points associated with the quadrature rule.
Definition at line 116 of file quadrature.h.
References libMesh::QBase::_points.
Referenced by libMesh::FEGenericBase< OutputType >::coarsened_dof_values(), libMesh::InfFE< Dim, T_radial, T_map >::combine_base_radial(), libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), libMesh::InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), libMesh::InfFE< Dim, T_radial, T_map >::init_shape_functions(), libMesh::FE< Dim, T >::n_quadrature_points(), libMesh::ProjectFEMSolution::operator()(), and libMesh::QBase::print_info().
| 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
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().
| 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().
| QuadratureType libMesh::QGauss::type | ( | ) | const [inline, virtual] |
- Returns:
QGAUSS
Implements libMesh::QBase.
Definition at line 61 of file quadrature_gauss.h.
References libMeshEnums::QGAUSS.
00061 { return QGAUSS; }
| Real libMesh::QBase::w | ( | const unsigned int | i | ) | const [inline, inherited] |
- Returns:
- the 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().
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 libMesh::ReferenceCounter::_counts [static, protected, inherited] |
Actually holds the data.
Definition at line 118 of file reference_counter.h.
Referenced by libMesh::ReferenceCounter::get_info(), libMesh::ReferenceCounter::increment_constructor_count(), and libMesh::ReferenceCounter::increment_destructor_count().
bool libMesh::ReferenceCounter::_enable_print_counter = true [static, protected, inherited] |
Flag to control whether reference count information is printed when print_info is called.
Definition at line 137 of file reference_counter.h.
Referenced by libMesh::ReferenceCounter::disable_print_counter_info(), libMesh::ReferenceCounter::enable_print_counter_info(), and libMesh::ReferenceCounter::print_info().
Threads::spin_mutex libMesh::ReferenceCounter::_mutex [static, protected, inherited] |
Mutual exclusion object to enable thread-safe reference counting.
Definition at line 131 of file reference_counter.h.
Threads::atomic< unsigned int > libMesh::ReferenceCounter::_n_objects [static, protected, inherited] |
The number of objects. Print the reference count information when the number returns to 0.
Definition at line 126 of file reference_counter.h.
Referenced by libMesh::ReferenceCounter::n_objects(), libMesh::ReferenceCounter::ReferenceCounter(), and libMesh::ReferenceCounter::~ReferenceCounter().
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] |
Definition at line 332 of file quadrature.h.
Referenced by libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), dunavant_rule(), dunavant_rule2(), libMesh::QBase::get_points(), libMesh::QGrundmann_Moller::gm_rule(), libMesh::QBase::init_0D(), libMesh::QTrap::init_1D(), libMesh::QSimpson::init_1D(), libMesh::QJacobi::init_1D(), libMesh::QGrid::init_1D(), init_1D(), libMesh::QClough::init_1D(), libMesh::QTrap::init_2D(), libMesh::QSimpson::init_2D(), libMesh::QMonomial::init_2D(), libMesh::QGrid::init_2D(), init_2D(), libMesh::QClough::init_2D(), libMesh::QTrap::init_3D(), libMesh::QSimpson::init_3D(), libMesh::QMonomial::init_3D(), libMesh::QGrid::init_3D(), init_3D(), keast_rule(), libMesh::QMonomial::kim_rule(), libMesh::QBase::n_points(), libMesh::QBase::print_info(), libMesh::QBase::qp(), libMesh::QBase::scale(), libMesh::QMonomial::stroud_rule(), and libMesh::QMonomial::wissmann_rule().
std::vector<Real> libMesh::QBase::_weights [protected, inherited] |
The value of the quadrature weights.
Definition at line 337 of file quadrature.h.
Referenced by libMesh::QConical::conical_product_pyramid(), libMesh::QConical::conical_product_tet(), libMesh::QConical::conical_product_tri(), dunavant_rule(), dunavant_rule2(), libMesh::QBase::get_weights(), libMesh::QGrundmann_Moller::gm_rule(), libMesh::QBase::init_0D(), libMesh::QTrap::init_1D(), libMesh::QSimpson::init_1D(), libMesh::QJacobi::init_1D(), libMesh::QGrid::init_1D(), init_1D(), libMesh::QClough::init_1D(), libMesh::QTrap::init_2D(), libMesh::QSimpson::init_2D(), libMesh::QMonomial::init_2D(), libMesh::QGrid::init_2D(), init_2D(), libMesh::QClough::init_2D(), libMesh::QTrap::init_3D(), libMesh::QSimpson::init_3D(), libMesh::QMonomial::init_3D(), libMesh::QGrid::init_3D(), init_3D(), keast_rule(), libMesh::QMonomial::kim_rule(), libMesh::QBase::print_info(), libMesh::QBase::scale(), libMesh::QMonomial::stroud_rule(), libMesh::QBase::w(), and libMesh::QMonomial::wissmann_rule().
bool libMesh::QBase::allow_rules_with_negative_weights [inherited] |
Flag (default true) controlling the use of quadrature rules with negative weights. Set this to false to ONLY use (potentially) safer but more expensive rules with all positive weights.
Negative weights typically appear in Gaussian quadrature rules over three-dimensional elements. Rules with negative weights can be unsuitable for some problems. For example, it is possible for a rule with negative weights to obtain a negative result when integrating a positive function.
A particular example: if rules with negative weights are not allowed, a request for TET,THIRD (5 points) will return the TET,FIFTH (14 points) rule instead, nearly tripling the computational effort required!
Definition at line 215 of file quadrature.h.
Referenced by libMesh::QMonomial::init_3D(), libMesh::QGrundmann_Moller::init_3D(), and 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: