Dune Core Modules (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
17#include <dune/common/stdthread.hh>
19
20#include <dune/geometry/type.hh>
22
28namespace Dune {
29
35
41 template<typename ct, int dim>
43 public:
45 enum { dimension = dim };
46
48 typedef ct Field;
49
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:
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:
159 QuadratureRule() : delivered_order(-1) {}
160
161 protected:
163 QuadratureRule(GeometryType t) : geometry_type(t), delivered_order(-1) {}
164
166 QuadratureRule(GeometryType t, int order) : geometry_type(t), delivered_order(order) {}
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:
186 GeometryType geometry_type;
187 int delivered_order;
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
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
315namespace 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 }
331 static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
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 {
347 }
348 DUNE_THROW(Exception, "Unknown GeometryType");
349 }
350 static QuadratureRule<ctype, dim> rule(const GeometryType& t, int , QuadratureType::Enum)
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 }
390 static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
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 }
431 static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
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
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:123
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 unsigned int id() const
Return the topology id of the type.
Definition: type.hh:374
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
Default exception for dummy implementations.
Definition: exceptions.hh:261
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
const Vector & position() const
return local coordinates of integration point i
Definition: quadraturerules.hh:60
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 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
Factory class for creation of quadrature rules, depending on GeometryType, order and QuadratureType.
Definition: quadraturerules.hh:324
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:152
virtual GeometryType type() const
return type of element
Definition: quadraturerules.hh:178
QuadratureRule(GeometryType t, int order)
Constructor for a given geometry type and a given quadrature order.
Definition: quadraturerules.hh:166
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
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
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::BasicType t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:286
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
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:156
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
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
@ 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
Dune namespace.
Definition: alignedallocator.hh:11
Standard Dune debug streams.
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 of macros controlling symbol visibility at the ABI level.
#define DUNE_EXPORT
Export a symbol as part of the public ABI.
Definition: visibility.hh:18
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)