dune-geometry  2.8.0
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 <algorithm>
8 #include <iostream>
9 #include <limits>
10 #include <mutex>
11 #include <utility>
12 #include <vector>
13 
14 #include <dune/common/fvector.hh>
15 #include <dune/common/exceptions.hh>
16 #include <dune/common/stdstreams.hh>
17 #include <dune/common/stdthread.hh>
18 #include <dune/common/visibility.hh>
19 
20 #include <dune/geometry/type.hh>
22 
28 namespace Dune {
29 
34  class QuadratureOrderOutOfRange : public NotImplemented {};
35 
41  template<typename ct, int dim>
43  public:
45  enum { dimension = dim };
46 
48  typedef ct Field;
49 
51  typedef Dune::FieldVector<ct,dim> Vector;
52 
54  QuadraturePoint (const Vector& x, ct w) : local(x)
55  {
56  weight_ = w;
57  }
58 
60  const Vector& position () const
61  {
62  return local;
63  }
64 
66  const ct &weight () const
67  {
68  return weight_;
69  }
70 
71  protected:
72  FieldVector<ct, dim> local;
73  ct weight_;
74  };
75 
79  namespace QuadratureType {
80  enum Enum {
91 
98 
105 
118 
126 
134 
143  size
144  };
145  }
146 
150  template<typename ct, int dim>
151  class QuadratureRule : public std::vector<QuadraturePoint<ct,dim> >
152  {
153  public:
160 
161  protected:
164 
167  public:
169  enum { d=dim };
170 
172  typedef ct CoordType;
173 
175  virtual int order () const { return delivered_order; }
176 
178  virtual GeometryType type () const { return geometry_type; }
179  virtual ~QuadratureRule(){}
180 
183  typedef typename std::vector<QuadraturePoint<ct,dim> >::const_iterator iterator;
184 
185  protected:
188  };
189 
190  // Forward declaration of the factory class,
191  // needed internally by the QuadratureRules container class.
192  template<typename ctype, int dim> class QuadratureRuleFactory;
193 
197  template<typename ctype, int dim>
199 
204  static void initQuadratureRule(QuadratureRule *qr, QuadratureType::Enum qt,
205  const GeometryType &t, int p)
206  {
208  }
209 
210  typedef std::vector<std::pair<std::once_flag, QuadratureRule> >
211  QuadratureOrderVector; // indexed by quadrature order
214  static void initQuadratureOrderVector(QuadratureOrderVector *qov,
216  const GeometryType &t)
217  {
218  if(dim == 0)
219  // we only need one quadrature rule for points, not maxint
220  *qov = QuadratureOrderVector(1);
221  else
222  *qov = QuadratureOrderVector(QuadratureRuleFactory<ctype,dim>::maxOrder(t,qt)+1);
223  }
224 
225  typedef std::vector<std::pair<std::once_flag, QuadratureOrderVector> >
226  GeometryTypeVector; // indexed by geometry type
229  static void initGeometryTypeVector(GeometryTypeVector *gtv)
230  {
231  *gtv = GeometryTypeVector(LocalGeometryTypeIndex::size(dim));
232  }
233 
235  DUNE_EXPORT const QuadratureRule& _rule(const GeometryType& t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
236  {
237  assert(t.dim()==dim);
238 
239  DUNE_ASSERT_CALL_ONCE();
240 
241  static std::vector<std::pair< // indexed by quadrature type
242  std::once_flag,
243  GeometryTypeVector
244  > > quadratureCache(QuadratureType::size);
245 
246  auto & quadratureTypeLevel = quadratureCache[qt];
247  std::call_once(quadratureTypeLevel.first, initGeometryTypeVector,
248  &quadratureTypeLevel.second);
249 
250  auto & geometryTypeLevel =
251  quadratureTypeLevel.second[LocalGeometryTypeIndex::index(t)];
252  std::call_once(geometryTypeLevel.first, initQuadratureOrderVector,
253  &geometryTypeLevel.second, qt, t);
254 
255  // we only have one quadrature rule for points
256  auto & quadratureOrderLevel = geometryTypeLevel.second[dim == 0 ? 0 : p];
257  std::call_once(quadratureOrderLevel.first, initQuadratureRule,
258  &quadratureOrderLevel.second, qt, t, p);
259 
260  return quadratureOrderLevel.second;
261  }
263  DUNE_EXPORT static QuadratureRules& instance()
264  {
265  static QuadratureRules instance;
266  return instance;
267  }
269  QuadratureRules () {}
270  public:
272  static unsigned
275  {
277  }
278 
281  {
282  return instance()._rule(t,p,qt);
283  }
284 
287  {
288  GeometryType gt(t,dim);
289  return instance()._rule(gt,p,qt);
290  }
291  };
292 
293 } // end namespace Dune
294 
295 #define DUNE_INCLUDING_IMPLEMENTATION
296 
297 // 0d rules
298 #include "quadraturerules/pointquadrature.hh"
299 // 1d rules
300 #include "quadraturerules/gausslobattoquadrature.hh"
301 #include "quadraturerules/gaussquadrature.hh"
302 #include "quadraturerules/gaussradauleftquadrature.hh"
303 #include "quadraturerules/gaussradaurightquadrature.hh"
304 #include "quadraturerules/jacobi1quadrature.hh"
305 #include "quadraturerules/jacobi2quadrature.hh"
306 #include "quadraturerules/jacobiNquadrature.hh"
307 // 3d rules
308 #include "quadraturerules/prismquadrature.hh"
309 // general rules
310 #include "quadraturerules/simplexquadrature.hh"
311 #include "quadraturerules/tensorproductquadrature.hh"
312 
313 #undef DUNE_INCLUDING_IMPLEMENTATION
314 
315 namespace Dune {
316 
323  template<typename ctype, int dim>
325  private:
326  friend class QuadratureRules<ctype, dim>;
327  static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
328  {
329  return TensorProductQuadratureRule<ctype,dim>::maxOrder(t.id(), qt);
330  }
332  {
333  return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
334  }
335  };
336 
337  template<typename ctype>
338  class QuadratureRuleFactory<ctype, 0> {
339  private:
340  enum { dim = 0 };
341  friend class QuadratureRules<ctype, dim>;
342  static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum)
343  {
344  if (t.isVertex())
345  {
346  return std::numeric_limits<int>::max();
347  }
348  DUNE_THROW(Exception, "Unknown GeometryType");
349  }
351  {
352  if (t.isVertex())
353  {
354  return PointQuadratureRule<ctype>();
355  }
356  DUNE_THROW(Exception, "Unknown GeometryType");
357  }
358  };
359 
360  template<typename ctype>
361  class QuadratureRuleFactory<ctype, 1> {
362  private:
363  enum { dim = 1 };
364  friend class QuadratureRules<ctype, dim>;
365  static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
366  {
367  if (t.isLine())
368  {
369  switch (qt) {
371  return GaussQuadratureRule1D<ctype>::highest_order;
373  return Jacobi1QuadratureRule1D<ctype>::highest_order;
375  return Jacobi2QuadratureRule1D<ctype>::highest_order;
377  return GaussLobattoQuadratureRule1D<ctype>::highest_order;
379  return JacobiNQuadratureRule1D<ctype>::maxOrder();
381  return GaussRadauLeftQuadratureRule1D<ctype>::highest_order;
383  return GaussRadauRightQuadratureRule1D<ctype>::highest_order;
384  default :
385  DUNE_THROW(Exception, "Unknown QuadratureType");
386  }
387  }
388  DUNE_THROW(Exception, "Unknown GeometryType");
389  }
391  {
392  if (t.isLine())
393  {
394  switch (qt) {
396  return GaussQuadratureRule1D<ctype>(p);
398  return Jacobi1QuadratureRule1D<ctype>(p);
400  return Jacobi2QuadratureRule1D<ctype>(p);
402  return GaussLobattoQuadratureRule1D<ctype>(p);
404  return JacobiNQuadratureRule1D<ctype>(p);
406  return GaussRadauLeftQuadratureRule1D<ctype>(p);
408  return GaussRadauRightQuadratureRule1D<ctype>(p);
409  default :
410  DUNE_THROW(Exception, "Unknown QuadratureType");
411  }
412  }
413  DUNE_THROW(Exception, "Unknown GeometryType");
414  }
415  };
416 
417  template<typename ctype>
418  class QuadratureRuleFactory<ctype, 2> {
419  private:
420  enum { dim = 2 };
421  friend class QuadratureRules<ctype, dim>;
422  static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
423  {
424  unsigned order =
425  TensorProductQuadratureRule<ctype,dim>::maxOrder(t.id(), qt);
426  if (t.isSimplex())
427  order = std::max
428  (order, unsigned(SimplexQuadratureRule<ctype,dim>::highest_order));
429  return order;
430  }
432  {
433  if (t.isSimplex()
435  && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
436  {
437  return SimplexQuadratureRule<ctype,dim>(p);
438  }
439  return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
440  }
441  };
442 
443  template<typename ctype>
444  class QuadratureRuleFactory<ctype, 3> {
445  private:
446  enum { dim = 3 };
447  friend class QuadratureRules<ctype, dim>;
448  static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
449  {
450  unsigned order =
451  TensorProductQuadratureRule<ctype,dim>::maxOrder(t.id(), qt);
452  if (t.isSimplex())
453  order = std::max
454  (order, unsigned(SimplexQuadratureRule<ctype,dim>::highest_order));
455  if (t.isPrism())
456  order = std::max
457  (order, unsigned(PrismQuadratureRule<ctype,dim>::highest_order));
458  return order;
459  }
460  static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
461  {
462 
463  if (t.isSimplex()
465  && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
466  {
467  return SimplexQuadratureRule<ctype,dim>(p);
468  }
469  if (t.isPrism()
471  && p <= PrismQuadratureRule<ctype,dim>::highest_order)
472  {
473  return PrismQuadratureRule<ctype,dim>(p);
474  }
475  return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
476  }
477  };
478 
479 #ifndef DUNE_NO_EXTERN_QUADRATURERULES
480  extern template class GaussLobattoQuadratureRule<double, 1>;
481  extern template class GaussQuadratureRule<double, 1>;
482  extern template class GaussRadauLeftQuadratureRule<double, 1>;
483  extern template class GaussRadauRightQuadratureRule<double, 1>;
484  extern template class Jacobi1QuadratureRule<double, 1>;
485  extern template class Jacobi2QuadratureRule<double, 1>;
486  extern template class JacobiNQuadratureRule<double, 1>;
487  extern template class PrismQuadratureRule<double, 3>;
488  extern template class SimplexQuadratureRule<double, 2>;
489  extern template class SimplexQuadratureRule<double, 3>;
490 #endif // !DUNE_NO_EXTERN_QUADRATURERULES
491 
492 } // end namespace
493 
494 #endif // DUNE_GEOMETRY_QUADRATURERULES_HH
A unique label for each type of element that can occur in a grid.
Helper classes to provide indices for geometrytypes for use in a vector.
Definition: affinegeometry.hh:19
Enum
Definition: quadraturerules.hh:80
@ GaussJacobi_n_0
Gauss-Legendre rules with .
Definition: quadraturerules.hh:117
@ GaussJacobi_2_0
Gauss-Legendre rules with .
Definition: quadraturerules.hh:104
@ GaussRadauRight
Gauss-Radau rules including the right endpoint.
Definition: quadraturerules.hh:142
@ GaussJacobi_1_0
Gauss-Jacobi rules with .
Definition: quadraturerules.hh:97
@ size
Definition: quadraturerules.hh:143
@ GaussLobatto
Gauss-Lobatto rules.
Definition: quadraturerules.hh:125
@ GaussRadauLeft
Gauss-Radau rules including the left endpoint.
Definition: quadraturerules.hh:133
@ GaussLegendre
Gauss-Legendre rules (default)
Definition: quadraturerules.hh:90
Exception thrown if a desired QuadratureRule is not available, because the requested order is to high...
Definition: quadraturerules.hh:34
Single evaluation point in a quadrature rule.
Definition: quadraturerules.hh:42
Dune::FieldVector< ct, dim > Vector
Type used for the position of a quadrature point.
Definition: quadraturerules.hh:51
ct Field
Number type used for coordinates and quadrature weights.
Definition: quadraturerules.hh:48
const Vector & position() const
return local coordinates of integration point i
Definition: quadraturerules.hh:60
@ dimension
Definition: quadraturerules.hh:45
ct weight_
Definition: quadraturerules.hh:73
const ct & weight() const
return weight associated with integration point i
Definition: quadraturerules.hh:66
QuadraturePoint(const Vector &x, ct w)
set up quadrature of given order in d dimensions
Definition: quadraturerules.hh:54
FieldVector< ct, dim > local
Definition: quadraturerules.hh:72
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:152
virtual ~QuadratureRule()
Definition: quadraturerules.hh:179
virtual GeometryType type() const
return type of element
Definition: quadraturerules.hh:178
int delivered_order
Definition: quadraturerules.hh:187
QuadratureRule(GeometryType t, int order)
Constructor for a given geometry type and a given quadrature order.
Definition: quadraturerules.hh:166
GeometryType geometry_type
Definition: quadraturerules.hh:186
ct CoordType
The type used for coordinates.
Definition: quadraturerules.hh:172
QuadratureRule()
Default constructor.
Definition: quadraturerules.hh:159
virtual int order() const
return order
Definition: quadraturerules.hh:175
@ d
Definition: quadraturerules.hh:169
QuadratureRule(GeometryType t)
Constructor for a given geometry type. Leaves the quadrature order invalid
Definition: quadraturerules.hh:163
std::vector< QuadraturePoint< ct, dim > >::const_iterator iterator
Definition: quadraturerules.hh:183
Factory class for creation of quadrature rules, depending on GeometryType, order and QuadratureType.
Definition: quadraturerules.hh:324
A container for all quadrature rules of dimension dim
Definition: quadraturerules.hh:198
static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
maximum quadrature order for given geometry type and quadrature type
Definition: quadraturerules.hh:273
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:280
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:286
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:123
constexpr bool isPrism() const
Return true if entity is a prism.
Definition: type.hh:318
constexpr bool isVertex() const
Return true if entity is a vertex.
Definition: type.hh:288
constexpr unsigned int dim() const
Return dimension of the type.
Definition: type.hh:369
BasicType
Each entity can be tagged by one of these basic types plus its space dimension.
Definition: type.hh:129
constexpr bool isLine() const
Return true if entity is a line segment.
Definition: type.hh:293
constexpr unsigned int id() const
Return the topology id of the type.
Definition: type.hh:374
constexpr bool isSimplex() const
Return true if entity is a simplex of any dimension.
Definition: type.hh:328
static constexpr std::size_t size(std::size_t dim)
Compute total number of geometry types for the given dimension.
Definition: typeindex.hh:59
static constexpr std::size_t index(const GeometryType &gt)
Compute the index for the given geometry type within its dimension.
Definition: typeindex.hh:71