dune-geometry  2.3.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 #include <dune/common/visibility.hh>
16 #include <dune/common/deprecated.hh>
17 
18 #include <dune/geometry/type.hh>
19 
25 namespace Dune {
26 
31  class QuadratureOrderOutOfRange : public NotImplemented {};
32 
38  template<typename ct, int dim>
40  public:
41  // compile time parameters
42  enum { d=dim } DUNE_DEPRECATED_MSG("Use 'dimension' instead");
43  typedef ct CoordType DUNE_DEPRECATED_MSG("Use type 'Field' instead");
44 
45  enum { dimension = dim };
46  typedef ct Field;
47  typedef Dune::FieldVector<ct,dim> Vector;
48 
50  QuadraturePoint (const Vector& x, ct w) : local(x)
51  {
52  weight_ = w;
53  }
54 
56  const Vector& position () const
57  {
58  return local;
59  }
60 
62  const ct &weight () const
63  {
64  return weight_;
65  }
66 
67  protected:
68  FieldVector<ct, dim> local;
69  ct weight_;
70  };
71 
75  namespace QuadratureType {
76  enum Enum {
77  Gauss = 0,
79 
84 
86  };
87  }
88 
92  template<typename ct, int dim>
93  class QuadratureRule : public std::vector<QuadraturePoint<ct,dim> >
94  {
95  protected:
96 
99 
102 
105  public:
107  enum { d=dim };
108 
110  typedef ct CoordType;
111 
113  virtual int order () const { return delivered_order; }
114 
116  virtual GeometryType type () const { return geometry_type; }
117  virtual ~QuadratureRule(){}
118 
121  typedef typename std::vector<QuadraturePoint<ct,dim> >::const_iterator iterator;
122 
123  protected:
126  };
127 
128  // Forward declaration of the factory class,
129  // needed internally by the QuadratureRules container class.
130  template<typename ctype, int dim> class QuadratureRuleFactory;
131 
135  template<typename ctype, int dim>
137 
141  typedef std::pair<GeometryType,int> QuadratureRuleKey;
142 
145 
147  DUNE_EXPORT const QuadratureRule& _rule(const GeometryType& t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
148  {
149  static std::map<QuadratureRuleKey, QuadratureRule> _quadratureMap;
150  QuadratureRuleKey key(t,p);
151  if (_quadratureMap.find(key) == _quadratureMap.end()) {
152  /*
153  The rule must be acquired before we can store it.
154  If we write this in one command, an invalid rule
155  would get stored in case of an exception.
156  */
159  _quadratureMap.insert(std::make_pair(key, rule));
160  }
161  return _quadratureMap.find(key)->second;
162  }
164  DUNE_EXPORT static QuadratureRules& instance()
165  {
166  static QuadratureRules instance;
167  return instance;
168  }
170  QuadratureRules () {}
171  public:
174  {
175  return instance()._rule(t,p,qt);
176  }
179  {
180  GeometryType gt(t,dim);
181  return instance()._rule(gt,p,qt);
182  }
183  };
184 
185 } // end namespace Dune
186 
187 #include "quadraturerules/pointquadrature.hh"
188 
189 namespace Dune {
190 
192  template<typename ct, bool fundamental = std::numeric_limits<ct>::is_specialized>
194  template<typename ct>
195  struct GaussQuadratureInitHelper<ct, true> {
196  static void init(int p,
197  std::vector< FieldVector<ct, 1> > & _points,
198  std::vector< ct > & _weight,
199  int & delivered_order);
200  };
201  template<typename ct>
202  struct GaussQuadratureInitHelper<ct, false> {
203  static void init(int p,
204  std::vector< FieldVector<ct, 1> > & _points,
205  std::vector< ct > & _weight,
206  int & delivered_order);
207  };
208 
210  template<typename ct>
212  public QuadratureRule<ct,1>
213  {
214  public:
215  // compile time parameters
216  enum { dim=1 };
217  enum { highest_order=61 };
218 
220  private:
221  friend class QuadratureRuleFactory<ct,dim>;
222  GaussQuadratureRule1D (int p)
223  : QuadratureRule<ct,1>(GeometryType(GeometryType::cube, 1))
224  {
226  std::vector< FieldVector<ct, dim> > _points;
227  std::vector< ct > _weight;
228 
230  (p, _points, _weight, this->delivered_order);
231 
232  assert(_points.size() == _weight.size());
233  for (size_t i = 0; i < _points.size(); i++)
234  this->push_back(QuadraturePoint<ct,dim>(_points[i], _weight[i]));
235  }
236  };
237 
238  extern template GaussQuadratureRule1D<float>::GaussQuadratureRule1D(int);
239  extern template GaussQuadratureRule1D<double>::GaussQuadratureRule1D(int);
240 
241 } // namespace Dune
242 
243 #define DUNE_INCLUDING_IMPLEMENTATION
244 #include "quadraturerules/gauss_imp.hh"
245 
246 namespace Dune {
247 
249  template<typename ct,
250  bool fundamental = std::numeric_limits<ct>::is_specialized>
252  template<typename ct>
253  struct Jacobi1QuadratureInitHelper<ct, true> {
254  static void init(int p,
255  std::vector< FieldVector<ct, 1> > & _points,
256  std::vector< ct > & _weight,
257  int & delivered_order);
258  };
259  template<typename ct>
260  struct Jacobi1QuadratureInitHelper<ct, false> {
261  static void init(int p,
262  std::vector< FieldVector<ct, 1> > & _points,
263  std::vector< ct > & _weight,
264  int & delivered_order);
265  };
266 
270  template<typename ct>
272  public QuadratureRule<ct,1>
273  {
274  public:
276  enum { dim=1 };
277 
279  enum { highest_order=61 };
280 
282  private:
283  friend class QuadratureRuleFactory<ct,dim>;
285  : QuadratureRule<ct,1>(GeometryType(GeometryType::cube, 1))
286  {
288  std::vector< FieldVector<ct, dim> > _points;
289  std::vector< ct > _weight;
290 
291  int delivered_order;
292 
294  (p, _points, _weight, delivered_order);
295  this->delivered_order = delivered_order;
296  assert(_points.size() == _weight.size());
297  for (size_t i = 0; i < _points.size(); i++)
298  this->push_back(QuadraturePoint<ct,dim>(_points[i], _weight[i]));
299  }
300  };
301 
302 #ifndef DOXYGEN
303  extern template Jacobi1QuadratureRule1D<float>::Jacobi1QuadratureRule1D(int);
304  extern template Jacobi1QuadratureRule1D<double>::Jacobi1QuadratureRule1D(int);
305 #endif // !DOXYGEN
306 
307 } // namespace Dune
308 
309 #define DUNE_INCLUDING_IMPLEMENTATION
310 #include "quadraturerules/jacobi_1_0_imp.hh"
311 
312 namespace Dune {
313 
315  template<typename ct,
316  bool fundamental = std::numeric_limits<ct>::is_specialized>
318  template<typename ct>
319  struct Jacobi2QuadratureInitHelper<ct, true> {
320  static void init(int p,
321  std::vector< FieldVector<ct, 1> > & _points,
322  std::vector< ct > & _weight,
323  int & delivered_order);
324  };
325  template<typename ct>
326  struct Jacobi2QuadratureInitHelper<ct, false> {
327  static void init(int p,
328  std::vector< FieldVector<ct, 1> > & _points,
329  std::vector< ct > & _weight,
330  int & delivered_order);
331  };
332 
336  template<typename ct>
338  public QuadratureRule<ct,1>
339  {
340  public:
342  enum { dim=1 };
343 
345  enum { highest_order=61 };
346 
348  private:
349  friend class QuadratureRuleFactory<ct,dim>;
351  : QuadratureRule<ct,1>(GeometryType(GeometryType::cube, 1))
352  {
354  std::vector< FieldVector<ct, dim> > _points;
355  std::vector< ct > _weight;
356 
357  int delivered_order;
358 
360  (p, _points, _weight, delivered_order);
361 
362  this->delivered_order = delivered_order;
363  assert(_points.size() == _weight.size());
364  for (size_t i = 0; i < _points.size(); i++)
365  this->push_back(QuadraturePoint<ct,dim>(_points[i], _weight[i]));
366  }
367  };
368 
369 #ifndef DOXYGEN
370  extern template Jacobi2QuadratureRule1D<float>::Jacobi2QuadratureRule1D(int);
371  extern template Jacobi2QuadratureRule1D<double>::Jacobi2QuadratureRule1D(int);
372 #endif // !DOXYGEN
373 
374 } // namespace Dune
375 
376 #define DUNE_INCLUDING_IMPLEMENTATION
377 #include "quadraturerules/jacobi_2_0_imp.hh"
378 
379 namespace Dune {
380 
382  template<typename ct,
383  bool fundamental = std::numeric_limits<ct>::is_specialized>
385  template<typename ct>
387  static void init(int p,
388  std::vector< FieldVector<ct, 1> > & _points,
389  std::vector< ct > & _weight,
390  int & delivered_order);
391  };
392  template<typename ct>
394  static void init(int p,
395  std::vector< FieldVector<ct, 1> > & _points,
396  std::vector< ct > & _weight,
397  int & delivered_order);
398  };
399 
403  template<typename ct>
405  public QuadratureRule<ct,1>
406  {
407  public:
409  enum { dim=1 };
410 
412  enum { highest_order=61 };
413 
415  private:
416  friend class QuadratureRuleFactory<ct,dim>;
418  : QuadratureRule<ct,1>(GeometryType(GeometryType::cube, 1))
419  {
421  std::vector< FieldVector<ct, dim> > _points;
422  std::vector< ct > _weight;
423 
424  int delivered_order;
425 
427  (p, _points, _weight, delivered_order);
428 
429  this->delivered_order = delivered_order;
430  assert(_points.size() == _weight.size());
431  for (size_t i = 0; i < _points.size(); i++)
432  this->push_back(QuadraturePoint<ct,dim>(_points[i], _weight[i]));
433  }
434  };
435 
436 #ifndef DOXYGEN
437  extern template GaussLobattoQuadratureRule1D<float>::GaussLobattoQuadratureRule1D(int);
438  extern template GaussLobattoQuadratureRule1D<double>::GaussLobattoQuadratureRule1D(int);
439 #endif // !DOXYGEN
440 
441 } // namespace Dune
442 
443 #define DUNE_INCLUDING_IMPLEMENTATION
444 #include "quadraturerules/gausslobatto_imp.hh"
445 
446 #include "quadraturerules/tensorproductquadrature.hh"
447 
448 #include "quadraturerules/simplexquadrature.hh"
449 
450 namespace Dune {
451 
452  /***********************************
453  * quadrature for Prism
454  **********************************/
455 
457  template<int dim>
459 
461  template<>
463  {
464  public:
465  enum { MAXP=6};
466  enum { highest_order=2 };
467 
470  {
471  int m = 0;
472  O[m] = 0;
473 
474  // polynom degree 0 ???
475  m = 6;
476  G[m][0][0] = 0.0;
477  G[m][0][1] = 0.0;
478  G[m][0][2] = 0.0;
479 
480  G[m][1][0] = 1.0;
481  G[m][1][1] = 0.0;
482  G[m][1][2] = 0.0;
483 
484  G[m][2][0] = 0.0;
485  G[m][2][1] = 1.0;
486  G[m][2][2] = 0.0;
487 
488  G[m][3][0] = 0.0;
489  G[m][3][1] = 0.0;
490  G[m][3][2] = 1.0;
491 
492  G[m][4][0] = 1.0;
493  G[m][4][1] = 0.0;
494  G[m][4][2] = 1.0;
495 
496  G[m][5][0] = 0.0;
497  G[m][5][1] = 0.1;
498  G[m][5][2] = 1.0;
499 
500  W[m][0] = 0.16666666666666666 / 2.0;
501  W[m][1] = 0.16666666666666666 / 2.0;
502  W[m][2] = 0.16666666666666666 / 2.0;
503  W[m][3] = 0.16666666666666666 / 2.0;
504  W[m][4] = 0.16666666666666666 / 2.0;
505  W[m][5] = 0.16666666666666666 / 2.0;
506 
507  O[m] = 0; // verify ????????
508 
509 
510  // polynom degree 2 ???
511  m = 6;
512  G[m][0][0] =0.66666666666666666 ;
513  G[m][0][1] =0.16666666666666666 ;
514  G[m][0][2] =0.211324865405187 ;
515 
516  G[m][1][0] = 0.16666666666666666;
517  G[m][1][1] =0.66666666666666666 ;
518  G[m][1][2] = 0.211324865405187;
519 
520  G[m][2][0] = 0.16666666666666666;
521  G[m][2][1] = 0.16666666666666666;
522  G[m][2][2] = 0.211324865405187;
523 
524  G[m][3][0] = 0.66666666666666666;
525  G[m][3][1] = 0.16666666666666666;
526  G[m][3][2] = 0.788675134594813;
527 
528  G[m][4][0] = 0.16666666666666666;
529  G[m][4][1] = 0.66666666666666666;
530  G[m][4][2] = 0.788675134594813;
531 
532  G[m][5][0] = 0.16666666666666666;
533  G[m][5][1] = 0.16666666666666666;
534  G[m][5][2] = 0.788675134594813;
535 
536  W[m][0] = 0.16666666666666666 / 2.0;
537  W[m][1] = 0.16666666666666666 / 2.0;
538  W[m][2] = 0.16666666666666666 / 2.0;
539  W[m][3] = 0.16666666666666666 / 2.0;
540  W[m][4] = 0.16666666666666666 / 2.0;
541  W[m][5] = 0.16666666666666666 / 2.0;
542 
543  O[m] = 2; // verify ????????
544 
545  }
546 
548  FieldVector<double, 3> point(int m, int i)
549  {
550  return G[m][i];
551  }
552 
554  double weight (int m, int i)
555  {
556  return W[m][i];
557  }
558 
560  int order (int m)
561  {
562  return O[m];
563  }
564 
565  private:
566  FieldVector<double, 3> G[MAXP+1][MAXP]; //positions
567 
568  double W[MAXP+1][MAXP]; // weights associated with points
569  int O[MAXP+1]; // order of the rule
570  };
571 
572 
576  template<int dim>
579  };
580 
584  template<>
587  };
588 
592  template<typename ct, int dim>
594 
598  template<typename ct>
600  {
601  public:
602 
604  enum { d = 3 };
605 
607  enum { highest_order = 2 };
608 
610  private:
613  {
614  int m=6;
615  this->delivered_order = PrismQuadraturePointsSingleton<3>::prqp.order(m);
616  for(int i=0; i<m; ++i)
617  {
618  FieldVector<ct,3> local;
619  for (int k=0; k<d; k++)
620  local[k] = PrismQuadraturePointsSingleton<3>::prqp.point(m,i)[k];
621  double weight =
623  // put in container
624  this->push_back(QuadraturePoint<ct,d>(local,weight));
625  }
626  }
627  };
628 
635  template<typename ctype, int dim>
636  class QuadratureRuleFactory {
637  private:
639  static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
640  {
641  return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
642  }
643  };
644 
645  template<typename ctype>
647  private:
648  enum { dim = 0 };
650  static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
651  {
652  if (t.isVertex())
653  {
654  return PointQuadratureRule<ctype>();
655  }
656  DUNE_THROW(Exception, "Unknown GeometryType");
657  }
658  };
659 
660  template<typename ctype>
662  private:
663  enum { dim = 1 };
665  static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
666  {
667  if (t.isLine())
668  {
669  switch (qt) {
678  default :
679  DUNE_THROW(Exception, "Unknown QuadratureType");
680  }
681  }
682  DUNE_THROW(Exception, "Unknown GeometryType");
683  }
684  };
685 
686  template<typename ctype>
688  private:
689  enum { dim = 2 };
691  static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
692  {
693  if (t.isSimplex() && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
694  {
695  return SimplexQuadratureRule<ctype,dim>(p);
696  }
697  return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
698  }
699  };
700 
701  template<typename ctype>
703  private:
704  enum { dim = 3 };
706  static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
707  {
708  if (t.isSimplex() && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
709  {
710  return SimplexQuadratureRule<ctype,dim>(p);
711  }
712  if (t.isPrism() && p <= PrismQuadratureRule<ctype,dim>::highest_order)
713  {
715  }
716  return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
717  }
718  };
719 
720 } // end namespace
721 
722 #endif // DUNE_GEOMETRY_QUADRATURERULES_HH
Definition: quadraturerules.hh:78
virtual GeometryType type() const
return type of element
Definition: quadraturerules.hh:116
ct CoordType
The type used for coordinates.
Definition: quadraturerules.hh:110
Definition: quadraturerules.hh:193
Quadrature rules for prisms.
Definition: quadraturerules.hh:593
double weight(int m, int i)
Definition: quadraturerules.hh:554
Exception thrown if a desired QuadratureRule is not available, because the requested order is to high...
Definition: quadraturerules.hh:31
Definition: quadraturerules.hh:45
Definition: quadraturerules.hh:77
A unique label for each type of element that can occur in a grid.
Definition: quadraturerules.hh:83
Dune::FieldVector< ct, dim > Vector
Definition: quadraturerules.hh:47
Factory class for creation of quadrature rules, depending on GeometryType, order and QuadratureType...
Definition: quadraturerules.hh:130
~Jacobi2QuadratureRule1D()
Definition: quadraturerules.hh:347
Definition: quadraturerules.hh:384
bool isSimplex() const
Return true if entity is a simplex of any dimension.
Definition: type.hh:307
Definition: quadraturerules.hh:458
bool isVertex() const
Return true if entity is a vertex.
Definition: type.hh:267
QuadratureRule(GeometryType t, int order)
Constructor for a given geometry type and a given quadrature order.
Definition: quadraturerules.hh:104
Definition: quadraturerules.hh:687
int delivered_order
Definition: quadraturerules.hh:125
~PrismQuadratureRule()
Definition: quadraturerules.hh:609
Definition: quadraturerules.hh:217
Jacobi-Gauss quadrature for alpha=2, beta=0.
Definition: quadraturerules.hh:404
QuadratureRule()
Default constructor.
Definition: quadraturerules.hh:98
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:173
std::vector< QuadraturePoint< ct, dim > >::const_iterator iterator
Definition: quadraturerules.hh:121
static const QuadratureRule & rule(const GeometryType::BasicType t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:178
Definition: quadraturerules.hh:107
const Vector & position() const
return local coordinates of integration point i
Definition: quadraturerules.hh:56
Definition: quadraturerules.hh:42
Enum
Definition: quadraturerules.hh:76
const ct & weight() const
return weight associated with integration point i
Definition: quadraturerules.hh:62
static PrismQuadraturePoints< 3 > prqp
Definition: quadraturerules.hh:586
Singleton holding the Prism Quadrature points.
Definition: quadraturerules.hh:585
Definition: quadraturerules.hh:81
Definition: quadraturerules.hh:317
static PrismQuadraturePoints< 3 > prqp
Definition: quadraturerules.hh:578
bool isLine() const
Return true if entity is a line segment.
Definition: type.hh:272
Jacobi-Gauss quadrature for alpha=2, beta=0.
Definition: quadraturerules.hh:337
unsigned int id() const
Return the topology id the type.
Definition: type.hh:327
Definition: quadraturerules.hh:216
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:93
QuadratureRule(GeometryType t)
Constructor for a given geometry type. Leaves the quadrature order invalid.
Definition: quadraturerules.hh:101
Definition: quadraturerules.hh:251
bool isPrism() const
Return true if entity is a prism.
Definition: type.hh:297
~GaussLobattoQuadratureRule1D()
Definition: quadraturerules.hh:414
Definition: quadraturerules.hh:661
GeometryType geometry_type
Definition: quadraturerules.hh:124
Definition: quadraturerules.hh:85
Singleton holding the Prism Quadrature points.
Definition: quadraturerules.hh:577
virtual ~QuadratureRule()
Definition: quadraturerules.hh:117
ct CoordType
Definition: quadraturerules.hh:43
ct Field
Definition: quadraturerules.hh:46
A container for all quadrature rules of dimension dim
Definition: quadraturerules.hh:136
Quadrature rules for prisms.
Definition: quadraturerules.hh:599
FieldVector< double, 3 > point(int m, int i)
Definition: quadraturerules.hh:548
virtual int order() const
return order
Definition: quadraturerules.hh:113
FieldVector< ct, dim > local
Definition: quadraturerules.hh:68
PrismQuadraturePoints()
initialize quadrature points on the interval for all orders
Definition: quadraturerules.hh:469
~Jacobi1QuadratureRule1D()
Definition: quadraturerules.hh:281
int order(int m)
Definition: quadraturerules.hh:560
Single evaluation point in a quadrature rule.
Definition: quadraturerules.hh:39
Jacobi-Gauss quadrature for alpha=1, beta=0.
Definition: quadraturerules.hh:271
~GaussQuadratureRule1D()
Definition: quadraturerules.hh:219
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:24
Definition: quadraturerules.hh:80
ct weight_
Definition: quadraturerules.hh:69
Definition: quadraturerules.hh:702
QuadraturePoint(const Vector &x, ct w)
set up quadrature of given order in d dimensions
Definition: quadraturerules.hh:50
BasicType
Each entity can be tagged by one of these basic types plus its space dimension.
Definition: type.hh:29
Definition: quadraturerules.hh:646
Gauss quadrature rule in 1D.
Definition: quadraturerules.hh:211
Definition: quadraturerules.hh:82
Definition: quadraturerules.hh:462