4 #ifndef DUNE_GEOMETRY_QUADRATURERULES_HH
5 #define DUNE_GEOMETRY_QUADRATURERULES_HH
12 #include <dune/common/fvector.hh>
13 #include <dune/common/exceptions.hh>
14 #include <dune/common/stdstreams.hh>
34 template<
typename ct,
int dim>
43 typedef Dune::FieldVector<ct,dim>
Vector;
71 namespace QuadratureType {
91 template<
typename ct,
int dim>
120 typedef typename std::vector<QuadraturePoint<ct,dim> >::const_iterator
iterator;
130 int n =
power(m,dim);
131 for (
int i=0; i<n; i++)
136 for (
int k=0; k<dim; ++k)
144 FieldVector<ct, dim> local;
145 for (
int k=0; k<dim; k++)
147 local[k] = q1d[x[k]].position()[0];
148 weight *= q1d[x[k]].weight();
159 for (
int i=0; i<
d; i++) m *= y;
166 template<
typename ctype,
int dim>
class QuadratureRuleFactory;
171 template<
typename ctype,
int dim>
175 typedef std::pair<GeometryType,int> QuadratureRuleKey;
183 static std::map<QuadratureRuleKey, QuadratureRule> _quadratureMap;
184 QuadratureRuleKey key(t,p);
185 if (_quadratureMap.find(key) == _quadratureMap.end()) {
193 _quadratureMap[key] =
rule;
195 return _quadratureMap[key];
209 return instance()._rule(t,p,qt);
215 return instance()._rule(gt,p,qt);
224 template<
typename ct,
int dim>
256 template<
typename ct>
274 FieldVector<ct, dim> point(0.0);
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);
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);
305 template<
typename ct>
306 class CubeQuadratureRule<ct,1> :
320 CubeQuadratureRule (
int p)
324 std::vector< FieldVector<ct, dim> > _points;
325 std::vector< ct > _weight;
327 CubeQuadratureInitHelper<ct>::init
330 assert(_points.size() == _weight.size());
331 for (
size_t i = 0; i < _points.size(); i++)
336 extern template CubeQuadratureRule<float, 1>::CubeQuadratureRule(
int);
337 extern template CubeQuadratureRule<double, 1>::CubeQuadratureRule(
int);
341 #define DUNE_INCLUDING_IMPLEMENTATION
342 #include "quadraturerules/cube_imp.hh"
349 template<
typename ct,
int dim>
350 class Jacobi1QuadratureRule;
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);
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);
374 template<
typename ct>
375 class Jacobi1QuadratureRule<ct,1> :
386 enum { highest_order=61 };
397 Jacobi1QuadratureRule (
int p)
401 std::vector< FieldVector<ct, dim> > _points;
402 std::vector< ct > _weight;
406 Jacobi1QuadratureInitHelper<ct>::init
407 (p, _points, _weight, delivered_order);
409 assert(_points.size() == _weight.size());
410 for (
size_t i = 0; i < _points.size(); i++)
416 extern template Jacobi1QuadratureRule<float, 1>::Jacobi1QuadratureRule(
int);
417 extern template Jacobi1QuadratureRule<double, 1>::Jacobi1QuadratureRule(
int);
422 #define DUNE_INCLUDING_IMPLEMENTATION
423 #include "quadraturerules/jacobi_1_0_imp.hh"
430 template<
typename ct,
int dim>
431 class Jacobi2QuadratureRule;
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);
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);
455 template<
typename ct>
456 class Jacobi2QuadratureRule<ct,1> :
468 enum { highest_order=61 };
479 Jacobi2QuadratureRule (
int p)
483 std::vector< FieldVector<ct, dim> > _points;
484 std::vector< ct > _weight;
488 Jacobi2QuadratureInitHelper<ct>::init
489 (p, _points, _weight, delivered_order);
492 assert(_points.size() == _weight.size());
493 for (
size_t i = 0; i < _points.size(); i++)
499 extern template Jacobi2QuadratureRule<float, 1>::Jacobi2QuadratureRule(
int);
500 extern template Jacobi2QuadratureRule<double, 1>::Jacobi2QuadratureRule(
int);
505 #define DUNE_INCLUDING_IMPLEMENTATION
506 #include "quadraturerules/jacobi_2_0_imp.hh"
517 template<
typename ct,
int dim>
518 class SimplexQuadratureRule;
523 template<
typename ct>
542 SimplexQuadratureRule (
int p);
548 template<
typename ct>
567 SimplexQuadratureRule (
int p);
576 class PrismQuadraturePoints;
584 enum { highest_order=2 };
587 PrismQuadraturePoints ()
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;
630 G[m][0][0] =0.66666666666666666 ;
631 G[m][0][1] =0.16666666666666666 ;
632 G[m][0][2] =0.211324865405187 ;
634 G[m][1][0] = 0.16666666666666666;
635 G[m][1][1] =0.66666666666666666 ;
636 G[m][1][2] = 0.211324865405187;
638 G[m][2][0] = 0.16666666666666666;
639 G[m][2][1] = 0.16666666666666666;
640 G[m][2][2] = 0.211324865405187;
642 G[m][3][0] = 0.66666666666666666;
643 G[m][3][1] = 0.16666666666666666;
644 G[m][3][2] = 0.788675134594813;
646 G[m][4][0] = 0.16666666666666666;
647 G[m][4][1] = 0.66666666666666666;
648 G[m][4][2] = 0.788675134594813;
650 G[m][5][0] = 0.16666666666666666;
651 G[m][5][1] = 0.16666666666666666;
652 G[m][5][2] = 0.788675134594813;
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;
666 FieldVector<double, 3> point(
int m,
int i)
672 double weight (
int m,
int i)
684 FieldVector<double, 3> G[MAXP+1][MAXP];
686 double W[MAXP+1][MAXP];
710 template<
typename ct,
int dim>
711 class PrismQuadratureRule;
716 template<
typename ct>
729 < (
int)SimplexQuadratureRule<ct,2>::highest_order
731 : (
int)SimplexQuadratureRule<ct,2>::highest_order
740 ~PrismQuadratureRule(){}
747 "QuadratureRule for order " << p <<
" and GeometryType "
748 << this->type() <<
" not available");
755 FieldVector<ct,3> local;
756 for (
int k=0; k<d; k++)
768 this->delivered_order = std::min(triangle.
order(),line.
order());
771 lit = line.begin(); lit != line.end(); ++lit)
774 tit = triangle.begin(); tit != triangle.end(); ++tit)
776 FieldVector<ct, d> local;
777 local[0] = tit->position()[0];
778 local[1] = tit->position()[1];
779 local[2] = lit->position()[0];
781 double weight = tit->weight() * lit->weight();
796 enum { highest_order=2 };
807 G[m][0][0] =0.58541020;
808 G[m][0][1] =0.72819660;
809 G[m][0][2] =0.13819660;
811 G[m][1][0] =0.13819660;
812 G[m][1][1] =0.72819660;
813 G[m][1][2] =0.13819660;
815 G[m][2][0] =0.13819660;
816 G[m][2][1] =0.27630920;
817 G[m][2][2] =0.58541020;
819 G[m][3][0] =0.13819660;
820 G[m][3][1] =0.27630920;
821 G[m][3][2] =0.13819660;
823 G[m][4][0] =0.72819660;
824 G[m][4][1] =0.13819660;
825 G[m][4][2] =0.13819660;
827 G[m][5][0] =0.72819660;
828 G[m][5][1] =0.58541020;
829 G[m][5][2] =0.13819660;
831 G[m][6][0] =0.27630920;
832 G[m][6][1] =0.13819660;
833 G[m][6][2] =0.58541020;
835 G[m][7][0] =0.27630920;
836 G[m][7][1] =0.13819660;
837 G[m][7][2] =0.13819660;
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;
853 FieldVector<double, 3> point(
int m,
int i)
859 double weight (
int m,
int i)
871 FieldVector<double, 3> G[MAXP+1][MAXP];
872 double W[MAXP+1][MAXP];
893 template<
typename ct,
int dim>
894 class PyramidQuadratureRule;
899 template<
typename ct>
925 "QuadratureRule for order " << p <<
" and GeometryType "
926 << this->type() <<
" not available");
932 FieldVector<ct, d> local;
950 it=simplex.begin(); it != simplex.end(); ++it)
952 FieldVector<ct,3> local = it->position();
953 ct weight = it->weight();
956 local[0] = local[0]+local[1];
960 local[0] = it->position()[0];
961 local[1] = local[0]+local[1];
965 this->delivered_order = simplex.
order();
976 template<
typename ctype,
int dim>
988 return SimplexQuadratureRule<ctype,dim>(p);
990 DUNE_THROW(Exception,
"Unknown GeometryType");
994 template<
typename ctype>
1005 DUNE_THROW(Exception,
"Unknown GeometryType");
1009 template<
typename ctype>
1022 return Jacobi1QuadratureRule<ctype,dim>(p);
1024 return Jacobi2QuadratureRule<ctype,dim>(p);
1026 DUNE_THROW(Exception,
"Unknown QuadratureType");
1029 DUNE_THROW(Exception,
"Unknown GeometryType");
1033 template<
typename ctype>
1046 return SimplexQuadratureRule<ctype,dim>(p);
1050 return PrismQuadratureRule<ctype,dim>(p);
1054 return PyramidQuadratureRule<ctype,dim>(p);
1056 DUNE_THROW(Exception,
"Unknown GeometryType");
1062 #endif // DUNE_GEOMETRY_QUADRATURERULES_HH