From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sat, 21 Dec 2013 04:47:05 +0000 (+0000)
Subject: New class InterpolatedUniformGridData.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=479bf9a83588e2371a8364626025d37ee3463167;p=dealii-svn.git

New class InterpolatedUniformGridData.

git-svn-id: https://svn.dealii.org/trunk@32072 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 33b62f9f15..a9a29235b4 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -64,7 +64,9 @@ inconvenience this causes.
   mesh of $x$ (and $y$ and $z$) values, for example to evaluate experimentally
   determined coefficients, or to assess the accuracy of a solution by
   comparing with a solution generated by a different code and written in
-  gridded data.
+  gridded data. There is also a new class InterpolatedUniformGridData that
+  can perform the same task more efficiently if the data is stored on meshes
+  that are uniform in each coordinate direction.
   <br>
   (Wolfgang Bangerth, 2013/12/20)
   </li>
diff --git a/deal.II/include/deal.II/base/function_lib.h b/deal.II/include/deal.II/base/function_lib.h
index 8ec2a3c11c..19f3824c64 100644
--- a/deal.II/include/deal.II/base/function_lib.h
+++ b/deal.II/include/deal.II/base/function_lib.h
@@ -23,7 +23,7 @@
 #include <deal.II/base/point.h>
 #include <deal.II/base/table.h>
 
-#include <boost/array.hpp>
+#include <deal.II/base/std_cxx1x/array.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -1229,8 +1229,8 @@ namespace Functions
        *   converting other data types into a table where you specify this
        *   argument.
        */
-      InterpolatedTensorProductGridData (const boost::array<std::vector<double>,dim> &coordinate_values,
-                                         const Table<dim,double>                     &data_values);
+      InterpolatedTensorProductGridData (const std_cxx1x::array<std::vector<double>,dim> &coordinate_values,
+                                         const Table<dim,double>                         &data_values);
 
       /**
        * Compute the value of the function set by bilinear interpolation of the
@@ -1251,7 +1251,96 @@ namespace Functions
       /**
        * The set of coordinate values in each of the coordinate directions.
        */
