dune-geometry
2.2.0
|
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set et ts=4 sw=2 sts=2: 00003 00004 #ifndef DUNE_GEOMETRY_QUADRATURERULES_HH 00005 #define DUNE_GEOMETRY_QUADRATURERULES_HH 00006 00007 #include <iostream> 00008 #include <limits> 00009 #include <vector> 00010 #include <map> 00011 00012 #include <dune/common/fvector.hh> 00013 #include <dune/common/exceptions.hh> 00014 #include <dune/common/stdstreams.hh> 00015 00016 #include <dune/geometry/type.hh> 00017 00023 namespace Dune { 00024 00029 class QuadratureOrderOutOfRange : public NotImplemented {}; 00030 00034 template<typename ct, int dim> 00035 class QuadraturePoint { 00036 public: 00037 // compile time parameters 00038 enum { d=dim }; 00039 typedef ct CoordType; 00040 00041 static const unsigned int dimension = d; 00042 typedef ct Field; 00043 typedef Dune::FieldVector<ct,dim> Vector; 00044 00046 QuadraturePoint (const Vector& x, ct w) : local(x) 00047 { 00048 wght = w; 00049 } 00050 00052 const Vector& position () const 00053 { 00054 return local; 00055 } 00056 00058 const ct &weight () const 00059 { 00060 return wght; 00061 } 00062 00063 protected: 00064 FieldVector<ct, dim> local; 00065 ct wght; 00066 }; 00067 00071 namespace QuadratureType { 00072 enum Enum { 00073 Gauss = 0, 00074 00075 Jacobian_1_0 = 1, 00076 Jacobian_2_0 = 2, 00077 00078 Simpson = 3, 00079 Trap = 4, 00080 Grid = 5, 00081 00082 Clough = 21, 00083 00084 Invalid_Rule = 127 00085 }; 00086 } 00087 00091 template<typename ct, int dim> 00092 class QuadratureRule : public std::vector<QuadraturePoint<ct,dim> > 00093 { 00094 public: 00095 00097 QuadratureRule() : delivered_order(-1) {} 00098 00100 QuadratureRule(GeometryType t) : geometry_type(t), delivered_order(-1) {} 00101 00103 QuadratureRule(GeometryType t, int order) : geometry_type(t), delivered_order(order) {} 00104 00106 enum { d=dim }; 00107 00109 typedef ct CoordType; 00110 00112 virtual int order () const { return delivered_order; } 00113 00115 virtual GeometryType type () const { return geometry_type; } 00116 virtual ~QuadratureRule(){} 00117 00120 typedef typename std::vector<QuadraturePoint<ct,dim> >::const_iterator iterator; 00121 00122 protected: 00123 GeometryType geometry_type; 00124 int delivered_order; 00125 00126 void tensor_product (const QuadratureRule<ct,1> & q1d) 00127 { 00128 // fill in all the gauss points 00129 int m = q1d.size(); 00130 int n = power(m,dim); 00131 for (int i=0; i<n; i++) 00132 { 00133 // compute multi index for Gauss point 00134 int x[dim]; 00135 int z = i; 00136 for (int k=0; k<dim; ++k) 00137 { 00138 x[k] = z%m; 00139 z = z/m; 00140 } 00141 00142 // compute coordinates and weight 00143 double weight = 1.0; 00144 FieldVector<ct, dim> local; 00145 for (int k=0; k<dim; k++) 00146 { 00147 local[k] = q1d[x[k]].position()[0]; 00148 weight *= q1d[x[k]].weight(); 00149 } 00150 00151 // put in container 00152 this->push_back(QuadraturePoint<ct,dim>(local,weight)); 00153 } 00154 } 00155 00156 int power (int y, int d) 00157 { 00158 int m=1; 00159 for (int i=0; i<d; i++) m *= y; 00160 return m; 00161 } 00162 }; 00163 00164 // Forward declaration of the factory class, 00165 // needed internally by the QuadratureRules container class. 00166 template<typename ctype, int dim> class QuadratureRuleFactory; 00167 00171 template<typename ctype, int dim> 00172 class QuadratureRules { 00173 00175 typedef std::pair<GeometryType,int> QuadratureRuleKey; 00176 00178 typedef Dune::QuadratureRule<ctype, dim> QuadratureRule; 00179 00181 const QuadratureRule& _rule(const GeometryType& t, int p, QuadratureType::Enum qt=QuadratureType::Gauss) 00182 { 00183 static std::map<QuadratureRuleKey, QuadratureRule> _quadratureMap; 00184 QuadratureRuleKey key(t,p); 00185 if (_quadratureMap.find(key) == _quadratureMap.end()) { 00186 /* 00187 The rule must be acquired before we can store it. 00188 If we write this in one command, an invalid rule 00189 would get stored in case of an exception. 00190 */ 00191 QuadratureRule rule = 00192 QuadratureRuleFactory<ctype,dim>::rule(t,p,qt); 00193 _quadratureMap[key] = rule; 00194 } 00195 return _quadratureMap[key]; 00196 } 00198 static QuadratureRules& instance() 00199 { 00200 static QuadratureRules instance; 00201 return instance; 00202 } 00204 QuadratureRules () {}; 00205 public: 00207 static const QuadratureRule& rule(const GeometryType& t, int p, QuadratureType::Enum qt=QuadratureType::Gauss) 00208 { 00209 return instance()._rule(t,p,qt); 00210 } 00212 static const QuadratureRule& rule(const GeometryType::BasicType t, int p, QuadratureType::Enum qt=QuadratureType::Gauss) 00213 { 00214 GeometryType gt(t,dim); 00215 return instance()._rule(gt,p,qt); 00216 } 00217 }; 00218 00219 /***********************************************************/ 00220 00224 template<typename ct, int dim> 00225 class CubeQuadratureRule : 00226 public QuadratureRule<ct,dim> 00227 { 00228 public: 00230 enum { d=dim }; 00231 00233 enum { highest_order=CubeQuadratureRule<ct,1>::highest_order }; 00234 00236 typedef ct CoordType; 00237 00239 typedef CubeQuadratureRule value_type; 00240 00241 ~CubeQuadratureRule(){} 00242 private: 00243 friend class QuadratureRuleFactory<ct,dim>; 00245 CubeQuadratureRule (int p) : QuadratureRule<ct,dim>(GeometryType(GeometryType::cube, d)) 00246 { 00247 QuadratureRule<ct,1> q1D = QuadratureRules<ct,1>::rule(GeometryType::cube, p); 00248 this->tensor_product( q1D ); 00249 this->delivered_order = q1D.order(); 00250 } 00251 00252 }; 00253 00256 template<typename ct> 00257 class CubeQuadratureRule<ct,0> : 00258 public QuadratureRule<ct,0> 00259 { 00260 public: 00261 // compile time parameters 00262 enum { d=0 }; 00263 enum { dim=0 }; 00264 enum { highest_order=61 }; 00265 typedef ct CoordType; 00266 typedef CubeQuadratureRule value_type; 00267 00268 ~CubeQuadratureRule(){} 00269 private: 00270 friend class QuadratureRuleFactory<ct,dim>; 00271 CubeQuadratureRule (int p): 00272 QuadratureRule<ct,0>(GeometryType(GeometryType::cube, 0)) 00273 { 00274 FieldVector<ct, dim> point(0.0); 00275 00276 if (p > highest_order) 00277 DUNE_THROW(QuadratureOrderOutOfRange, "Quadrature rule " << p << " not supported!"); 00278 00279 this->delivered_order = highest_order; 00280 this->push_back(QuadraturePoint<ct,dim>(point, 1.0)); 00281 } 00282 }; 00283 00285 template<typename ct, 00286 bool fundamental = std::numeric_limits<ct>::is_specialized> 00287 struct CubeQuadratureInitHelper; 00288 template<typename ct> 00289 struct CubeQuadratureInitHelper<ct, true> { 00290 static void init(int p, 00291 std::vector< FieldVector<ct, 1> > & _points, 00292 std::vector< ct > & _weight, 00293 int & delivered_order); 00294 }; 00295 template<typename ct> 00296 struct CubeQuadratureInitHelper<ct, false> { 00297 static void init(int p, 00298 std::vector< FieldVector<ct, 1> > & _points, 00299 std::vector< ct > & _weight, 00300 int & delivered_order); 00301 }; 00302 00305 template<typename ct> 00306 class CubeQuadratureRule<ct,1> : 00307 public QuadratureRule<ct,1> 00308 { 00309 public: 00310 // compile time parameters 00311 enum { d=1 }; 00312 enum { dim=1 }; 00313 enum { highest_order=61 }; 00314 typedef ct CoordType; 00315 typedef CubeQuadratureRule value_type; 00316 00317 ~CubeQuadratureRule(){} 00318 private: 00319 friend class QuadratureRuleFactory<ct,dim>; 00320 CubeQuadratureRule (int p) 00321 : QuadratureRule<ct,1>(GeometryType(GeometryType::cube, 1)) 00322 { 00324 std::vector< FieldVector<ct, dim> > _points; 00325 std::vector< ct > _weight; 00326 00327 CubeQuadratureInitHelper<ct>::init 00328 (p, _points, _weight, this->delivered_order); 00329 00330 assert(_points.size() == _weight.size()); 00331 for (size_t i = 0; i < _points.size(); i++) 00332 this->push_back(QuadraturePoint<ct,dim>(_points[i], _weight[i])); 00333 } 00334 }; 00335 00336 extern template CubeQuadratureRule<float, 1>::CubeQuadratureRule(int); 00337 extern template CubeQuadratureRule<double, 1>::CubeQuadratureRule(int); 00338 00339 } // namespace Dune 00340 00341 #define DUNE_INCLUDING_IMPLEMENTATION 00342 #include "quadraturerules/cube_imp.hh" 00343 00344 namespace Dune { 00345 00349 template<typename ct, int dim> 00350 class Jacobi1QuadratureRule; 00351 00353 template<typename ct, 00354 bool fundamental = std::numeric_limits<ct>::is_specialized> 00355 struct Jacobi1QuadratureInitHelper; 00356 template<typename ct> 00357 struct Jacobi1QuadratureInitHelper<ct, true> { 00358 static void init(int p, 00359 std::vector< FieldVector<ct, 1> > & _points, 00360 std::vector< ct > & _weight, 00361 int & delivered_order); 00362 }; 00363 template<typename ct> 00364 struct Jacobi1QuadratureInitHelper<ct, false> { 00365 static void init(int p, 00366 std::vector< FieldVector<ct, 1> > & _points, 00367 std::vector< ct > & _weight, 00368 int & delivered_order); 00369 }; 00370 00374 template<typename ct> 00375 class Jacobi1QuadratureRule<ct,1> : 00376 public QuadratureRule<ct,1> 00377 { 00378 public: 00380 enum { d=1 }; 00383 enum { dim=1 }; 00384 00386 enum { highest_order=61 }; 00387 00389 typedef ct CoordType; 00390 00392 typedef Jacobi1QuadratureRule value_type; 00393 00394 ~Jacobi1QuadratureRule(){} 00395 private: 00396 friend class QuadratureRuleFactory<ct,dim>; 00397 Jacobi1QuadratureRule (int p) 00398 : QuadratureRule<ct,1>(GeometryType(GeometryType::cube, 1)) 00399 { 00401 std::vector< FieldVector<ct, dim> > _points; 00402 std::vector< ct > _weight; 00403 00404 int delivered_order; 00405 00406 Jacobi1QuadratureInitHelper<ct>::init 00407 (p, _points, _weight, delivered_order); 00408 this->delivered_order = delivered_order; 00409 assert(_points.size() == _weight.size()); 00410 for (size_t i = 0; i < _points.size(); i++) 00411 this->push_back(QuadraturePoint<ct,dim>(_points[i], _weight[i])); 00412 } 00413 }; 00414 00415 #ifndef DOXYGEN 00416 extern template Jacobi1QuadratureRule<float, 1>::Jacobi1QuadratureRule(int); 00417 extern template Jacobi1QuadratureRule<double, 1>::Jacobi1QuadratureRule(int); 00418 #endif // !DOXYGEN 00419 00420 } // namespace Dune 00421 00422 #define DUNE_INCLUDING_IMPLEMENTATION 00423 #include "quadraturerules/jacobi_1_0_imp.hh" 00424 00425 namespace Dune { 00426 00430 template<typename ct, int dim> 00431 class Jacobi2QuadratureRule; 00432 00434 template<typename ct, 00435 bool fundamental = std::numeric_limits<ct>::is_specialized> 00436 struct Jacobi2QuadratureInitHelper; 00437 template<typename ct> 00438 struct Jacobi2QuadratureInitHelper<ct, true> { 00439 static void init(int p, 00440 std::vector< FieldVector<ct, 1> > & _points, 00441 std::vector< ct > & _weight, 00442 int & delivered_order); 00443 }; 00444 template<typename ct> 00445 struct Jacobi2QuadratureInitHelper<ct, false> { 00446 static void init(int p, 00447 std::vector< FieldVector<ct, 1> > & _points, 00448 std::vector< ct > & _weight, 00449 int & delivered_order); 00450 }; 00451 00455 template<typename ct> 00456 class Jacobi2QuadratureRule<ct,1> : 00457 public QuadratureRule<ct,1> 00458 { 00459 public: 00461 enum { d=1 }; 00462 00465 enum { dim=1 }; 00466 00468 enum { highest_order=61 }; 00469 00471 typedef ct CoordType; 00472 00474 typedef Jacobi2QuadratureRule value_type; 00475 00476 ~Jacobi2QuadratureRule(){} 00477 private: 00478 friend class QuadratureRuleFactory<ct,dim>; 00479 Jacobi2QuadratureRule (int p) 00480 : QuadratureRule<ct,1>(GeometryType(GeometryType::cube, 1)) 00481 { 00483 std::vector< FieldVector<ct, dim> > _points; 00484 std::vector< ct > _weight; 00485 00486 int delivered_order; 00487 00488 Jacobi2QuadratureInitHelper<ct>::init 00489 (p, _points, _weight, delivered_order); 00490 00491 this->delivered_order = delivered_order; 00492 assert(_points.size() == _weight.size()); 00493 for (size_t i = 0; i < _points.size(); i++) 00494 this->push_back(QuadraturePoint<ct,dim>(_points[i], _weight[i])); 00495 } 00496 }; 00497 00498 #ifndef DOXYGEN 00499 extern template Jacobi2QuadratureRule<float, 1>::Jacobi2QuadratureRule(int); 00500 extern template Jacobi2QuadratureRule<double, 1>::Jacobi2QuadratureRule(int); 00501 #endif // !DOXYGEN 00502 00503 } // namespace Dune 00504 00505 #define DUNE_INCLUDING_IMPLEMENTATION 00506 #include "quadraturerules/jacobi_2_0_imp.hh" 00507 00508 namespace Dune { 00509 00510 /************************************************ 00511 * Quadraturerule for Simplices/Triangle 00512 *************************************************/ 00513 00517 template<typename ct, int dim> 00518 class SimplexQuadratureRule; 00519 00523 template<typename ct> 00524 class SimplexQuadratureRule<ct,2> : public QuadratureRule<ct,2> 00525 { 00526 public: 00527 00529 enum{d=2}; 00530 00532 enum { highest_order=CubeQuadratureRule<ct,1>::highest_order -1 }; 00533 00535 typedef ct CoordType; 00536 00538 typedef SimplexQuadratureRule value_type; 00539 ~SimplexQuadratureRule(){} 00540 private: 00541 friend class QuadratureRuleFactory<ct,d>; 00542 SimplexQuadratureRule (int p); 00543 }; 00544 00548 template<typename ct> 00549 class SimplexQuadratureRule<ct,3> : public QuadratureRule<ct,3> 00550 { 00551 public: 00552 00554 enum{d=3}; 00555 00557 enum { highest_order=CubeQuadratureRule<ct,1>::highest_order -2 }; 00558 00560 typedef ct CoordType; 00561 00563 typedef SimplexQuadratureRule<ct,3> value_type; 00564 ~SimplexQuadratureRule(){} 00565 private: 00566 friend class QuadratureRuleFactory<ct,d>; 00567 SimplexQuadratureRule (int p); 00568 }; 00569 00570 /*********************************** 00571 * quadrature for Prism 00572 **********************************/ 00573 00575 template<int dim> 00576 class PrismQuadraturePoints; 00577 00579 template<> 00580 class PrismQuadraturePoints<3> 00581 { 00582 public: 00583 enum { MAXP=6}; 00584 enum { highest_order=2 }; 00585 00587 PrismQuadraturePoints () 00588 { 00589 int m = 0; 00590 O[m] = 0; 00591 00592 // polynom degree 0 ??? 00593 m = 6; 00594 G[m][0][0] = 0.0; 00595 G[m][0][1] = 0.0; 00596 G[m][0][2] = 0.0; 00597 00598 G[m][1][0] = 1.0; 00599 G[m][1][1] = 0.0; 00600 G[m][1][2] = 0.0; 00601 00602 G[m][2][0] = 0.0; 00603 G[m][2][1] = 1.0; 00604 G[m][2][2] = 0.0; 00605 00606 G[m][3][0] = 0.0; 00607 G[m][3][1] = 0.0; 00608 G[m][3][2] = 1.0; 00609 00610 G[m][4][0] = 1.0; 00611 G[m][4][1] = 0.0; 00612 G[m][4][2] = 1.0; 00613 00614 G[m][5][0] = 0.0; 00615 G[m][5][1] = 0.1; 00616 G[m][5][2] = 1.0; 00617 00618 W[m][0] = 0.16666666666666666 / 2.0; 00619 W[m][1] = 0.16666666666666666 / 2.0; 00620 W[m][2] = 0.16666666666666666 / 2.0; 00621 W[m][3] = 0.16666666666666666 / 2.0; 00622 W[m][4] = 0.16666666666666666 / 2.0; 00623 W[m][5] = 0.16666666666666666 / 2.0; 00624 00625 O[m] = 0;// verify ???????? 00626 00627 00628 // polynom degree 2 ??? 00629 m = 6; 00630 G[m][0][0] =0.66666666666666666 ; 00631 G[m][0][1] =0.16666666666666666 ; 00632 G[m][0][2] =0.211324865405187 ; 00633 00634 G[m][1][0] = 0.16666666666666666; 00635 G[m][1][1] =0.66666666666666666 ; 00636 G[m][1][2] = 0.211324865405187; 00637 00638 G[m][2][0] = 0.16666666666666666; 00639 G[m][2][1] = 0.16666666666666666; 00640 G[m][2][2] = 0.211324865405187; 00641 00642 G[m][3][0] = 0.66666666666666666; 00643 G[m][3][1] = 0.16666666666666666; 00644 G[m][3][2] = 0.788675134594813; 00645 00646 G[m][4][0] = 0.16666666666666666; 00647 G[m][4][1] = 0.66666666666666666; 00648 G[m][4][2] = 0.788675134594813; 00649 00650 G[m][5][0] = 0.16666666666666666; 00651 G[m][5][1] = 0.16666666666666666; 00652 G[m][5][2] = 0.788675134594813; 00653 00654 W[m][0] = 0.16666666666666666 / 2.0; 00655 W[m][1] = 0.16666666666666666 / 2.0; 00656 W[m][2] = 0.16666666666666666 / 2.0; 00657 W[m][3] = 0.16666666666666666 / 2.0; 00658 W[m][4] = 0.16666666666666666 / 2.0; 00659 W[m][5] = 0.16666666666666666 / 2.0; 00660 00661 O[m] = 2;// verify ???????? 00662 00663 } 00664 00666 FieldVector<double, 3> point(int m, int i) 00667 { 00668 return G[m][i]; 00669 } 00670 00672 double weight (int m, int i) 00673 { 00674 return W[m][i]; 00675 } 00676 00678 int order (int m) 00679 { 00680 return O[m]; 00681 } 00682 00683 private: 00684 FieldVector<double, 3> G[MAXP+1][MAXP];//positions 00685 00686 double W[MAXP+1][MAXP]; // weights associated with points 00687 int O[MAXP+1]; // order of the rule 00688 }; 00689 00690 00694 template<int dim> 00695 struct PrismQuadraturePointsSingleton { 00696 static PrismQuadraturePoints<3> prqp; 00697 }; 00698 00702 template<> 00703 struct PrismQuadraturePointsSingleton<3> { 00704 static PrismQuadraturePoints<3> prqp; 00705 }; 00706 00710 template<typename ct, int dim> 00711 class PrismQuadratureRule; 00712 00716 template<typename ct> 00717 class PrismQuadratureRule<ct,3> : public QuadratureRule<ct,3> 00718 { 00719 public: 00720 00722 enum{ d=3 }; 00723 00725 enum{ 00726 /* min(Line::order, Triangle::order) */ 00727 highest_order = 00728 (int)CubeQuadratureRule<ct,1>::highest_order 00729 < (int)SimplexQuadratureRule<ct,2>::highest_order 00730 ? (int)CubeQuadratureRule<ct,1>::highest_order 00731 : (int)SimplexQuadratureRule<ct,2>::highest_order 00732 }; 00733 00735 typedef ct CoordType; 00736 00738 typedef PrismQuadratureRule<ct,3> value_type; 00739 00740 ~PrismQuadratureRule(){} 00741 private: 00742 friend class QuadratureRuleFactory<ct,d>; 00743 PrismQuadratureRule(int p) : QuadratureRule<ct,3>(GeometryType(GeometryType::prism, d)) 00744 { 00745 if (p>highest_order) 00746 DUNE_THROW(QuadratureOrderOutOfRange, 00747 "QuadratureRule for order " << p << " and GeometryType " 00748 << this->type() << " not available"); 00749 00750 if (p<=2) { 00751 int m=6; 00752 this->delivered_order = PrismQuadraturePointsSingleton<3>::prqp.order(m); 00753 for(int i=0;i<m;++i) 00754 { 00755 FieldVector<ct,3> local; 00756 for (int k=0; k<d; k++) 00757 local[k] = PrismQuadraturePointsSingleton<3>::prqp.point(m,i)[k]; 00758 double weight = 00759 PrismQuadraturePointsSingleton<3>::prqp.weight(m,i); 00760 // put in container 00761 this->push_back(QuadraturePoint<ct,d>(local,weight)); 00762 } 00763 } 00764 else { 00765 const QuadratureRule<ct,2> & triangle = QuadratureRules<ct,2>::rule(GeometryType::simplex, p); 00766 const QuadratureRule<ct,1> & line = QuadratureRules<ct,1>::rule(GeometryType::cube, p); 00767 00768 this->delivered_order = std::min(triangle.order(),line.order()); 00769 00770 for (typename QuadratureRule<ct,1>::const_iterator 00771 lit = line.begin(); lit != line.end(); ++lit) 00772 { 00773 for (typename QuadratureRule<ct,2>::const_iterator 00774 tit = triangle.begin(); tit != triangle.end(); ++tit) 00775 { 00776 FieldVector<ct, d> local; 00777 local[0] = tit->position()[0]; 00778 local[1] = tit->position()[1]; 00779 local[2] = lit->position()[0]; 00780 00781 double weight = tit->weight() * lit->weight(); 00782 00783 // put in container 00784 this->push_back(QuadraturePoint<ct,d>(local,weight)); 00785 } 00786 } 00787 } 00788 } 00789 }; 00790 00792 class PyramidQuadraturePoints 00793 { 00794 public: 00795 enum { MAXP=8}; 00796 enum { highest_order=2 }; 00797 00799 PyramidQuadraturePoints() 00800 { 00801 int m = 0; 00802 O[m] = 0; 00803 00804 00805 // polynom degree 2 ??? 00806 m = 8; 00807 G[m][0][0] =0.58541020; 00808 G[m][0][1] =0.72819660; 00809 G[m][0][2] =0.13819660; 00810 00811 G[m][1][0] =0.13819660; 00812 G[m][1][1] =0.72819660; 00813 G[m][1][2] =0.13819660; 00814 00815 G[m][2][0] =0.13819660; 00816 G[m][2][1] =0.27630920; 00817 G[m][2][2] =0.58541020; 00818 00819 G[m][3][0] =0.13819660; 00820 G[m][3][1] =0.27630920; 00821 G[m][3][2] =0.13819660; 00822 00823 G[m][4][0] =0.72819660; 00824 G[m][4][1] =0.13819660; 00825 G[m][4][2] =0.13819660; 00826 00827 G[m][5][0] =0.72819660; 00828 G[m][5][1] =0.58541020; 00829 G[m][5][2] =0.13819660; 00830 00831 G[m][6][0] =0.27630920; 00832 G[m][6][1] =0.13819660; 00833 G[m][6][2] =0.58541020; 00834 00835 G[m][7][0] =0.27630920; 00836 G[m][7][1] =0.13819660; 00837 G[m][7][2] =0.13819660; 00838 00839 W[m][0] = 0.125 / 3.0; 00840 W[m][1] = 0.125 / 3.0; 00841 W[m][2] = 0.125 / 3.0; 00842 W[m][3] = 0.125 / 3.0; 00843 W[m][4] = 0.125 / 3.0; 00844 W[m][5] = 0.125 / 3.0; 00845 W[m][6] = 0.125 / 3.0; 00846 W[m][7] = 0.125 / 3.0; 00847 00848 O[m] = 2;// verify ???????? 00849 00850 } 00851 00853 FieldVector<double, 3> point(int m, int i) 00854 { 00855 return G[m][i]; 00856 } 00857 00859 double weight (int m, int i) 00860 { 00861 return W[m][i]; 00862 } 00863 00865 int order (int m) 00866 { 00867 return O[m]; 00868 } 00869 00870 private: 00871 FieldVector<double, 3> G[MAXP+1][MAXP]; 00872 double W[MAXP+1][MAXP]; // weights associated with points 00873 int O[MAXP+1]; // order of the rule 00874 }; 00875 00879 template<int dim> 00880 struct PyramidQuadraturePointsSingleton {}; 00881 00885 template<> 00886 struct PyramidQuadraturePointsSingleton<3> { 00887 static PyramidQuadraturePoints pyqp; 00888 }; 00889 00893 template<typename ct, int dim> 00894 class PyramidQuadratureRule; 00895 00899 template<typename ct> 00900 class PyramidQuadratureRule<ct,3> : public QuadratureRule<ct,3> 00901 { 00902 public: 00903 00905 enum{d=3}; 00906 00908 enum{highest_order=CubeQuadratureRule<ct,1>::highest_order}; 00909 00911 typedef ct CoordType; 00912 00914 typedef PyramidQuadratureRule<ct,3> value_type; 00915 00916 ~PyramidQuadratureRule(){} 00917 private: 00918 friend class QuadratureRuleFactory<ct,d>; 00919 PyramidQuadratureRule(int p) : QuadratureRule<ct,3>(GeometryType(GeometryType::pyramid, d)) 00920 { 00921 int m; 00922 00923 if (p>highest_order) 00924 DUNE_THROW(QuadratureOrderOutOfRange, 00925 "QuadratureRule for order " << p << " and GeometryType " 00926 << this->type() << " not available"); 00927 00928 if(false) { 00929 // if(p<=2) { 00930 m=8; 00931 this->delivered_order = PyramidQuadraturePointsSingleton<3>::pyqp.order(m); 00932 FieldVector<ct, d> local; 00933 double weight; 00934 for(int i=0;i<m;++i) 00935 { 00936 for(int k=0;k<d;++k) 00937 local[k]=PyramidQuadraturePointsSingleton<3>::pyqp.point(m,i)[k]; 00938 weight=PyramidQuadraturePointsSingleton<3>::pyqp.weight(m,i); 00939 // put in container 00940 this->push_back(QuadraturePoint<ct,d>(local,weight)); 00941 } 00942 } 00943 else 00944 { 00945 // Define the quadrature rules... 00946 QuadratureRule<ct,3> simplex = 00947 QuadratureRules<ct,3>::rule(GeometryType::simplex,p); 00948 00949 for (typename QuadratureRule<ct,3>::const_iterator 00950 it=simplex.begin(); it != simplex.end(); ++it) 00951 { 00952 FieldVector<ct,3> local = it->position(); 00953 ct weight = it->weight(); 00954 // Simplex 1 00955 // x:=x+y 00956 local[0] = local[0]+local[1]; 00957 this->push_back(QuadraturePoint<ct,d>(local,weight)); 00958 // Simplex 2 00959 // y:=x+y 00960 local[0] = it->position()[0]; 00961 local[1] = local[0]+local[1]; 00962 this->push_back(QuadraturePoint<ct,d>(local,weight)); 00963 } 00964 00965 this->delivered_order = simplex.order(); 00966 } 00967 } 00968 }; 00969 00976 template<typename ctype, int dim> 00977 class QuadratureRuleFactory { 00978 private: 00979 friend class QuadratureRules<ctype, dim>; 00980 static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt) 00981 { 00982 if (t.isCube()) 00983 { 00984 return CubeQuadratureRule<ctype,dim>(p); 00985 } 00986 if (t.isSimplex()) 00987 { 00988 return SimplexQuadratureRule<ctype,dim>(p); 00989 } 00990 DUNE_THROW(Exception, "Unknown GeometryType"); 00991 } 00992 }; 00993 00994 template<typename ctype> 00995 class QuadratureRuleFactory<ctype, 0> { 00996 private: 00997 enum { dim = 0 }; 00998 friend class QuadratureRules<ctype, dim>; 00999 static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt) 01000 { 01001 if (t.isVertex()) 01002 { 01003 return CubeQuadratureRule<ctype,dim>(p); 01004 } 01005 DUNE_THROW(Exception, "Unknown GeometryType"); 01006 } 01007 }; 01008 01009 template<typename ctype> 01010 class QuadratureRuleFactory<ctype, 1> { 01011 private: 01012 enum { dim = 1 }; 01013 friend class QuadratureRules<ctype, dim>; 01014 static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt) 01015 { 01016 if (t.isLine()) 01017 { 01018 switch (qt) { 01019 case QuadratureType::Gauss: 01020 return CubeQuadratureRule<ctype,dim>(p); 01021 case QuadratureType::Jacobian_1_0: 01022 return Jacobi1QuadratureRule<ctype,dim>(p); 01023 case QuadratureType::Jacobian_2_0: 01024 return Jacobi2QuadratureRule<ctype,dim>(p); 01025 default: 01026 DUNE_THROW(Exception, "Unknown QuadratureType"); 01027 } 01028 } 01029 DUNE_THROW(Exception, "Unknown GeometryType"); 01030 } 01031 }; 01032 01033 template<typename ctype> 01034 class QuadratureRuleFactory<ctype, 3> { 01035 private: 01036 enum { dim = 3 }; 01037 friend class QuadratureRules<ctype, dim>; 01038 static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt) 01039 { 01040 if (t.isCube()) 01041 { 01042 return CubeQuadratureRule<ctype,dim>(p); 01043 } 01044 if (t.isSimplex()) 01045 { 01046 return SimplexQuadratureRule<ctype,dim>(p); 01047 } 01048 if (t.isPrism()) 01049 { 01050 return PrismQuadratureRule<ctype,dim>(p); 01051 } 01052 if (t.isPyramid()) 01053 { 01054 return PyramidQuadratureRule<ctype,dim>(p); 01055 } 01056 DUNE_THROW(Exception, "Unknown GeometryType"); 01057 } 01058 }; 01059 01060 } // end namespace 01061 01062 #endif // DUNE_GEOMETRY_QUADRATURERULES_HH