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