* the Gauss-Lobatto-Legendre quadrature (<i>alpha = beta = 0</i>)
* is a special case.
*
- * @sa http://en.wikipedia.org/wiki/Handbook_of_Mathematical_Functions
+ * @sa http://en.wikipedia.org/wiki/Handbook_of_Mathematical_Functions
* @sa Karniadakis, G.E. and Sherwin, S.J.:
- * Spectral/hp element methods for computational fluid dynamics.
- * Oxford: Oxford University Press, 2005
+ * Spectral/hp element methods for computational fluid dynamics.
+ * Oxford: Oxford University Press, 2005
*
* @author Guido Kanschat, 2005, 2006; F. Prill, 2006
*/
compute_quadrature_weights (const std::vector<long double> &x,
const int alpha,
const int beta) const;
-
+
/**
* Evaluate a Jacobi polynomial
* $ P^{\alpha, \beta}_n(x) $
/**
* Evaluate the Gamma function
- * $ \Gamma(n) = (n-1)! $.
+ * $ \Gamma(n) = (n-1)! $.
* @param n point of evaluation (integer).
*/
long double gamma(const unsigned int n) const;
};
-
+
/**
* @deprecated Use QGauss for arbitrary order Gauss formulae instead!
/**
- * Simpson quadrature rule, exact for polynomials of degree 3.
+ * Simpson quadrature rule, exact for polynomials of degree 3.
*/
template <int dim>
class QSimpson : public Quadrature<dim>
public:
/**
* Generate a formula with
- * <tt>n</tt> quadrature points
+ * <tt>n</tt> quadrature points
*/
QGaussLog(const unsigned int n,
const bool revert=false);
-
- protected:
- /**
+
+ protected:
+ /**
* Sets the points of the
* quadrature formula.
*/
std::vector<double>
set_quadrature_points(const unsigned int n) const;
- /**
+ /**
* Sets the weights of the
* quadrature formula.
*/
* alpha is 1, then this quadrature is the same as QGaussLog.
*
* The last argument from the constructor allows you to use this
- * quadrature rule in one of two possible ways:
+ * quadrature rule in one of two possible ways:
* \f[
* \int_0^1 g(x) dx =
* \int_0^1 f(x) \ln\left(\frac{|x-x_0|}{\alpha}\right) dx
* it is factored out, to be included in
* the integrand.
*/
- QGaussLogR(const unsigned int n,
- const Point<dim> x0 = Point<dim>(),
+ QGaussLogR(const unsigned int n,
+ const Point<dim> x0 = Point<dim>(),
const double alpha = 1,
const bool factor_out_singular_weight=false);
* element $[0,1]^2$, where $f$ is a smooth function without
* singularities, and $R$ is the distance from the point $x$ to the vertex
* $\xi$, given at construction time by specifying its index. Notice that
- * this distance is evaluated in the reference element.
+ * this distance is evaluated in the reference element.
*
* This quadrature formula is obtained from two QGauss quadrature
* formulas, upon transforming them into polar coordinate system
* Upon construction it is possible to specify wether we want the
* singularity removed, or not. In other words, this quadrature can be
* used to integrate $g(x) = 1/R\ f(x)$, or simply $f(x)$, with the $1/R$
- * factor already included in the quadrature weights.
+ * factor already included in the quadrature weights.
*/
template<int dim>
class QGaussOneOverR : public Quadrature<dim>
* }
* @endcode
*/
- QGaussOneOverR(const unsigned int n,
+ QGaussOneOverR(const unsigned int n,
const Point<dim> singularity,
const bool factor_out_singular_weight=false);
/**
* }
* @endcode
*/
- QGaussOneOverR(const unsigned int n,
+ QGaussOneOverR(const unsigned int n,
const unsigned int vertex_index,
const bool factor_out_singular_weight=false);
private:
long double QGaussLobatto<1>::
JacobiP(const long double, const int, const int, const unsigned int) const;
template <>
-long double
+long double
QGaussLobatto<1>::gamma(const unsigned int n) const;
template <> std::vector<double> QGaussLog<1>::set_quadrature_points(const unsigned int) const;
// double exists and is described
// by std::numeric_limits, we may
// not actually get the additional
- // precission. One case where this
+ // precision. One case where this
// happens is on x86, where one can
// set hardware flags that disable
// long double precision even for
// a similar situation exists, btw,
// when running programs under
// valgrind up to and including at
- // least version 3.1: valgrind's
+ // least version 3.3: valgrind's
// emulator only supports 64 bit
// arithmetic, even for 80 bit long
// doubles.
origins.push_back(Point<2>(singularity[0],0.));
origins.push_back(Point<2>(0.,singularity[1]));
origins.push_back(singularity);
-
+
// Lexycographical ordering.
-
+
double eps = 1e-8;
unsigned int q_id = 0; // Current quad point index.
double area = 0;
Point<2> dist;
- for(unsigned int box=0; box<4; ++box)
+ for(unsigned int box=0; box<4; ++box)
{
dist = (singularity-GeometryInfo<2>::unit_cell_vertex(box));
- dist = Point<2>(std::abs(dist[0]), std::abs(dist[1]));
+ dist = Point<2>(std::abs(dist[0]), std::abs(dist[1]));
area = dist[0]*dist[1];
- if(area > eps)
- for(unsigned int q=0; q<quads[box].size(); ++q, ++q_id)
+ if(area > eps)
+ for(unsigned int q=0; q<quads[box].size(); ++q, ++q_id)
{
const Point<2> &qp = quads[box].point(q);
this->quadrature_points[q_id] =