dune-geometry  2.2.1
quadraturerules.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #ifndef DUNE_GEOMETRY_QUADRATURERULES_HH
5 #define DUNE_GEOMETRY_QUADRATURERULES_HH
6 
7 #include <iostream>
8 #include <limits>
9 #include <vector>
10 #include <map>
11 
12 #include <dune/common/fvector.hh>
13 #include <dune/common/exceptions.hh>
14 #include <dune/common/stdstreams.hh>
15 
16 #include <dune/geometry/type.hh>
17 
23 namespace Dune {
24 
29  class QuadratureOrderOutOfRange : public NotImplemented {};
30 
34  template<typename ct, int dim>
36  public:
37  // compile time parameters
38  enum { d=dim };
39  typedef ct CoordType;
40 
41  static const unsigned int dimension = d;
42  typedef ct Field;
43  typedef Dune::FieldVector<ct,dim> Vector;
44 
46  QuadraturePoint (const Vector& x, ct w) : local(x)
47  {
48  wght = w;
49  }
50 
52  const Vector& position () const
53  {
54  return local;
55  }
56 
58  const ct &weight () const
59  {
60  return wght;
61  }
62 
63  protected:
64  FieldVector<ct, dim> local;
65  ct wght;
66  };
67 
71  namespace QuadratureType {
72  enum Enum {
73  Gauss = 0,
74 
77 
78  Simpson = 3,
79  Trap = 4,
80  Grid = 5,
81 
82  Clough = 21,
83 
85  };
86  }
87 
91  template<typename ct, int dim>
92  class QuadratureRule : public std::vector<QuadraturePoint<ct,dim> >
93  {
94  public:
95 
98 
101 
104 
106  enum { d=dim };
107 
109  typedef ct CoordType;
110 
112  virtual int order () const { return delivered_order; }
113 
115  virtual GeometryType type () const { return geometry_type; }
116  virtual ~QuadratureRule(){}
117 
120  typedef typename std::vector<QuadraturePoint<ct,dim> >::const_iterator iterator;
121 
122  protected:
125 
127  {
128  // fill in all the gauss points
129  int m = q1d.size();
130  int n = power(m,dim);
131  for (int i=0; i<n; i++)
132  {
133  // compute multi index for Gauss point
134  int x[dim];
135  int z = i;
136  for (int k=0; k<dim; ++k)
137  {
138  x[k] = z%m;
139  z = z/m;
140  }
141 
142  // compute coordinates and weight
143  double weight = 1.0;
144  FieldVector<ct, dim> local;
145  for (int k=0; k<dim; k++)
146  {
147  local[k] = q1d[x[k]].position()[0];
148  weight *= q1d[x[k]].weight();
149  }
150 
151  // put in container
152  this->push_back(QuadraturePoint<ct,dim>(local,weight));
153  }
154  }
155 
156  int power (int y, int d)
157  {
158  int m=1;
159  for (int i=0; i<d; i++) m *= y;
160  return m;
161  }
162  };
163 
164  // Forward declaration of the factory class,
165  // needed internally by the QuadratureRules container class.
166  template<typename ctype, int dim> class QuadratureRuleFactory;
167 
171  template<typename ctype, int dim>
173 
175  typedef std::pair<GeometryType,int> QuadratureRuleKey;
176 
179 
182  {
183  static std::map<QuadratureRuleKey, QuadratureRule> _quadratureMap;
184  QuadratureRuleKey key(t,p);
185  if (_quadratureMap.find(key) == _quadratureMap.end()) {
186  /*
187  The rule must be acquired before we can store it.
188  If we write this in one command, an invalid rule
189  would get stored in case of an exception.
190  */
193  _quadratureMap[key] = rule;
194  }
195  return _quadratureMap[key];
196  }
198  static QuadratureRules& instance()
199  {
200  static QuadratureRules instance;
201  return instance;
202  }
204  QuadratureRules () {}
205  public:
208  {
209  return instance()._rule(t,p,qt);
210  }
213  {
214  GeometryType gt(t,dim);
215  return instance()._rule(gt,p,qt);
216  }
217  };
218 
219  /***********************************************************/
220 
224  template<typename ct, int dim>
226  public QuadratureRule<ct,dim>
227  {
228  public:
230  enum { d=dim };
231 
234 
236  typedef ct CoordType;
237 
240 
242  private:
243  friend class QuadratureRuleFactory<ct,dim>;
245  CubeQuadratureRule (int p) : QuadratureRule<ct,dim>(GeometryType(GeometryType::cube, d))
246  {
248  this->tensor_product( q1D );
249  this->delivered_order = q1D.order();
250  }
251 
252  };
253 
256  template<typename ct>
257  class CubeQuadratureRule<ct,0> :
258  public QuadratureRule<ct,0>
259  {
260  public:
261  // compile time parameters
262  enum { d=0 };
263  enum { dim=0 };
264  enum { highest_order=61 };
265  typedef ct CoordType;
267 
269  private:
270  friend class QuadratureRuleFactory<ct,dim>;
271  CubeQuadratureRule (int p):
272  QuadratureRule<ct,0>(GeometryType(GeometryType::cube, 0))
273  {
274  FieldVector<ct, dim> point(0.0);
275 
276  if (p > highest_order)
277  DUNE_THROW(QuadratureOrderOutOfRange, "Quadrature rule " << p << " not supported!");
278 
280  this->push_back(QuadraturePoint<ct,dim>(point, 1.0));
281  }
282  };
283 
285  template<typename ct,
286  bool fundamental = std::numeric_limits<ct>::is_specialized>
287  struct CubeQuadratureInitHelper;
288  template<typename ct>
289  struct CubeQuadratureInitHelper<ct, true> {
290  static void init(int p,
291  std::vector< FieldVector<ct, 1> > & _points,
292  std::vector< ct > & _weight,
293  int & delivered_order);
294  };
295  template<typename ct>
296  struct CubeQuadratureInitHelper<ct, false> {
297  static void init(int p,
298  std::vector< FieldVector<ct, 1> > & _points,
299  std::vector< ct > & _weight,
300  int & delivered_order);
301  };
302 
305  template<typename ct>
306  class CubeQuadratureRule<ct,1> :
307  public QuadratureRule<ct,1>
308  {
309  public:
310  // compile time parameters
311  enum { d=1 };
312  enum { dim=1 };
313  enum { highest_order=61 };
314  typedef ct CoordType;
315  typedef CubeQuadratureRule value_type;
316 
318  private:
319  friend class QuadratureRuleFactory<ct,dim>;
320  CubeQuadratureRule (int p)
321  : QuadratureRule<ct,1>(GeometryType(GeometryType::cube, 1))
322  {
324  std::vector< FieldVector<ct, dim> > _points;
325  std::vector< ct > _weight;
326 
327  CubeQuadratureInitHelper<ct>::init
328  (p, _points, _weight, this->delivered_order);
329 
330  assert(_points.size() == _weight.size());
331  for (size_t i = 0; i < _points.size(); i++)
332  this->push_back(QuadraturePoint<ct,dim>(_points[i], _weight[i]));
333  }
334  };
335 
336  extern template CubeQuadratureRule<float, 1>::CubeQuadratureRule(int);
337  extern template CubeQuadratureRule<double, 1>::CubeQuadratureRule(int);
338 
339 } // namespace Dune
340 
341 #define DUNE_INCLUDING_IMPLEMENTATION
342 #include "quadraturerules/cube_imp.hh"
343 
344 namespace Dune {
345 
349  template<typename ct, int dim>
350  class Jacobi1QuadratureRule;
351 
353  template<typename ct,
354  bool fundamental = std::numeric_limits<ct>::is_specialized>
355  struct Jacobi1QuadratureInitHelper;
356  template<typename ct>
357  struct Jacobi1QuadratureInitHelper<ct, true> {
358  static void init(int p,
359  std::vector< FieldVector<ct, 1> > & _points,
360  std::vector< ct > & _weight,
361  int & delivered_order);
362  };
363  template<typename ct>
364  struct Jacobi1QuadratureInitHelper<ct, false> {
365  static void init(int p,
366  std::vector< FieldVector<ct, 1> > & _points,
367  std::vector< ct > & _weight,
368  int & delivered_order);
369  };
370 
374  template<typename ct>
375  class Jacobi1QuadratureRule<ct,1> :
376  public QuadratureRule<ct,1>
377  {
378  public:
380  enum { d=1 };
383  enum { dim=1 };
384 
386  enum { highest_order=61 };
387 
389  typedef ct CoordType;
390 
392  typedef Jacobi1QuadratureRule value_type;
393 
395  private:
396  friend class QuadratureRuleFactory<ct,dim>;
397  Jacobi1QuadratureRule (int p)
398  : QuadratureRule<ct,1>(GeometryType(GeometryType::cube, 1))
399  {
401  std::vector< FieldVector<ct, dim> > _points;
402  std::vector< ct > _weight;
403 
404  int delivered_order;
405 
406  Jacobi1QuadratureInitHelper<ct>::init
407  (p, _points, _weight, delivered_order);
408  this->delivered_order = delivered_order;
409  assert(_points.size() == _weight.size());
410  for (size_t i = 0; i < _points.size(); i++)
411  this->push_back(QuadraturePoint<ct,dim>(_points[i], _weight[i]));
412  }
413  };
414 
415 #ifndef DOXYGEN
416  extern template Jacobi1QuadratureRule<float, 1>::Jacobi1QuadratureRule(int);
417  extern template Jacobi1QuadratureRule<double, 1>::Jacobi1QuadratureRule(int);
418 #endif // !DOXYGEN
419 
420 } // namespace Dune
421 
422 #define DUNE_INCLUDING_IMPLEMENTATION
423 #include "quadraturerules/jacobi_1_0_imp.hh"
424 
425 namespace Dune {
426 
430  template<typename ct, int dim>
431  class Jacobi2QuadratureRule;
432 
434  template<typename ct,
435  bool fundamental = std::numeric_limits<ct>::is_specialized>
436  struct Jacobi2QuadratureInitHelper;
437  template<typename ct>
438  struct Jacobi2QuadratureInitHelper<ct, true> {
439  static void init(int p,
440  std::vector< FieldVector<ct, 1> > & _points,
441  std::vector< ct > & _weight,
442  int & delivered_order);
443  };
444  template<typename ct>
445  struct Jacobi2QuadratureInitHelper<ct, false> {
446  static void init(int p,
447  std::vector< FieldVector<ct, 1> > & _points,
448  std::vector< ct > & _weight,
449  int & delivered_order);
450  };
451 
455  template<typename ct>
456  class Jacobi2QuadratureRule<ct,1> :
457  public QuadratureRule<ct,1>
458  {
459  public:
461  enum { d=1 };
462 
465  enum { dim=1 };
466 
468  enum { highest_order=61 };
469 
471  typedef ct CoordType;
472 
474  typedef Jacobi2QuadratureRule value_type;
475 
477  private:
478  friend class QuadratureRuleFactory<ct,dim>;
479  Jacobi2QuadratureRule (int p)
480  : QuadratureRule<ct,1>(GeometryType(GeometryType::cube, 1))
481  {
483  std::vector< FieldVector<ct, dim> > _points;
484  std::vector< ct > _weight;
485 
486  int delivered_order;
487 
488  Jacobi2QuadratureInitHelper<ct>::init
489  (p, _points, _weight, delivered_order);
490 
491  this->delivered_order = delivered_order;
492  assert(_points.size() == _weight.size());
493  for (size_t i = 0; i < _points.size(); i++)
494  this->push_back(QuadraturePoint<ct,dim>(_points[i], _weight[i]));
495  }
496  };
497 
498 #ifndef DOXYGEN
499  extern template Jacobi2QuadratureRule<float, 1>::Jacobi2QuadratureRule(int);
500  extern template Jacobi2QuadratureRule<double, 1>::Jacobi2QuadratureRule(int);
501 #endif // !DOXYGEN
502 
503 } // namespace Dune
504 
505 #define DUNE_INCLUDING_IMPLEMENTATION
506 #include "quadraturerules/jacobi_2_0_imp.hh"
507 
508 namespace Dune {
509 
510  /************************************************
511  * Quadraturerule for Simplices/Triangle
512  *************************************************/
513 
517  template<typename ct, int dim>
518  class SimplexQuadratureRule;
519 
523  template<typename ct>
525  {
526  public:
527 
529  enum{d=2};
530 
532  enum { highest_order=CubeQuadratureRule<ct,1>::highest_order -1 };
533 
535  typedef ct CoordType;
536 
538  typedef SimplexQuadratureRule value_type;
540  private:
542  SimplexQuadratureRule (int p);
543  };
544 
548  template<typename ct>
550  {
551  public:
552 
554  enum{d=3};
555 
557  enum { highest_order=CubeQuadratureRule<ct,1>::highest_order -2 };
558 
560  typedef ct CoordType;
561 
565  private:
567  SimplexQuadratureRule (int p);
568  };
569 
570 /***********************************
571  * quadrature for Prism
572  **********************************/
573 
575  template<int dim>
576  class PrismQuadraturePoints;
577 
579  template<>
581  {
582  public:
583  enum { MAXP=6};
584  enum { highest_order=2 };
585 
587  PrismQuadraturePoints ()
588  {
589  int m = 0;
590  O[m] = 0;
591 
592  // polynom degree 0 ???
593  m = 6;
594  G[m][0][0] = 0.0;
595  G[m][0][1] = 0.0;
596  G[m][0][2] = 0.0;
597 
598  G[m][1][0] = 1.0;
599  G[m][1][1] = 0.0;
600  G[m][1][2] = 0.0;
601 
602  G[m][2][0] = 0.0;
603  G[m][2][1] = 1.0;
604  G[m][2][2] = 0.0;
605 
606  G[m][3][0] = 0.0;
607  G[m][3][1] = 0.0;
608  G[m][3][2] = 1.0;
609 
610  G[m][4][0] = 1.0;
611  G[m][4][1] = 0.0;
612  G[m][4][2] = 1.0;
613 
614  G[m][5][0] = 0.0;
615  G[m][5][1] = 0.1;
616  G[m][5][2] = 1.0;
617 
618  W[m][0] = 0.16666666666666666 / 2.0;
619  W[m][1] = 0.16666666666666666 / 2.0;
620  W[m][2] = 0.16666666666666666 / 2.0;
621  W[m][3] = 0.16666666666666666 / 2.0;
622  W[m][4] = 0.16666666666666666 / 2.0;
623  W[m][5] = 0.16666666666666666 / 2.0;
624 
625  O[m] = 0;// verify ????????
626 
627 
628  // polynom degree 2 ???
629  m = 6;
630  G[m][0][0] =0.66666666666666666 ;
631  G[m][0][1] =0.16666666666666666 ;
632  G[m][0][2] =0.211324865405187 ;
633 
634  G[m][1][0] = 0.16666666666666666;
635  G[m][1][1] =0.66666666666666666 ;
636  G[m][1][2] = 0.211324865405187;
637 
638  G[m][2][0] = 0.16666666666666666;
639  G[m][2][1] = 0.16666666666666666;
640  G[m][2][2] = 0.211324865405187;
641 
642  G[m][3][0] = 0.66666666666666666;
643  G[m][3][1] = 0.16666666666666666;
644  G[m][3][2] = 0.788675134594813;
645 
646  G[m][4][0] = 0.16666666666666666;
647  G[m][4][1] = 0.66666666666666666;
648  G[m][4][2] = 0.788675134594813;
649 
650  G[m][5][0] = 0.16666666666666666;
651  G[m][5][1] = 0.16666666666666666;
652  G[m][5][2] = 0.788675134594813;
653 
654  W[m][0] = 0.16666666666666666 / 2.0;
655  W[m][1] = 0.16666666666666666 / 2.0;
656  W[m][2] = 0.16666666666666666 / 2.0;
657  W[m][3] = 0.16666666666666666 / 2.0;
658  W[m][4] = 0.16666666666666666 / 2.0;
659  W[m][5] = 0.16666666666666666 / 2.0;
660 
661  O[m] = 2;// verify ????????
662 
663  }
664 
666  FieldVector<double, 3> point(int m, int i)
667  {
668  return G[m][i];
669  }
670 
672  double weight (int m, int i)
673  {
674  return W[m][i];
675  }
676 
678  int order (int m)
679  {
680  return O[m];
681  }
682 
683  private:
684  FieldVector<double, 3> G[MAXP+1][MAXP];//positions
685 
686  double W[MAXP+1][MAXP]; // weights associated with points
687  int O[MAXP+1]; // order of the rule
688  };
689 
690 
694  template<int dim>
697  };
698 
702  template<>
705  };
706 
710  template<typename ct, int dim>
711  class PrismQuadratureRule;
712 
716  template<typename ct>
718  {
719  public:
720 
722  enum{ d=3 };
723 
725  enum{
726  /* min(Line::order, Triangle::order) */
727  highest_order =
729  < (int)SimplexQuadratureRule<ct,2>::highest_order
731  : (int)SimplexQuadratureRule<ct,2>::highest_order
732  };
733 
735  typedef ct CoordType;
736 
738  typedef PrismQuadratureRule<ct,3> value_type;
739 
740  ~PrismQuadratureRule(){}
741  private:
743  PrismQuadratureRule(int p) : QuadratureRule<ct,3>(GeometryType(GeometryType::prism, d))
744  {
745  if (p>highest_order)
746  DUNE_THROW(QuadratureOrderOutOfRange,
747  "QuadratureRule for order " << p << " and GeometryType "
748  << this->type() << " not available");
749 
750  if (p<=2) {
751  int m=6;
752  this->delivered_order = PrismQuadraturePointsSingleton<3>::prqp.order(m);
753  for(int i=0;i<m;++i)
754  {
755  FieldVector<ct,3> local;
756  for (int k=0; k<d; k++)
757  local[k] = PrismQuadraturePointsSingleton<3>::prqp.point(m,i)[k];
758  double weight =
760  // put in container
761  this->push_back(QuadraturePoint<ct,d>(local,weight));
762  }
763  }
764  else {
765  const QuadratureRule<ct,2> & triangle = QuadratureRules<ct,2>::rule(GeometryType::simplex, p);
766  const QuadratureRule<ct,1> & line = QuadratureRules<ct,1>::rule(GeometryType::cube, p);
767 
768  this->delivered_order = std::min(triangle.order(),line.order());
769 
771  lit = line.begin(); lit != line.end(); ++lit)
772  {
774  tit = triangle.begin(); tit != triangle.end(); ++tit)
775  {
776  FieldVector<ct, d> local;
777  local[0] = tit->position()[0];
778  local[1] = tit->position()[1];
779  local[2] = lit->position()[0];
780 
781  double weight = tit->weight() * lit->weight();
782 
783  // put in container
784  this->push_back(QuadraturePoint<ct,d>(local,weight));
785  }
786  }
787  }
788  }
789  };
790 
793  {
794  public:
795  enum { MAXP=8};
796  enum { highest_order=2 };
797 
800  {
801  int m = 0;
802  O[m] = 0;
803 
804 
805  // polynom degree 2 ???
806  m = 8;
807  G[m][0][0] =0.58541020;
808  G[m][0][1] =0.72819660;
809  G[m][0][2] =0.13819660;
810 
811  G[m][1][0] =0.13819660;
812  G[m][1][1] =0.72819660;
813  G[m][1][2] =0.13819660;
814 
815  G[m][2][0] =0.13819660;
816  G[m][2][1] =0.27630920;
817  G[m][2][2] =0.58541020;
818 
819  G[m][3][0] =0.13819660;
820  G[m][3][1] =0.27630920;
821  G[m][3][2] =0.13819660;
822 
823  G[m][4][0] =0.72819660;
824  G[m][4][1] =0.13819660;
825  G[m][4][2] =0.13819660;
826 
827  G[m][5][0] =0.72819660;
828  G[m][5][1] =0.58541020;
829  G[m][5][2] =0.13819660;
830 
831  G[m][6][0] =0.27630920;
832  G[m][6][1] =0.13819660;
833  G[m][6][2] =0.58541020;
834 
835  G[m][7][0] =0.27630920;
836  G[m][7][1] =0.13819660;
837  G[m][7][2] =0.13819660;
838 
839  W[m][0] = 0.125 / 3.0;
840  W[m][1] = 0.125 / 3.0;
841  W[m][2] = 0.125 / 3.0;
842  W[m][3] = 0.125 / 3.0;
843  W[m][4] = 0.125 / 3.0;
844  W[m][5] = 0.125 / 3.0;
845  W[m][6] = 0.125 / 3.0;
846  W[m][7] = 0.125 / 3.0;
847 
848  O[m] = 2;// verify ????????
849 
850  }
851 
853  FieldVector<double, 3> point(int m, int i)
854  {
855  return G[m][i];
856  }
857 
859  double weight (int m, int i)
860  {
861  return W[m][i];
862  }
863 
865  int order (int m)
866  {
867  return O[m];
868  }
869 
870  private:
871  FieldVector<double, 3> G[MAXP+1][MAXP];
872  double W[MAXP+1][MAXP]; // weights associated with points
873  int O[MAXP+1]; // order of the rule
874  };
875 
879  template<int dim>
881 
885  template<>
888  };
889 
893  template<typename ct, int dim>
894  class PyramidQuadratureRule;
895 
899  template<typename ct>
901  {
902  public:
903 
905  enum{d=3};
906 
909 
911  typedef ct CoordType;
912 
915 
917  private:
919  PyramidQuadratureRule(int p) : QuadratureRule<ct,3>(GeometryType(GeometryType::pyramid, d))
920  {
921  int m;
922 
923  if (p>highest_order)
924  DUNE_THROW(QuadratureOrderOutOfRange,
925  "QuadratureRule for order " << p << " and GeometryType "
926  << this->type() << " not available");
927 
928  if(false) {
929 // if(p<=2) {
930  m=8;
931  this->delivered_order = PyramidQuadraturePointsSingleton<3>::pyqp.order(m);
932  FieldVector<ct, d> local;
933  double weight;
934  for(int i=0;i<m;++i)
935  {
936  for(int k=0;k<d;++k)
937  local[k]=PyramidQuadraturePointsSingleton<3>::pyqp.point(m,i)[k];
939  // put in container
940  this->push_back(QuadraturePoint<ct,d>(local,weight));
941  }
942  }
943  else
944  {
945  // Define the quadrature rules...
946  QuadratureRule<ct,3> simplex =
947  QuadratureRules<ct,3>::rule(GeometryType::simplex,p);
948 
950  it=simplex.begin(); it != simplex.end(); ++it)
951  {
952  FieldVector<ct,3> local = it->position();
953  ct weight = it->weight();
954  // Simplex 1
955  // x:=x+y
956  local[0] = local[0]+local[1];
957  this->push_back(QuadraturePoint<ct,d>(local,weight));
958  // Simplex 2
959  // y:=x+y
960  local[0] = it->position()[0];
961  local[1] = local[0]+local[1];
962  this->push_back(QuadraturePoint<ct,d>(local,weight));
963  }
964 
965  this->delivered_order = simplex.order();
966  }
967  }
968  };
969 
976  template<typename ctype, int dim>
978  private:
980  static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
981  {
982  if (t.isCube())
983  {
985  }
986  if (t.isSimplex())
987  {
988  return SimplexQuadratureRule<ctype,dim>(p);
989  }
990  DUNE_THROW(Exception, "Unknown GeometryType");
991  }
992  };
993 
994  template<typename ctype>
996  private:
997  enum { dim = 0 };
999  static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
1000  {
1001  if (t.isVertex())
1002  {
1004  }
1005  DUNE_THROW(Exception, "Unknown GeometryType");
1006  }
1007  };
1008 
1009  template<typename ctype>
1011  private:
1012  enum { dim = 1 };
1014  static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
1015  {
1016  if (t.isLine())
1017  {
1018  switch (qt) {
1019  case QuadratureType::Gauss:
1022  return Jacobi1QuadratureRule<ctype,dim>(p);
1024  return Jacobi2QuadratureRule<ctype,dim>(p);
1025  default:
1026  DUNE_THROW(Exception, "Unknown QuadratureType");
1027  }
1028  }
1029  DUNE_THROW(Exception, "Unknown GeometryType");
1030  }
1031  };
1032 
1033  template<typename ctype>
1035  private:
1036  enum { dim = 3 };
1038  static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
1039  {
1040  if (t.isCube())
1041  {
1043  }
1044  if (t.isSimplex())
1045  {
1046  return SimplexQuadratureRule<ctype,dim>(p);
1047  }
1048  if (t.isPrism())
1049  {
1050  return PrismQuadratureRule<ctype,dim>(p);
1051  }
1052  if (t.isPyramid())
1053  {
1054  return PyramidQuadratureRule<ctype,dim>(p);
1055  }
1056  DUNE_THROW(Exception, "Unknown GeometryType");
1057  }
1058  };
1059 
1060 } // end namespace
1061 
1062 #endif // DUNE_GEOMETRY_QUADRATURERULES_HH