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>
193  struct GaussQuadratureInitHelper;
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 
229  GaussQuadratureInitHelper<ct>::init
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>
251  struct Jacobi1QuadratureInitHelper;
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 
293  Jacobi1QuadratureInitHelper<ct>::init
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>
317  struct Jacobi2QuadratureInitHelper;
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 
359  Jacobi2QuadratureInitHelper<ct>::init
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>
384  struct GaussLobattoQuadratureInitHelper;
385  template<typename ct>
386  struct GaussLobattoQuadratureInitHelper<ct, true> {
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>
393  struct GaussLobattoQuadratureInitHelper<ct, false> {
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 
426  GaussLobattoQuadratureInitHelper<ct>::init
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>
458  class PrismQuadraturePoints;
459 
461  template<>
463  {
464  public:
465  enum { MAXP=6};
466  enum { highest_order=2 };
467 
469  PrismQuadraturePoints ()
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>
593  class PrismQuadratureRule;
594 
598  template<typename ct>
600  {
601  public:
602 
604  enum { d = 3 };
605 
607  enum { highest_order = 2 };
608 
610  private:
612  PrismQuadratureRule(int p) : QuadratureRule<ct,3>(GeometryType(GeometryType::prism, d))
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>
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  {
714  return PrismQuadratureRule<ctype,dim>(p);
715  }
716  return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
717  }
718  };
719 
720 } // end namespace
721 
722 #endif // DUNE_GEOMETRY_QUADRATURERULES_HH