dune-fem 2.8.0
Loading...
Searching...
No Matches
femquadratures_inline.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_FEMQUADRATURES_INLINE_HH
2#define DUNE_FEM_FEMQUADRATURES_INLINE_HH
3
4#include "femquadratures.hh"
5
6namespace Dune
7{
8
9 namespace Fem
10 {
11
12#define SimplexPointsAdapter ParDGSimplexQuadrature::Quadrature
13
14 // only if we use dune-fem quadratures
15 template <class ct, int dim>
16 SimplexQuadrature<ct, dim>::SimplexQuadrature(const GeometryType&, int order, size_t id) :
17 QuadratureImp<ct, dim>(id),
18 order_((order <= 0) ? 1 : order)
19 {
20 // make sure that we only use orders that are available
21 assert( order_ <= SimplexMaxOrder :: maxOrder( dim ) );
22
23 SimplexPointsAdapter<dim> points(order);
24
25 order_ = points.order();
26
27 for (int i = 0; i < points.numPoints(); ++i) {
28 this->addQuadraturePoint(points.point(i), points.weight(i));
29 }
30 }
31
32 template <class ct, int dim>
33 CubeQuadrature<ct, dim>::CubeQuadrature(const GeometryType&, int order, size_t id) :
34 QuadratureImp<ct, dim>(id),
35 order_((order <= 0) ? 1 : order)
36 {
37 assert( order_ <= GaussPts::highestOrder );
38
39 typedef FieldVector< ct, dim > CoordinateType;
40
41 const GaussPts gp;
42
43 // find the right Gauss Rule from given order
44 int m = 0;
45 for (int i = 0; i <= GaussPts::MAXP; i++) {
46 if (gp.order(i)>=order_) {
47 m = i;
48 break;
49 }
50 }
51
52 if (m==0) DUNE_THROW(NotImplemented, "order " << order_ << " not implemented");
53 order_ = gp.order(m);
54
55 // fill in all the gauss points
56 int n = gp.power(m,dim);
57 for (int i = 0; i < n; i++) {
58 // compute multidimensional coordinates of Gauss point
59 int x[dim];
60 int z = i;
61 for (int k=0; k<dim; ++k) {
62 x[k] = z%m;
63 z = z/m;
64 }
65
66 // compute coordinates and weight
67 double weight = 1.0;
68 CoordinateType local;
69 for (int k = 0; k < dim; k++)
70 {
71 local[k] = gp.point(m,x[k]);
72 weight *= gp.weight(m,x[k]);
73 }
74
75 // put in container
76 this->addQuadraturePoint(local, weight);
77 }
78
79 }
80
81 template <>
82 inline CubeQuadrature<double, 0>::CubeQuadrature(const GeometryType&, int order, size_t id) :
83 QuadratureImp<double, 0>(id),
84 order_((order <= 0) ? 1 : order)
85 {
86 typedef FieldVector< double, 0 > CoordinateType;
87
88 order_ = 20;
89
90 // fill in all the gauss points
91 // compute coordinates and weight
92 double weight = 1.0;
93 CoordinateType local( 0 );
94
95 // put in container
96 this->addQuadraturePoint(local, weight);
97
98 }
99
100
101 template <class ct>
102 PrismQuadrature<ct>::PrismQuadrature(const GeometryType&, int order, size_t id) :
103 QuadratureImp<ct, 3>(id),
104 order_((order <= 0) ? 1 : order)
105 {
106 SimplexPointsAdapter<2> simplexPoints(order);
107 int simplexOrder = simplexPoints.order();
108
109 const GaussPts gp;
110
111 // find the right Gauss Rule from given order
112 int m = 0;
113 for (int i = 0; i <= GaussPts::MAXP; i++) {
114 if (gp.order(i)>=order_) {
115 m = i;
116 break;
117 }
118 }
119 if (m==0) DUNE_THROW(NotImplemented, "order not implemented");
120
121 int gaussOrder = gp.order(m);
122 int minOrder = ((simplexOrder < gaussOrder) ? simplexOrder : gaussOrder);
123 order_ = minOrder;
124
125 int numSimplexPoints = simplexPoints.numPoints();
126 int numGaussPoints = gp.power(m,1);
127
128 FieldVector<ct, 3> local;
129 double weight, simplexWeight;
130
131 for (int i = 0; i < numSimplexPoints; ++i) {
132 local[0] = simplexPoints.point(i)[0];
133 local[1] = simplexPoints.point(i)[1];
134 simplexWeight = simplexPoints.weight(i);
135 for (int j = 0; j < numGaussPoints; ++j) {
136 local[2] = gp.point(m,j);
137 weight = simplexWeight;
138 weight *= gp.weight(m,j);
139 this->addQuadraturePoint(local, weight);
140 }
141 }
142 }
143
144 template <class ct>
145 PyramidQuadrature<ct>::PyramidQuadrature(const GeometryType&, int order, size_t id) :
146 QuadratureImp<ct, 3>(id),
147 order_((order <= 0) ? 1 : order)
148 {
149 const PyramidPoints points;
150
151 int m = 0;
152 for (int i = 0; i < PyramidPoints::numQuads; ++i) {
153 if (points.order(i) >= order_) {
154 m = i;
155 break;
156 }
157 }
158
159 if (m==0) DUNE_THROW(NotImplemented, "order not implemented");
160 order_ = points.order(m);
161
162 // fill in the points
163 for (int i = 0; i < points.numPoints(m); ++i) {
164 this->addQuadraturePoint(points.point(m, i), points.weight(m, i));
165 }
166 }
167
168
169 template <class ct, int dim>
170 PolyhedronQuadrature<ct, dim>::PolyhedronQuadrature(const GeometryType& geomType, int order, size_t id) :
171 QuadratureImp<ct, dim>(id),
172 geometryType_( geomType ),
173 order_((order <= 0) ? 1 : order)
174 {
175 }
176
177
178 } // namespace Fem
179
180} // namespace Dune
181
182#endif // #ifndef DUNE_FEM_FEMQUADRATURES_INLINE_HH
Definition: bindguard.hh:11
int order_
Definition: femquadratures.hh:66
virtual int order() const
obtain order of the integration point list
Definition: femquadratures.hh:86
SimplexQuadrature(const GeometryType &geometry, int order, size_t id)
constructor filling the list of points and weights
Definition: femquadratures_inline.hh:16
CubeQuadrature(const GeometryType &geometry, int order, size_t id)
constructor filling the list of points and weights
Definition: femquadratures_inline.hh:33
int order_
Definition: femquadratures.hh:132
BaseType::CoordinateType CoordinateType
type of local coordinates
Definition: femquadratures.hh:129
PrismQuadrature(const GeometryType &geometry, int order, size_t id)
constructor filling the list of points and weights
Definition: femquadratures_inline.hh:102
virtual int order() const
obtain order of the integration point list
Definition: femquadratures.hh:208
PyramidQuadrature(const GeometryType &geometry, int order, size_t id)
constructor filling the list of points and weights
Definition: femquadratures_inline.hh:145
PolyhedronQuadrature(const GeometryType &geometry, int order, size_t id)
constructor filling the list of points and weights
Definition: femquadratures_inline.hh:170
one-dimensional Gauss points and their weights
Definition: gausspoints.hh:26
double point(int m, int i) const
obtain the i-th point of the m-th quadratre
Definition: gausspoints.hh:51
@ highestOrder
Definition: gausspoints.hh:32
int order(int m) const
obtain the order of the m-th quadratre
Definition: gausspoints.hh:76
@ MAXP
Definition: gausspoints.hh:29
int power(int y, int d) const
a simple power method
Definition: gausspoints.hh:90
double weight(int m, int i) const
obtain the i-th weight of the m-th quadratre
Definition: gausspoints.hh:64
Definition: pyramidpoints.hh:17
const FieldVector< double, 3 > & point(int m, int i) const
Access to the ith point of quadrature rule m.
Definition: pyramidpoints.hh:24
int numPoints(int m) const
Number of points in the quadrature rule m.
Definition: pyramidpoints.hh:42
int order(int m) const
Actual order of quadrature rule m.
Definition: pyramidpoints.hh:36
double weight(int m, int i) const
Access to the ith weight of quadrature rule m.
Definition: pyramidpoints.hh:30
@ numQuads
Definition: pyramidpoints.hh:19
Generic implementation of a Dune quadrature.
Definition: quadratureimp.hh:196
void addQuadraturePoint(const CoordinateType &point, const FieldType weight)
Adds a point-weight pair to the quadrature.
Definition: quadratureimp.hh:270
const FieldType & weight(size_t i) const
obtain weight of i-th integration point
Definition: quadratureimp.hh:251