Change weight in Hermite-like polynomials to make constant represented by (1,1,1...

diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h
index 3fdee6c474d312fcb961342cb736c240ec951567..129eb330502262f9187965f1407f308625752bb8 100644 (file)
--- a/include/deal.II/base/polynomial.h
+++ b/include/deal.II/base/polynomial.h
@@ -613,7 +613,10 @@ namespace Polynomials
    * property. Then, these polynomials are constructed in the usual way as
    * 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)$.
+   * through the factor $(x-0.5)$. In summary, this basis is dominated by
+   * nodal contributions, but it is not a nodal one because the second and
+   * second to last polynomials that are non-nodal, and due to the presence of
+   * double nodes in $x=0$ and $x=1$.
    *
    * 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
@@ -639,42 +642,42 @@ namespace Polynomials
    *   <tr>
    *    <th>n=3</th>
    *    <th>1057</th>
-   *    <th>21.40</th>
+   *    <th>17.18</th>
    *   </tr>
    *   <tr>
    *    <th>n=4</th>
    *    <th>6580</th>
-   *    <th>15.52</th>
+   *    <th>16.83</th>
    *   </tr>
    *   <tr>
    *    <th>n=5</th>
    *    <th>1.875e+04</th>
-   *    <th>18.52</th>
+   *    <th>19.37</th>
    *   </tr>
    *   <tr>
    *    <th>n=6</th>
    *    <th>6.033e+04</th>
-   *    <th>19.42</th>
+   *    <th>18.99</th>
    *   </tr>
    *   <tr>
    *    <th>n=10</th>
    *    <th>9.756e+05</th>
-   *    <th>27.85</th>
+   *    <th>25.65</th>
    *   </tr>
    *   <tr>
    *    <th>n=15</th>
    *    <th>9.431e+06</th>
-   *    <th>40.48</th>
+   *    <th>36.47</th>
    *   </tr>
    *   <tr>
    *    <th>n=25</th>
    *    <th>2.220e+08</th>
-   *    <th>68.30</th>
+   *    <th>62.28</th>
    *   </tr>
    *   <tr>
    *    <th>n=35</th>
    *    <th>2.109e+09</th>
-   *    <th>98.06</th>
+   *    <th>91.50</th>
    *   </tr>
    * </table>
    *

diff --git a/include/deal.II/fe/fe_dgq.h b/include/deal.II/fe/fe_dgq.h
index cf7ff20b679258503f566f27fc72ab997e34033d..3fb2441995a51f771b7b26556ac44a12c03e3306 100644 (file)
--- a/include/deal.II/fe/fe_dgq.h
+++ b/include/deal.II/fe/fe_dgq.h
@@ -506,16 +506,6 @@ public:
    */
   FE_DGQHermite (const unsigned int degree);
 