-      const boost::array<std::vector<double>,dim> coordinate_values;
+      const std_cxx1x::array<std::vector<double>,dim> coordinate_values;
+
+      /**
+       * The data that is to be interpolated.
+       */
+      const Table<dim,double>                     data_values;
+  };
+
+
+  /**
+   * A scalar function that computes its values by (bi-, tri-)linear interpolation
+   * from a set of point data that are arranged on a uniformly spaced
+   * tensor product mesh. In other words, considering the three-dimensional case,
+   * let there be points $x_0,\ldotx, x_{K-1}$ that result from a uniform subdivision
+   * of the interval $[x_0,x_{K-1}]$ into $K-1$ sub-intervals of size $\Delta x
+   * = (x_{K-1}-x_0)/(K-1)$, and similarly
+   * $y_0,\ldots,y_{L-1}$, $z_1,\ldots,z_{M-1}$. Also consider data $d_{klm}$ defined at
+   * point $(x_k,y_l,z_m)^T$, then evaluating the function at a point
+   * $\mathbf x=(x,y,z)$ will find the box so that
+   * $x_k\le x\le x_{k+1}, y_l\le x\le y_{l+1}, z_m\le z\le z_{m+1}$, and do a
+   * trilinear interpolation of the data on this cell. Similar operations are
+   * done in lower dimensions.
+   *
+   * This class is most often used for either evaluating coefficients or right
+   * hand sides that are provided experimentally at a number of points inside the
+   * domain, or for comparing outputs of a solution on a finite element mesh
+   * against previously obtained data defined on a grid.
+   *
+   * @note If you have a problem where the points $x_i$ are not equally spaced
+   * (e.g., they result from a computation on a graded mesh that is denser
+   * closer to one boundary), then use the InterpolatedTensorProductGridData
+   * class instead.
+   *
+   * If a point is requested outside the box defined by the end points of the
+   * coordinate arrays, then the function is assumed to simply extend by
+   * constant values beyond the last data point in each coordinate
+   * direction. (The class does not throw an error if a point lies outside the
+   * box since it frequently happens that a point lies just outside the box
+   * by an amount on the order of numerical roundoff.)
+   *
+   * @author Wolfgang Bangerth, 2013
+   */
+  template <int dim>
+  class InterpolatedUniformGridData : public Function<dim>
+  {
+    public:
+      /**
+       * Constructor
+       * @param interval_endpoints The left and right end points of the (uniformly
+       *   subdivided) intervals in each of the coordinate directions.
+       * @param n_subdivisions The number of subintervals of the subintervals
+       *   in each coordinate direction. A value of one for a coordinate
+       *   means that the interval is considered as one subinterval consisting
+       *   of the entire range. A value of two means that there are two subintervals
+       *   each with one half of the range, etc.
+       * @param data_values A dim-dimensional table of data at each of the
+       *   mesh points defined by the coordinate arrays above. Note that the
+       *   Table class has a number of conversion constructors that allow
+       *   converting other data types into a table where you specify this
+       *   argument.
+       */
+      InterpolatedUniformGridData (const std_cxx1x::array<std::pair<double,double>,dim> &interval_endpoints,
+                                   const std_cxx1x::array<unsigned int,dim>             &n_subintervals,
+                                   const Table<dim,double>                              &data_values);
+
+      /**
+       * Compute the value 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 interpolated value at this point. If the point lies outside
+       *   the set of coordinates, the function is extended by a constant.
+       */
+      virtual
+      double
+      value (const Point<dim> &p,
+             const unsigned int component = 0) const;
+
+    private:
+      /**
+       * The set of interval endpoints in each of the coordinate directions.
+       */
+      const std_cxx1x::array<std::pair<double,double>,dim> interval_endpoints;
+
+      /**
+       * The number of subintervals in each of the coordinate directions.
+       */
+      const std_cxx1x::array<unsigned int,dim>             n_subintervals;
 
       /**
        * The data that is to be interpolated.
diff --git a/deal.II/source/base/function_lib.cc b/deal.II/source/base/function_lib.cc
index 710ace921d..8cdbc261f0 100644
--- a/deal.II/source/base/function_lib.cc
+++ b/deal.II/source/base/function_lib.cc
@@ -2266,29 +2266,6 @@ namespace Functions
 
 
 
-  template <int dim>
-  InterpolatedTensorProductGridData<dim>::
-  InterpolatedTensorProductGridData(const boost::array<std::vector<double>,dim> &coordinate_values,
-                                    const Table<dim,double>                     &data_values)
-  :
-  coordinate_values (coordinate_values),
-  data_values (data_values)
-  {
-    for (unsigned int d=0; d<dim; ++d)
-      {
-        Assert (coordinate_values[d].size() >= 2,
-                ExcMessage ("Coordinate arrays must have at least two coordinate values!"));
-        for (unsigned int i=0; i<coordinate_values[d].size()-1; ++i)
-          Assert (coordinate_values[d][i] < coordinate_values[d][i+1],
-              ExcMessage ("Coordinate arrays must be sorted in strictly ascending order."));
-
-        Assert (data_values.size()[d] == coordinate_values[d].size(),
-                ExcMessage ("Data and coordinate tables do not have the same size."));
-      }
-  }
-
-
-
   namespace
   {
     // interpolate a data value from a table where ix denotes
@@ -2308,12 +2285,12 @@ namespace Functions
                         const Point<2>        &p_unit)
     {
       return (((1-p_unit[0])*data_values[ix[0]][ix[1]]
-	       +
-	       p_unit[0]*data_values[ix[0]+1][ix[1]])*(1-p_unit[1])
-	      +
-	      ((1-p_unit[0])*data_values[ix[0]][ix[1]+1]
-	       +
-	       p_unit[0]*data_values[ix[0]+1][ix[1]+1])*p_unit[1]);
+               +
+               p_unit[0]*data_values[ix[0]+1][ix[1]])*(1-p_unit[1])
+              +
+              ((1-p_unit[0])*data_values[ix[0]][ix[1]+1]
+               +
+               p_unit[0]*data_values[ix[0]+1][ix[1]+1])*p_unit[1]);
     }
 
     double interpolate (const Table<3,double> &data_values,
@@ -2323,22 +2300,45 @@ namespace Functions
       return ((((1-p_unit[0])*data_values[ix[0]][ix[1]][ix[2]]
                 +
                 p_unit[0]*data_values[ix[0]+1][ix[1]][ix[2]])*(1-p_unit[1])
-	       +
-	       ((1-p_unit[0])*data_values[ix[0]][ix[1]+1][ix[2]]
-		+
-		p_unit[0]*data_values[ix[0]+1][ix[1]+1][ix[2]])*p_unit[1]) * (1-p_unit[2])
-	      +
-	      (((1-p_unit[0])*data_values[ix[0]][ix[1]][ix[2]+1]
-		+
-		p_unit[0]*data_values[ix[0]+1][ix[1]][ix[2]+1])*(1-p_unit[1])
-	       +
-	       ((1-p_unit[0])*data_values[ix[0]][ix[1]+1][ix[2]+1]
-		+
-		p_unit[0]*data_values[ix[0]+1][ix[1]+1][ix[2]+1])*p_unit[1]) * p_unit[2]);
+               +
+               ((1-p_unit[0])*data_values[ix[0]][ix[1]+1][ix[2]]
+                +
+                p_unit[0]*data_values[ix[0]+1][ix[1]+1][ix[2]])*p_unit[1]) * (1-p_unit[2])
+              +
+              (((1-p_unit[0])*data_values[ix[0]][ix[1]][ix[2]+1]
+                +
+                p_unit[0]*data_values[ix[0]+1][ix[1]][ix[2]+1])*(1-p_unit[1])
+               +
+               ((1-p_unit[0])*data_values[ix[0]][ix[1]+1][ix[2]+1]
+                +
+                p_unit[0]*data_values[ix[0]+1][ix[1]+1][ix[2]+1])*p_unit[1]) * p_unit[2]);
     }
   }
 
 
+  template <int dim>
+  InterpolatedTensorProductGridData<dim>::
+  InterpolatedTensorProductGridData(const std_cxx1x::array<std::vector<double>,dim> &coordinate_values,
+                                    const Table<dim,double>                         &data_values)
+  :
+  coordinate_values (coordinate_values),
+  data_values (data_values)
+  {
+    for (unsigned int d=0; d<dim; ++d)
+      {
+        Assert (coordinate_values[d].size() >= 2,
+                ExcMessage ("Coordinate arrays must have at least two coordinate values!"));
+        for (unsigned int i=0; i<coordinate_values[d].size()-1; ++i)
+          Assert (coordinate_values[d][i] < coordinate_values[d][i+1],
+              ExcMessage ("Coordinate arrays must be sorted in strictly ascending order."));
+
+        Assert (data_values.size()[d] == coordinate_values[d].size(),
+                ExcMessage ("Data and coordinate tables do not have the same size."));
+      }
+  }
+
+
+
   template <int dim>
   double
   InterpolatedTensorProductGridData<dim>::value(const Point<dim> &p,
@@ -2388,6 +2388,70 @@ namespace Functions
 
 
 
+  template <int dim>
+  InterpolatedUniformGridData<dim>::
+  InterpolatedUniformGridData(const std_cxx1x::array<std::pair<double,double>,dim> &interval_endpoints,
+                              const std_cxx1x::array<unsigned int,dim>             &n_subintervals,
+                              const Table<dim,double>                              &data_values)
+  :
+  interval_endpoints (interval_endpoints),
+  n_subintervals (n_subintervals),
+  data_values (data_values)
+  {
+    for (unsigned int d=0; d<dim; ++d)
+      {
+        Assert (n_subintervals[d] >= 1,
+                ExcMessage ("There needs to be at least one subinterval in each "
+                            "coordinate direction."));
+        Assert (interval_endpoints[d].first < interval_endpoints[d].second,
+                ExcMessage ("The interval in each coordinate direction needs "
+                            "to have positive size"));
+        Assert (data_values.size()[d] == n_subintervals[d]+1,
+                ExcMessage ("The data table does not have the correct size."));
+      }
+  }
+
+
+  template <int dim>
+  double
+  InterpolatedUniformGridData<dim>::value(const Point<dim> &p,
+                                          const unsigned int component) const
+  {
+    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 = ((interval_endpoints[d].second - interval_endpoints[d].first) /
+                                n_subintervals[d]);
+        if (p[d] <= interval_endpoints[d].first)
+          ix[d] = 0;
+        else if (p[d] >= interval_endpoints[d].second-delta_x)
+          ix[d] = n_subintervals[d]-1;
+        else
+          ix[d] = (unsigned int)((p[d]-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;
+    for (unsigned int d=0; d<dim; ++d)
+      {
+        const double delta_x = ((interval_endpoints[d].second - interval_endpoints[d].first) /
+                                n_subintervals[d]);	
+	p_unit[d] = std::max(std::min((p[d]-interval_endpoints[d].first-ix[d]*delta_x)/delta_x,
+				      1.),
+			     0.);
+      }
+    
+    return interpolate (data_values, ix, p_unit);
+  }
+
+
 
 // explicit instantiations
   template class SquareFunction<1>;
@@ -2433,6 +2497,9 @@ namespace Functions
   template class InterpolatedTensorProductGridData<1>;
   template class InterpolatedTensorProductGridData<2>;
   template class InterpolatedTensorProductGridData<3>;
+  template class InterpolatedUniformGridData<1>;
+  template class InterpolatedUniformGridData<2>;
+  template class InterpolatedUniformGridData<3>;
 }
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/base/functions_10.cc b/tests/base/functions_10.cc
index 6c3555a45e..0f6714c0a5 100644
--- a/tests/base/functions_10.cc
+++ b/tests/base/functions_10.cc
@@ -23,7 +23,7 @@
 // now interpolate the function x*y*z onto points. note that this function is
 // (bi/tri)linear and so we can later know what the correct value is that the
 // function should provide
-Table<1,double> fill (const boost::array<std::vector<double>,1> &coordinates)
+Table<1,double> fill (const std_cxx1x::array<std::vector<double>,1> &coordinates)
 {
   Table<1,double> data(coordinates[0].size());
   for (unsigned int i=0; i<coordinates[0].size(); ++i)
@@ -31,7 +31,7 @@ Table<1,double> fill (const boost::array<std::vector<double>,1> &coordinates)
   return data;
 }
 
-Table<2,double> fill (const boost::array<std::vector<double>,2> &coordinates)
+Table<2,double> fill (const std_cxx1x::array<std::vector<double>,2> &coordinates)
 {
   Table<2,double> data(coordinates[0].size(),
 		coordinates[1].size());
@@ -41,7 +41,7 @@ Table<2,double> fill (const boost::array<std::vector<double>,2> &coordinates)
   return data;
 }
 
-Table<3,double> fill (const boost::array<std::vector<double>,3> &coordinates)
+Table<3,double> fill (const std_cxx1x::array<std::vector<double>,3> &coordinates)
 {
   Table<3,double> data(coordinates[0].size(),
 		coordinates[1].size(),
@@ -61,7 +61,7 @@ void check ()
 {
   // have coordinate arrays that span an interval starting at d+1
   // d+5 nonuniform intervals
-  boost::array<std::vector<double>,dim> coordinates;
+  std_cxx1x::array<std::vector<double>,dim> coordinates;
   for (unsigned int d=0; d<dim; ++d)
     for (unsigned int i=0; i<d+5; ++i)
       coordinates[d].push_back (d+1 + 1.*i*i);
diff --git a/tests/base/functions_11.cc b/tests/base/functions_11.cc
new file mode 100644
index 0000000000..00dffb5301
--- /dev/null
+++ b/tests/base/functions_11.cc
@@ -0,0 +1,140 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test InterpolatedUniformGridData
+
+#include "../tests.h"
+#include <deal.II/base/function_lib.h>
+
+// now interpolate the function x*y*z onto points. note that this function is
+// (bi/tri)linear and so we can later know what the correct value is that the
+// function should provide
+Table<1,double> fill (const std_cxx1x::array<std::vector<double>,1> &coordinates)
+{
+  Table<1,double> data(coordinates[0].size());
+  for (unsigned int i=0; i<coordinates[0].size(); ++i)
+    data[i] = coordinates[0][i];
+  return data;
+}
+
+Table<2,double> fill (const std_cxx1x::array<std::vector<double>,2> &coordinates)
+{
+  Table<2,double> data(coordinates[0].size(),
+		coordinates[1].size());
+  for (unsigned int i=0; i<coordinates[0].size(); ++i)
+    for (unsigned int j=0; j<coordinates[1].size(); ++j)
+      data[i][j] = coordinates[0][i] * coordinates[1][j];
+  return data;
+}
+
+Table<3,double> fill (const std_cxx1x::array<std::vector<double>,3> &coordinates)
+{
+  Table<3,double> data(coordinates[0].size(),
+		coordinates[1].size(),
+		coordinates[2].size());
+  for (unsigned int i=0; i<coordinates[0].size(); ++i)
+    for (unsigned int j=0; j<coordinates[1].size(); ++j)
+      for (unsigned int k=0; k<coordinates[2].size(); ++k)
+	data[i][j][k] = coordinates[0][i] *
+			coordinates[1][j] *
+			coordinates[2][k];
+  return data;
+}
+
+
+template <int dim>
+void check ()
+{
+  // have coordinate arrays that span an interval starting at d+1
+  // d+5 nonuniform intervals
+  std_cxx1x::array<std::pair<double,double>,dim> intervals;
+  std_cxx1x::array<unsigned int,dim> n_subintervals;
+  for (unsigned int d=0; d<dim; ++d)
+    {
+      intervals[d] = std::make_pair(d+2., 2*d+5.);
+      n_subintervals[d] = d+1 + d*d;
+    }
+
+  std_cxx1x::array<std::vector<double>,dim> coordinates;
+  for (unsigned int d=0; d<dim; ++d)
+    {
+      const double x = intervals[d].first;
+      const double dx = (intervals[d].second-intervals[d].first)/n_subintervals[d];
+      
+      for (unsigned int i=0; i<n_subintervals[d]+1; ++i)
+	coordinates[d].push_back (x+dx*i);
+    }
+  
+  const Table<dim,double> data = fill(coordinates);
+
+  Functions::InterpolatedUniformGridData<dim> f(intervals, n_subintervals, data);
+
+  // now choose a number of randomly chosen points inside the box and
+  // verify that the functions returned are correct
+  for (unsigned int i=0; i<10; ++i)
+    {
+      Point<dim> p;
+      for (unsigned int d=0; d<dim; ++d)
+	p[d] = coordinates[d][0] +
+	       (1. * Testing::rand() / RAND_MAX) * (coordinates[d].back() -
+						   coordinates[d][0]);
+
+      double exact_value = 1;
+      for (unsigned int d=0; d<dim; ++d)
+	exact_value *= p[d];
+      
+      Assert (std::fabs (exact_value - f.value(p)) < 1e-12,
+	      ExcInternalError());
+    }
+
+  // now also verify that it computes values outside the box correctly, as
+  // documented
+  double value_at_bottom_left = 1;
+  for (unsigned int d=0; d<dim; ++d)
+    value_at_bottom_left *= coordinates[d][0];
+  
+  Assert (std::fabs(f.value(Point<dim>()) - value_at_bottom_left) < 1e-12,
+	  ExcInternalError());
+
+  Point<dim> top_right;
+  double value_at_top_right = 1;
+  for (unsigned int d=0; d<dim; ++d)
+    {
+      top_right[d] = 1000;
+      value_at_top_right *= coordinates[d].back();
+    }
+  Assert (std::fabs(f.value(top_right) - value_at_top_right) < 1e-12,
+	  ExcInternalError());
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main()
+{
+  std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  check<1>();
+  check<2>();
+  check<3>();
+}
+
diff --git a/tests/base/functions_11.output b/tests/base/functions_11.output
new file mode 100644
index 0000000000..fb71de2867
--- /dev/null
+++ b/tests/base/functions_11.output
@@ -0,0 +1,4 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK