/**
- * Gauss Quadrature Formula to integrate <tt>ln[(x-a)/(b-a)]*f(x)</tt> on the
- * interval <tt>[a,b]</tt>, where f is a smooth function without
- * singularities. See <tt>Numerical Recipes</tt> for more information. The
- * quadrature formula weights are chosen in order to integrate
- * <tt>f(x)*ln[(x-a)/(b-a)]</tt> (that is, the argument of the logarithm goes
- * linearly from 0 to 1 on the interval of integration). So, the only
- * thing one can specify is the function <tt>f(x)</tt>. Using logarithm
- * properties, one can add the integral of <tt>f(x)*ln(b-a)</tt> (calculated
- * with QGauss class formulas) to the integral of <tt>f(x)*ln[(x-a)/(b-a)]</tt>
- * to finally obtain the integral of <tt>f(x)*ln(x-a)</tt>. It is also
- * possible to integrate with weighting function <tt>ln[(b-x)/(b-a)]</tt> if
- * revert is set to true.
- *
- * It works up to <tt>n=12</tt>.
+ * Gauss Quadrature Formula with logarithmic weighting function. This
+ * formula is used to to integrate <tt>ln|x|*f(x)</tt> on the interval
+ * <tt>[0,1]</tt>, where f is a smooth function without
+ * singularities. The collection of quadrature points and weights has
+ * been obtained using <tt>Numerical Recipes</tt>.
*
-//TODO: implement the calculation of points and weights as it is
-//described in Numerical Recipes, instead of passing them to set_*
-//functions, so that n can also be grater than 12.
+ * Notice that only the function <tt>f(x)</tt> should be provided,
+ * i.e., $\int_0^1 f(x) ln|x| dx = \sum_{i=0}^N w_i f(q_i)$. Setting
+ * the @p revert flag to true at construction time swithces the weight
+ * from <tt>ln|x|</tt> to <tt>ln|1-x|</tt>.
+ *
+ * The weights and functions have been tabulated up to order 12.
*
*/
template <int dim>
std::vector<double>
set_quadrature_weights(const unsigned int n) const;
- /**
- * If the singularity is at the rightmost point
- * of the interval just revert the vectors
- * of points and weights and everything
- * works.
- */
- void
- revert_points_and_weights(std::vector<double> &p,
- std::vector<double> &w) const;
-
};
+/**
+ * Gauss Quadrature Formula with arbitrary logarithmic weighting
+ * function. This formula is used to to integrate
+ * <tt>ln(|x-x0|/alpha)*f(x)</tt> on the interval <tt>[0,1]</tt>,
+ * where f is a smooth function without singularities, and x0 and
+ * alpha are given at construction time, and are the location of the
+ * singularity <tt>x0</tt> and an arbitrary scaling factor in the
+ * singularity.
+ *
+ * This quadrature formula is rather expensive, since it uses
+ * internally two Gauss quadrature formulas of order n to integrate
+ * the nonsingular part of the factor, and two GaussLog quadrature
+ * formulas to integrate on the separate segments [0,x0] and
+ * [x0,1]. If the singularity is one of the extremes and the factor
+ * alpha is 1, then this quadrature is the same as QGaussLog.
+ *
+ * Notice again that only the function <tt>f(x)</tt> should be
+ * provided, i.e., $\int_0^1 f(x) ln(|x-x0|/alpha) dx = \sum_{i=0}^N
+ * w_i f(q_i)$.
+ *
+ * The weights and functions have been tabulated up to order 12.
+ *
+ */
+template<int dim>
+class QGaussLogR : public Quadrature<dim> {
+public:
+ /** The constructor takes three arguments arguments: the order of
+ * the gauss formula on each of the segments [0,x0] and [x0,1],
+ * the actual location of the singularity and the scale factor for
+ * the logarithmic function. */
+ QGaussLogR(const unsigned int n,
+ const Point<dim> x0 = Point<dim>(),
+ const double alpha = 1);
+
+ double weight_correction(const Point<dim> &q);
+
+protected:
+ /** This is the length of interval (0,origin), or 1 if either of
+ * the two extremes have been selected. */
+ const double fraction;
+};
+
+
/*@}*/
unsigned int
QGaussLobatto<1>::gamma(const unsigned int n) const;
-template <> void QGaussLog<1>::revert_points_and_weights(std::vector<double> &, std::vector<double> &) const;
template <> std::vector<double> QGaussLog<1>::set_quadrature_points(const unsigned int) const;
template <> std::vector<double> QGaussLog<1>::set_quadrature_weights(const unsigned int) const;
template <> QMilne<1>::QMilne ();
template <> QWeddle<1>::QWeddle ();
template <> QGaussLog<1>::QGaussLog (const unsigned int n, const bool revert);
+template <> QGaussLogR<1>::QGaussLogR (const unsigned int n, const Point<1> x0, const double alpha);
std::vector<double> p=set_quadrature_points(n);
std::vector<double> w=set_quadrature_weights(n);
-
- if (revert == true)
- revert_points_and_weights(p,w);
for (unsigned int i=0; i<this->size(); ++i)
{
- this->quadrature_points[i] = Point<1>(p[i]);
- this->weights[i] = w[i];
+ // Using the change of variables x=1-t, it's possible to show
+ // that int f(x)ln|1-x| = int f(1-t) ln|t|, which implies that
+ // we can use this quadrature formula also with weight ln|1-x|.
+ this->quadrature_points[i] = revert ? Point<1>(1-p[n-1-i]) : Point<1>(p[i]);
+ this->weights[i] = revert ? w[n-1-i] : w[i];
}
-
+
}
template <>
}
-template <>
-void
-QGaussLog<1>::revert_points_and_weights(std::vector<double> &p,
- std::vector<double> &w) const
+template<>
+QGaussLogR<1>::QGaussLogR(const unsigned int n,
+ const Point<1> origin,
+ const double alpha) :
+ Quadrature<1>( ( (origin[0] == 0) || (origin[0] == 1) ) ?
+ (alpha == 1 ? n : 2*n ) : 4*n ),
+ fraction( ( (origin[0] == 0) || (origin[0] == 1.) ) ? 1. : origin[0] )
{
+ // The three quadrature formulas that make this one up. There are
+ // at most two when the origin is one of the extremes, and there is
+ // only one if the origin is one of the extremes and alpha is
+ // equal to one.
+ //
+ // If alpha is different from one, then we need a correction which
+ // is performed with a standard Gauss quadrature rule on each
+ // segment. This is not needed in the standard case where alpha is
+ // equal to one and the origin is on one of the extremes. We
+ // integrate with weight ln|(x-o)/alpha|. In the easy cases, we
+ // only need n quadrature points. In the most difficult one, we
+ // need 2*n points for the first segment, and 2*n points for the
+ // second segment.
+ QGaussLog<1> quad1(n, origin[0] != 0);
+ QGaussLog<1> quad2(n);
+ QGauss<1> quad(n);
+
+ // Check that the origin is inside 0,1
+ Assert( (fraction >= 0) && (fraction <= 1),
+ ExcMessage("Origin is outside [0,1]."));
- double temp;
- const unsigned int sz = this->size();
+ // Non singular offset. This is the start of non singular quad
+ // points.
+ unsigned int ns_offset = (fraction == 1) ? n : 2*n;
+
+ for(unsigned int i=0, j=ns_offset; i<n; ++i, ++j) {
+ // The first i quadrature points are the same as quad1, and
+ // are by default singular.
+ this->quadrature_points[i] = quad1.point(i)*fraction;
+ this->weights[i] = quad1.weight(i)*fraction;
+
+ // We need to scale with -log|fraction*alpha|
+ if( (alpha != 1) || (fraction != 1) ) {
+ this->quadrature_points[j] = quad.point(i)*fraction;
+ this->weights[j] = -log(alpha/fraction)*quad.weight(i)*fraction;
+ }
- for (unsigned int i=0; i<sz/2; i++)
- {
- temp = p[sz-1-i];
- p[sz-1-i] = p[i];
- p[i]=temp;
-
- temp = w[sz-1-i];
- w[sz-1-i]=w[i];
- w[i]=temp;
+ // In case we need the second quadrature as well, do it now.
+ if(fraction != 1) {
+ this->quadrature_points[i+n] = quad2.point(i)*(1-fraction)+Point<1>(fraction);
+ this->weights[i+n] = quad2.weight(i)*(1-fraction);
+
+ // We need to scale with -log|fraction*alpha|
+ this->quadrature_points[j+n] = quad.point(i)*(1-fraction)+Point<1>(fraction);
+ this->weights[j+n] = -log(alpha/(1-fraction))*quad.weight(i)*(1-fraction);
+ }
}
-
}
-
-
+
// construct the quadrature formulae in higher dimensions by