4 #ifndef DUNE_GEOMETRY_QUADRATURERULES_HH
5 #define DUNE_GEOMETRY_QUADRATURERULES_HH
12 #include <dune/common/fvector.hh>
13 #include <dune/common/exceptions.hh>
14 #include <dune/common/stdstreams.hh>
15 #include <dune/common/visibility.hh>
16 #include <dune/common/deprecated.hh>
38 template<
typename ct,
int dim>
42 enum {
d=dim } DUNE_DEPRECATED_MSG(
"Use 'dimension' instead");
43 typedef ct
CoordType DUNE_DEPRECATED_MSG(
"Use type 'Field' instead");
47 typedef Dune::FieldVector<ct,dim>
Vector;
75 namespace QuadratureType {
92 template<
typename ct,
int dim>
121 typedef typename std::vector<QuadraturePoint<ct,dim> >::const_iterator
iterator;
135 template<
typename ctype,
int dim>
141 typedef std::pair<GeometryType,int> QuadratureRuleKey;
149 static std::map<QuadratureRuleKey, QuadratureRule> _quadratureMap;
150 QuadratureRuleKey key(t,p);
151 if (_quadratureMap.find(key) == _quadratureMap.end()) {
159 _quadratureMap.insert(std::make_pair(key, rule));
161 return _quadratureMap.find(key)->second;
175 return instance()._rule(t,p,qt);
181 return instance()._rule(gt,p,qt);
187 #include "quadraturerules/pointquadrature.hh"
192 template<typename ct, bool fundamental = std::numeric_limits<ct>::is_specialized>
193 struct GaussQuadratureInitHelper;
194 template<
typename ct>
195 struct GaussQuadratureInitHelper<ct, true> {
196 static void init(
int p,
197 std::vector< FieldVector<ct, 1> > & _points,
198 std::vector< ct > & _weight,
199 int & delivered_order);
201 template<
typename ct>
202 struct GaussQuadratureInitHelper<ct, false> {
203 static void init(
int p,
204 std::vector< FieldVector<ct, 1> > & _points,
205 std::vector< ct > & _weight,
206 int & delivered_order);
210 template<
typename ct>
226 std::vector< FieldVector<ct, dim> > _points;
227 std::vector< ct > _weight;
229 GaussQuadratureInitHelper<ct>::init
232 assert(_points.size() == _weight.size());
233 for (
size_t i = 0; i < _points.size(); i++)
238 extern template GaussQuadratureRule1D<float>::GaussQuadratureRule1D(
int);
239 extern template GaussQuadratureRule1D<double>::GaussQuadratureRule1D(
int);
243 #define DUNE_INCLUDING_IMPLEMENTATION
244 #include "quadraturerules/gauss_imp.hh"
249 template<
typename ct,
250 bool fundamental = std::numeric_limits<ct>::is_specialized>
251 struct Jacobi1QuadratureInitHelper;
252 template<
typename ct>
253 struct Jacobi1QuadratureInitHelper<ct, true> {
254 static void init(
int p,
255 std::vector< FieldVector<ct, 1> > & _points,
256 std::vector< ct > & _weight,
257 int & delivered_order);
259 template<
typename ct>
260 struct Jacobi1QuadratureInitHelper<ct, false> {
261 static void init(
int p,
262 std::vector< FieldVector<ct, 1> > & _points,
263 std::vector< ct > & _weight,
264 int & delivered_order);
270 template<
typename ct>
279 enum { highest_order=61 };
288 std::vector< FieldVector<ct, dim> > _points;
289 std::vector< ct > _weight;
293 Jacobi1QuadratureInitHelper<ct>::init
294 (p, _points, _weight, delivered_order);
295 this->delivered_order = delivered_order;
296 assert(_points.size() == _weight.size());
297 for (
size_t i = 0; i < _points.size(); i++)
303 extern template Jacobi1QuadratureRule1D<float>::Jacobi1QuadratureRule1D(
int);
304 extern template Jacobi1QuadratureRule1D<double>::Jacobi1QuadratureRule1D(
int);
309 #define DUNE_INCLUDING_IMPLEMENTATION
310 #include "quadraturerules/jacobi_1_0_imp.hh"
315 template<
typename ct,
316 bool fundamental = std::numeric_limits<ct>::is_specialized>
317 struct Jacobi2QuadratureInitHelper;
318 template<
typename ct>
319 struct Jacobi2QuadratureInitHelper<ct, true> {
320 static void init(
int p,
321 std::vector< FieldVector<ct, 1> > & _points,
322 std::vector< ct > & _weight,
323 int & delivered_order);
325 template<
typename ct>
326 struct Jacobi2QuadratureInitHelper<ct, false> {
327 static void init(
int p,
328 std::vector< FieldVector<ct, 1> > & _points,
329 std::vector< ct > & _weight,
330 int & delivered_order);
336 template<
typename ct>
345 enum { highest_order=61 };
354 std::vector< FieldVector<ct, dim> > _points;
355 std::vector< ct > _weight;
359 Jacobi2QuadratureInitHelper<ct>::init
360 (p, _points, _weight, delivered_order);
362 this->delivered_order = delivered_order;
363 assert(_points.size() == _weight.size());
364 for (
size_t i = 0; i < _points.size(); i++)
370 extern template Jacobi2QuadratureRule1D<float>::Jacobi2QuadratureRule1D(
int);
371 extern template Jacobi2QuadratureRule1D<double>::Jacobi2QuadratureRule1D(
int);
376 #define DUNE_INCLUDING_IMPLEMENTATION
377 #include "quadraturerules/jacobi_2_0_imp.hh"
382 template<
typename ct,
383 bool fundamental = std::numeric_limits<ct>::is_specialized>
384 struct GaussLobattoQuadratureInitHelper;
385 template<
typename ct>
386 struct GaussLobattoQuadratureInitHelper<ct, true> {
387 static void init(
int p,
388 std::vector< FieldVector<ct, 1> > & _points,
389 std::vector< ct > & _weight,
390 int & delivered_order);
392 template<
typename ct>
393 struct GaussLobattoQuadratureInitHelper<ct, false> {
394 static void init(
int p,
395 std::vector< FieldVector<ct, 1> > & _points,
396 std::vector< ct > & _weight,
397 int & delivered_order);
403 template<
typename ct>
412 enum { highest_order=61 };
421 std::vector< FieldVector<ct, dim> > _points;
422 std::vector< ct > _weight;
426 GaussLobattoQuadratureInitHelper<ct>::init
427 (p, _points, _weight, delivered_order);
429 this->delivered_order = delivered_order;
430 assert(_points.size() == _weight.size());
431 for (
size_t i = 0; i < _points.size(); i++)
437 extern template GaussLobattoQuadratureRule1D<float>::GaussLobattoQuadratureRule1D(
int);
438 extern template GaussLobattoQuadratureRule1D<double>::GaussLobattoQuadratureRule1D(
int);
443 #define DUNE_INCLUDING_IMPLEMENTATION
444 #include "quadraturerules/gausslobatto_imp.hh"
446 #include "quadraturerules/tensorproductquadrature.hh"
448 #include "quadraturerules/simplexquadrature.hh"
458 class PrismQuadraturePoints;
466 enum { highest_order=2 };
469 PrismQuadraturePoints ()
500 W[m][0] = 0.16666666666666666 / 2.0;
501 W[m][1] = 0.16666666666666666 / 2.0;
502 W[m][2] = 0.16666666666666666 / 2.0;
503 W[m][3] = 0.16666666666666666 / 2.0;
504 W[m][4] = 0.16666666666666666 / 2.0;
505 W[m][5] = 0.16666666666666666 / 2.0;
512 G[m][0][0] =0.66666666666666666 ;
513 G[m][0][1] =0.16666666666666666 ;
514 G[m][0][2] =0.211324865405187 ;
516 G[m][1][0] = 0.16666666666666666;
517 G[m][1][1] =0.66666666666666666 ;
518 G[m][1][2] = 0.211324865405187;
520 G[m][2][0] = 0.16666666666666666;
521 G[m][2][1] = 0.16666666666666666;
522 G[m][2][2] = 0.211324865405187;
524 G[m][3][0] = 0.66666666666666666;
525 G[m][3][1] = 0.16666666666666666;
526 G[m][3][2] = 0.788675134594813;
528 G[m][4][0] = 0.16666666666666666;
529 G[m][4][1] = 0.66666666666666666;
530 G[m][4][2] = 0.788675134594813;
532 G[m][5][0] = 0.16666666666666666;
533 G[m][5][1] = 0.16666666666666666;
534 G[m][5][2] = 0.788675134594813;
536 W[m][0] = 0.16666666666666666 / 2.0;
537 W[m][1] = 0.16666666666666666 / 2.0;
538 W[m][2] = 0.16666666666666666 / 2.0;
539 W[m][3] = 0.16666666666666666 / 2.0;
540 W[m][4] = 0.16666666666666666 / 2.0;
541 W[m][5] = 0.16666666666666666 / 2.0;
548 FieldVector<double, 3> point(
int m,
int i)
554 double weight (
int m,
int i)
566 FieldVector<double, 3> G[MAXP+1][MAXP];
568 double W[MAXP+1][MAXP];
592 template<
typename ct,
int dim>
593 class PrismQuadratureRule;
598 template<
typename ct>
607 enum { highest_order = 2 };
616 for(
int i=0; i<m; ++i)
618 FieldVector<ct,3> local;
619 for (
int k=0; k<d; k++)
635 template<
typename ctype,
int dim>
641 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
645 template<
typename ctype>
654 return PointQuadratureRule<ctype>();
656 DUNE_THROW(Exception,
"Unknown GeometryType");
660 template<
typename ctype>
679 DUNE_THROW(Exception,
"Unknown QuadratureType");
682 DUNE_THROW(Exception,
"Unknown GeometryType");
686 template<
typename ctype>
693 if (t.
isSimplex() && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
695 return SimplexQuadratureRule<ctype,dim>(p);
697 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
701 template<
typename ctype>
708 if (t.
isSimplex() && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
710 return SimplexQuadratureRule<ctype,dim>(p);
712 if (t.
isPrism() && p <= PrismQuadratureRule<ctype,dim>::highest_order)
714 return PrismQuadratureRule<ctype,dim>(p);
716 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
722 #endif // DUNE_GEOMETRY_QUADRATURERULES_HH