FE_Q_iso_Q1: varying subdivisions

include/deal.II/base/polynomials_piecewise.h

@@ -80,6 +80,15 @@ namespace Polynomials
                         const unsigned int        interval,
                         const bool                spans_next_interval);
 
+    /**
+     * Constructor for linear Lagrange polynomial on an interval that is a
+     * subset of the unit interval. It uses a polynomial description that is
+     * scaled to the size of the subinterval compared to the unit interval.
+     * The subintervals are bounded by the adjacent points in @p points.
+     */
+    PiecewisePolynomial(const std::vector<Point<1, number>> &points,
+                        const unsigned int                   index);
+
     /**
      * Return the value of this polynomial at the given point, evaluating the
      * underlying polynomial. The polynomial evaluates to zero when outside of
@@ -172,6 +181,24 @@ namespace Polynomials
      * one given in subinterval and the next one.
      */
     bool spans_two_intervals;
+
+    /**
+     * Points bounding the subintervals in the case that piecewise linear
+     * polynomial on varying subintervals is requested.
+     */
+    std::vector<number> points;
+
+    /**
+     * Precomputed inverses of the lengths of the subintervals, i.e.,
+     * `one_over_lengths[i] = 1.0 / (points[i + 1] - points[i]`.
+     */
+    std::vector<number> one_over_lengths;
+
+    /**
+     * A variable storing the index of the current polynomial in the case that
+     * piecewise linear polynomial on varying subintervals is requested.
+     */
+    unsigned int index;
   };
 
 
@@ -186,6 +213,14 @@ namespace Polynomials
     const unsigned int n_subdivisions,
     const unsigned int base_degree);
 
+  /**
+   * Generates a complete linear basis on a subdivision of the unit interval
+   * in smaller intervals for a given vector of points.
+   */
+  std::vector<PiecewisePolynomial<double>>
+  generate_complete_linear_basis_on_subdivisions(
+    const std::vector<Point<1>> &points);
+
 } // namespace Polynomials
 
 
@@ -199,6 +234,8 @@ namespace Polynomials
   inline unsigned int
   PiecewisePolynomial<number>::degree() const
   {
+    if (points.size() > 0)
+      return 1;
     return polynomial.degree();
   }
 
@@ -208,6 +245,20 @@ namespace Polynomials
   inline number
   PiecewisePolynomial<number>::value(const number x) const
   {
+    if (points.size() > 0)
+      {
+        if (x > points[index])
+          return std::max<number>(0.0,
+                                  1.0 - (x - points[index]) *
+                                          one_over_lengths[index]);
+        else if (x < points[index])
+          return std::max<number>(0.0,
+                                  0.0 + (x - points[index - 1]) *
+                                          one_over_lengths[index - 1]);
+        else
+          return 1.0;
+      }
+
     AssertIndexRange(interval, n_intervals);
     number y = x;
     // shift polynomial if necessary
@@ -256,8 +307,10 @@ namespace Polynomials
     ar &n_intervals;
     ar &interval;
     ar &spans_two_intervals;
+    ar &points;
+    ar &one_over_lengths;
+    ar &index;
   }
-
 } // namespace Polynomials
 
 DEAL_II_NAMESPACE_CLOSE

include/deal.II/fe/fe_q_iso_q1.h

@@ -130,6 +130,11 @@ class FE_Q_iso_Q1 : public FE_Q_Base<dim, spacedim>
    */
   FE_Q_iso_Q1(const unsigned int n_subdivisions);
 
