quadrature: Add pfasst::augment_nodes and handle is_proper flags appropriately.

Author: Matthew Emmett

include/pfasst/encap/encap_sweeper.hpp

@@ -25,6 +25,7 @@ namespace pfasst
     {
         //! @{
         vector<time> nodes;
+        vector<bool> is_proper;
         shared_ptr<EncapFactory<time>> factory;
         //! @}
 
@@ -69,14 +70,21 @@ namespace pfasst
         //! @{
         void set_nodes(vector<time> nodes)
         {
-          this->nodes = nodes;
+          auto augmented = pfasst::augment_nodes(nodes);
+          this->nodes = get<0>(augmented);
+          this->is_proper = get<1>(augmented);
         }
 
         const vector<time> get_nodes() const
         {
           return nodes;
         }
 
+        const vector<bool> get_is_proper() const
+        {
+          return is_proper;
+        }
+
         void set_factory(shared_ptr<EncapFactory<time>> factory)
         {
           this->factory = shared_ptr<EncapFactory<time>>(factory);

include/pfasst/encap/imex_sweeper.hpp

@@ -30,9 +30,9 @@ namespace pfasst
      * explicitly and the stiff part implicitly.
      * Therefore, we define the splitting \\( F(t,u) = F_{expl}(t,u) + F_{impl}(t,u) \\).
      *
-     * This sweeper requires three interfaces to be implement for such ODEs: two routines to 
-     * evaluate the explicit \\( F_{\\rm expl} \\) and implicit \\( F_{\\rm impl} \\) pieces for a 
-     * given state, and one that solves (perhaps with an external solver) the backward-Euler 
+     * This sweeper requires three interfaces to be implement for such ODEs: two routines to
+     * evaluate the explicit \\( F_{\\rm expl} \\) and implicit \\( F_{\\rm impl} \\) pieces for a
+     * given state, and one that solves (perhaps with an external solver) the backward-Euler
      * equation \\( U^{n+1} - \\Delta t F_{\\rm impl}(U^{n+1}) = RHS \\) for \\( U^{n+1} \\).
      *
      * @tparam time precision type of the time dimension
@@ -173,9 +173,10 @@ namespace pfasst
         virtual void setup(bool coarse) override
         {
           auto nodes = this->get_nodes();
+          auto is_proper = this->get_is_proper();
           assert(nodes.size() >= 1);
 
-          this->s_mat = compute_quadrature(nodes, nodes, 's');
+          this->s_mat = compute_quadrature(nodes, nodes, is_proper, 's');
 
           this->s_mat_expl = this->s_mat;
           this->s_mat_impl = this->s_mat;

include/pfasst/quadrature.hpp

@@ -210,8 +210,27 @@ namespace pfasst
     return nodes;
   }
 
+  template<typename node>
+  auto augment_nodes(vector<node> const orig) -> pair<vector<node>,vector<bool>>
+  {
+    vector<node> nodes = orig;
+
+    bool left = nodes.front() == node(0.0);
+    bool right = nodes.back() == node(1.0);
+
+    if (!left)  { nodes.insert(nodes.begin(), node(0.0)); }
+    if (!right) { nodes.insert(nodes.end(),   node(1.0)); }
+
+    vector<bool> is_proper(nodes.size(), true);
+    is_proper.front() = left;
+    is_proper.back() = right;
+
+    return pair<vector<node>,vector<bool>>(nodes, is_proper);
+  }
+
+
   template<typename node = time_precision>
-  matrix<node> compute_quadrature(vector<node> dst, vector<node> src, char type)
+  matrix<node> compute_quadrature(vector<node> dst, vector<node> src, vector<bool> is_proper, char type)
   {
     const size_t ndst = dst.size();
     const size_t nsrc = src.size();
@@ -222,11 +241,13 @@ namespace pfasst
     Polynomial<node> p(nsrc + 1), p1(nsrc + 1);
 
     for (size_t i = 0; i < nsrc; i++) {
+      if (!is_proper[i]) { continue; }
+
       // construct interpolating polynomial coefficients
       p[0] = 1.0;
       for (size_t j = 1; j < nsrc + 1; j++) { p[j] = 0.0; }
       for (size_t m = 0; m < nsrc; m++) {
-        if (m == i) { continue; }
+        if ((!is_proper[m]) || (m == i)) { continue; }
 
         // p_{m+1}(x) = (x - x_j) * p_m(x)
         p1[0] = 0.0;

tests/test_quadrature.cpp

@@ -222,7 +222,8 @@ TEST(NodesTest, UniformNodes)
 TEST(QuadratureTest, GaussLobattoNodes)
 {
   auto l3 = pfasst::compute_nodes<long double>(3, "gauss-lobatto");
-  auto s3 = pfasst::compute_quadrature(l3, l3, 's');
+  auto a3 = pfasst::augment_nodes(l3);
+  auto s3 = pfasst::compute_quadrature(get<0>(a3), get<0>(a3), get<1>(a3), 's');
   const long double s3e[6] = { 0.20833333333333333,
                                0.33333333333333333,
                                -0.04166666666666666,
@@ -234,7 +235,8 @@ TEST(QuadratureTest, GaussLobattoNodes)
   EXPECT_THAT(s3.data(), testing::Pointwise(DoubleNear(), s3e));
 
   auto l5 = pfasst::compute_nodes<long double>(5, "gauss-lobatto");
-  auto s5 = pfasst::compute_quadrature(l5, l5, 's');
+  auto a5 = pfasst::augment_nodes(l5);
+  auto s5 = pfasst::compute_quadrature(get<0>(a5), get<0>(a5), get<1>(a5), 's');
   const long double s5e[] = { 0.067728432186156897969267419174073482,
                               0.11974476934341168251615379970493965,
                               -0.021735721866558113665511351745074292,
@@ -262,7 +264,8 @@ TEST(QuadratureTest, GaussLobattoNodes)
 TEST(QuadratureTest, ClenshawCurtisNodes)
 {
   auto c4 = pfasst::compute_nodes<long double>(4, "clenshaw-curtis");
-  auto s4 = pfasst::compute_quadrature(c4, c4, 's');
+  auto a4 = pfasst::augment_nodes(c4);
+  auto s4 = pfasst::compute_quadrature(get<0>(a4), get<0>(a4), get<1>(a4), 's');
   const long double s4e[] = { 0.10243055555555555555555555555555556,
                                0.16319444444444444444444444444444444,
                               -0.024305555555555555555555555555555556,