From: Jake Harmon <jake.harmon@ieee.org>
Date: Thu, 22 Apr 2021 04:49:38 +0000 (-0600)
Subject: Corrects 2-D FE_Nedelec Interpolation
X-Git-Tag: v9.3.0-rc1~191^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F12076%2Fhead;p=dealii.git

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.
---

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.
+<br>
+(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<dim>::convert_generalized_support_point_values_to_dof_values(
         {
           // Let us begin with the
           // interpolation part.
-          const QGauss<dim - 1> 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<dim>::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 <deal.II/base/function.h>
+#include <deal.II/base/function_lib.h>
 #include <deal.II/base/quadrature_lib.h>
 
 #include <deal.II/fe/fe.h>
@@ -146,3 +147,68 @@ public:
       }
   }
 };
+
+
+// Local implementation of an alternative test function
+
+
+template <int dim, int degree, int COMP = 1>
+class PolynomialFunction : public Function<dim>
+{
+public:
+  PolynomialFunction()
+    : Function<dim>(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<double> coeffs(degree);
+    for (unsigned int i = 0; i < degree; ++i)
+      coeffs[i] = std::pow(-1.0, static_cast<double>(i)) * (i + 1);
+    poly = std::make_unique<Functions::Polynomial<dim>>(
+      Functions::Polynomial<dim>(exponents, coeffs));
+  }
+
+  double
+  value(const Point<dim> &p, const unsigned int c) const
+  {
+    return poly->value(p, 0) + c;
+  }
+
+  void
+  value_list(const std::vector<Point<dim>> &points,
+             std::vector<double> &          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<dim> &p = points[i];
+
+        values[i] = poly->value(p, 0) + c;
+      }
+  }
+
+  void
+  vector_value_list(const std::vector<Point<dim>> &points,
+                    std::vector<Vector<double>> &  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<dim> &p = points[i];
+        for (unsigned int c = 0; c < COMP; ++c)
+          values[i](c) = poly->value(p, 0);
+      }
+  }
+
+private:
+  std::unique_ptr<Functions::Polynomial<dim>> 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<dim>::interpolate(...)
 
 template <int dim>
@@ -37,6 +36,7 @@ check1(const Function<dim> &f, const unsigned int degree)
                                      Vector<double>(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