From 4ceadbbaffa17c812a7c3b426204c49fa7a89fa4 Mon Sep 17 00:00:00 2001
From: Bob Myhill <myhill.bob@gmail.com>
Date: Tue, 28 May 2019 20:39:51 -0600
Subject: [PATCH] Add InterpolatedUniformGridData::gradient().

---
 include/deal.II/base/function_lib.h | 15 ++++++++
 source/base/function_lib.cc         | 53 +++++++++++++++++++++++++++++
 2 files changed, 68 insertions(+)

diff --git a/include/deal.II/base/function_lib.h b/include/deal.II/base/function_lib.h
index 6b59bc11e9..2c539ab93e 100644
--- a/include/deal.II/base/function_lib.h
+++ b/include/deal.II/base/function_lib.h
@@ -1562,6 +1562,21 @@ namespace Functions
     virtual double
     value(const Point<dim> &p, const unsigned int component = 0) const override;
 
+    /**
+     * Compute the gradient of the function set by bilinear interpolation of the
+     * given data set.
+     *
+     * @param p The point at which the function is to be evaluated.
+     * @param component The vector component. Since this function is scalar,
+     *   only zero is a valid argument here.
+     * @return The gradient of the interpolated function at this point. If the
+     *   point lies outside the set of coordinates, the function is extended
+     *   by a constant whose gradient is then of course zero.
+     */
+    virtual Tensor<1, dim>
+    gradient(const Point<dim> & p,
+             const unsigned int component = 0) const override;
+
   private:
     /**
      * The set of interval endpoints in each of the coordinate directions.
diff --git a/source/base/function_lib.cc b/source/base/function_lib.cc
index 2dac6bbb64..55ea5ccd4a 100644
--- a/source/base/function_lib.cc
+++ b/source/base/function_lib.cc
@@ -2670,6 +2670,59 @@ namespace Functions
     return interpolate(data_values, ix, p_unit);
   }
 
+
+
+  template <int dim>
+  Tensor<1, dim>
+  InterpolatedUniformGridData<dim>::gradient(const Point<dim> & p,
+                                             const unsigned int component) const
+  {
+    (void)component;
+    Assert(
+      component == 0,
+      ExcMessage(
+        "This is a scalar function object, the component can only be zero."));
+
+    // find out where this data point lies, relative to the given
+    // subdivision points
+    TableIndices<dim> ix;
+    for (unsigned int d = 0; d < dim; ++d)
+      {
+        const double delta_x = ((this->interval_endpoints[d].second -
+                                 this->interval_endpoints[d].first) /
+                                this->n_subintervals[d]);
+        if (p[d] <= this->interval_endpoints[d].first)
+          ix[d] = 0;
+        else if (p[d] >= this->interval_endpoints[d].second - delta_x)
+          ix[d] = this->n_subintervals[d] - 1;
+        else
+          ix[d] = static_cast<unsigned int>(
+            (p[d] - this->interval_endpoints[d].first) / delta_x);
+      }
+
+    // now compute the relative point within the interval/rectangle/box
+    // defined by the point coordinates found above. truncate below and
+    // above to accommodate points that may lie outside the range
+    Point<dim> p_unit;
+    Point<dim> delta_x;
+    for (unsigned int d = 0; d < dim; ++d)
+      {
+        delta_x[d] = ((this->interval_endpoints[d].second -
+                       this->interval_endpoints[d].first) /
+                      this->n_subintervals[d]);
+        p_unit[d] =
+          std::max(std::min((p[d] - this->interval_endpoints[d].first -
+                             ix[d] * delta_x[d]) /
+                              delta_x[d],
+                            1.),
+                   0.);
+      }
+
+    return gradient_interpolate(this->data_values, ix, p_unit, delta_x);
+  }
+
+
+
   /* ---------------------- Polynomial ----------------------- */
 
 
-- 
2.39.5