/usr/include/dune/geometry/quadraturerules/gaussquadrature.hh is in libdune-geometry-dev 2.2.1-2ubuntu2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 | #ifndef DUNE_GEOMETRY_QUADRATURERULES_GAUSSQUADRATURE_HH
#define DUNE_GEOMETRY_QUADRATURERULES_GAUSSQUADRATURE_HH
#if HAVE_ALGLIB
#include <alglib/gqgengauss.h>
#warning ALGLIB support is deprecated, thus higher precision Gauss points will be dropped after DUNE 2.2 (cf. FS#931)
#endif
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/quadraturerules/genericquadrature.hh>
namespace Dune
{
namespace GenericGeometry
{
// GaussPoints
// -----------
/**
* @brief Gauss quadrature points and weights in 1d.
*
* if ALGLib is found higher precision Gauss points can be used by
* prescribing the amp::ampf field type; otherwise
* the dune-geometry quadrature is used.
**/
template< class F>
class GaussPoints
: public std::vector< QuadraturePoint<F,1> >
{
typedef std::vector< QuadraturePoint<F,1> > Base;
typedef GaussPoints< F > This;
public:
typedef F Field;
// n is number of points required
explicit GaussPoints ( unsigned int n )
{
Base::reserve( n );
const QuadratureRule<Field,1>& points =
QuadratureRules<Field,1>::rule(GeometryType(GeometryType::cube,1),
2*n-2, QuadratureType::Gauss);
for( unsigned int i = 0; i < n; ++i )
{
QuadraturePoint<Field,1> q( points[i].position()[0], points[i].weight() );
Base::push_back( q );
}
}
// order is the maximal order of polynomials which can be exactly integrated
static unsigned int minPoints( unsigned int order )
{
return (order+2)/2;
}
};
#if HAVE_ALGLIB
template< unsigned int precision >
class GaussPoints< amp::ampf< precision > >
: public std::vector< QuadraturePoint<amp::ampf< precision>,1> >
{
typedef amp::ampf< precision > F;
typedef std::vector< QuadraturePoint<F,1> > Base;
typedef GaussPoints< F > This;
public:
typedef F Field;
explicit GaussPoints ( unsigned int n )
{
Base::reserve( n );
typedef ap::template_1d_array< Field > AlgLibVector;
AlgLibVector p,w;
p.setbounds( 0, n-1 );
w.setbounds( 0, n-1 );
AlgLibVector alpha,beta;
alpha.setbounds( 0, n-1 );
beta.setbounds( 0, n-1 );
for( unsigned int i = 0; i < n; ++i )
{
alpha( i ) = 0;
beta( i ) = Field( i*i ) / Field( 4*(i*i)-1 );
}
gqgengauss::generategaussquadrature< precision >( alpha, beta, Field( 2 ), n, p, w );
const Field half = Field( 1 ) / Field( 2 );
for( unsigned int i = 0; i < n; ++i )
{
QuadraturePoint<Field,1> q( (p( i ) + Field( 1 )) * half, w( i )*half );
Base::push_back( q );
}
}
// order is the maximal order of polynomials which can be exactly integrated
static unsigned int minPoints( unsigned int order )
{
return (order+2)/2;
}
};
#endif
/**
* @brief Singleton provider for Gauss quadratures
*
* \tparam dim dimension of the reference elements contained in the factory
* \tparam F field in which weight and point of the quadrature are stored
* \tparam CF the compute field for the points and weights
**/
template< int dim, class F, class CF=F >
struct GaussQuadratureProvider
: public TopologySingletonFactory< GenericQuadratureFactory< dim, F, GaussPoints<CF> > >
{};
}
}
#endif // #ifndef DUNE_GEOMETRY_QUADRATURERULES_GAUSSQUADRATURE_HH
|