Extend Hermite-like interpolation to case degree=2.

diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h
index c78532899f6b6e4e973e5e3365d62f4eb9749e36..3fdee6c474d312fcb961342cb736c240ec951567 100644 (file)
--- a/include/deal.II/base/polynomial.h
+++ b/include/deal.II/base/polynomial.h
@@ -591,9 +591,7 @@ namespace Polynomials
 
   /**
    * Polynomials for a variant of Hermite polynomials with better condition
-   * number in the interpolation than the basis from
-   * HermiteInterpolation. This class is only implemented for degree at least
-   * three, $n\geq 3$.
+   * number in the interpolation than the basis from HermiteInterpolation.
    *
    * In analogy to the actual Hermite polynomials this basis evaluates the
    * first polynomial $p_0$ to 1 at $x=0$ and has both a zero value and zero
@@ -616,7 +614,14 @@ namespace Polynomials
    * Lagrange polynomials with double roots at $x=0$ and $x=1$. For example at
    * $n=4$, all of $p_0, p_1, p_3, p_4$ get an additional root at $x=0.5$
    * through the factor $(x-0.5)$.
-
+   *
+   * The basis only contains Hermite information at <code>degree>=3</code>,
+   * but it is also implemented for degrees between 0 and two. For the linear
+   * case, the usual hat functions are implemented, whereas the polynomials
+   * for <code>degree=2</code> are $p_0(x)=(1-x)^2$, $p_1(x)=4x(x-1)$, and
+   * $p_2(x)=x^2$, in accordance with the construction principle for degree 3
+   * that allows a non-zero of $p_0$ and $p_2$.
+   *
    * These two relaxations improve the condition number of the mass matrix
    * (i.e., interpolation) significantly, as can be seen from the following
    * table:
@@ -696,8 +701,7 @@ namespace Polynomials
 
     /**
      * Return the polynomials with index <tt>0</tt> up to <tt>degree+1</tt> in
-     * a space of degree up to <tt>degree</tt>. Here, <tt>degree</tt> has to
-     * be at least 3.
+     * a space of degree up to <tt>degree</tt>.
      */
     static std::vector<Polynomial<double> >
     generate_complete_basis (const unsigned int degree);

