From: Martin Kronbichler Date: Tue, 5 Apr 2016 12:49:54 +0000 (+0200) Subject: Make Gauss-Lobatto points the default node distribution in FE_Q/FE_DGQ. X-Git-Tag: v8.5.0-rc1~1119^2~7 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2034c8d7886febd6dfb1f906f5cf1048330cf17e;p=dealii.git Make Gauss-Lobatto points the default node distribution in FE_Q/FE_DGQ. --- diff --git a/include/deal.II/fe/fe_dgq.h b/include/deal.II/fe/fe_dgq.h index 2c44352cae..e8beccbfc2 100644 --- a/include/deal.II/fe/fe_dgq.h +++ b/include/deal.II/fe/fe_dgq.h @@ -75,6 +75,30 @@ template class Quadrature; * element are exactly the same as those of the FE_Q element where they are * shown visually. * + *

Unit support point distribution and conditioning of interpolation

+ * + * When constructing an FE_DGQ element at polynomial degrees one or two, + * equidistant support points at 0 and 1 (linear case) or 0, 0.5, and 1 + * (quadratic case) are used. The unit support or nodal points + * xi are those points where the jth Lagrange + * polynomial satisfies the $\delta_{ij}$ property, i.e., where one polynomial + * is one and all the others are zero. For higher polynomial degrees, the + * support points are non-equidistant by default, and chosen to be the support + * points of the (degree+1)-order Gauss-Lobatto quadrature rule. This + * point distribution yields well-conditioned Lagrange interpolation at + * arbitrary polynomial degrees. By contrast, polynomials based on equidistant + * points get increasingly ill-conditioned as the polynomial degree + * increases. In interpolation, this effect is known as the Runge + * phenomenon. For Galerkin methods, the Runge phenomenon is typically not + * visible in the solution quality but rather in the condition number of the + * associated system matrices. For example, the elemental mass matrix of + * equidistant points at degree 10 has condition number 2.6e6, whereas the + * condition number for Gauss-Lobatto points is around 400. + * + * The Gauss-Lobatto points in 1D include the end points 0 and +1 of the unit + * interval. The interior points are shifted towards the end points, which + * gives a denser point distribution close to the element boundary. + * * @author Ralf Hartmann, Guido Kanschat 2001, 2004 */ template @@ -84,7 +108,7 @@ public: /** * Constructor for tensor product polynomials of degree p. The * shape functions created using this constructor correspond to Lagrange - * interpolation polynomials for equidistantly spaced support points in each + * interpolation polynomials for Gauss-Lobatto support (node) points in each * coordinate direction. */ FE_DGQ (const unsigned int p); @@ -360,7 +384,8 @@ private: * quadrature points. If this set of quadrature points is then also used in * integrating the mass matrix, then it will be diagonal. The number of * quadrature points automatically determines the polynomial degree chosen for - * this element. + * this element. The typical applications are the Gauss quadrature or the + * Gauss-Lobatto quadrature (provided through the base class). * * See the base class documentation in FE_DGQ for details. * diff --git a/include/deal.II/fe/fe_q.h b/include/deal.II/fe/fe_q.h index f47cba2e7c..c821dc8f72 100644 --- a/include/deal.II/fe/fe_q.h +++ b/include/deal.II/fe/fe_q.h @@ -30,7 +30,8 @@ DEAL_II_NAMESPACE_OPEN * Implementation of a scalar Lagrange finite element @p Qp that yields the * finite element space of continuous, piecewise polynomials of degree @p p in * each coordinate direction. This class is realized using tensor product - * polynomials based on equidistant or given support points. + * polynomials based on 1D Lagrange polynomials with equidistant (degree up to + * 2), Gauss-Lobatto (starting from degree 3), or given support points. * * The standard constructor of this class takes the degree @p p of this finite * element. Alternatively, it can take a quadrature formula @p points defining @@ -54,6 +55,35 @@ DEAL_II_NAMESPACE_OPEN * implemented only up to a certain degree and may not be available for very * high polynomial degree. * + *

Unit support point distribution and conditioning of interpolation

