* element are exactly the same as those of the FE_Q element where they are
* shown visually.
*
+ * <h3>Unit support point distribution and conditioning of interpolation</h3>
+ *
+ * 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
+ * <i>x<sub>i</sub></i> are those points where the <i>j</i>th 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 <tt>(degree+1)</tt>-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 <int dim, int spacedim=dim>
/**
* Constructor for tensor product polynomials of degree <tt>p</tt>. 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);
* 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.
*
* 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
* implemented only up to a certain degree and may not be available for very
* high polynomial degree.
*
+ * <h3>Unit support point distribution and conditioning of interpolation</h3>
+ *
+ * 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
+ * <i>x<sub>i</sub></i> are those points where the <i>j</i>th 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 <tt>(degree+1)</tt>-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.
*
* <h3>Numbering of the degrees of freedom (DoFs)</h3>
*
{
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 <tt>points.size()-1</tt>. Note that the first point has to be
- * 0 and the last one 1. If
- * <tt>FE_Q<dim>(QGaussLobatto<1>(fe_degree+1))</tt> 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 <tt>points.size()-1</tt>. Note that the first point has to be
+ * 0 and the last one 1. Constructing
+ * <tt>FE_Q<dim>(QGaussLobatto<1>(fe_degree+1))</tt> is equivalent to the
+ * constructor that specifies the polynomial degree only. For selecting
+ * equidistant nodes at <tt>fe_degree > 2</tt>, construct
+ * <tt>FE_Q<dim>(QIterated<1>(QTrapez<1>(),fe_degree))</tt>.
*/
FE_Q (const Quadrature<1> &points);
#include <deal.II/base/quadrature.h>
#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/template_constraints.h>
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_tools.h>
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<Point<1> >
+ 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<N, ExcInternalError());
- if (N==1)
- return Point<1> (.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<N, ExcInternalError());
-
- if (N==1)
- return Point<2> (.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<N, ExcInternalError());
- if (N==1)
- return Point<3> (.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<Point<1> >(1, Point<1>(0.5));
}
}
-
template <int dim, int spacedim>
FE_DGQ<dim, spacedim>::FE_DGQ (const unsigned int degree)
:
- FE_Poly<TensorProductPolynomials<dim>, dim, spacedim> (
- TensorProductPolynomials<dim>(Polynomials::LagrangeEquidistant::generate_complete_basis(degree)),
- FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
- std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(degree),1, degree).dofs_per_cell, true),
- std::vector<ComponentMask>(FiniteElementData<dim>(
- get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector<bool>(1,true)))
+ FE_Poly<TensorProductPolynomials<dim>, dim, spacedim>
+ (TensorProductPolynomials<dim>(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))),
+ FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
+ std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(degree),1, degree).dofs_per_cell, true),
+ std::vector<ComponentMask>(FiniteElementData<dim>(get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector<bool>(1,true)))
{
- // fill in support points
- if (degree == 0)
- {
- // constant elements, take
- // midpoint
- this->unit_support_points.resize(1);
- for (unsigned int i=0; i<dim; ++i)
- this->unit_support_points[0](i) = 0.5;
- }
- else
- {
- // number of points: (degree+1)^dim
- unsigned int n = degree+1;
- for (unsigned int i=1; i<dim; ++i)
- n *= degree+1;
-
- this->unit_support_points.resize(n);
-
- const double step = 1./degree;
- Point<dim> 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<dim> support_quadrature(get_QGaussLobatto_points(degree));
+ Assert (support_quadrature.get_points().size() > 0,
+ (typename FiniteElement<dim, spacedim>::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
TensorProductPolynomials<dim>(Polynomials::generate_complete_Lagrange_basis(points.get_points())),
FiniteElementData<dim>(get_dpo_vector(points.size()-1), 1, points.size()-1, FiniteElementData<dim>::L2),
std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(points.size()-1),1, points.size()-1).dofs_per_cell, true),
- std::vector<ComponentMask>(FiniteElementData<dim>(
- get_dpo_vector(points.size()-1),1, points.size()-1).dofs_per_cell, std::vector<bool>(1,true)))
+ std::vector<ComponentMask>(FiniteElementData<dim>(get_dpo_vector(points.size()-1),1, points.size()-1).dofs_per_cell, std::vector<bool>(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<dim, spacedim>::ExcFEHasNoSupportPoints ()));
Quadrature<dim> support_quadrature(points);
std::string
FE_DGQ<dim, spacedim>::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<"
// generate a point on this
// cell and evaluate the
// shape functions there
- const Point<dim> p = generate_unit_point (j, this->dofs_per_cell,
- dealii::internal::int2type<dim>());
+ const Point<dim> p = this->unit_support_points[j];
for (unsigned int i=0; i<this->dofs_per_cell; ++i)
cell_interpolation(j,i)
= this->poly_space.compute_value (i, p);
// 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();
}
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();
}
DEAL_II_NAMESPACE_OPEN
+namespace
+{
+ std::vector<Point<1> >
+ 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<dim>(0)."));
+ return std::vector<Point<1> >();
+ }
+}
+
+
template <int dim, int spacedim>
FE_Q<dim,spacedim>::FE_Q (const unsigned int degree)
:
- FE_Q_Base<TensorProductPolynomials<dim>, dim, spacedim> (
- TensorProductPolynomials<dim>(Polynomials::LagrangeEquidistant::generate_complete_basis(degree)),
- FiniteElementData<dim>(this->get_dpo_vector(degree),
- 1, degree,
- FiniteElementData<dim>::H1),
- std::vector<bool> (1, false))
+ FE_Q_Base<TensorProductPolynomials<dim>, dim, spacedim>
+ (TensorProductPolynomials<dim>(Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))),
+ FiniteElementData<dim>(this->get_dpo_vector(degree),
+ 1, degree,
+ FiniteElementData<dim>::H1),
+ std::vector<bool> (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<dim>(0)."));
- std::vector<Point<1> > support_points_1d(degree+1);
- for (unsigned int i=0; i<=degree; ++i)
- support_points_1d[i][0] = static_cast<double>(i)/degree;
-
- this->initialize(support_points_1d);
+ this->initialize(get_QGaussLobatto_points(degree));
}
}
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.
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)
{
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; j<q_dofs_per_cell; ++j)
if (j!=i)
Assert (std::fabs (this->poly_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
double sum = 0;
for (unsigned int col=0; col<this->dofs_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
// distinguish q/q_dg0 case
const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(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<unsigned int> &index_map_inverse =
}
FE_Q_Helper::increment_indices<dim> (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
FE_Q_Bubbles<dim,spacedim>::FE_Q_Bubbles (const unsigned int q_degree)
:
FE_Q_Base<TensorProductPolynomialsBubbles<dim>, dim, spacedim> (
- TensorProductPolynomialsBubbles<dim>(Polynomials::LagrangeEquidistant::generate_complete_basis(q_degree)),
+ TensorProductPolynomialsBubbles<dim>(Polynomials::generate_complete_Lagrange_basis(QGaussLobatto<1>(q_degree+1).get_points())),
FiniteElementData<dim>(get_dpo_vector(q_degree),
1, q_degree+1,
FiniteElementData<dim>::H1),
ExcMessage ("This element can only be used for polynomial degrees "
"greater than zero"));
- std::vector<Point<1> > support_points_1d(q_degree+1);
- for (unsigned int i=0; i<=q_degree; ++i)
- support_points_1d[i][0] = static_cast<double>(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<dim> point;
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"));
}
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
{
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();
}
FE_Q_DG0<dim,spacedim>::FE_Q_DG0 (const unsigned int degree)
:
FE_Q_Base<TensorProductPolynomialsConst<dim>, dim, spacedim> (
- TensorProductPolynomialsConst<dim>(Polynomials::LagrangeEquidistant::generate_complete_basis(degree)),
+ TensorProductPolynomialsConst<dim>(Polynomials::generate_complete_Lagrange_basis(QGaussLobatto<1>(degree+1).get_points())),
FiniteElementData<dim>(get_dpo_vector(degree),
1, degree,
FiniteElementData<dim>::L2),
ExcMessage ("This element can only be used for polynomial degrees "
"greater than zero"));
- std::vector<Point<1> > support_points_1d(degree+1);
- for (unsigned int i=0; i<=degree; ++i)
- support_points_1d[i][0] = static_cast<double>(i)/degree;
-
- this->initialize(support_points_1d);
+ this->initialize(QGaussLobatto<1>(degree+1).get_points());
// adjust unit support point for discontinuous node
Point<dim> point;
}
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
{
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)
const FEFactoryBase<dim,spacedim> *fef=dynamic_cast<const FEFactoryBase<dim,spacedim>*>(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<int,unsigned int> 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<dim,spacedim> *fef=dynamic_cast<const FEFactoryBase<dim,spacedim>*>(ptr);
+ return fef->get(QIterated<1>(QTrapez<1>(),tmp.first));
+ }
else
{
AssertThrow (false,ExcNotImplemented());
projects / dealii.git / commitdiff