]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add InterpolatedUniformGridData::gradient().
authorBob Myhill <myhill.bob@gmail.com>
Wed, 29 May 2019 02:39:51 +0000 (20:39 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Wed, 29 May 2019 02:39:51 +0000 (20:39 -0600)
include/deal.II/base/function_lib.h
source/base/function_lib.cc

index 6b59bc11e94bb2f4a1858bcc939a6670138b4ae8..2c539ab93e5414c2661e39d97845e6eca5177199 100644 (file)
@@ -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.
index 2dac6bbb6481535c569a226c7551e4e65faf1909..55ea5ccd4ab322edfff984b7a9675c0c7606f0fa 100644 (file)
@@ -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 ----------------------- */
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.