From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Sat, 10 Mar 2018 10:41:10 +0000 (+0100)
Subject: Change weight in Hermite-like polynomials to make constant represented by (1,1,1... 
X-Git-Tag: v9.0.0-rc1~331^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F6028%2Fhead;p=dealii.git

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

diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h
index 3fdee6c474..129eb33050 100644
--- 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 cf7ff20b67..3fb2441995 100644
--- 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 3fb971418b..693d6b1d32 100644
--- 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 321944f493..d5d92a3f5c 100644
--- 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 a547bc58dd..a21770105d 100644
--- 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 8a0594ffd0..6e48f5b9a5 100644
--- 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 4fe22916e4..14be190e6f 100644
--- 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 ee07412c5f..7c0f38548f 100644
--- 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::