--- /dev/null
+New: Added QGaussRadau quadrature up to and including 8 quadrature points.
+<br>
+(Jan Philipp Thiele, 2023/09/15)
\ No newline at end of file
QGauss(const unsigned int n);
};
+/**
+ * The Gauss-Radau family of quadrature rules for numerical integration.
+ *
+ * This modification of the Gauss quadrature uses one of the two interval end
+ * points as well. Being exact for polynomials of degree $2n-2$, this
+ * formula is suboptimal by one degree.
+ *
+ * This formula is often used in the context of discontinuous Galerkin
+ * discretizations of ODEs and the temporal part of PDEs.
+ *
+ * The quadrature points are the left interval end point plus the $n-1$
+ * roots of the polynomial
+ * \f[
+ * \frac{P_{n-1}(x)+P_n(x)}{1+x}
+ * \f]
+ * where $P_{n-1}$ and $P_n$ are Legendre polynomials.
+ * The quadrature weights are
+ * \f[
+ * w_0=\frac{2}{n^2}\quad\text{and}
+ * \quad w_i=\frac{1-x_i}{n^2(P_{n-1}(x_i))^2}\text{ for }i>0
+ * \f]
+ *
+ * For the right Gauss-Radau formula the quadrature points are
+ * $\tilde{x}_i=1-x_{n-i-1}$ and the weights are $\tilde{w}_i=w_{n-i-1}$,
+ * with $(x_i,w_i)$ as quadrature points
+ * and weights of the left Gauss-Radau formula.
+ *
+ * @see https://mathworld.wolfram.com/RadauQuadrature.html
+ */
+template <int dim>
+class QGaussRadau : public Quadrature<dim>
+{
+public:
+ /**
+ * EndPoint is used to specify which of the two endpoints of the unit interval
+ * is used also as quadrature point.
+ */
+ enum class EndPoint
+ {
+ /**
+ * Left end point.
+ */
+ left,
+ /**
+ * Right end point.
+ */
+ right
+ };
+ /**
+ * Generate a formula with <tt>n</tt> quadrature points (in each space
+ * direction).
+ * <tt>ep</tt> defines whether the left/lower/front endpoint(s) (default)
+ * or the right/upper/back endpoint(s) are part of the quadrature points.
+ */
+ QGaussRadau(const unsigned int n,
+ const EndPoint end_point = QGaussRadau::EndPoint::left);
+
+ /**
+ * Move constructor.
+ */
+ QGaussRadau(QGaussRadau<dim> &&) noexcept = default;
+
+private:
+ const EndPoint end_point;
+};
+
/**
* The Gauss-Lobatto family of quadrature rules for numerical integration.
template <>
QGauss<1>::QGauss(const unsigned int n);
template <>
+QGaussRadau<1>::QGaussRadau(const unsigned int n,
+ const QGaussRadau<1>::EndPoint end_poin);
+template <>
QGaussLobatto<1>::QGaussLobatto(const unsigned int n);
template <>
} // namespace internal
+namespace internal
+{
+ namespace QGaussRadau
+ {
+
+ // Implements lookup table after affine transformation to [0,1].
+ //
+ // Analytical values for [-1,1] and n < 4 listed on
+ // https://mathworld.wolfram.com/RadauQuadrature.html
+ // Values for n > 3 calculated with the Julia Package
+ // FastGaussQuadrature.jl
+ // https://github.com/JuliaApproximation/FastGaussQuadrature.jl
+ //
+ std::vector<double>
+ get_left_quadrature_points(const unsigned int n)
+ {
+ std::vector<double> q_points(n);
+ switch (n)
+ {
+ case 1:
+ q_points[0] = 0.;
+ break;
+ case 2:
+ q_points[0] = 0.;
+ q_points[1] = 2. / 3.;
+ break;
+ case 3:
+ q_points[0] = 0.;
+ q_points[1] = (6. - std::sqrt(6)) * 0.1;
+ q_points[2] = (6. + std::sqrt(6)) * 0.1;
+ break;
+
+ case 4:
+ q_points[0] = 0.000000000000000000;
+ q_points[1] = 0.212340538239152943;
+ q_points[2] = 0.590533135559265343;
+ q_points[3] = 0.911412040487296071;
+ break;
+ case 5:
+ q_points[0] = 0.000000000000000000;
+ q_points[1] = 0.139759864343780571;
+ q_points[2] = 0.416409567631083166;
+ q_points[3] = 0.723156986361876197;
+ q_points[4] = 0.942895803885482331;
+ break;
+ case 6:
+ q_points[0] = 0.000000000000000000;
+ q_points[1] = 0.098535085798826416;
+ q_points[2] = 0.304535726646363913;
+ q_points[3] = 0.562025189752613841;
+ q_points[4] = 0.801986582126391845;
+ q_points[5] = 0.960190142948531222;
+ break;
+ case 7:
+ q_points[0] = 0.000000000000000000;
+ q_points[1] = 0.073054328680258851;
+ q_points[2] = 0.230766137969945495;
+ q_points[3] = 0.441328481228449865;
+ q_points[4] = 0.663015309718845702;
+ q_points[5] = 0.851921400331515644;
+ q_points[6] = 0.970683572840215114;
+ break;
+ case 8:
+ q_points[0] = 0.000000000000000000;
+ q_points[1] = 0.056262560536922135;
+ q_points[2] = 0.180240691736892389;
+ q_points[3] = 0.352624717113169672;
+ q_points[4] = 0.547153626330555420;
+ q_points[5] = 0.734210177215410598;
+ q_points[6] = 0.885320946839095790;
+ q_points[7] = 0.977520613561287499;
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ break;
+ }
+ return q_points;
+ }
+
+ std::vector<double>
+ get_quadrature_points(const unsigned int n,
+ const ::dealii::QGaussRadau<1>::EndPoint end_point)
+ {
+ std::vector<double> left_points = get_left_quadrature_points(n);
+ switch (end_point)
+ {
+ case ::dealii::QGaussRadau<1>::EndPoint::left:
+ return left_points;
+ case ::dealii::QGaussRadau<1>::EndPoint::right:
+ {
+ std::vector<double> points(n);
+ for (unsigned int i = 0; i < n; ++i)
+ {
+ points[n - i - 1] = 1. - left_points[i];
+ }
+ return points;
+ }
+ default:
+ Assert(
+ false,
+ ExcMessage(
+ "This constructor can only be called with either "
+ "QGaussRadau::left or QGaussRadau::right as second argument."));
+ return {};
+ }
+ }
+
+ // Implements lookup table after affine transformation to [0,1].
+ //
+ // Analytical values for [-1,1] and n < 4 listed on
+ // https://mathworld.wolfram.com/RadauQuadrature.html
+ // Values for n > 3 calculated with the Julia Package
+ // FastGaussQuadrature.jl
+ // https://github.com/JuliaApproximation/FastGaussQuadrature.jl
+ //
+ std::vector<double>
+ get_left_quadrature_weights(const unsigned int n)
+ {
+ std::vector<double> weights(n);
+ switch (n)
+ {
+ case 1:
+ weights[0] = 1.;
+ break;
+ case 2:
+ weights[0] = 0.25;
+ weights[1] = 0.75;
+ break;
+ case 3:
+ weights[0] = 1. / 9.;
+ weights[1] = (16. + std::sqrt(6)) / 36.;
+ weights[2] = (16. - std::sqrt(6)) / 36.;
+ break;
+ case 4:
+ weights[0] = 0.062500000000000000;
+ weights[1] = 0.328844319980059696;
+ weights[2] = 0.388193468843171852;
+ weights[3] = 0.220462211176768369;
+ break;
+ case 5:
+ weights[0] = 0.040000000000000001;
+ weights[1] = 0.223103901083570894;
+ weights[2] = 0.311826522975741427;
+ weights[3] = 0.281356015149462124;
+ weights[4] = 0.143713560791225797;
+ break;
+ case 6:
+ weights[0] = 0.027777777777777776;
+ weights[1] = 0.159820376610255471;
+ weights[2] = 0.242693594234484888;
+ weights[3] = 0.260463391594787597;
+ weights[4] = 0.208450667155953895;
+ weights[5] = 0.100794192626740456;
+ break;
+ case 7:
+ weights[0] = 0.020408163265306121;
+ weights[1] = 0.119613744612656100;
+ weights[2] = 0.190474936822115581;
+ weights[3] = 0.223554914507283209;
+ weights[4] = 0.212351889502977870;
+ weights[5] = 0.159102115733650767;
+ weights[6] = 0.074494235556010341;
+ break;
+ case 8:
+ weights[0] = 0.015625000000000000;
+ weights[1] = 0.092679077401489660;
+ weights[2] = 0.152065310323392683;
+ weights[3] = 0.188258772694559262;
+ weights[4] = 0.195786083726246729;
+ weights[5] = 0.173507397817250691;
+ weights[6] = 0.124823950664932445;
+ weights[7] = 0.057254407372128648;
+ break;
+
+ default:
+ Assert(false, dealii::StandardExceptions::ExcNotImplemented());
+ break;
+ }
+ return weights;
+ }
+
+ std::vector<double>
+ get_quadrature_weights(const unsigned int n,
+ const ::dealii::QGaussRadau<1>::EndPoint end_point)
+ {
+ std::vector<double> left_weights = get_left_quadrature_weights(n);
+ switch (end_point)
+ {
+ case ::dealii::QGaussRadau<1>::EndPoint::left:
+ return left_weights;
+ case ::dealii::QGaussRadau<1>::EndPoint::right:
+ {
+ std::vector<double> weights(n);
+ for (unsigned int i = 0; i < n; ++i)
+ {
+ weights[n - i - 1] = left_weights[i];
+ }
+ return weights;
+ }
+ default:
+ Assert(false,
+ ExcMessage(
+ "This constructor can only be called with either "
+ "QGaussRadau::EndPoint::left or "
+ "QGaussRadau::EndPoint::right as second argument."));
+ return {};
+ }
+ }
+ } // namespace QGaussRadau
+} // namespace internal
+
+#ifndef DOXYGEN
+template <>
+QGaussRadau<1>::QGaussRadau(const unsigned int n, const EndPoint end_point)
+ : Quadrature<1>(n)
+ , end_point(end_point)
+{
+ Assert(n > 0, ExcMessage("Need at least one point for quadrature rules."));
+ std::vector<double> p =
+ internal::QGaussRadau::get_quadrature_points(n, end_point);
+ std::vector<double> w =
+ internal::QGaussRadau::get_quadrature_weights(n, end_point);
+
+ for (unsigned int i = 0; i < this->size(); ++i)
+ {
+ this->quadrature_points[i] = dealii::Point<1>(p[i]);
+ this->weights[i] = w[i];
+ }
+}
+#endif
+
+
#ifndef DOXYGEN
template <>
QGaussLobatto<1>::QGaussLobatto(const unsigned int n)
+template <int dim>
+QGaussRadau<dim>::QGaussRadau(const unsigned int n, EndPoint end_point)
+ : Quadrature<dim>(
+ QGaussRadau<1>(n, static_cast<QGaussRadau<1>::EndPoint>(end_point)))
+ , end_point(end_point)
+{}
+
+
+
template <int dim>
QGaussLobatto<dim>::QGaussLobatto(const unsigned int n)
: Quadrature<dim>(QGaussLobatto<dim - 1>(n), QGaussLobatto<1>(n))
// explicit specialization
// note that 1d formulae are specialized by implementation above
template class QGauss<2>;
+template class QGaussRadau<2>;
template class QGaussLobatto<2>;
template class QMidpoint<2>;
template class QTrapezoid<2>;
template class QWeddle<2>;
template class QGauss<3>;
+template class QGaussRadau<3>;
template class QGaussLobatto<3>;
template class QMidpoint<3>;
template class QTrapezoid<3>;
{
quadratures.push_back(new QGaussLobatto<dim>(i));
}
+ for (unsigned int i = 1; i < 9; ++i)
+ {
+ quadratures.push_back(new QGaussRadau<dim>(i));
+ quadratures.push_back(
+ new QGaussRadau<dim>(i, QGaussRadau<dim>::EndPoint::right));
+ }
}
template <int dim>
DEAL:1d::Quadrature no.18 is exact for polynomials of degree 7
DEAL:1d::Quadrature no.19 is exact for polynomials of degree 9
DEAL:1d::Quadrature no.20 is exact for polynomials of degree 11
+DEAL:1d::Quadrature no.21 is exact for polynomials of degree 0
+DEAL:1d::Quadrature no.22 is exact for polynomials of degree 0
+DEAL:1d::Quadrature no.23 is exact for polynomials of degree 2
+DEAL:1d::Quadrature no.24 is exact for polynomials of degree 2
+DEAL:1d::Quadrature no.25 is exact for polynomials of degree 4
+DEAL:1d::Quadrature no.26 is exact for polynomials of degree 4
+DEAL:1d::Quadrature no.27 is exact for polynomials of degree 6
+DEAL:1d::Quadrature no.28 is exact for polynomials of degree 6
+DEAL:1d::Quadrature no.29 is exact for polynomials of degree 8
+DEAL:1d::Quadrature no.30 is exact for polynomials of degree 8
+DEAL:1d::Quadrature no.31 is exact for polynomials of degree 10
+DEAL:1d::Quadrature no.32 is exact for polynomials of degree 10
+DEAL:1d::Quadrature no.33 is exact for polynomials of degree 12
+DEAL:1d::Quadrature no.34 is exact for polynomials of degree 12
+DEAL:1d::Quadrature no.35 is exact for polynomials of degree 14
+DEAL:1d::Quadrature no.36 is exact for polynomials of degree 14
DEAL:2d::Quadrature no.0 is exact for polynomials of degree 1
DEAL:2d::Quadrature no.1 is exact for polynomials of degree 1
DEAL:2d::Quadrature no.2 is exact for polynomials of degree 3
DEAL:2d::Quadrature no.18 is exact for polynomials of degree 7
DEAL:2d::Quadrature no.19 is exact for polynomials of degree 9
DEAL:2d::Quadrature no.20 is exact for polynomials of degree 11
+DEAL:2d::Quadrature no.21 is exact for polynomials of degree 0
+DEAL:2d::Quadrature no.22 is exact for polynomials of degree 0
+DEAL:2d::Quadrature no.23 is exact for polynomials of degree 2
+DEAL:2d::Quadrature no.24 is exact for polynomials of degree 2
+DEAL:2d::Quadrature no.25 is exact for polynomials of degree 4
+DEAL:2d::Quadrature no.26 is exact for polynomials of degree 4
+DEAL:2d::Quadrature no.27 is exact for polynomials of degree 6
+DEAL:2d::Quadrature no.28 is exact for polynomials of degree 6
+DEAL:2d::Quadrature no.29 is exact for polynomials of degree 8
+DEAL:2d::Quadrature no.30 is exact for polynomials of degree 8
+DEAL:2d::Quadrature no.31 is exact for polynomials of degree 10
+DEAL:2d::Quadrature no.32 is exact for polynomials of degree 10
+DEAL:2d::Quadrature no.33 is exact for polynomials of degree 12
+DEAL:2d::Quadrature no.34 is exact for polynomials of degree 12
+DEAL:2d::Quadrature no.35 is exact for polynomials of degree 14
+DEAL:2d::Quadrature no.36 is exact for polynomials of degree 14
DEAL:2d:faces::Quadrature no.0 is exact for polynomials of degree 1
DEAL:2d:faces::Quadrature no.1 is exact for polynomials of degree 1
DEAL:2d:faces::Quadrature no.2 is exact for polynomials of degree 3
DEAL:2d:faces::Quadrature no.18 is exact for polynomials of degree 7
DEAL:2d:faces::Quadrature no.19 is exact for polynomials of degree 9
DEAL:2d:faces::Quadrature no.20 is exact for polynomials of degree 11
+DEAL:2d:faces::Quadrature no.21 is exact for polynomials of degree 0
+DEAL:2d:faces::Quadrature no.22 is exact for polynomials of degree 0
+DEAL:2d:faces::Quadrature no.23 is exact for polynomials of degree 2
+DEAL:2d:faces::Quadrature no.24 is exact for polynomials of degree 2
+DEAL:2d:faces::Quadrature no.25 is exact for polynomials of degree 4
+DEAL:2d:faces::Quadrature no.26 is exact for polynomials of degree 4
+DEAL:2d:faces::Quadrature no.27 is exact for polynomials of degree 6
+DEAL:2d:faces::Quadrature no.28 is exact for polynomials of degree 6
+DEAL:2d:faces::Quadrature no.29 is exact for polynomials of degree 8
+DEAL:2d:faces::Quadrature no.30 is exact for polynomials of degree 8
+DEAL:2d:faces::Quadrature no.31 is exact for polynomials of degree 10
+DEAL:2d:faces::Quadrature no.32 is exact for polynomials of degree 10
+DEAL:2d:faces::Quadrature no.33 is exact for polynomials of degree 12
+DEAL:2d:faces::Quadrature no.34 is exact for polynomials of degree 12
+DEAL:2d:faces::Quadrature no.35 is exact for polynomials of degree 14
+DEAL:2d:faces::Quadrature no.36 is exact for polynomials of degree 14
DEAL:2d:subfaces::Quadrature no.0 is exact for polynomials of degree 1
DEAL:2d:subfaces::Quadrature no.1 is exact for polynomials of degree 1
DEAL:2d:subfaces::Quadrature no.2 is exact for polynomials of degree 3
DEAL:2d:subfaces::Quadrature no.18 is exact for polynomials of degree 7
DEAL:2d:subfaces::Quadrature no.19 is exact for polynomials of degree 9
DEAL:2d:subfaces::Quadrature no.20 is exact for polynomials of degree 11
+DEAL:2d:subfaces::Quadrature no.21 is exact for polynomials of degree 0
+DEAL:2d:subfaces::Quadrature no.22 is exact for polynomials of degree 0
+DEAL:2d:subfaces::Quadrature no.23 is exact for polynomials of degree 2
+DEAL:2d:subfaces::Quadrature no.24 is exact for polynomials of degree 2
+DEAL:2d:subfaces::Quadrature no.25 is exact for polynomials of degree 4
+DEAL:2d:subfaces::Quadrature no.26 is exact for polynomials of degree 4
+DEAL:2d:subfaces::Quadrature no.27 is exact for polynomials of degree 6
+DEAL:2d:subfaces::Quadrature no.28 is exact for polynomials of degree 6
+DEAL:2d:subfaces::Quadrature no.29 is exact for polynomials of degree 8
+DEAL:2d:subfaces::Quadrature no.30 is exact for polynomials of degree 8
+DEAL:2d:subfaces::Quadrature no.31 is exact for polynomials of degree 10
+DEAL:2d:subfaces::Quadrature no.32 is exact for polynomials of degree 10
+DEAL:2d:subfaces::Quadrature no.33 is exact for polynomials of degree 12
+DEAL:2d:subfaces::Quadrature no.34 is exact for polynomials of degree 12
+DEAL:2d:subfaces::Quadrature no.35 is exact for polynomials of degree 14
+DEAL:2d:subfaces::Quadrature no.36 is exact for polynomials of degree 14
DEAL:3d::Quadrature no.0 is exact for polynomials of degree 1
DEAL:3d::Quadrature no.1 is exact for polynomials of degree 1
DEAL:3d::Quadrature no.2 is exact for polynomials of degree 3
DEAL:3d::Quadrature no.18 is exact for polynomials of degree 7
DEAL:3d::Quadrature no.19 is exact for polynomials of degree 9
DEAL:3d::Quadrature no.20 is exact for polynomials of degree 11
+DEAL:3d::Quadrature no.21 is exact for polynomials of degree 0
+DEAL:3d::Quadrature no.22 is exact for polynomials of degree 0
+DEAL:3d::Quadrature no.23 is exact for polynomials of degree 2
+DEAL:3d::Quadrature no.24 is exact for polynomials of degree 2
+DEAL:3d::Quadrature no.25 is exact for polynomials of degree 4
+DEAL:3d::Quadrature no.26 is exact for polynomials of degree 4
+DEAL:3d::Quadrature no.27 is exact for polynomials of degree 6
+DEAL:3d::Quadrature no.28 is exact for polynomials of degree 6
+DEAL:3d::Quadrature no.29 is exact for polynomials of degree 8
+DEAL:3d::Quadrature no.30 is exact for polynomials of degree 8
+DEAL:3d::Quadrature no.31 is exact for polynomials of degree 10
+DEAL:3d::Quadrature no.32 is exact for polynomials of degree 10
+DEAL:3d::Quadrature no.33 is exact for polynomials of degree 12
+DEAL:3d::Quadrature no.34 is exact for polynomials of degree 12
+DEAL:3d::Quadrature no.35 is exact for polynomials of degree 14
+DEAL:3d::Quadrature no.36 is exact for polynomials of degree 14
DEAL:3d:faces::Quadrature no.0 is exact for polynomials of degree 1
DEAL:3d:faces::Quadrature no.1 is exact for polynomials of degree 1
DEAL:3d:faces::Quadrature no.2 is exact for polynomials of degree 3
DEAL:3d:faces::Quadrature no.18 is exact for polynomials of degree 7
DEAL:3d:faces::Quadrature no.19 is exact for polynomials of degree 9
DEAL:3d:faces::Quadrature no.20 is exact for polynomials of degree 11
+DEAL:3d:faces::Quadrature no.21 is exact for polynomials of degree 1
+DEAL:3d:faces::Quadrature no.22 is exact for polynomials of degree 1
+DEAL:3d:faces::Quadrature no.23 is exact for polynomials of degree 3
+DEAL:3d:faces::Quadrature no.24 is exact for polynomials of degree 3
+DEAL:3d:faces::Quadrature no.25 is exact for polynomials of degree 5
+DEAL:3d:faces::Quadrature no.26 is exact for polynomials of degree 5
+DEAL:3d:faces::Quadrature no.27 is exact for polynomials of degree 7
+DEAL:3d:faces::Quadrature no.28 is exact for polynomials of degree 7
+DEAL:3d:faces::Quadrature no.29 is exact for polynomials of degree 9
+DEAL:3d:faces::Quadrature no.30 is exact for polynomials of degree 9
+DEAL:3d:faces::Quadrature no.31 is exact for polynomials of degree 11
+DEAL:3d:faces::Quadrature no.32 is exact for polynomials of degree 11
+DEAL:3d:faces::Quadrature no.33 is exact for polynomials of degree 13
+DEAL:3d:faces::Quadrature no.34 is exact for polynomials of degree 13
+DEAL:3d:faces::Quadrature no.35 is exact for polynomials of degree 15
+DEAL:3d:faces::Quadrature no.36 is exact for polynomials of degree 15
DEAL:3d:subfaces::Quadrature no.0 is exact for polynomials of degree 1
DEAL:3d:subfaces::Quadrature no.1 is exact for polynomials of degree 1
DEAL:3d:subfaces::Quadrature no.2 is exact for polynomials of degree 3
DEAL:3d:subfaces::Quadrature no.18 is exact for polynomials of degree 7
DEAL:3d:subfaces::Quadrature no.19 is exact for polynomials of degree 9
DEAL:3d:subfaces::Quadrature no.20 is exact for polynomials of degree 11
+DEAL:3d:subfaces::Quadrature no.21 is exact for polynomials of degree 1
+DEAL:3d:subfaces::Quadrature no.22 is exact for polynomials of degree 1
+DEAL:3d:subfaces::Quadrature no.23 is exact for polynomials of degree 3
+DEAL:3d:subfaces::Quadrature no.24 is exact for polynomials of degree 3
+DEAL:3d:subfaces::Quadrature no.25 is exact for polynomials of degree 5
+DEAL:3d:subfaces::Quadrature no.26 is exact for polynomials of degree 5
+DEAL:3d:subfaces::Quadrature no.27 is exact for polynomials of degree 7
+DEAL:3d:subfaces::Quadrature no.28 is exact for polynomials of degree 7
+DEAL:3d:subfaces::Quadrature no.29 is exact for polynomials of degree 9
+DEAL:3d:subfaces::Quadrature no.30 is exact for polynomials of degree 9
+DEAL:3d:subfaces::Quadrature no.31 is exact for polynomials of degree 11
+DEAL:3d:subfaces::Quadrature no.32 is exact for polynomials of degree 11
+DEAL:3d:subfaces::Quadrature no.33 is exact for polynomials of degree 13
+DEAL:3d:subfaces::Quadrature no.34 is exact for polynomials of degree 13
+DEAL:3d:subfaces::Quadrature no.35 is exact for polynomials of degree 15
+DEAL:3d:subfaces::Quadrature no.36 is exact for polynomials of degree 15
projects / dealii.git / commitdiff