From 5e67dd32ad6dc34b3779f7a0a66a610d33a92054 Mon Sep 17 00:00:00 2001 From: gpitton Date: Wed, 22 Apr 2015 11:26:15 +0200 Subject: [PATCH] Chebyshev quadratures. Some issues on the test output. acceptable Chebyshev quadrature added fixed test added documentation and indentation added q_chebyshev in changes.h --- doc/news/changes.h | 9 +- include/deal.II/base/quadrature_lib.h | 106 +- source/base/quadrature_lib.cc | 224 +++ tests/base/quadrature_chebyshev.cc | 174 ++ tests/base/quadrature_chebyshev.output | 2144 ++++++++++++++++++++++++ 5 files changed, 2655 insertions(+), 2 deletions(-) create mode 100644 tests/base/quadrature_chebyshev.cc create mode 100644 tests/base/quadrature_chebyshev.output diff --git a/doc/news/changes.h b/doc/news/changes.h index 8afab0efed..e459754899 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -506,11 +506,18 @@ inconvenience this causes.

Specific improvements

    +
  1. New: Three new quadrature formulas in quadrature_lib, based on + Chebyshev quadrature rules. See functions QGaussChebyshev, + QGaussRadauChebyshev and QGaussLobattoChebyshev. +
    + (Giuseppe Pitton, Luca Heltai 2015/05/11) +
    2. +
  2. Fixed: MatrixOut now also works with Trilinos and PETSc matrices.
    (Wolfgang Bangerth, 2015/05/11)
    3. - +
  3. Changed: TrilinosWrappers::Vector, TrilinosWrappers::BlockVector, PETScWrappers::Vector, and PETScWrappers::BlockVector are deprecated. Either use the MPI or the deal.II version of the Vector/BlockVector. diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index da8cf11af5..33a633563a 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -511,6 +511,111 @@ public: /*@}*/ +/** +* Gauss-Chebyshev quadrature rules integrate the weighted product +* $\int_{-1}^1 f(x) w(x) dx$ with weight given by: +* $w(x) = 1/\sqrt{1-x^2}$. The nodes and weights are known analytically, +* and are exact for monomials up to the order $2n-1$, where $n$ is the number +* of quadrature points. +* Here we rescale the quadrature formula so that it is defined on t +* he interval [0,1] instead of [-1,1]. So the quadrature formulas inte +* grate exactly the integral $\int_0^1 f(x) w(x) dx$ with the weight: +* $w(x) = 1/sqrt{x(1-x)}$. +* For details see: +* M. Abramowitz & I.A. Stegun: Handbook of Mathematical Functions, par. 25.4.38 +* +* @author Giuseppe Pitton, Luca Heltai 2015 +**/ +template +class QGaussChebyshev : public Quadrature +{ +public: + //Generate a formula with n quadrature points + QGaussChebyshev(const unsigned int n); + +protected: + // Sets the points of the quadrature formula. + std::vector + set_quadrature_points(const unsigned int n) const; + + // Sets the weights of the quadrature formula. + std::vector + set_quadrature_weights(const unsigned int n) const; + +}; + + +/** +* Gauss-Radau-Chebyshev quadrature rules integrate the weighted product +* $\int_{-1}^1 f(x) w(x) dx$ with weight given by: +* $w(x) = 1/\sqrt{1-x^2}$ with the additional constraint that a quadrature point +* lies at one of the two extrema of the interval. +* The nodes and weights are known analytically, +* and are exact for monomials up to the order $2n-2$, where $n$ is the number +* of quadrature points. +* Here we rescale the quadrature formula so that it is defined on t +* he interval [0,1] instead of [-1,1]. So the quadrature formulas inte +* grate exactly the integral $\int_0^1 f(x) w(x) dx$ with the weight: +* $w(x) = 1/sqrt{x(1-x)}$. By default the quadrature is constructed with the +* left endpoint as quadrature node, but the quadrature node can be imposed at the +* right endpoint through the variable ep that can assume the values left or right. +* +* @author Giuseppe Pitton, Luca Heltai 2015 +**/ +template +class QGaussRadauChebyshev : public Quadrature +{ +public: + enum EndPoint { left,right }; + EndPoint ep; + //Generate a formula with n quadrature points + QGaussRadauChebyshev(const unsigned int n, + QGaussRadauChebyshev::EndPoint ep=QGaussRadauChebyshev::left); + +protected: + // Sets the points of the quadrature formula. + std::vector + set_quadrature_points(const unsigned int n) const; + + // Sets the weights of the quadrature formula. + std::vector + set_quadrature_weights(const unsigned int n) const; + +}; + +/** +* Gauss-Lobatto-Chebyshev quadrature rules integrate the weighted product +* $\int_{-1}^1 f(x) w(x) dx$ with weight given by: +* $w(x) = 1/\sqrt{1-x^2}$. The nodes and weights are known analytically, +* and are exact for monomials up to the order $2n-1$, where $n$ is the number +* of quadrature points. +* Here we rescale the quadrature formula so that it is defined on t +* he interval [0,1] instead of [-1,1]. So the quadrature formulas inte +* grate exactly the integral $\int_0^1 f(x) w(x) dx$ with the weight: +* $w(x) = 1/sqrt{x(1-x)}$. +* For details see: +* M. Abramowitz & I.A. Stegun: Handbook of Mathematical Functions, par. 25.4.40 +* +* @author Giuseppe Pitton, Luca Heltai 2015 +**/ +template +class QGaussLobattoChebyshev : public Quadrature +{ +public: + //Generate a formula with n quadrature points + QGaussLobattoChebyshev(const unsigned int n); + +protected: + // Sets the points of the quadrature formula. + std::vector + set_quadrature_points(const unsigned int n) const; + + // Sets the weights of the quadrature formula. + std::vector + set_quadrature_weights(const unsigned int n) const; + +}; + /* -------------- declaration of explicit specializations ------------- */ template <> QGauss<1>::QGauss (const unsigned int n); @@ -543,7 +648,6 @@ template <> QTelles<1>::QTelles(const Quadrature<1> &base_quad, const Point<1> & - DEAL_II_NAMESPACE_CLOSE #endif diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index b2df7848a9..0abe19989e 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -16,6 +16,7 @@ #include #include + #include #include #include @@ -1127,6 +1128,217 @@ QTelles<1>::QTelles ( } } + + +template <> +std::vector +QGaussChebyshev<1>::set_quadrature_points(const unsigned int n) const +{ + + std::vector points(n); + // n point quadrature: index from 0 to n-1 + for (unsigned short i=0; i +std::vector +QGaussChebyshev<1>::set_quadrature_weights(const unsigned int n) const +{ + + std::vector weights(n); + + for (unsigned short i=0; i +QGaussChebyshev<1>::QGaussChebyshev(const unsigned int n) + : + Quadrature<1> (n) +{ + + Assert(n>0,ExcMessage("Need at least one point for the quadrature rule")); + std::vector p=set_quadrature_points(n); + std::vector w=set_quadrature_weights(n); + + for (unsigned int i=0; isize(); ++i) + { + this->quadrature_points[i] = Point<1>(p[i]); + this->weights[i] = w[i]; + } + +} + + +template +QGaussChebyshev::QGaussChebyshev (const unsigned int n) + : + Quadrature (QGaussChebyshev(n), QGaussChebyshev<1>(n)) +{} + + + + + +template <> +std::vector +QGaussRadauChebyshev<1>::set_quadrature_points(const unsigned int n) const +{ + + std::vector points(n); + // n point quadrature: index from 0 to n-1 + for (unsigned short i=0; i +std::vector +QGaussRadauChebyshev<1>::set_quadrature_weights(const unsigned int n) const +{ + + std::vector weights(n); + + for (unsigned short i=0; i +QGaussRadauChebyshev<1>::QGaussRadauChebyshev(const unsigned int n, + QGaussRadauChebyshev<1>::EndPoint ep) + : + Quadrature<1> (n), + ep (ep) +{ + + Assert(n>0,ExcMessage("Need at least one point for quadrature rules")); + std::vector p=set_quadrature_points(n); + std::vector w=set_quadrature_weights(n); + + for (unsigned int i=0; isize(); ++i) + { + this->quadrature_points[i] = Point<1>(p[i]); + this->weights[i] = w[i]; + } +} + + +template <> +QGaussRadauChebyshev<2>::QGaussRadauChebyshev (const unsigned int n, + QGaussRadauChebyshev::EndPoint ep) + : + Quadrature<2> (QGaussRadauChebyshev<1>(n, static_cast::EndPoint>(ep)), + QGaussRadauChebyshev<1>(n, static_cast::EndPoint>(ep))) +{} + + +template +QGaussRadauChebyshev::QGaussRadauChebyshev (const unsigned int n, + QGaussRadauChebyshev::EndPoint ep) + : + Quadrature (QGaussRadauChebyshev(n,static_cast::EndPoint>(ep)), + QGaussRadauChebyshev<1>(n,static_cast::EndPoint>(ep))) +{} + + +template <> +std::vector +QGaussLobattoChebyshev<1>::set_quadrature_points(const unsigned int n) const +{ + + std::vector points(n); + // n point quadrature: index from 0 to n-1 + for (unsigned short i=0; i +std::vector +QGaussLobattoChebyshev<1>::set_quadrature_weights(const unsigned int n) const +{ + + std::vector weights(n); + + for (unsigned short i=0; i +QGaussLobattoChebyshev<1>::QGaussLobattoChebyshev(const unsigned int n) + : + Quadrature<1> (n) +{ + + Assert(n>1,ExcMessage("Need at least two points for Gauss-Lobatto quadrature rule")); + std::vector p=set_quadrature_points(n); + std::vector w=set_quadrature_weights(n); + + for (unsigned int i=0; isize(); ++i) + { + this->quadrature_points[i] = Point<1>(p[i]); + this->weights[i] = w[i]; + } + +} + + +template +QGaussLobattoChebyshev::QGaussLobattoChebyshev (const unsigned int n) + : + Quadrature (QGaussLobattoChebyshev(n), QGaussLobattoChebyshev<1>(n)) +{} + // explicit specialization // note that 1d formulae are specialized by implementation above template class QGauss<2>; @@ -1153,4 +1365,16 @@ template class QTelles<1> ; template class QTelles<2> ; template class QTelles<3> ; +template class QGaussChebyshev<1>; +template class QGaussChebyshev<2>; +template class QGaussChebyshev<3>; + +template class QGaussRadauChebyshev<1>; +template class QGaussRadauChebyshev<2>; +template class QGaussRadauChebyshev<3>; + +template class QGaussLobattoChebyshev<1>; +template class QGaussLobattoChebyshev<2>; +template class QGaussLobattoChebyshev<3>; + DEAL_II_NAMESPACE_CLOSE diff --git a/tests/base/quadrature_chebyshev.cc b/tests/base/quadrature_chebyshev.cc new file mode 100644 index 0000000000..5314b558cf --- /dev/null +++ b/tests/base/quadrature_chebyshev.cc @@ -0,0 +1,174 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// check accuracy of the Chebyshev quadrature formulas by using them to +// integrate polynomials of increasing degree, and finding the degree +// until which they integrate exactly + + +#include "../tests.h" +#include +#include + +#include +#include +#include + +using namespace dealii; + + +template +void check_quadrature (double*); +void check_GRC_right(double*); + + +int main() +{ + // this stores the exact values of \int_0^1 x^i/sqrt(x(1-x)) dx + static double exact_monomials[32]; + + exact_monomials[0] = 3.141592653589793; + exact_monomials[1] = 1.570796326794897; + exact_monomials[2] = 1.178097245096172; + exact_monomials[3] = 0.9817477042468104; + exact_monomials[4] = 0.8590292412159591; + exact_monomials[5] = 0.7731263170943632; + exact_monomials[6] = 0.7086991240031662; + exact_monomials[7] = 0.6580777580029401; + exact_monomials[8] = 0.6169478981277563; + exact_monomials[9] = 0.5826730148984365; + exact_monomials[10] = 0.5535393641535147; + exact_monomials[11] = 0.5283784839647186; + exact_monomials[12] = 0.5063627137995220; + exact_monomials[13] = 0.4868872248072327; + exact_monomials[14] = 0.4694983953498315; + exact_monomials[15] = 0.4538484488381705; + exact_monomials[16] = 0.4396656848119776; + exact_monomials[17] = 0.4267343411410371; + exact_monomials[18] = 0.4148806094426750; + exact_monomials[19] = 0.4039626986678677; + exact_monomials[20] = 0.3938636312011710; + exact_monomials[21] = 0.3844859256963813; + exact_monomials[22] = 0.3757476092032817; + exact_monomials[23] = 0.3675791829162538; + exact_monomials[24] = 0.3599212832721652; + exact_monomials[25] = 0.3527228576067219; + exact_monomials[26] = 0.3459397257296695; + exact_monomials[27] = 0.3395334345124534; + exact_monomials[28] = 0.3334703374675882; + exact_monomials[29] = 0.3277208488905608; + exact_monomials[30] = 0.3222588347423848; + exact_monomials[31] = 0.3170611116013786; + + + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + deallog.push("1d Gauss-Chebyshev\n"); + deallog.pop(); + check_quadrature,1>(&exact_monomials[0]); + + deallog.push("1d Gauss-Radau-Chebyshev, left endpoint\n"); + deallog.pop(); + check_quadrature,1>(&exact_monomials[0]); + + deallog.push("1d Gauss-Radau-Chebyshev, right endpoint\n"); + deallog.pop(); + check_GRC_right(&exact_monomials[0]); + + deallog.push("1d Gauss-Lobatto-Chebyshev\n"); + deallog.pop(); + check_quadrature,2>(&exact_monomials[0]); + + return 0; +} + + +template +void check_quadrature(double *exact_monomials) +{ + + for (unsigned int n=startn; n<18; ++n) + { + quadrature_type quadrature(n); + const std::vector > &points =quadrature.get_points(); + const std::vector &weights=quadrature.get_weights(); + + + for (unsigned int i=0; i<32; ++i) + { + double quadrature_int=0; + double err = 0; + + // Check the integral + // x^i/sqrt(x(1-x)) + long double f=1.; + for (unsigned int x=0; x quadrature(n,QGaussRadauChebyshev<1>::right); + const std::vector > &points=quadrature.get_points(); + const std::vector &weights=quadrature.get_weights(); + + + for (unsigned int i=0; i<32; ++i) + { + double quadrature_int=0; + double err = 0; + + // Check the integral + // x^i/sqrt(x(1-x)) + long double f=1.; + for (unsigned int x=0; x