+ * + * When constructing an FE_Q element at polynomial degrees one or two, + * equidistant support points at 0 and 1 (linear case) or 0, 0.5, and 1 + * (quadratic case) are used. The unit support or nodal points + * xi are those points where the jth Lagrange + * polynomial satisfies the $\delta_{ij}$ property, i.e., where one polynomial + * is one and all the others are zero. For higher polynomial degrees, the + * support points are non-equidistant by default, and chosen to be the support + * points of the (degree+1)-order Gauss-Lobatto quadrature rule. This + * point distribution yields well-conditioned Lagrange interpolation at + * arbitrary polynomial degrees. By contrast, polynomials based on equidistant + * points get increasingly ill-conditioned as the polynomial degree + * increases. In interpolation, this effect is known as the Runge + * phenomenon. For Galerkin methods, the Runge phenomenon is typically not + * visible in the solution quality but rather in the condition number of the + * associated system matrices. For example, the elemental mass matrix of + * equidistant points at degree 10 has condition number 2.6e6, whereas the + * condition number for Gauss-Lobatto points is around 400. + * + * The Gauss-Lobatto points in 1D include the end points 0 and +1 of the unit + * interval. The interior points are shifted towards the end points, which + * gives a denser point distribution close to the element boundary. + * + * If combined with Gauss-Lobatto quadrature, FE_Q based on the default + * support points gives diagonal mass matrices. This case is demonstrated in + * step-48. However, this element can be combined with arbitrary quadrature + * rules through the usual FEValues approach, including full Gauss + * quadrature. In the general case, the mass matrix is non-diagonal. * *

Numbering of the degrees of freedom (DoFs)

