From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 20 Dec 2013 11:23:53 +0000 (+0000)
Subject: New class InterpolatedTensorProductData.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cbb21676522914d8d6b8b58198596e252d68b55d;p=dealii-svn.git

New class InterpolatedTensorProductData.

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

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index bdde5bfdd4..a69b0d4ead 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -59,6 +59,16 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> New: There is now a new class InterpolatedTensorProductData that can
+  be used to (bi-/tri-)linearly interpolate data given on a tensor product
+  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.
+  <br>
+  (Wolfgang Bangerth, 2013/12/20)
+  </li>
+
   <li> Fixed: ParameterHandler::get_double() and ParameterHandler::get_integer()
   had bugs in that they didn't detect if they were asked to return a number
   for a parameter whose value was in fact not a number but some general
diff --git a/deal.II/include/deal.II/base/function_lib.h b/deal.II/include/deal.II/base/function_lib.h
index 5e49212269..267eaacd34 100644
--- a/deal.II/include/deal.II/base/function_lib.h
+++ b/deal.II/include/deal.II/base/function_lib.h
@@ -21,6 +21,9 @@
 #include <deal.II/base/config.h>
 #include <deal.II/base/function.h>
 #include <deal.II/base/point.h>
+#include <deal.II/base/table.h>
+
+#include <boost/array.hpp>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -1175,6 +1178,86 @@ namespace Functions
     const Tensor<1,dim> exponents;
   };
 
+
+
+  /**
+   * A scalar function that computes its values by (bi-, tri-)linear interpolation
+   * from a set of point data that are arranged on a possibly non-uniform
+   * tensor product mesh. In other words, considering the three-dimensional case,
+   * let there be points $x_0,\ldotx, x_{K-1}$,
+   * $y_0,\ldots,y_{L-1}$, $z_1,\ldots,z_{M-1}$, and 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 the points $x_i$ are actually equally spaced on an interval $[x_0,x_1]$
+   * and the same is true for the other data points in higher dimensions, you should
+   * use the InterpolatedUniformGridData 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 InterpolatedTensorProductData : public Function<dim>
+  {
+    public:
+      /**
+       * Constructor.
+       * @param coordinate_values An array of dim arrays. Each of the inner
+       *   arrays contains the coordinate values $x_0,\ldotx, x_{K-1}$ and
+       *   similarly for the other coordinate directions. These arrays
+       *   need not have the same size. Obviously, we need dim such arrays
+       *   for a dim-dimensional function object. The coordinate values
+       *   within this array are assumed to be strictly ascending to allow
+       *   for efficient lookup.
+       * @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.
+       */
+      InterpolatedTensorProductData (const boost::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
+       * 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 coordinate values in each of the coordinate directions.
+       */
+      const boost::array<std::vector<double>,dim> coordinate_values;
+
+      /**
+       * The data that is to be interpolated.
+       */
+      const Table<dim,double>                     data_values;
+  };
 }
 DEAL_II_NAMESPACE_CLOSE
 
diff --git a/deal.II/source/base/function_lib.cc b/deal.II/source/base/function_lib.cc
index 92d7097165..8ca33e67ba 100644
--- a/deal.II/source/base/function_lib.cc
+++ b/deal.II/source/base/function_lib.cc
@@ -2265,6 +2265,130 @@ namespace Functions
   }
 
 
+
+  template <int dim>
+  InterpolatedTensorProductData<dim>::
+  InterpolatedTensorProductData(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
+    // the (lower) left endpoint of the interval to interpolate
+    // in, and p_unit denotes the point in unit coordinates to do so.
+    double interpolate (const Table<1,double> &data_values,
+                        const TableIndices<1> &ix,
+                        const Point<1>        &xi)
+    {
+      return ((1-xi[0])*data_values[ix[0]]
+                                    +
+                                    xi[0]*data_values[ix[0]+1]);
+    }
+
+    double interpolate (const Table<2,double> &data_values,
+                        const TableIndices<2> &ix,
+                        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]);
+    }
+
+    double interpolate (const Table<3,double> &data_values,
+                        const TableIndices<3> &ix,
+                        const Point<3>        &p_unit)
+    {
+      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]);
+    }
+  }
+
+
+  template <int dim>
+  double
+  InterpolatedTensorProductData<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
+    // 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
+    // intervals, starting at x.size()-2 and going to x.size()-1.
+    TableIndices<dim> ix;
+    for (unsigned int d=0; d<dim; ++d)
+      {
+        // 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(),
+                                   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
+        // us to extend the function by a constant value)
+        //
+        // to make this work, if we got coordinate_values[d].end(), we actually have
+        // to consider the last box which has index size()-2
+        if (ix[d] == coordinate_values[d].size())
+          ix[d] = coordinate_values[d].size()-2;
+        else
+          if (ix[d] > 0)
+            --ix[d];
+      }
+
+    // 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)
+      p_unit[d] =  std::max(std::min((p[d]-coordinate_values[d][ix[d]]) /
+                                     (coordinate_values[d][ix[d]+1]-coordinate_values[d][ix[d]]),
+                                     1.),
+                            0.);
+
+    return interpolate (data_values, ix, p_unit);
+  }
+
+
+
+
 // explicit instantiations
   template class SquareFunction<1>;
   template class SquareFunction<2>;
@@ -2306,6 +2430,9 @@ namespace Functions
   template class Monomial<3>;
   template class Bessel1<2>;
   template class Bessel1<3>;
+  template class InterpolatedTensorProductData<1>;
+  template class InterpolatedTensorProductData<2>;
+  template class InterpolatedTensorProductData<3>;
 }
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/base/functions_10.cc b/tests/base/functions_10.cc
new file mode 100644
index 0000000000..3fa85d08ea
--- /dev/null
+++ b/tests/base/functions_10.cc
@@ -0,0 +1,127 @@
+// ---------------------------------------------------------------------
+// $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 InterpolatedTensorProductData
+
+#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 boost::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 boost::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 boost::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
+  boost::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);
+
+  const Table<dim,double> data = fill(coordinates);
+
+  Functions::InterpolatedTensorProductData<dim> f(coordinates, 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_10.output b/tests/base/functions_10.output
new file mode 100644
index 0000000000..fb71de2867
--- /dev/null
+++ b/tests/base/functions_10.output
@@ -0,0 +1,4 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK