From: bangerth 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.

Specific improvements

    +
  1. 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. +
    + (Wolfgang Bangerth, 2013/12/20) +
  2. +
  3. 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 #include #include +#include + +#include 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 + class InterpolatedTensorProductData : public Function + { + 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,dim> &coordinate_values, + 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 coordinate values in each of the coordinate directions. + */ + const boost::array,dim> coordinate_values; + + /** + * The data that is to be interpolated. + */ + const Table 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 + InterpolatedTensorProductData:: + InterpolatedTensorProductData(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 &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 + double + InterpolatedTensorProductData::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 + // 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 ix; + for (unsigned int d=0; 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 p_unit; + for (unsigned int d=0; d; 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 + +// 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> data(coordinates[0].size()); + for (unsigned int i=0; i fill (const boost::array,2> &coordinates) +{ + Table<2,double> data(coordinates[0].size(), + coordinates[1].size()); + for (unsigned int i=0; i fill (const boost::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 + boost::array,dim> coordinates; + for (unsigned int d=0; d data = fill(coordinates); + + Functions::InterpolatedTensorProductData 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 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_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