From 479bf9a83588e2371a8364626025d37ee3463167 Mon Sep 17 00:00:00 2001 From: bangerth Date: Sat, 21 Dec 2013 04:47:05 +0000 Subject: [PATCH] New class InterpolatedUniformGridData. git-svn-id: https://svn.dealii.org/trunk@32072 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 4 +- deal.II/include/deal.II/base/function_lib.h | 97 ++++++++++++- deal.II/source/base/function_lib.cc | 149 ++++++++++++++------ tests/base/functions_10.cc | 8 +- tests/base/functions_11.cc | 140 ++++++++++++++++++ tests/base/functions_11.output | 4 + 6 files changed, 352 insertions(+), 50 deletions(-) create mode 100644 tests/base/functions_11.cc create mode 100644 tests/base/functions_11.output 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.
(Wolfgang Bangerth, 2013/12/20) 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 #include -#include +#include 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,dim> &coordinate_values, - const Table &data_values); + InterpolatedTensorProductGridData (const std_cxx1x::array,dim> &coordinate_values, + const Table &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,dim> coordinate_values; + const std_cxx1x::array,dim> coordinate_values; + + /** + * The data that is to be interpolated. + */ + const Table 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 + class InterpolatedUniformGridData : public Function + { + 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,dim> &interval_endpoints, + const std_cxx1x::array &n_subintervals, + const Table &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 &p, + const unsigned int component = 0) const; + + private: + /** + * The set of interval endpoints in each of the coordinate directions. + */ + const std_cxx1x::array,dim> interval_endpoints; + + /** + * The number of subintervals in each of the coordinate directions. + */ + const std_cxx1x::array 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 - InterpolatedTensorProductGridData:: - InterpolatedTensorProductGridData(const boost::array,dim> &coordinate_values, - const Table &data_values) - : - coordinate_values (coordinate_values), - data_values (data_values) - { - for (unsigned int d=0; d= 2, - ExcMessage ("Coordinate arrays must have at least two coordinate values!")); - for (unsigned int i=0; i &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 + InterpolatedTensorProductGridData:: + InterpolatedTensorProductGridData(const std_cxx1x::array,dim> &coordinate_values, + const Table &data_values) + : + coordinate_values (coordinate_values), + data_values (data_values) + { + for (unsigned int d=0; d= 2, + ExcMessage ("Coordinate arrays must have at least two coordinate values!")); + for (unsigned int i=0; i double InterpolatedTensorProductGridData::value(const Point &p, @@ -2388,6 +2388,70 @@ namespace Functions + template + InterpolatedUniformGridData:: + InterpolatedUniformGridData(const std_cxx1x::array,dim> &interval_endpoints, + const std_cxx1x::array &n_subintervals, + const Table &data_values) + : + interval_endpoints (interval_endpoints), + n_subintervals (n_subintervals), + data_values (data_values) + { + for (unsigned int d=0; 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 + double + InterpolatedUniformGridData::value(const Point &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 ix; + for (unsigned int d=0; 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 p_unit; + for (unsigned int d=0; d; @@ -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,1> &coordinates) +Table<1,double> fill (const std_cxx1x::array,1> &coordinates) { Table<1,double> data(coordinates[0].size()); for (unsigned int i=0; i fill (const boost::array,1> &coordinates) return data; } -Table<2,double> fill (const boost::array,2> &coordinates) +Table<2,double> fill (const std_cxx1x::array,2> &coordinates) { Table<2,double> data(coordinates[0].size(), coordinates[1].size()); @@ -41,7 +41,7 @@ Table<2,double> fill (const boost::array,2> &coordinates) return data; } -Table<3,double> fill (const boost::array,3> &coordinates) +Table<3,double> fill (const std_cxx1x::array,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,dim> coordinates; + std_cxx1x::array,dim> coordinates; for (unsigned int d=0; d + +// 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,1> &coordinates) +{ + Table<1,double> data(coordinates[0].size()); + for (unsigned int i=0; i fill (const std_cxx1x::array,2> &coordinates) +{ + Table<2,double> data(coordinates[0].size(), + coordinates[1].size()); + for (unsigned int i=0; i fill (const std_cxx1x::array,3> &coordinates) +{ + Table<3,double> data(coordinates[0].size(), + coordinates[1].size(), + coordinates[2].size()); + for (unsigned int i=0; i +void check () +{ + // have coordinate arrays that span an interval starting at d+1 + // d+5 nonuniform intervals + std_cxx1x::array,dim> intervals; + std_cxx1x::array n_subintervals; + for (unsigned int d=0; d,dim> coordinates; + for (unsigned int d=0; d data = fill(coordinates); + + Functions::InterpolatedUniformGridData 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 p; + for (unsigned int d=0; d()) - value_at_bottom_left) < 1e-12, + ExcInternalError()); + + Point top_right; + double value_at_top_right = 1; + for (unsigned int d=0; d(); + 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 -- 2.39.5