]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make Gauss-Lobatto points the default node distribution in FE_Q/FE_DGQ.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 5 Apr 2016 12:49:54 +0000 (14:49 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 9 Apr 2016 13:49:09 +0000 (15:49 +0200)
include/deal.II/fe/fe_dgq.h
include/deal.II/fe/fe_q.h
source/fe/fe_dgq.cc
source/fe/fe_q.cc
source/fe/fe_q_base.cc
source/fe/fe_q_bubbles.cc
source/fe/fe_q_dg0.cc
source/fe/fe_tools.cc

index 2c44352caec415f1e7261dc3f231ecd77fdf3ad7..e8beccbfc24644cf2a6a8869fce0c07b371dfad8 100644 (file)
@@ -75,6 +75,30 @@ template <int dim> class Quadrature;
  * 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>
@@ -84,7 +108,7 @@ public:
   /**
    * 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);
@@ -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.
  *
index f47cba2e7cb39c0c4d1f1f1a245fccc3792b585f..c821dc8f720d9f119fd558468fe1cc4c447f882a 100644 (file)
@@ -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.
  *
+ * <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>
  *
@@ -523,19 +553,21 @@ class FE_Q : public FE_Q_Base<TensorProductPolynomials<dim>,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 <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);
 
index e8206eefd77bdfa7966428532dc0c6485e885e7f..12c1973c4c5aea7bd49e9071b37edd7ffd5892de 100644 (file)
@@ -16,7 +16,6 @@
 
 #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
@@ -194,13 +73,10 @@ FE_DGQ<dim, spacedim>::FE_DGQ (const Quadrature<1> &points)
     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);
@@ -218,12 +94,9 @@ template <int dim, int spacedim>
 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<"
@@ -378,8 +251,7 @@ get_interpolation_matrix (const FiniteElement<dim, spacedim> &x_source_fe,
       // 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);
@@ -807,15 +679,20 @@ FE_DGQArbitraryNodes<dim,spacedim>::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<dim,spacedim>::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();
     }
 
index d170222b30fd6b04517ab459414e941f033d03fd..a4c969dee57945e4198668a1c4670df3d3df1911 100644 (file)
 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));
 }
 
 
@@ -95,9 +104,16 @@ FE_Q<dim,spacedim>::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<dim,spacedim>::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)
index 9137d03e3703141bb6d413f0646cf19440ee71aa..a521c6c5b8a748b0820edc4dbb8b2a8ff557d79c 100644 (file)
@@ -1180,12 +1180,22 @@ FE_Q_Base<PolynomialType,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
 
@@ -1301,7 +1311,14 @@ FE_Q_Base<PolynomialType,dim,spacedim>
           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
 
@@ -1343,27 +1360,23 @@ FE_Q_Base<PolynomialType,dim,spacedim>
       // 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 =
@@ -1427,8 +1440,14 @@ FE_Q_Base<PolynomialType,dim,spacedim>
                     }
                   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
index 44b7e5fc6ce05b945dcbc910fbefef7527f62d30..9868ae39abb8526e0f9829e52909810b6dc33555 100644 (file)
@@ -172,7 +172,7 @@ template <int dim, int spacedim>
 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),
@@ -183,11 +183,7 @@ FE_Q_Bubbles<dim,spacedim>::FE_Q_Bubbles (const unsigned int q_degree)
           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;
@@ -221,9 +217,7 @@ FE_Q_Bubbles<dim,spacedim>::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<dim,spacedim>::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<dim,spacedim>::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();
 }
index d8cb1c808a0aad222f8430f2a3474252e83fb43b..d2218eadc9bf38b4e284bee6a21ad0ce277e38ca 100644 (file)
@@ -33,7 +33,7 @@ template <int dim, int spacedim>
 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),
@@ -43,11 +43,7 @@ FE_Q_DG0<dim,spacedim>::FE_Q_DG0 (const unsigned int degree)
           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;
@@ -133,9 +129,16 @@ FE_Q_DG0<dim,spacedim>::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<dim,spacedim>::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)
index eee128e3679c05f9a2dc71459328a675455953a4..dc82e00ba6f37dbd3b2c7bf1d8f617f1e3daf1a1 100644 (file)
@@ -1540,6 +1540,23 @@ namespace FETools
                     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());

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.