-  /**
-   * Return a list of constant modes of the element. Due to the special
-   * construction of the Hermite basis, there is no simple representation of
-   * the constant value 1 in terms of a 0/1 selection, so get_constant_modes()
-   * throws an exception in case it is called for degrees larger or equal to
-   * three.
-   */
-  virtual std::pair<Table<2,bool>, std::vector<unsigned int> >
-  get_constant_modes () const;
-
   /**
    * Return a string that uniquely identifies a finite element. This class
    * returns <tt>FE_DGQHermite<dim>(degree)</tt>, with <tt>dim</tt> and

diff --git a/source/base/polynomial.cc b/source/base/polynomial.cc
index 3fb971418b27e404bc8767a96499e756430afedd..693d6b1d32b4b9d85d23508cbcdfbfff0bff7dfc 100644 (file)
--- a/source/base/polynomial.cc
+++ b/source/base/polynomial.cc
@@ -1320,7 +1320,7 @@ namespace Polynomials
           {
             this->lagrange_support_points[0] = 0;
             this->lagrange_support_points[1] = 1;
-            this->lagrange_weight = 4.;
+            this->lagrange_weight = -2.;
           }
         else
           {
@@ -1351,14 +1351,17 @@ namespace Polynomials
             this->lagrange_support_points[0] = 0.;
             this->lagrange_support_points[1] = 1.;
             this->lagrange_support_points[2] = 1.;
-            this->lagrange_weight = 6.75;
+
+            // this magic value 5.5 is obtained when evaluating the general
+            // formula below for the degree=3 case
+            this->lagrange_weight = 5.5;
           }
         else if (index==2)
           {
             this->lagrange_support_points[0] = 0.;
             this->lagrange_support_points[1] = 0.;
             this->lagrange_support_points[2] = 1.;
-            this->lagrange_weight = -6.75;
+            this->lagrange_weight = -5.5;
           }
         else if (index==3)
           {
@@ -1383,7 +1386,7 @@ namespace Polynomials
         //     | 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
@@ -1445,8 +1448,30 @@ namespace Polynomials
             this->lagrange_support_points[degree-2] = 1.;
             this->lagrange_support_points[degree-1] = 1.;
 
-            // scale close to approximate maximum
-            this->lagrange_weight = 1./this->value(0.4*jacobi_roots(0));
+            // Select the weight to make the derivative of the sum of P_0 and
+            // P_1 in zero to be 0. The derivative in x=0 is simply given by
+            // p~(0)/auxiliary_zero+p~'(0) + a*p~(0), where p~(x) is the
+            // Lagrange polynomial in all points except the first one which is
+            // the same for P_0 and P_1, and a is the weight we seek here. If
+            // we solve this for a, we obtain the desired property. Since the
+            // basis is nodal for all interior points, this property ensures
+            // that the sum of all polynomials with weight 1 is one.
+            std::vector<Point<1>> points(degree);
+            double ratio = 1.;
+            for (unsigned int i=0; i<degree; ++i)
+              {
+                points[i][0] = this->lagrange_support_points[i];
+                if (i>0)
+                  ratio *= -this->lagrange_support_points[i];
+              }
+            Polynomial<double> helper(points, 0);
+            std::vector<double> value_and_grad(2);
+            helper.value(0., value_and_grad);
+            Assert(std::abs(value_and_grad[0]) > 1e-10,
+                   ExcInternalError("There should not be a zero at x=0."));
+
+            const double auxiliary_zero = find_support_point_x_star(jacobi_roots);
+            this->lagrange_weight = (1./auxiliary_zero - value_and_grad[1]/value_and_grad[0])/ratio;
           }
         else if (index>=2 && index<degree-1)
           {
@@ -1470,7 +1495,22 @@ namespace Polynomials
               this->lagrange_support_points[m+2] = jacobi_roots(m);
             this->lagrange_support_points[degree-1] = 1.;
 
-            this->lagrange_weight = 1./this->value(1.0-0.4*jacobi_roots(0));
+            std::vector<Point<1>> points(degree);
+            double ratio = 1.;
+            for (unsigned int i=0; i<degree; ++i)
+              {
+                points[i][0] = this->lagrange_support_points[i];
+                if (i<degree-1)
+                  ratio *= 1.-this->lagrange_support_points[i];
+              }
+            Polynomial<double> helper(points, degree-1);
+            std::vector<double> value_and_grad(2);
+            helper.value(1., value_and_grad);
+            Assert(std::abs(value_and_grad[0]) > 1e-10,
+                   ExcInternalError("There should not be a zero at x=1."));
+
+            const double auxiliary_zero = find_support_point_x_star(jacobi_roots);
+            this->lagrange_weight = (-1./auxiliary_zero - value_and_grad[1]/value_and_grad[0])/ratio;
           }
         else if (index==degree)
           {

diff --git a/source/fe/fe_dgq.cc b/source/fe/fe_dgq.cc
index 321944f493ff822c97262ceb00a888a7b0ed661d..d5d92a3f5c9614117efd907e99eb0b4f5a64b03d 100644 (file)
--- a/source/fe/fe_dgq.cc
+++ b/source/fe/fe_dgq.cc
@@ -880,24 +880,6 @@ FE_DGQHermite<dim,spacedim>::FE_DGQHermite (const unsigned int degree)
 
 
 
-template <int dim, int spacedim>
-std::pair<Table<2,bool>, std::vector<unsigned int> >
-FE_DGQHermite<dim,spacedim>::get_constant_modes () const
-{
-  if (this->degree < 3)
-    return this->FE_DGQ<dim,spacedim>::get_constant_modes();
-  else
-    {
-      AssertThrow(false,
-                  ExcMessage("Constant mode cannot be represented by 0/1 vector"));
-      return std::pair<Table<2,bool>, std::vector<unsigned int> >
-             (Table<2,bool>(1, this->dofs_per_cell),
-              std::vector<unsigned int>(1, 0));
-    }
-}
-
-
-
 template <int dim, int spacedim>
 std::string
 FE_DGQHermite<dim,spacedim>::get_name () const

diff --git a/tests/base/polynomial_hermite_like.cc b/tests/base/polynomial_hermite_like.cc
index a547bc58dd3ec207de82b4f932ebc88c804c6191..a21770105d900fb5509cf987579c6c93becf69ed 100644 (file)
--- a/tests/base/polynomial_hermite_like.cc
+++ b/tests/base/polynomial_hermite_like.cc
@@ -51,6 +51,20 @@ int main()
 {
   initlog();
 
+  deallog << "degree 1" << std::endl;
+  plot(HermiteLikeInterpolation::generate_complete_basis(1));
+  deallog << std::endl << "degree 2" << std::endl;
+  plot(HermiteLikeInterpolation::generate_complete_basis(2));
+  deallog << std::endl << "degree 3" << std::endl;
+  interpolation_conditions(HermiteLikeInterpolation::generate_complete_basis(3));
+  plot(HermiteLikeInterpolation::generate_complete_basis(3));
+  deallog << std::endl << "degree 4" << std::endl;
+  interpolation_conditions(HermiteLikeInterpolation::generate_complete_basis(4));
+  plot(HermiteLikeInterpolation::generate_complete_basis(4));
+  deallog << std::endl << "degree 6" << std::endl;
   interpolation_conditions(HermiteLikeInterpolation::generate_complete_basis(6));
   plot(HermiteLikeInterpolation::generate_complete_basis(6));
+  deallog << std::endl << "degree 9" << std::endl;
+  interpolation_conditions(HermiteLikeInterpolation::generate_complete_basis(9));
+  plot(HermiteLikeInterpolation::generate_complete_basis(9));
 }

diff --git a/tests/base/polynomial_hermite_like.with_lapack=true.output b/tests/base/polynomial_hermite_like.with_lapack=true.output
index 8a0594ffd0db7cacd52a6ed925270b2e6b779287..6e48f5b9a594747d75892ea1e86c573852b0f8ed 100644 (file)
--- a/tests/base/polynomial_hermite_like.with_lapack=true.output
+++ b/tests/base/polynomial_hermite_like.with_lapack=true.output
@@ -1,7 +1,70 @@
 
+DEAL::degree 1
+DEAL::0.00000  1.00000 0.00000
+DEAL::0.125000 0.875000        0.125000
+DEAL::0.250000 0.750000        0.250000
+DEAL::0.375000 0.625000        0.375000
+DEAL::0.500000 0.500000        0.500000
+DEAL::0.625000 0.375000        0.625000
+DEAL::0.750000 0.250000        0.750000
+DEAL::0.875000 0.125000        0.875000
+DEAL::1.00000  0.00000 1.00000
+DEAL::
+DEAL::degree 2
+DEAL::0.00000  1.00000 0.00000 0.00000
+DEAL::0.125000 0.765625        0.218750        0.0156250
+DEAL::0.250000 0.562500        0.375000        0.0625000
+DEAL::0.375000 0.390625        0.468750        0.140625
+DEAL::0.500000 0.250000        0.500000        0.250000
+DEAL::0.625000 0.140625        0.468750        0.390625
+DEAL::0.750000 0.0625000       0.375000        0.562500
+DEAL::0.875000 0.0156250       0.218750        0.765625
+DEAL::1.00000  0.00000 0.00000 1.00000
+DEAL::
+DEAL::degree 3
+DEAL::0         0:     1.00000 -5.50000
+DEAL::0         1:     0.00000 0.00000
+DEAL::1         0:     0.00000 5.50000
+DEAL::1         1:     0.00000 0.00000
+DEAL::2         0:     0.00000 0.00000
+DEAL::2         1:     0.00000 -5.50000
+DEAL::3         0:     0.00000 0.00000
+DEAL::3         1:     1.00000 5.50000
+DEAL::0.00000  1.00000 0.00000 0.00000 0.00000
+DEAL::0.125000 0.430664        0.526367        0.0751953       -0.0322266
+DEAL::0.250000 0.0703125       0.773438        0.257812        -0.101562
+DEAL::0.375000 -0.122070       0.805664        0.483398        -0.166992
+DEAL::0.500000 -0.187500       0.687500        0.687500        -0.187500
+DEAL::0.625000 -0.166992       0.483398        0.805664        -0.122070
+DEAL::0.750000 -0.101562       0.257812        0.773438        0.0703125
+DEAL::0.875000 -0.0322266      0.0751953       0.526367        0.430664
+DEAL::1.00000  0.00000 0.00000 0.00000 1.00000
+DEAL::
+DEAL::degree 4
+DEAL::0         0:     1.00000 -10.0000
+DEAL::0         1:     0.00000 0.00000
+DEAL::1         0:     0.00000 10.0000
+DEAL::1         1:     0.00000 0.00000
+DEAL::2         0:     0.00000 0.00000
+DEAL::2         1:     0.00000 0.00000
+DEAL::3         0:     0.00000 0.00000
+DEAL::3         1:     0.00000 -10.0000
+DEAL::4         0:     0.00000 0.00000
+DEAL::4         1:     1.00000 10.0000
+DEAL::0.00000  1.00000 0.00000 0.00000 0.00000 0.00000
+DEAL::0.125000 0.143555        0.717773        0.191406        -0.102539       0.0498047
+DEAL::0.250000 -0.140625       0.703125        0.562500        -0.234375       0.109375
+DEAL::0.375000 -0.122070       0.366211        0.878906        -0.219727       0.0966797
+DEAL::0.500000 0.00000 0.00000 1.00000 0.00000 0.00000
+DEAL::0.625000 0.0966797       -0.219727       0.878906        0.366211        -0.122070
+DEAL::0.750000 0.109375        -0.234375       0.562500        0.703125        -0.140625
+DEAL::0.875000 0.0498047       -0.102539       0.191406        0.717773        0.143555
+DEAL::1.00000  0.00000 0.00000 0.00000 0.00000 1.00000
+DEAL::
+DEAL::degree 6
 DEAL::0         0:     1.00000 -21.1429
 DEAL::0         1:     0.00000 0.00000
-DEAL::1         0:     0.00000 31.7114
+DEAL::1         0:     0.00000 21.1429
 DEAL::1         1:     0.00000 0.00000
 DEAL::2         0:     0.00000 0.00000
 DEAL::2         1:     0.00000 0.00000
@@ -10,15 +73,46 @@ DEAL::3      1:     0.00000 0.00000
 DEAL::4         0:     0.00000 0.00000
 DEAL::4         1:     0.00000 0.00000
 DEAL::5         0:     0.00000 0.00000
-DEAL::5         1:     0.00000 -31.7114
+DEAL::5         1:     0.00000 -21.1429
 DEAL::6         0:     0.00000 0.00000
 DEAL::6         1:     1.00000 21.1429
 DEAL::0.00000  1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
-DEAL::0.125000 -0.0775452      0.782431        0.643097        -0.131592       0.0836483       -0.111776       0.0352478
-DEAL::0.250000 0.0627790       -0.278714       1.02264 0.140625        -0.0734223      0.0929046       -0.0287388
-DEAL::0.375000 0.126103        -0.471781       0.613544        0.714111        -0.242755       0.283069        -0.0851833
+DEAL::0.125000 -0.0775452      0.521667        0.643097        -0.131592       0.0836483       -0.0745239      0.0352478
+DEAL::0.250000 0.0627790       -0.185826       1.02264 0.140625        -0.0734223      0.0619420       -0.0287388
+DEAL::0.375000 0.126103        -0.314549       0.613544        0.714111        -0.242755       0.188729        -0.0851833
 DEAL::0.500000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000
-DEAL::0.625000 -0.0851833      0.283069        -0.242755       0.714111        0.613544        -0.471781       0.126103
-DEAL::0.750000 -0.0287388      0.0929046       -0.0734223      0.140625        1.02264 -0.278714       0.0627790
-DEAL::0.875000 0.0352478       -0.111776       0.0836483       -0.131592       0.643097        0.782431        -0.0775452
+DEAL::0.625000 -0.0851833      0.188729        -0.242755       0.714111        0.613544        -0.314549       0.126103
+DEAL::0.750000 -0.0287388      0.0619420       -0.0734223      0.140625        1.02264 -0.185826       0.0627790
+DEAL::0.875000 0.0352478       -0.0745239      0.0836483       -0.131592       0.643097        0.521667        -0.0775452
 DEAL::1.00000  0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000
+DEAL::
+DEAL::degree 9
+DEAL::0         0:     1.00000 -43.0000
+DEAL::0         1:     0.00000 0.00000
+DEAL::1         0:     0.00000 43.0000
+DEAL::1         1:     0.00000 0.00000
+DEAL::2         0:     0.00000 0.00000
+DEAL::2         1:     0.00000 0.00000
+DEAL::3         0:     0.00000 0.00000
+DEAL::3         1:     0.00000 0.00000
+DEAL::4         0:     0.00000 0.00000
+DEAL::4         1:     0.00000 0.00000
+DEAL::5         0:     0.00000 0.00000
+DEAL::5         1:     0.00000 0.00000
+DEAL::6         0:     0.00000 0.00000
+DEAL::6         1:     0.00000 0.00000
+DEAL::7         0:     0.00000 0.00000
+DEAL::7         1:     0.00000 0.00000
+DEAL::8         0:     0.00000 0.00000
+DEAL::8         1:     0.00000 -43.0000
+DEAL::9         0:     0.00000 0.00000
+DEAL::9         1:     1.00000 43.0000
+DEAL::0.00000  1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL::0.125000 0.0734501       -0.287123       1.00390 0.270272        -0.0928510      0.0555242       -0.0438505      0.0446656       -0.0410176      0.0170339
+DEAL::0.250000 -0.0360489      0.103340        -0.157805       1.00328 0.120509        -0.0544771      0.0391483       -0.0382418      0.0344467       -0.0141525
+DEAL::0.375000 -0.0402864      0.106060        -0.136357       0.215809        0.941854        -0.131662       0.0793256       -0.0726297      0.0636361       -0.0257504
+DEAL::0.500000 0.0664062       -0.167969       0.200133        -0.244967       0.646397        0.646397        -0.244967       0.200133        -0.167969       0.0664062
+DEAL::0.625000 -0.0257504      0.0636361       -0.0726297      0.0793256       -0.131662       0.941854        0.215809        -0.136357       0.106060        -0.0402864
+DEAL::0.750000 -0.0141525      0.0344467       -0.0382418      0.0391483       -0.0544771      0.120509        1.00328 -0.157805       0.103340        -0.0360489
+DEAL::0.875000 0.0170339       -0.0410176      0.0446656       -0.0438505      0.0555242       -0.0928510      0.270272        1.00390 -0.287123       0.0734501
+DEAL::1.00000  0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000

diff --git a/tests/fe/element_constant_modes.cc b/tests/fe/element_constant_modes.cc
index 4fe22916e4700f2e0ffc6372598fef45b614c2e4..14be190e6fff2b8675e872374b29dec6fa3e57d2 100644 (file)
--- a/tests/fe/element_constant_modes.cc
+++ b/tests/fe/element_constant_modes.cc
@@ -50,6 +50,7 @@ void test()
   print_constant_modes(FE_DGQ<dim>(1));
   print_constant_modes(FE_DGQLegendre<dim>(2));
   print_constant_modes(FE_DGQHermite<dim>(2));
+  print_constant_modes(FE_DGQHermite<dim>(3));
   print_constant_modes(FE_DGP<dim>(2));
   print_constant_modes(FE_Q_Hierarchical<dim>(1));
   print_constant_modes(FE_Q_Hierarchical<dim>(2));

diff --git a/tests/fe/element_constant_modes.output b/tests/fe/element_constant_modes.output
index ee07412c5fc81fc71ee5e40b5f9ab004f344493d..7c0f38548fb96d26e7b9ac33429c88ea73f771e2 100644 (file)
--- a/tests/fe/element_constant_modes.output
+++ b/tests/fe/element_constant_modes.output
@@ -26,6 +26,9 @@ DEAL::
 DEAL::Testing FE_DGQHermite<2>(2)
 DEAL::1 1 1 1 1 1 1 1 1 
 DEAL::
+DEAL::Testing FE_DGQHermite<2>(3)
+DEAL::1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
+DEAL::
 DEAL::Testing FE_DGP<2>(2)
 DEAL::1 0 0 0 0 0 
 DEAL::
@@ -83,6 +86,9 @@ DEAL::
 DEAL::Testing FE_DGQHermite<3>(2)
 DEAL::1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
 DEAL::
+DEAL::Testing FE_DGQHermite<3>(3)
+DEAL::1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
+DEAL::
 DEAL::Testing FE_DGP<3>(2)
 DEAL::1 0 0 0 0 0 0 0 0 0 
 DEAL::