From 4ceadbbaffa17c812a7c3b426204c49fa7a89fa4 Mon Sep 17 00:00:00 2001 From: Bob Myhill 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 &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 & 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 + Tensor<1, dim> + InterpolatedUniformGridData::gradient(const Point & 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 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( + (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 p_unit; + Point 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