+  /**
+   * Construct a FE_Q_iso_Q1 element with a given vector of support points.
+   */
+  FE_Q_iso_Q1(const std::vector<Point<1>> &support_points);
+
   /**
    * Return a string that uniquely identifies a finite element. This class
    * returns <tt>FE_Q_iso_q1<dim>(equivalent_degree)</tt>, with @p dim and @p

source/base/polynomials_piecewise.cc

@@ -33,13 +33,38 @@ namespace Polynomials
     , n_intervals(n_intervals)
     , interval(interval)
     , spans_two_intervals(spans_next_interval)
+    , index(numbers::invalid_unsigned_int)
   {
     Assert(n_intervals > 0, ExcMessage("No intervals given"));
     AssertIndexRange(interval, n_intervals);
   }
 
 
 
+  template <typename number>
+  PiecewisePolynomial<number>::PiecewisePolynomial(
+    const std::vector<Point<1, number>> &points,
+    const unsigned int                   index)
+    : n_intervals(numbers::invalid_unsigned_int)
+    , interval(numbers::invalid_unsigned_int)
+    , spans_two_intervals(false)
+    , index(index)
+  {
+    Assert(points.size() > 1, ExcMessage("No enough points given!"));
+    AssertIndexRange(index, points.size());
+
+    this->points.resize(points.size());
+    for (unsigned int i = 0; i < points.size(); ++i)
+      this->points[i] = points[i][0];
+
+    this->one_over_lengths.resize(points.size() - 1);
+    for (unsigned int i = 0; i < points.size() - 1; ++i)
+      this->one_over_lengths[i] =
+        number(1.0) / (points[i + 1][0] - points[i][0]);
+  }
+
+
+
   template <typename number>
   void
   PiecewisePolynomial<number>::value(const number         x,
@@ -58,6 +83,36 @@ namespace Polynomials
                                      const unsigned int n_derivatives,
                                      number *           values) const
   {
+    if (points.size() > 0)
+      {
+        if (x > points[index])
+          values[0] = std::max<number>(0.0,
+                                       1.0 - (x - points[index]) *
+                                               one_over_lengths[index]);
+        else if (x < points[index])
+          values[0] = std::max<number>(0.0,
+                                       0.0 + (x - points[index - 1]) *
+                                               one_over_lengths[index - 1]);
+        else
+          values[0] = 1.0;
+
+        if (n_derivatives >= 1)
+          {
+            if ((x > points[index]) && (points[index + 1] >= x))
+              values[1] = -1.0 * one_over_lengths[index];
+            else if ((x < points[index]) && (points[index - 1] <= x))
+              values[1] = +1.0 * one_over_lengths[index - 1];
+            else
+              values[1] = 0.0;
+          }
+
+        // all other derivatives are zero
+        for (unsigned int i = 2; i <= n_derivatives; ++i)
+          values[i] = 0.0;
+
+        return;
+      }
+
     // shift polynomial if necessary
     number y                      = x;
     double derivative_change_sign = 1.;
@@ -125,7 +180,9 @@ namespace Polynomials
     return (polynomial.memory_consumption() +
             MemoryConsumption::memory_consumption(n_intervals) +
             MemoryConsumption::memory_consumption(interval) +
-            MemoryConsumption::memory_consumption(spans_two_intervals));
+            MemoryConsumption::memory_consumption(spans_two_intervals) +
+            MemoryConsumption::memory_consumption(points) +
+            MemoryConsumption::memory_consumption(index));
   }
 
 
@@ -153,6 +210,21 @@ namespace Polynomials
     return p;
   }
 
+
+
+  std::vector<PiecewisePolynomial<double>>
+  generate_complete_linear_basis_on_subdivisions(
+    const std::vector<Point<1>> &points)
+  {
+    std::vector<PiecewisePolynomial<double>> p;
+    p.reserve(points.size());
+
+    for (unsigned int s = 0; s < points.size(); ++s)
+      p.emplace_back(points, s);
+
+    return p;
+  }
+
 } // namespace Polynomials
 
 // ------------------ explicit instantiations --------------- //

source/fe/fe_q_iso_q1.cc

@@ -55,6 +55,28 @@ FE_Q_iso_Q1<dim, spacedim>::FE_Q_iso_Q1(const unsigned int subdivisions)
 
 
 
+template <int dim, int spacedim>
+FE_Q_iso_Q1<dim, spacedim>::FE_Q_iso_Q1(
+  const std::vector<Point<1>> &support_points)
+  : FE_Q_Base<dim, spacedim>(
+      TensorProductPolynomials<dim, Polynomials::PiecewisePolynomial<double>>(
+        Polynomials::generate_complete_linear_basis_on_subdivisions(
+          support_points)),
+      FiniteElementData<dim>(this->get_dpo_vector(support_points.size() - 1),
+                             1,
+                             support_points.size() - 1,
+                             FiniteElementData<dim>::H1),
+      std::vector<bool>(1, false))
+{
+  Assert(support_points.size() > 1,
+         ExcMessage("This element can only be used with a positive number of "
+                    "subelements"));
+
+  this->initialize(support_points);
+}
+
+
+
 template <int dim, int spacedim>
 std::string
 FE_Q_iso_Q1<dim, spacedim>::get_name() const

tests/fe/shapes_q_iso_q1.cc

@@ -42,6 +42,16 @@ plot_FE_Q_shape_functions()
   plot_face_shape_functions(m, q2, "Q2_iso_Q1");
   test_compute_functions(m, q2, "Q2_iso_Q1");
 
+  FE_Q_iso_Q1<dim> q1_gl(QGaussLobatto<1>(2).get_points());
+  plot_shape_functions(m, q1_gl, "Q1_gl");
+  plot_face_shape_functions(m, q1_gl, "Q1_gl");
+  test_compute_functions(m, q1_gl, "Q1_gl");
+
+  FE_Q_iso_Q1<dim> q2_gl(QGaussLobatto<1>(3).get_points());
+  plot_shape_functions(m, q2_gl, "Q2_iso_Q1_gl");
+  plot_face_shape_functions(m, q2_gl, "Q2_iso_Q1_gl");
+  test_compute_functions(m, q2_gl, "Q2_iso_Q1_gl");
+
   // skip the following tests to
   // reduce run-time
   if (dim < 3)
@@ -55,6 +65,16 @@ plot_FE_Q_shape_functions()
       plot_shape_functions(m, q4, "Q4_iso_Q1");
       plot_face_shape_functions(m, q4, "Q4_iso_Q1");
       test_compute_functions(m, q4, "Q4_iso_Q1");
+
+      FE_Q_iso_Q1<dim> q3_gl(QGaussLobatto<1>(4).get_points());
+      plot_shape_functions(m, q3_gl, "Q3_iso_Q1_gl");
+      plot_face_shape_functions(m, q3_gl, "Q3_iso_Q1_gl");
+      test_compute_functions(m, q3_gl, "Q3_iso_Q1_gl");
+
+      FE_Q_iso_Q1<dim> q4_gl(QGaussLobatto<1>(5).get_points());
+      plot_shape_functions(m, q4_gl, "Q4_iso_Q1_gl");
+      plot_face_shape_functions(m, q4_gl, "Q4_iso_Q1_gl");
+      test_compute_functions(m, q4_gl, "Q4_iso_Q1_gl");
     };
 }