From e2d43186e4253782901efbaacee8b09d335bbffe Mon Sep 17 00:00:00 2001 From: Jake Harmon Date: Wed, 21 Apr 2021 22:49:38 -0600 Subject: [PATCH] Corrects 2-D FE_Nedelec Interpolation Fixes a conditional statement in FE_Nedelec::convert_...to_dof_values. Adds a new test function in tests/fe/interpolation_common.h for use in tests/fe/interpolation_nedelec.cc as the previous test case indicated (incorrectly) that the function returned correct DoF values. --- doc/news/changes/minor/20210422Harmon | 4 ++ source/fe/fe_nedelec.cc | 20 ++++---- tests/fe/interpolate_common.h | 66 +++++++++++++++++++++++++++ tests/fe/interpolate_nedelec.cc | 18 +++++++- tests/fe/interpolate_nedelec.output | 15 +++--- 5 files changed, 106 insertions(+), 17 deletions(-) create mode 100644 doc/news/changes/minor/20210422Harmon diff --git a/doc/news/changes/minor/20210422Harmon b/doc/news/changes/minor/20210422Harmon new file mode 100644 index 0000000000..ef9308a049 --- /dev/null +++ b/doc/news/changes/minor/20210422Harmon @@ -0,0 +1,4 @@ +Bugfix: FE_Nedelec<2>::convert_generalized_support_point_values_to_dof_values() +now works correctly for every degree. +
+(Jake Harmon, 2021/04/22) diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc index 1b301fe080..95b01f6511 100644 --- a/source/fe/fe_nedelec.cc +++ b/source/fe/fe_nedelec.cc @@ -3171,7 +3171,7 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( { // Let us begin with the // interpolation part. - const QGauss reference_edge_quadrature(this->degree); + const QGauss<1> reference_edge_quadrature(this->degree); const unsigned int n_edge_points = reference_edge_quadrature.size(); for (unsigned int i = 0; i < 2; ++i) @@ -3191,17 +3191,17 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( nodal_values[(i + 2 * j) * this->degree] = 0.0; } - // If the degree is greater - // than 0, then we have still - // some higher order edge - // shape functions to - // consider. - // Here the projection part - // starts. The dof support_point_values - // are obtained by solving + // If the Nedelec element degree is greater + // than 0 (i.e., the polynomial degree is greater than 1), + // then we have still some higher order edge + // shape functions to consider. + // Note that this->degree returns the polynomial + // degree. + // Here the projection part starts. + // The dof support_point_values are obtained by solving // a linear system of // equations. - if (this->degree - 1 > 1) + if (this->degree > 1) { // We start with projection // on the higher order edge diff --git a/tests/fe/interpolate_common.h b/tests/fe/interpolate_common.h index 97a4541a81..6f4fbfe171 100644 --- a/tests/fe/interpolate_common.h +++ b/tests/fe/interpolate_common.h @@ -15,6 +15,7 @@ #include +#include #include #include @@ -146,3 +147,68 @@ public: } } }; + + +// Local implementation of an alternative test function + + +template +class PolynomialFunction : public Function +{ +public: + PolynomialFunction() + : Function(COMP) + { + Table<2, double> exponents(degree, dim); + for (unsigned int i = 0; i < degree; ++i) + for (unsigned int d = 0; d < dim; ++d) + exponents[i][d] = i + d; + std::vector coeffs(degree); + for (unsigned int i = 0; i < degree; ++i) + coeffs[i] = std::pow(-1.0, static_cast(i)) * (i + 1); + poly = std::make_unique>( + Functions::Polynomial(exponents, coeffs)); + } + + double + value(const Point &p, const unsigned int c) const + { + return poly->value(p, 0) + c; + } + + void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int c) const + { + Assert(values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int i = 0; i < points.size(); ++i) + { + const Point &p = points[i]; + + values[i] = poly->value(p, 0) + c; + } + } + + void + vector_value_list(const std::vector> &points, + std::vector> & values) const + { + Assert(values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + Assert(values[0].size() == this->n_components, + ExcDimensionMismatch(values.size(), this->n_components)); + + for (unsigned int i = 0; i < points.size(); ++i) + { + const Point &p = points[i]; + for (unsigned int c = 0; c < COMP; ++c) + values[i](c) = poly->value(p, 0); + } + } + +private: + std::unique_ptr> poly; +}; diff --git a/tests/fe/interpolate_nedelec.cc b/tests/fe/interpolate_nedelec.cc index 5ce2b14c97..085885cc05 100644 --- a/tests/fe/interpolate_nedelec.cc +++ b/tests/fe/interpolate_nedelec.cc @@ -20,7 +20,6 @@ #include "interpolate_common.h" - // FE_Nedelec::interpolate(...) template @@ -37,6 +36,7 @@ check1(const Function &f, const unsigned int degree) Vector(dim)); f.vector_value_list(fe.get_generalized_support_points(), values); fe.convert_generalized_support_point_values_to_dof_values(values, dofs); + deallog << " vector " << vector_difference(fe, dofs, f, 0) << std::endl; } @@ -56,16 +56,32 @@ main() check1(w22, 2); Q1WedgeFunction<2, 3, 2> w23; check1(w23, 3); + + PolynomialFunction<2, 1, 2> p21; + check1(p21, 1); + PolynomialFunction<2, 2, 2> p22; + check1(p22, 2); + PolynomialFunction<2, 3, 2> p23; + check1(p23, 3); } { Q1WedgeFunction<3, 1, 3> w21; check1(w21, 1); check1(w21, 2); + + + // FIXME - higher order interpolation is currently broken + // PolynomialFunction<3, 1, 3> p21; + // check1(p21, 1); // Q1WedgeFunction<3, 2, 3> w22; // check1(w22, 2); // Q1WedgeFunction<3, 3, 3> w23; // check1(w23, 3); + // PolynomialFunction<3, 2, 3> p22; + // check1(p22, 2); + // PolynomialFunction<3, 3, 3> p23; + // check1(p23, 3); } } diff --git a/tests/fe/interpolate_nedelec.output b/tests/fe/interpolate_nedelec.output index aeb8cd57d6..2c684ca007 100644 --- a/tests/fe/interpolate_nedelec.output +++ b/tests/fe/interpolate_nedelec.output @@ -1,9 +1,12 @@ -DEAL::FE_Nedelec<2>(0) 4 vector 0.00000 -DEAL::FE_Nedelec<2>(1) 12 vector 0.00000 +DEAL::FE_Nedelec<2>(0) 4 vector 2.77556e-17 +DEAL::FE_Nedelec<2>(1) 12 vector 1.38778e-17 DEAL::FE_Nedelec<2>(2) 21 vector 5.55112e-17 DEAL::FE_Nedelec<2>(2) 21 vector 1.11022e-16 -DEAL::FE_Nedelec<2>(2) 21 vector 1.58727e-16 -DEAL::FE_Nedelec<2>(3) 32 vector 1.17788e-15 -DEAL::FE_Nedelec<3>(1) 56 vector 5.55112e-17 -DEAL::FE_Nedelec<3>(2) 117 vector 1.11022e-16 +DEAL::FE_Nedelec<2>(2) 21 vector 2.81025e-16 +DEAL::FE_Nedelec<2>(3) 32 vector 7.56830e-16 +DEAL::FE_Nedelec<2>(1) 12 vector 5.55112e-17 +DEAL::FE_Nedelec<2>(2) 21 vector 4.44089e-16 +DEAL::FE_Nedelec<2>(3) 32 vector 1.60288e-15 +DEAL::FE_Nedelec<3>(1) 56 vector 2.77556e-17 +DEAL::FE_Nedelec<3>(2) 117 vector 1.66533e-16 -- 2.39.5