From ab46993556503f170afb9e16a55f890f78ef2d68 Mon Sep 17 00:00:00 2001 From: heltai Date: Wed, 11 Mar 2009 18:06:11 +0000 Subject: [PATCH] Added QGaussLogR class git-svn-id: https://svn.dealii.org/trunk@18481 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/quadrature_lib.h | 83 ++++++++++++++-------- deal.II/base/source/quadrature_lib.cc | 82 ++++++++++++++------- deal.II/doc/news/changes.h | 8 +++ 3 files changed, 121 insertions(+), 52 deletions(-) diff --git a/deal.II/base/include/base/quadrature_lib.h b/deal.II/base/include/base/quadrature_lib.h index 0606d258ea..fc4429dcb8 100644 --- a/deal.II/base/include/base/quadrature_lib.h +++ b/deal.II/base/include/base/quadrature_lib.h @@ -309,24 +309,18 @@ class QWeddle : public Quadrature /** - * Gauss Quadrature Formula to integrate ln[(x-a)/(b-a)]*f(x) on the - * interval [a,b], where f is a smooth function without - * singularities. See Numerical Recipes for more information. The - * quadrature formula weights are chosen in order to integrate - * f(x)*ln[(x-a)/(b-a)] (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 f(x). Using logarithm - * properties, one can add the integral of f(x)*ln(b-a) (calculated - * with QGauss class formulas) to the integral of f(x)*ln[(x-a)/(b-a)] - * to finally obtain the integral of f(x)*ln(x-a). It is also - * possible to integrate with weighting function ln[(b-x)/(b-a)] if - * revert is set to true. - * - * It works up to n=12. + * Gauss Quadrature Formula with logarithmic weighting function. This + * formula is used to to integrate ln|x|*f(x) on the interval + * [0,1], where f is a smooth function without + * singularities. The collection of quadrature points and weights has + * been obtained using Numerical Recipes. * -//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 f(x) 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 ln|x| to ln|1-x|. + * + * The weights and functions have been tabulated up to order 12. * */ template @@ -355,21 +349,54 @@ class QGaussLog : public Quadrature std::vector 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 &p, - std::vector &w) const; - }; +/** + * Gauss Quadrature Formula with arbitrary logarithmic weighting + * function. This formula is used to to integrate + * ln(|x-x0|/alpha)*f(x) on the interval [0,1], + * where f is a smooth function without singularities, and x0 and + * alpha are given at construction time, and are the location of the + * singularity x0 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 f(x) 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 +class QGaussLogR : public Quadrature { +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 x0 = Point(), + const double alpha = 1); + + double weight_correction(const Point &q); + +protected: + /** This is the length of interval (0,origin), or 1 if either of + * the two extremes have been selected. */ + const double fraction; +}; + + /*@}*/ @@ -391,7 +418,6 @@ template <> unsigned int QGaussLobatto<1>::gamma(const unsigned int n) const; -template <> void QGaussLog<1>::revert_points_and_weights(std::vector &, std::vector &) const; template <> std::vector QGaussLog<1>::set_quadrature_points(const unsigned int) const; template <> std::vector QGaussLog<1>::set_quadrature_weights(const unsigned int) const; @@ -407,6 +433,7 @@ template <> QSimpson<1>::QSimpson (); 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); diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc index d05d70bc67..5723a507c1 100644 --- a/deal.II/base/source/quadrature_lib.cc +++ b/deal.II/base/source/quadrature_lib.cc @@ -645,16 +645,16 @@ QGaussLog<1>::QGaussLog(const unsigned int n, std::vector p=set_quadrature_points(n); std::vector w=set_quadrature_weights(n); - - if (revert == true) - revert_points_and_weights(p,w); for (unsigned int i=0; isize(); ++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 <> @@ -925,29 +925,63 @@ QGaussLog<1>::set_quadrature_weights(const unsigned int n) const } -template <> -void -QGaussLog<1>::revert_points_and_weights(std::vector &p, - std::vector &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; iquadrature_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; iquadrature_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 diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 50dd32393a..edd544096b 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -400,6 +400,14 @@ inconvenience this causes.

base

    +
  1. +

    + New: There is now a new QGaussLogR class, that generalizes the QGaussLog class to + allow for arbitrary location of singularities, and singularity factors. +
    + (Luca Heltai 2009/03/11) +

    +

  2. New: The FunctionParser class now supports the fparser library's interface to use -- 2.39.5