From: danshapero Date: Thu, 11 Feb 2016 22:10:22 +0000 (-0800) Subject: Factored out repeated code in interpolated grid data X-Git-Tag: v8.5.0-rc1~1320^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F2173%2Fhead;p=dealii.git Factored out repeated code in interpolated grid data InterpolatedTensorProductGridData has some repeated code for finding the right table index for an input point when computed the value or gradient; this has been factored out into a protected member function and the member data have been made protected instead of private for ease of extension. --- diff --git a/include/deal.II/base/function_lib.h b/include/deal.II/base/function_lib.h index 67d34f3b40..9d459f14f9 100644 --- a/include/deal.II/base/function_lib.h +++ b/include/deal.II/base/function_lib.h @@ -1170,7 +1170,12 @@ namespace Functions gradient (const Point &p, const unsigned int component = 0) const; - private: + protected: + /** + * Find the index in the table of the rectangle containing an input point + */ + TableIndices table_index_of_point (const Point &p) const; + /** * The set of coordinate values in each of the coordinate directions. */ diff --git a/source/base/function_lib.cc b/source/base/function_lib.cc index b2763657c9..44aacc114b 100644 --- a/source/base/function_lib.cc +++ b/source/base/function_lib.cc @@ -2445,14 +2445,9 @@ namespace Functions template - double - InterpolatedTensorProductGridData::value(const Point &p, - const unsigned int component) const + TableIndices + InterpolatedTensorProductGridData::table_index_of_point(const Point &p) 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 // points. if we run all the way to the end of the range, // set the indices so that we will simply query the last of the @@ -2462,9 +2457,11 @@ namespace Functions { // get the index of the first element of the coordinate arrays that is // larger than p[d] - ix[d] = (std::lower_bound (coordinate_values[d].begin(), coordinate_values[d].end(), + ix[d] = (std::lower_bound (coordinate_values[d].begin(), + coordinate_values[d].end(), p[d]) - coordinate_values[d].begin()); + // the one we want is the index of the coordinate to the left, however, // so decrease it by one (unless we have a point to the left of all, in which // case we stay where we are; the formulas below are made in a way that allow @@ -2478,6 +2475,23 @@ namespace Functions --ix[d]; } + return ix; + } + + + + template + double + InterpolatedTensorProductGridData::value(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 the index in the data table of the cell containing the input point + const TableIndices ix = table_index_of_point (p); + // 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 @@ -2503,19 +2517,7 @@ namespace Functions ExcMessage ("This is a scalar function object, the component can only be zero.")); // find out where this data point lies - TableIndices ix; - for (unsigned int d=0; d 0) - --ix[d]; - } + const TableIndices ix = table_index_of_point (p); Point dx; for (unsigned int d=0; d