<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
#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
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
}
+
+ 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>;
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
--- /dev/null
+// ---------------------------------------------------------------------
+// $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>();
+}
+
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
projects / dealii-svn.git / commitdiff