* @@ -523,19 +553,21 @@ class FE_Q : public FE_Q_Base,dim,spacedim> { public: /** - * Constructor for tensor product polynomials of degree @p p. + * Constructor for tensor product polynomials of degree @p p based on + * Gauss-Lobatto support (node) points. For polynomial degrees of one and + * two, these are the usual equidistant points. */ FE_Q (const unsigned int p); /** * Constructor for tensor product polynomials with support points @p points * based on a one-dimensional quadrature formula. The degree of the finite - * element is points.size()-1. Note that the first point has to be - * 0 and the last one 1. If - * FE_Q(QGaussLobatto<1>(fe_degree+1)) is specified, so-called - * Gauss-Lobatto elements are obtained which can give a diagonal mass matrix - * if combined with Gauss-Lobatto quadrature on the same points. Their use - * is shown in step-48. + * element is points.size()-1. Note that the first point has to be + * 0 and the last one 1. Constructing + * FE_Q(QGaussLobatto<1>(fe_degree+1)) is equivalent to the + * constructor that specifies the polynomial degree only. For selecting + * equidistant nodes at fe_degree > 2, construct + * FE_Q(QIterated<1>(QTrapez<1>(),fe_degree)). */ FE_Q (const Quadrature<1> &points); diff --git a/source/fe/fe_dgq.cc b/source/fe/fe_dgq.cc index e8206eefd7..12c1973c4c 100644 --- a/source/fe/fe_dgq.cc +++ b/source/fe/fe_dgq.cc @@ -16,7 +16,6 @@ #include #include -#include #include #include #include @@ -28,155 +27,35 @@ DEAL_II_NAMESPACE_OPEN -// namespace for some functions that are used in this file. they are -// specific to numbering conventions used for the FE_DGQ element, and -// are thus not very interesting to the outside world namespace { - // given an integer N, compute its - // integer square root (if it - // exists, otherwise give up) - inline unsigned int int_sqrt (const unsigned int N) + std::vector > + get_QGaussLobatto_points (const unsigned int degree) { - for (unsigned int i=0; i<=N; ++i) - if (i*i == N) - return i; - Assert (false, ExcInternalError()); - return numbers::invalid_unsigned_int; - } - - - // given an integer N, compute its - // integer cube root (if it - // exists, otherwise give up) - inline unsigned int int_cuberoot (const unsigned int N) - { - for (unsigned int i=0; i<=N; ++i) - if (i*i*i == N) - return i; - Assert (false, ExcInternalError()); - return numbers::invalid_unsigned_int; - } - - - // given N, generate i=0...N-1 - // equidistant points in the - // interior of the interval [0,1] - inline Point<1> - generate_unit_point (const unsigned int i, - const unsigned int N, - const dealii::internal::int2type<1> ) - { - Assert (i (.5); - else - { - const double h = 1./(N-1); - return Point<1>(i*h); - } - } - - - // given N, generate i=0...N-1 - // equidistant points in the domain - // [0,1]^2 - inline Point<2> - generate_unit_point (const unsigned int i, - const unsigned int N, - const dealii::internal::int2type<2> ) - { - Assert (i (.5, .5); - else - { - Assert (N>=4, ExcInternalError()); - const unsigned int N1d = int_sqrt(N); - const double h = 1./(N1d-1); - - return Point<2> (i%N1d * h, - i/N1d * h); - } - } - - - - - // given N, generate i=0...N-1 - // equidistant points in the domain - // [0,1]^3 - inline Point<3> - generate_unit_point (const unsigned int i, - const unsigned int N, - const dealii::internal::int2type<3> ) - { - Assert (i (.5, .5, .5); + if (degree > 0) + return QGaussLobatto<1>(degree+1).get_points(); else - { - Assert (N>=8, ExcInternalError()); - - const unsigned int N1d = int_cuberoot(N); - const double h = 1./(N1d-1); - - return Point<3> (i%N1d * h, - (i/N1d)%N1d * h, - i/(N1d*N1d) * h); - } + return std::vector >(1, Point<1>(0.5)); } } - template FE_DGQ::FE_DGQ (const unsigned int degree) : - FE_Poly, dim, spacedim> ( - TensorProductPolynomials(Polynomials::LagrangeEquidistant::generate_complete_basis(degree)), - FiniteElementData(get_dpo_vector(degree), 1, degree, FiniteElementData::L2), - std::vector(FiniteElementData(get_dpo_vector(degree),1, degree).dofs_per_cell, true), - std::vector(FiniteElementData( - get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector(1,true))) + FE_Poly, dim, spacedim> + (TensorProductPolynomials(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))), + FiniteElementData(get_dpo_vector(degree), 1, degree, FiniteElementData::L2), + std::vector(FiniteElementData(get_dpo_vector(degree),1, degree).dofs_per_cell, true), + std::vector(FiniteElementData(get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector(1,true))) { - // fill in support points - if (degree == 0) - { - // constant elements, take - // midpoint - this->unit_support_points.resize(1); - for (unsigned int i=0; iunit_support_points[0](i) = 0.5; - } - else - { - // number of points: (degree+1)^dim - unsigned int n = degree+1; - for (unsigned int i=1; iunit_support_points.resize(n); - - const double step = 1./degree; - Point p; - - unsigned int k=0; - for (unsigned int iz=0; iz <= ((dim>2) ? degree : 0) ; ++iz) - for (unsigned int iy=0; iy <= ((dim>1) ? degree : 0) ; ++iy) - for (unsigned int ix=0; ix<=degree; ++ix) - { - p(0) = ix * step; - if (dim>1) - p(1) = iy * step; - if (dim>2) - p(2) = iz * step; - - this->unit_support_points[k++] = p; - }; - }; + // Compute support points, which are the tensor product of the Lagrange + // interpolation points in the constructor. + Quadrature support_quadrature(get_QGaussLobatto_points(degree)); + Assert (support_quadrature.get_points().size() > 0, + (typename FiniteElement::ExcFEHasNoSupportPoints ())); + this->unit_support_points = support_quadrature.get_points(); // do not initialize embedding and restriction here. these matrices are // initialized on demand in get_restriction_matrix and @@ -194,13 +73,10 @@ FE_DGQ::FE_DGQ (const Quadrature<1> &points) TensorProductPolynomials(Polynomials::generate_complete_Lagrange_basis(points.get_points())), FiniteElementData(get_dpo_vector(points.size()-1), 1, points.size()-1, FiniteElementData::L2), std::vector(FiniteElementData(get_dpo_vector(points.size()-1),1, points.size()-1).dofs_per_cell, true), - std::vector(FiniteElementData( - get_dpo_vector(points.size()-1),1, points.size()-1).dofs_per_cell, std::vector(1,true))) + std::vector(FiniteElementData(get_dpo_vector(points.size()-1),1, points.size()-1).dofs_per_cell, std::vector(1,true))) { - // Compute support points, which - // are the tensor product of the - // Lagrange interpolation points in - // the constructor. + // Compute support points, which are the tensor product of the Lagrange + // interpolation points in the constructor. Assert (points.size() > 0, (typename FiniteElement::ExcFEHasNoSupportPoints ())); Quadrature support_quadrature(points); @@ -218,12 +94,9 @@ template std::string FE_DGQ::get_name () const { - // note that the - // FETools::get_fe_from_name - // function depends on the - // particular format of the string - // this function returns, so they - // have to be kept in synch + // note that the FETools::get_fe_from_name function depends on the + // particular format of the string this function returns, so they have to be + // kept in sync std::ostringstream namebuf; namebuf << "FE_DGQ<" @@ -378,8 +251,7 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, // generate a point on this // cell and evaluate the // shape functions there - const Point p = generate_unit_point (j, this->dofs_per_cell, - dealii::internal::int2type()); + const Point p = this->unit_support_points[j]; for (unsigned int i=0; idofs_per_cell; ++i) cell_interpolation(j,i) = this->poly_space.compute_value (i, p); @@ -807,15 +679,20 @@ FE_DGQArbitraryNodes::get_name () const // Check whether the support points are equidistant. for (unsigned int j=0; j<=this->degree; j++) - if (std::fabs(points[j] - (double)j/this->degree) > 1e-15) + if (std::abs(points[j] - (double)j/this->degree) > 1e-15) { equidistant = false; break; } + if (this->degree == 0 && std::abs(points[0]-0.5) < 1e-15) + equidistant = true; if (equidistant == true) { - namebuf << "FE_DGQ<" << Utilities::dim_string(dim,spacedim) << ">(" << this->degree << ")"; + if (this->degree > 2) + namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim,spacedim) << ">(QIterated(QTrapez()," << this->degree << "))"; + else + namebuf << "FE_DGQ<" << Utilities::dim_string(dim,spacedim) << ">(" << this->degree << ")"; return namebuf.str(); } @@ -831,7 +708,7 @@ FE_DGQArbitraryNodes::get_name () const if (gauss_lobatto == true) { - namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim,spacedim) << ">(QGaussLobatto(" << this->degree+1 << "))"; + namebuf << "FE_DGQ<" << Utilities::dim_string(dim,spacedim) << ">(" << this->degree << ")"; return namebuf.str(); } diff --git a/source/fe/fe_q.cc b/source/fe/fe_q.cc index d170222b30..a4c969dee5 100644 --- a/source/fe/fe_q.cc +++ b/source/fe/fe_q.cc @@ -23,27 +23,36 @@ DEAL_II_NAMESPACE_OPEN +namespace +{ + std::vector > + get_QGaussLobatto_points (const unsigned int degree) + { + if (degree > 0) + return QGaussLobatto<1>(degree+1).get_points(); + else + AssertThrow(false, + ExcMessage ("FE_Q can only be used for polynomial degrees " + "greater than zero. If you want an element of polynomial " + "degree zero, then it cannot be continuous and you " + "will want to use FE_DGQ(0).")); + return std::vector >(); + } +} + + template FE_Q::FE_Q (const unsigned int degree) : - FE_Q_Base, dim, spacedim> ( - TensorProductPolynomials(Polynomials::LagrangeEquidistant::generate_complete_basis(degree)), - FiniteElementData(this->get_dpo_vector(degree), - 1, degree, - FiniteElementData::H1), - std::vector (1, false)) + FE_Q_Base, dim, spacedim> + (TensorProductPolynomials(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))), + FiniteElementData(this->get_dpo_vector(degree), + 1, degree, + FiniteElementData::H1), + std::vector (1, false)) { - Assert (degree > 0, - ExcMessage ("This element can only be used for polynomial degrees " - "greater than zero. If you want an element of polynomial " - "degree zero, then it cannot be continuous and you " - "will want to use FE_DGQ(0).")); - std::vector > support_points_1d(degree+1); - for (unsigned int i=0; i<=degree; ++i) - support_points_1d[i][0] = static_cast(i)/degree; - - this->initialize(support_points_1d); + this->initialize(get_QGaussLobatto_points(degree)); } @@ -95,9 +104,16 @@ FE_Q::get_name () const } if (equidistant == true) - namebuf << "FE_Q<" - << Utilities::dim_string(dim,spacedim) - << ">(" << this->degree << ")"; + { + if (this->degree > 2) + namebuf << "FE_Q<" + << Utilities::dim_string(dim,spacedim) + << ">(QIterated(QTrapez()," << this->degree << "))"; + else + namebuf << "FE_Q<" + << Utilities::dim_string(dim,spacedim) + << ">(" << this->degree << ")"; + } else { // Check whether the support points come from QGaussLobatto. @@ -112,7 +128,7 @@ FE_Q::get_name () const if (gauss_lobatto == true) namebuf << "FE_Q<" << Utilities::dim_string(dim,spacedim) - << ">(QGaussLobatto(" << this->degree+1 << "))"; + << ">(" << this->degree << ")"; else namebuf << "FE_Q<" << Utilities::dim_string(dim,spacedim) diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index 9137d03e37..a521c6c5b8 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -1180,12 +1180,22 @@ FE_Q_Base { Assert (std::fabs (1.-this->poly_space.compute_value (i, this->unit_support_points[i])) < eps, - ExcInternalError()); + ExcInternalError("The Lagrange polynomial does not evaluate " + "to one or zero in a nodal point. " + "This typically indicates that the " + "polynomial interpolation is " + "ill-conditioned such that round-off " + "prevents the sum to be one.")); for (unsigned int j=0; jpoly_space.compute_value (i, this->unit_support_points[j])) < eps, - ExcInternalError()); + ExcInternalError("The Lagrange polynomial does not evaluate " + "to one or zero in a nodal point. " + "This typically indicates that the " + "polynomial interpolation is " + "ill-conditioned such that round-off " + "prevents the sum to be one.")); } #endif @@ -1301,7 +1311,14 @@ FE_Q_Base double sum = 0; for (unsigned int col=0; coldofs_per_cell; ++col) sum += prolongate(row,col); - Assert (std::fabs(sum-1.) < eps, ExcInternalError()); + Assert (std::fabs(sum-1.) < + std::max(eps, 5e-16*std::sqrt(this->dofs_per_cell)), + ExcInternalError("The entries in a row of the local " + "prolongation matrix do not add to one. " + "This typically indicates that the " + "polynomial interpolation is " + "ill-conditioned such that round-off " + "prevents the sum to be one.")); } #endif @@ -1343,27 +1360,23 @@ FE_Q_Base // distinguish q/q_dg0 case const unsigned int q_dofs_per_cell = Utilities::fixed_power(q_degree+1); - // for these Lagrange interpolation polynomials, construction of the - // restriction matrices is relatively simple. the reason is that the - // interpolation points on the mother cell are (except for the case with - // arbitrary nonequidistant nodes) always also interpolation points for - // some shape function on one or the other child, because we have chosen - // equidistant Lagrange interpolation points for the polynomials + // for Lagrange interpolation polynomials based on equidistant points, + // construction of the restriction matrices is relatively simple. the + // reason is that in this case the interpolation points on the mother + // cell are always also interpolation points for some shape function on + // one or the other child. // - // so the only thing we have to find out is: for each shape function on - // the mother cell, which is the child cell (possibly more than one) on - // which it is located, and which is the corresponding shape function - // there. rather than doing it for the shape functions on the mother - // cell, we take the interpolation points there + // in the general case with non-equidistant points, we need to actually + // do an interpolation. thus, we take the interpolation points on the + // mother cell and evaluate the shape functions of the child cell on + // those points. it does not hurt in the equidistant case because then + // simple one shape function evaluates to one and the others to zero. // - // note that the interpolation point of a shape function can be on the - // boundary between subcells. in that case, restriction from children to - // mother may produce two or more entries for a dof on the mother - // cell. however, this doesn't hurt: since the element is continuous, - // the contribution from each child should yield the same result, and - // since the element is non-additive we just overwrite one value - // (compute on one child) by the same value (compute on a later child), - // so we don't have to care about this + // this element is non-additive in all its degrees of freedom by + // default, which requires care in downstream use. fortunately, even the + // interpolation on non-equidistant points is invariant under the + // assumption that whenever a row makes a non-zero contribution to the + // mother's residual, the correct value is interpolated. const double eps = 1e-15*q_degree*dim; const std::vector &index_map_inverse = @@ -1427,8 +1440,14 @@ FE_Q_Base } FE_Q_Helper::increment_indices (j_indices, dofs1d); } - Assert (std::fabs(sum_check-1) < eps, - ExcInternalError()); + Assert (std::fabs(sum_check-1) < + std::max(eps, 5e-16*std::sqrt(this->dofs_per_cell)), + ExcInternalError("The entries in a row of the local " + "restriction matrix do not add to one. " + "This typically indicates that the " + "polynomial interpolation is " + "ill-conditioned such that round-off " + "prevents the sum to be one.")); } // part for FE_Q_DG0 diff --git a/source/fe/fe_q_bubbles.cc b/source/fe/fe_q_bubbles.cc index 44b7e5fc6c..9868ae39ab 100644 --- a/source/fe/fe_q_bubbles.cc +++ b/source/fe/fe_q_bubbles.cc @@ -172,7 +172,7 @@ template FE_Q_Bubbles::FE_Q_Bubbles (const unsigned int q_degree) : FE_Q_Base, dim, spacedim> ( - TensorProductPolynomialsBubbles(Polynomials::LagrangeEquidistant::generate_complete_basis(q_degree)), + TensorProductPolynomialsBubbles(Polynomials::generate_complete_Lagrange_basis(QGaussLobatto<1>(q_degree+1).get_points())), FiniteElementData(get_dpo_vector(q_degree), 1, q_degree+1, FiniteElementData::H1), @@ -183,11 +183,7 @@ FE_Q_Bubbles::FE_Q_Bubbles (const unsigned int q_degree) ExcMessage ("This element can only be used for polynomial degrees " "greater than zero")); - std::vector > support_points_1d(q_degree+1); - for (unsigned int i=0; i<=q_degree; ++i) - support_points_1d[i][0] = static_cast(i)/q_degree; - - this->initialize(support_points_1d); + this->initialize(QGaussLobatto<1>(q_degree+1).get_points()); // adjust unit support point for discontinuous node Point point; @@ -221,9 +217,7 @@ FE_Q_Bubbles::FE_Q_Bubbles (const Quadrature<1> &points) get_riaf_vector(points.size()-1)), n_bubbles((points.size()-1<=1)?1:dim) { - const unsigned int q_degree = points.size()-1; - (void) q_degree; - Assert (q_degree > 0, + Assert (points.size() > 1, ExcMessage ("This element can only be used for polynomial degrees " "at least one")); @@ -294,7 +288,16 @@ FE_Q_Bubbles::get_name () const } if (type == true) - namebuf << "FE_Q_Bubbles<" << dim << ">(" << this->degree-1 << ")"; + { + if (this->degree > 3) + namebuf << "FE_Q_Bubbles<" + << Utilities::dim_string(dim,spacedim) + << ">(QIterated(QTrapez()," << this->degree-1 << "))"; + else + namebuf << "FE_Q_Bubbles<" + << Utilities::dim_string(dim,spacedim) + << ">(" << this->degree-1 << ")"; + } else { @@ -308,9 +311,9 @@ FE_Q_Bubbles::get_name () const break; } if (type == true) - namebuf << "FE_Q_Bubbles<" << dim << ">(QGaussLobatto(" << this->degree << "))"; + namebuf << "FE_Q_Bubbles<" << dim << ">(" << this->degree-1 << ")"; else - namebuf << "FE_Q_Bubbles<" << dim << ">(QUnknownNodes(" << this->degree-1 << "))"; + namebuf << "FE_Q_Bubbles<" << dim << ">(QUnknownNodes(" << this->degree << "))"; } return namebuf.str(); } diff --git a/source/fe/fe_q_dg0.cc b/source/fe/fe_q_dg0.cc index d8cb1c808a..d2218eadc9 100644 --- a/source/fe/fe_q_dg0.cc +++ b/source/fe/fe_q_dg0.cc @@ -33,7 +33,7 @@ template FE_Q_DG0::FE_Q_DG0 (const unsigned int degree) : FE_Q_Base, dim, spacedim> ( - TensorProductPolynomialsConst(Polynomials::LagrangeEquidistant::generate_complete_basis(degree)), + TensorProductPolynomialsConst(Polynomials::generate_complete_Lagrange_basis(QGaussLobatto<1>(degree+1).get_points())), FiniteElementData(get_dpo_vector(degree), 1, degree, FiniteElementData::L2), @@ -43,11 +43,7 @@ FE_Q_DG0::FE_Q_DG0 (const unsigned int degree) ExcMessage ("This element can only be used for polynomial degrees " "greater than zero")); - std::vector > support_points_1d(degree+1); - for (unsigned int i=0; i<=degree; ++i) - support_points_1d[i][0] = static_cast(i)/degree; - - this->initialize(support_points_1d); + this->initialize(QGaussLobatto<1>(degree+1).get_points()); // adjust unit support point for discontinuous node Point point; @@ -133,9 +129,16 @@ FE_Q_DG0::get_name () const } if (type == true) - namebuf << "FE_Q_DG0<" - << Utilities::dim_string(dim,spacedim) - << ">(" << this->degree << ")"; + { + if (this->degree > 2) + namebuf << "FE_Q_DG0<" + << Utilities::dim_string(dim,spacedim) + << ">(QIterated(QTrapez()," << this->degree << "))"; + else + namebuf << "FE_Q_DG0<" + << Utilities::dim_string(dim,spacedim) + << ">(" << this->degree << ")"; + } else { @@ -151,7 +154,7 @@ FE_Q_DG0::get_name () const if (type == true) namebuf << "FE_Q_DG0<" << Utilities::dim_string(dim,spacedim) - << ">(QGaussLobatto(" << this->degree+1 << "))"; + << ">(" << this->degree << ")"; else namebuf << "FE_Q_DG0<" << Utilities::dim_string(dim,spacedim) diff --git a/source/fe/fe_tools.cc b/source/fe/fe_tools.cc index eee128e367..dc82e00ba6 100644 --- a/source/fe/fe_tools.cc +++ b/source/fe/fe_tools.cc @@ -1540,6 +1540,23 @@ namespace FETools const FEFactoryBase *fef=dynamic_cast*>(ptr); return fef->get(QGauss<1>(tmp.first)); } + else if (quadrature_name.compare("QIterated") == 0) + { + // find sub-quadrature + position = name.find('('); + const std::string subquadrature_name(name, 0, position); + AssertThrow(subquadrature_name.compare("QTrapez") == 0, + ExcNotImplemented("Could not detect quadrature of name " + subquadrature_name)); + // delete "QTrapez()," + name.erase(0,position+3); + const std::pair tmp + = Utilities::get_integer_at_position (name, 0); + // delete "))" + name.erase(0, tmp.second+2); + const Subscriptor *ptr = fe_name_map.find(name_part)->second.get(); + const FEFactoryBase *fef=dynamic_cast*>(ptr); + return fef->get(QIterated<1>(QTrapez<1>(),tmp.first)); + } else { AssertThrow (false,ExcNotImplemented());