dune-geometry  2.2.0
quadraturerules.hh
Go to the documentation of this file.
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