dune-fem  2.8-git
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 
6 namespace 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
static int maxOrder(const int dim)
Definition: femquadratures.hh:24
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