From: Martin Kronbichler 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=24956a148fdf7de8510f4f0c217777d1702e9232;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 degree>=3, * but it is also implemented for degrees between 0 and two. For the linear @@ -639,42 +642,42 @@ namespace Polynomials * * n=3 * 1057 - * 21.40 + * 17.18 * * * n=4 * 6580 - * 15.52 + * 16.83 * * * n=5 * 1.875e+04 - * 18.52 + * 19.37 * * * n=6 * 6.033e+04 - * 19.42 + * 18.99 * * * n=10 * 9.756e+05 - * 27.85 + * 25.65 * * * n=15 * 9.431e+06 - * 40.48 + * 36.47 * * * n=25 * 2.220e+08 - * 68.30 + * 62.28 * * * n=35 * 2.109e+09 - * 98.06 + * 91.50 * * * 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, std::vector > - get_constant_modes () const; - /** * Return a string that uniquely identifies a finite element. This class * returns FE_DGQHermite(degree), with dim 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> points(degree); + double ratio = 1.; + for (unsigned int i=0; ilagrange_support_points[i]; + if (i>0) + ratio *= -this->lagrange_support_points[i]; + } + Polynomial helper(points, 0); + std::vector 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 && indexlagrange_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> points(degree); + double ratio = 1.; + for (unsigned int i=0; ilagrange_support_points[i]; + if (ilagrange_support_points[i]; + } + Polynomial helper(points, degree-1); + std::vector 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::FE_DGQHermite (const unsigned int degree) -template -std::pair, std::vector > -FE_DGQHermite::get_constant_modes () const -{ - if (this->degree < 3) - return this->FE_DGQ::get_constant_modes(); - else - { - AssertThrow(false, - ExcMessage("Constant mode cannot be represented by 0/1 vector")); - return std::pair, std::vector > - (Table<2,bool>(1, this->dofs_per_cell), - std::vector(1, 0)); - } -} - - - template std::string FE_DGQHermite::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(1)); print_constant_modes(FE_DGQLegendre(2)); print_constant_modes(FE_DGQHermite(2)); + print_constant_modes(FE_DGQHermite(3)); print_constant_modes(FE_DGP(2)); print_constant_modes(FE_Q_Hierarchical(1)); print_constant_modes(FE_Q_Hierarchical(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::