diff --git a/source/base/polynomial.cc b/source/base/polynomial.cc
index c495c34d9102a3617f7e08ed18ac97a624549526..3fb971418b27e404bc8767a96499e756430afedd 100644 (file)
--- a/source/base/polynomial.cc
+++ b/source/base/polynomial.cc
@@ -1286,9 +1286,6 @@ namespace Polynomials
     :
     Polynomial<double>(0)
   {
-    Assert(degree>=3,
-           ExcNotImplemented("Hermite interpolation makes no sense for "
-                             "degrees less than three"));
     AssertIndexRange(index, degree+1);
 
     this->coefficients.clear();
@@ -1296,16 +1293,52 @@ namespace Polynomials
 
     this->lagrange_support_points.resize(degree);
 
-    // 4 Polynomials with degree 3
-    // entries (1,0) and (3,2) of the mass matrix will be equal to 0
-    //
-    //     | x  0  x  x |
-    //     | 0  x  x  x |
-    // M = | x  x  x  0 |
-    //     | x  x  0  x |
-    //
-    if (degree==3)
+    if (degree == 0)
+      this->lagrange_weight = 1.;
+    else if (degree == 1)
       {
+        if (index == 0)
+          {
+            this->lagrange_support_points[0] = 1.;
+            this->lagrange_weight = -1.;
+          }
+        else
+          {
+            this->lagrange_support_points[0] = 0.;
+            this->lagrange_weight = 1.;
+          }
+      }
+    else if (degree == 2)
+      {
+        if (index == 0)
+          {
+            this->lagrange_support_points[0] = 1.;
+            this->lagrange_support_points[1] = 1.;
+            this->lagrange_weight = 1.;
+          }
+        else if (index == 1)
+          {
+            this->lagrange_support_points[0] = 0;
+            this->lagrange_support_points[1] = 1;
+            this->lagrange_weight = 4.;
+          }
+        else
+          {
+            this->lagrange_support_points[0] = 0.;
+            this->lagrange_support_points[1] = 0.;
+            this->lagrange_weight = 1.;
+          }
+      }
+    else if (degree==3)
+      {
+        // 4 Polynomials with degree 3
+        // entries (1,0) and (3,2) of the mass matrix will be equal to 0
+        //
+        //     | x  0  x  x |
+        //     | 0  x  x  x |
+        // M = | x  x  x  0 |
+        //     | x  x  0  x |
+        //
         if (index==0)
           {
             this->lagrange_support_points[0] = 2./7.;
@@ -1335,23 +1368,22 @@ namespace Polynomials
             this->lagrange_weight = 3.5;
           }
       }
-
-    // Higher order Polynomials degree>=4: the entries (1,0) and
-    // (degree,degree-1) of the mass matrix will be equal to 0
-    //
-    //     | x  0  x  x         x  x  x |
-    //     | 0  x  x  x  . . .  x  x  x |
-    //     | x  x  x  x         x  x  x |
-    //     | x  x  x  x         x  x  x |
-    //     |     .       .         .    |
-    // M = |     .         .       .    |
-    //     |     .           .     .    |
-    //     | x  x  x  x         x  x  x |
-    //     | x  x  x  x  . . .  x  x  0 |
-    //     | x  x  x  x         x  0  x |
-    //
-    if (degree >= 4)
+    else
       {
+        // Higher order Polynomials degree>=4: the entries (1,0) and
+        // (degree,degree-1) of the mass matrix will be equal to 0
+        //
+        //     | x  0  x  x         x  x  x |
+        //     | 0  x  x  x  . . .  x  x  x |
+        //     | x  x  x  x         x  x  x |
+        //     | x  x  x  x         x  x  x |
+        //     |     .       .         .    |
+        // M = |     .         .       .    |
+        //     |     .           .     .    |
+        //     | x  x  x  x         x  x  x |
+        //     | x  x  x  x  . . .  x  x  0 |
+        //     | x  x  x  x         x  0  x |
+
         // We find the inner points as the zeros of the Jacobi polynomials
         // with alpha = beta = 2 which is the polynomial with the kernel
         // (1-x)^2 (1+x)^2, the two polynomials achieving zero value and zero

diff --git a/source/fe/fe_dgq.cc b/source/fe/fe_dgq.cc
index 9abded35b875acf05b024c1ab8a2601507a03409..321944f493ff822c97262ceb00a888a7b0ed661d 100644 (file)
--- a/source/fe/fe_dgq.cc
+++ b/source/fe/fe_dgq.cc
@@ -62,10 +62,8 @@ FE_DGQ<dim, spacedim>::FE_DGQ (const unsigned int degree)
 {
   // Compute support points, which are the tensor product of the Lagrange
   // interpolation points in the constructor.
-  Quadrature<dim> support_quadrature(internal::FE_DGQ::get_QGaussLobatto_points(degree));
-  Assert (support_quadrature.get_points().size() > 0,
-          (typename FiniteElement<dim, spacedim>::ExcFEHasNoSupportPoints ()));
-  this->unit_support_points = support_quadrature.get_points();
+  this->unit_support_points =
+    Quadrature<dim>(internal::FE_DGQ::get_QGaussLobatto_points(degree)).get_points();
 
   // do not initialize embedding and restriction here. these matrices are
   // initialized on demand in get_restriction_matrix and
@@ -877,11 +875,7 @@ FE_DGQLegendre<dim,spacedim>::clone() const
 
 template <int dim, int spacedim>
 FE_DGQHermite<dim,spacedim>::FE_DGQHermite (const unsigned int degree)
-  : FE_DGQ<dim,spacedim>(degree < 3 ?
-                         Polynomials::generate_complete_Lagrange_basis
-                         (internal::FE_DGQ::get_QGaussLobatto_points(degree))
-                         :
-                         Polynomials::HermiteLikeInterpolation::generate_complete_basis(degree))
+  : FE_DGQ<dim,spacedim>(Polynomials::HermiteLikeInterpolation::generate_complete_basis(degree))
 {}