DUNE PDELab (git)

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// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5
6#ifndef DUNE_GEOMETRY_QUADRATURERULES_HH
7#define DUNE_GEOMETRY_QUADRATURERULES_HH
8
9#include <algorithm>
10#include <iostream>
11#include <limits>
12#include <mutex>
13#include <utility>
14#include <vector>
15
19#include <dune/common/stdthread.hh>
21
22#include <dune/geometry/type.hh>
24
31namespace Dune {
32 // forward declaration
33 template<typename ct, int dim>
34 class QuadraturePoint;
35}
36
37// class specialization of standard classes that allow to use structured bindings on QuadraturePoint
38namespace std {
39 template<typename ct, int dim>
40 struct tuple_size<Dune::QuadraturePoint<ct,dim>> : public std::integral_constant<std::size_t,2> {};
41
42 template<typename ct, int dim>
43 struct tuple_element<0, Dune::QuadraturePoint<ct,dim>> { using type = Dune::FieldVector<ct, dim>; };
44
45 template<typename ct, int dim>
46 struct tuple_element<1, Dune::QuadraturePoint<ct,dim>> { using type = ct; };
47}
48
49namespace Dune {
50
56
65 template<typename ct, int dim>
67 public:
69 constexpr static int dimension = dim;
70
72 typedef ct Field;
73
76
78 QuadraturePoint (const Vector& x, ct w) : local(x), weight_(w)
79 {}
80
82 const Vector& position () const
83 {
84 return local;
85 }
86
88 const ct &weight () const
89 {
90 return weight_;
91 }
92
111 template<std::size_t index, std::enable_if_t<(index<=1), int> = 0>
112 std::tuple_element_t<index, QuadraturePoint<ct, dim>> get() const
113 {
114 if constexpr (index == 0) {
115 return local;
116 }
117 else {
118 return weight_;
119 }
120 }
121
122 protected:
123 FieldVector<ct, dim> local;
124 ct weight_;
125 };
126
130 namespace QuadratureType {
131 enum Enum {
141 GaussLegendre = 0,
142
148 GaussJacobi_1_0 = 1,
149
155 GaussJacobi_2_0 = 2,
156
168 GaussJacobi_n_0 = 3,
169
176 GaussLobatto = 4,
177
184 GaussRadauLeft = 5,
185
193 GaussRadauRight = 6,
194 size
195 };
196 }
197
212 template<typename ct, int dim>
213 class QuadratureRule : public std::vector<QuadraturePoint<ct,dim> >
214 {
215 public:
221 QuadratureRule() : delivered_order(-1) {}
222
223 protected:
225 QuadratureRule(GeometryType t) : geometry_type(t), delivered_order(-1) {}
226
228 QuadratureRule(GeometryType t, int order) : geometry_type(t), delivered_order(order) {}
229 public:
231 constexpr static int d = dim;
232
234 typedef ct CoordType;
235
237 virtual int order () const { return delivered_order; }
238
240 virtual GeometryType type () const { return geometry_type; }
241 virtual ~QuadratureRule(){}
242
245 typedef typename std::vector<QuadraturePoint<ct,dim> >::const_iterator iterator;
246
247 protected:
248 GeometryType geometry_type;
249 int delivered_order;
250 };
251
252 // Forward declaration of the factory class,
253 // needed internally by the QuadratureRules container class.
254 template<typename ctype, int dim> class QuadratureRuleFactory;
255
259 template<typename ctype, int dim>
261
264
265 // indexed by quadrature order
266 using QuadratureOrderVector = std::vector<std::pair<std::once_flag, QuadratureRule> >;
267
268 // indexed by geometry type
269 using GeometryTypeVector = std::vector<std::pair<std::once_flag, QuadratureOrderVector> >;
270
271 // indexed by quadrature type enum
272 using QuadratureCacheVector = std::vector<std::pair<std::once_flag, GeometryTypeVector> >;
273
276 {
277 assert(t.dim()==dim);
278
279 DUNE_ASSERT_CALL_ONCE();
280
281 static QuadratureCacheVector quadratureCache(QuadratureType::size);
282
283 auto& [ onceFlagQuadratureType, geometryTypes ] = quadratureCache[qt];
284 // initialize geometry types for this quadrature type once
285 std::call_once(onceFlagQuadratureType, [&types = geometryTypes]{
286 types = GeometryTypeVector(LocalGeometryTypeIndex::size(dim));
287 });
288
289 auto& [ onceFlagGeometryType, quadratureOrders ] = geometryTypes[LocalGeometryTypeIndex::index(t)];
290 // initialize quadrature orders for this geometry type and quadrature type once
291 std::call_once(onceFlagGeometryType, [&, &orders = quadratureOrders]{
292 // we only need one quadrature rule for points, not maxint
293 const auto numRules = dim == 0 ? 1 : QuadratureRuleFactory<ctype,dim>::maxOrder(t, qt)+1;
294 orders = QuadratureOrderVector(numRules);
295 });
296
297 // we only have one quadrature rule for points
298 auto& [ onceFlagQuadratureOrder, quadratureRule ] = quadratureOrders[dim == 0 ? 0 : p];
299 // initialize quadrature rule once
300 std::call_once(onceFlagQuadratureOrder, [&, &rule = quadratureRule]{
302 });
303
304 return quadratureRule;
305 }
306
308 DUNE_EXPORT static QuadratureRules& instance()
309 {
310 static QuadratureRules instance;
311 return instance;
312 }
313
315 QuadratureRules () = default;
316 public:
318 static unsigned
321 {
323 }
324
327 {
328 return instance()._rule(t,p,qt);
329 }
330
333 {
334 GeometryType gt(t,dim);
335 return instance()._rule(gt,p,qt);
336 }
337 };
338
339} // end namespace Dune
340
341#define DUNE_INCLUDING_IMPLEMENTATION
342
343// 0d rules
344#include "quadraturerules/pointquadrature.hh"
345// 1d rules
346#include "quadraturerules/gausslobattoquadrature.hh"
347#include "quadraturerules/gaussquadrature.hh"
348#include "quadraturerules/gaussradauleftquadrature.hh"
349#include "quadraturerules/gaussradaurightquadrature.hh"
350#include "quadraturerules/jacobi1quadrature.hh"
351#include "quadraturerules/jacobi2quadrature.hh"
352#include "quadraturerules/jacobiNquadrature.hh"
353// 3d rules
354#include "quadraturerules/prismquadrature.hh"
355// general rules
356#include "quadraturerules/simplexquadrature.hh"
357#include "quadraturerules/tensorproductquadrature.hh"
358
359#undef DUNE_INCLUDING_IMPLEMENTATION
360
361namespace Dune {
362
369 template<typename ctype, int dim>
371 private:
372 friend class QuadratureRules<ctype, dim>;
373 static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
374 {
375 return TensorProductQuadratureRule<ctype,dim>::maxOrder(t.id(), qt);
376 }
377 static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
378 {
379 return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
380 }
381 };
382
383 template<typename ctype>
384 class QuadratureRuleFactory<ctype, 0> {
385 private:
386 constexpr static int dim = 0;
387 friend class QuadratureRules<ctype, dim>;
388 static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum)
389 {
390 if (t.isVertex())
391 {
393 }
394 DUNE_THROW(Exception, "Unknown GeometryType");
395 }
396 static QuadratureRule<ctype, dim> rule(const GeometryType& t, int , QuadratureType::Enum)
397 {
398 if (t.isVertex())
399 {
400 return PointQuadratureRule<ctype>();
401 }
402 DUNE_THROW(Exception, "Unknown GeometryType");
403 }
404 };
405
406 template<typename ctype>
407 class QuadratureRuleFactory<ctype, 1> {
408 private:
409 constexpr static int dim = 1;
410 friend class QuadratureRules<ctype, dim>;
411 static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
412 {
413 if (t.isLine())
414 {
415 switch (qt) {
417 return GaussQuadratureRule1D<ctype>::highest_order;
419 return Jacobi1QuadratureRule1D<ctype>::highest_order;
421 return Jacobi2QuadratureRule1D<ctype>::highest_order;
423 return GaussLobattoQuadratureRule1D<ctype>::highest_order;
425 return JacobiNQuadratureRule1D<ctype>::maxOrder();
427 return GaussRadauLeftQuadratureRule1D<ctype>::highest_order;
429 return GaussRadauRightQuadratureRule1D<ctype>::highest_order;
430 default :
431 DUNE_THROW(Exception, "Unknown QuadratureType");
432 }
433 }
434 DUNE_THROW(Exception, "Unknown GeometryType");
435 }
436 static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
437 {
438 if (t.isLine())
439 {
440 switch (qt) {
442 return GaussQuadratureRule1D<ctype>(p);
444 return Jacobi1QuadratureRule1D<ctype>(p);
446 return Jacobi2QuadratureRule1D<ctype>(p);
448 return GaussLobattoQuadratureRule1D<ctype>(p);
450 return JacobiNQuadratureRule1D<ctype>(p);
452 return GaussRadauLeftQuadratureRule1D<ctype>(p);
454 return GaussRadauRightQuadratureRule1D<ctype>(p);
455 default :
456 DUNE_THROW(Exception, "Unknown QuadratureType");
457 }
458 }
459 DUNE_THROW(Exception, "Unknown GeometryType");
460 }
461 };
462
463 template<typename ctype>
464 class QuadratureRuleFactory<ctype, 2> {
465 private:
466 constexpr static int dim = 2;
467 friend class QuadratureRules<ctype, dim>;
468 static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
469 {
470 unsigned order =
471 TensorProductQuadratureRule<ctype,dim>::maxOrder(t.id(), qt);
472 if (t.isSimplex())
473 order = std::max
474 (order, static_cast<unsigned>(SimplexQuadratureRule<ctype,dim>::highest_order));
475 return order;
476 }
477 static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
478 {
479 if (t.isSimplex()
481 && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
482 {
483 return SimplexQuadratureRule<ctype,dim>(p);
484 }
485 return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
486 }
487 };
488
489 template<typename ctype>
490 class QuadratureRuleFactory<ctype, 3> {
491 private:
492 constexpr static int dim = 3;
493 friend class QuadratureRules<ctype, dim>;
494 static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
495 {
496 unsigned order =
497 TensorProductQuadratureRule<ctype,dim>::maxOrder(t.id(), qt);
498 if (t.isSimplex())
499 order = std::max
500 (order, static_cast<unsigned>(SimplexQuadratureRule<ctype,dim>::highest_order));
501 if (t.isPrism())
502 order = std::max
503 (order, static_cast<unsigned>(PrismQuadratureRule<ctype,dim>::highest_order));
504 return order;
505 }
506 static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
507 {
508
509 if (t.isSimplex()
511 && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
512 {
513 return SimplexQuadratureRule<ctype,dim>(p);
514 }
515 if (t.isPrism()
517 && p <= PrismQuadratureRule<ctype,dim>::highest_order)
518 {
519 return PrismQuadratureRule<ctype,dim>(p);
520 }
521 return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
522 }
523 };
524
525#ifndef DUNE_NO_EXTERN_QUADRATURERULES
526 extern template class GaussLobattoQuadratureRule<double, 1>;
527 extern template class GaussQuadratureRule<double, 1>;
528 extern template class GaussRadauLeftQuadratureRule<double, 1>;
529 extern template class GaussRadauRightQuadratureRule<double, 1>;
530 extern template class Jacobi1QuadratureRule<double, 1>;
531 extern template class Jacobi2QuadratureRule<double, 1>;
532 extern template class JacobiNQuadratureRule<double, 1>;
533 extern template class PrismQuadratureRule<double, 3>;
534 extern template class SimplexQuadratureRule<double, 2>;
535 extern template class SimplexQuadratureRule<double, 3>;
536#endif // !DUNE_NO_EXTERN_QUADRATURERULES
537
538} // end namespace
539
540#endif // DUNE_GEOMETRY_QUADRATURERULES_HH
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
constexpr bool isVertex() const
Return true if entity is a vertex.
Definition: type.hh:279
constexpr unsigned int dim() const
Return dimension of the type.
Definition: type.hh:360
BasicType
Each entity can be tagged by one of these basic types plus its space dimension.
Definition: type.hh:120
constexpr unsigned int id() const
Return the topology id of the type.
Definition: type.hh:365
static constexpr std::size_t size(std::size_t dim)
Compute total number of geometry types for the given dimension.
Definition: typeindex.hh:61
static constexpr std::size_t index(const GeometryType &gt)
Compute the index for the given geometry type within its dimension.
Definition: typeindex.hh:73
Default exception for dummy implementations.
Definition: exceptions.hh:263
Exception thrown if a desired QuadratureRule is not available, because the requested order is to high...
Definition: quadraturerules.hh:55
Single evaluation point in a quadrature rule.
Definition: quadraturerules.hh:66
const Vector & position() const
return local coordinates of integration point i
Definition: quadraturerules.hh:82
Dune::FieldVector< ct, dim > Vector
Type used for the position of a quadrature point.
Definition: quadraturerules.hh:75
ct Field
Number type used for coordinates and quadrature weights.
Definition: quadraturerules.hh:72
const ct & weight() const
return weight associated with integration point i
Definition: quadraturerules.hh:88
static constexpr int dimension
Dimension of the integration domain.
Definition: quadraturerules.hh:69
QuadraturePoint(const Vector &x, ct w)
set up quadrature of given order in d dimensions
Definition: quadraturerules.hh:78
Factory class for creation of quadrature rules, depending on GeometryType, order and QuadratureType.
Definition: quadraturerules.hh:370
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:214
virtual GeometryType type() const
return type of element
Definition: quadraturerules.hh:240
QuadratureRule(GeometryType t, int order)
Constructor for a given geometry type and a given quadrature order.
Definition: quadraturerules.hh:228
ct CoordType
The type used for coordinates.
Definition: quadraturerules.hh:234
QuadratureRule()
Default constructor.
Definition: quadraturerules.hh:221
virtual int order() const
return order
Definition: quadraturerules.hh:237
QuadratureRule(GeometryType t)
Constructor for a given geometry type. Leaves the quadrature order invalid
Definition: quadraturerules.hh:225
std::vector< QuadraturePoint< ct, dim > >::const_iterator iterator
Definition: quadraturerules.hh:245
A container for all quadrature rules of dimension dim
Definition: quadraturerules.hh:260
static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
maximum quadrature order for given geometry type and quadrature type
Definition: quadraturerules.hh:319
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:332
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:326
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:218
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:158
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
Enum
Definition: quadraturerules.hh:131
@ GaussJacobi_n_0
Gauss-Legendre rules with .
Definition: quadraturerules.hh:168
@ GaussJacobi_2_0
Gauss-Legendre rules with .
Definition: quadraturerules.hh:155
@ GaussRadauRight
Gauss-Radau rules including the right endpoint.
Definition: quadraturerules.hh:193
@ GaussJacobi_1_0
Gauss-Jacobi rules with .
Definition: quadraturerules.hh:148
@ GaussLobatto
Gauss-Lobatto rules.
Definition: quadraturerules.hh:176
@ GaussRadauLeft
Gauss-Radau rules including the left endpoint.
Definition: quadraturerules.hh:184
@ GaussLegendre
Gauss-Legendre rules (default)
Definition: quadraturerules.hh:141
Dune namespace.
Definition: alignedallocator.hh:13
STL namespace.
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:20
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Oct 13, 22:30, 2024)