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 67d34f3b408da1f35c840481a41b6cadb6bf6e9f..9d459f14f9c6fd61de6fc367e9cf2a4afcc4729f 100644 (file)
--- 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<dim>    &p,
               const unsigned int component = 0) const;
 
-  private:
+  protected:
+    /**
+     * Find the index in the table of the rectangle containing an input point
+     */
+    TableIndices<dim> table_index_of_point (const Point<dim> &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 b2763657c9e254969ba17fedc45a7de5a2cfb99a..44aacc114b5e033ed06e2094341c2af64fa8e156 100644 (file)
--- a/source/base/function_lib.cc
+++ b/source/base/function_lib.cc
@@ -2445,14 +2445,9 @@ namespace Functions
 
 
   template <int dim>
-  double
-  InterpolatedTensorProductGridData<dim>::value(const Point<dim> &p,
-                                                const unsigned int component) const
+  TableIndices<dim>
+  InterpolatedTensorProductGridData<dim>::table_index_of_point(const Point<dim> &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 <int dim>
+  double
+  InterpolatedTensorProductGridData<dim>::value(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 the index in the data table of the cell containing the input point
+    const TableIndices<dim> 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<dim> ix;
-    for (unsigned int d=0; d<dim; ++d)
-      {
-        ix[d] = (std::lower_bound (coordinate_values[d].begin(),
-                                   coordinate_values[d].end(),
-                                   p[d])
-                 - coordinate_values[d].begin());
-
-        if (ix[d] == coordinate_values[d].size())
-          ix[d] = coordinate_values[d].size()-2;
-        else if (ix[d] > 0)
-          --ix[d];
-      }
+    const TableIndices<dim> ix = table_index_of_point (p);
 
     Point<dim> dx;
     for (unsigned int d=0; d<dim; ++d)