From: danshapero Date: Sun, 2 Aug 2015 23:10:33 +0000 (-0700) Subject: gradient of InterpolatedTensorProductGridData X-Git-Tag: v8.4.0-rc2~405^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b7cb965ac1bd850180de2f89df7dca774c87ec27;p=dealii.git gradient of InterpolatedTensorProductGridData --- diff --git a/doc/news/changes.h b/doc/news/changes.h index 6f5c18bc5e..29a4d4626e 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -150,6 +150,11 @@ inconvenience this causes.

General

+
  • New: implemented the gradient method for + InterpolatedTensorProductGridData +
    + (Daniel Shapero, 2015/08/12) +
    1. Changed: All doxygen-generated pages now contain a link to the diff --git a/include/deal.II/base/function_lib.h b/include/deal.II/base/function_lib.h index d05b71e1b8..67d34f3b40 100644 --- a/include/deal.II/base/function_lib.h +++ b/include/deal.II/base/function_lib.h @@ -1154,6 +1154,22 @@ namespace Functions value (const Point &p, const unsigned int component = 0) const; + /** + * Compute the gradient of the function defined by bilinear interpolation + * of the given data set. + * + * @param p The point at which the function gradient is to be evaluated. + * @param component The vector component. Since this function is scalar, + * only zero is a valid argument here. + * @return The value of the gradient of the interpolated function at this + * point. If the point lies outside the set of coordinates, the function + * is extended by a constant and so its gradient is extended by 0. + */ + virtual + Tensor<1, dim> + gradient (const Point &p, + const unsigned int component = 0) const; + private: /** * The set of coordinate values in each of the coordinate directions. diff --git a/source/base/function_lib.cc b/source/base/function_lib.cc index e86263e0f6..ab271f81b6 100644 --- a/source/base/function_lib.cc +++ b/source/base/function_lib.cc @@ -2340,9 +2340,87 @@ namespace Functions + p_unit[0]*data_values[ix[0]+1][ix[1]+1][ix[2]+1])*p_unit[1]) * p_unit[2]); } - } + // Interpolate the gradient of a data value from a table where ix + // denotes the lower left endpoint of the interval to interpolate + // in, p_unit denotes the point in unit coordinates, and dx + // denotes the width of the interval in each dimension. + Tensor<1,1> gradient_interpolate (const Table<1,double> &data_values, + const TableIndices<1> &ix, + const Point<1> &p_unit, + const Point<1> &dx) + { + (void)p_unit; + Tensor<1,1> grad; + grad[0] = (data_values[ix[0]+1] - data_values[ix[0]]) / dx[0]; + return grad; + } + + + Tensor<1,2> gradient_interpolate (const Table<2,double> &data_values, + const TableIndices<2> &ix, + const Point<2> &p_unit, + const Point<2> &dx) + { + Tensor<1,2> grad; + double + u00 = data_values[ix[0]][ix[1]], + u01 = data_values[ix[0]+1][ix[1]], + u10 = data_values[ix[0]][ix[1]+1], + u11 = data_values[ix[0]+1][ix[1]+1]; + + grad[0] = ((1-p_unit[1])*(u01-u00) + p_unit[1]*(u11-u10))/dx[0]; + grad[1] = ((1-p_unit[0])*(u10-u00) + p_unit[0]*(u11-u01))/dx[1]; + return grad; + } + + + Tensor<1,3> gradient_interpolate (const Table<3,double> &data_values, + const TableIndices<3> &ix, + const Point<3> &p_unit, + const Point<3> &dx) + { + Tensor<1,3> grad; + double + u000 = data_values[ix[0]][ix[1]][ix[2]], + u001 = data_values[ix[0]+1][ix[1]][ix[2]], + u010 = data_values[ix[0]][ix[1]+1][ix[2]], + u100 = data_values[ix[0]][ix[1]][ix[2]+1], + u011 = data_values[ix[0]+1][ix[1]+1][ix[2]], + u101 = data_values[ix[0]+1][ix[1]][ix[2]+1], + u110 = data_values[ix[0]][ix[1]+1][ix[2]+1], + u111 = data_values[ix[0]+1][ix[1]+1][ix[2]+1]; + + grad[0] = ((1-p_unit[2]) + * + ((1-p_unit[1])*(u001-u000) + p_unit[1]*(u011-u010)) + + + p_unit[2] + * + ((1-p_unit[1])*(u101-u100) + p_unit[1]*(u111-u110)) + )/dx[0]; + grad[1] = ((1-p_unit[2]) + * + ((1-p_unit[0])*(u010-u000) + p_unit[0]*(u011-u001)) + + + p_unit[2] + * + ((1-p_unit[0])*(u110-u100) + p_unit[0]*(u111-u101)) + )/dx[1]; + grad[2] = ((1-p_unit[1]) + * + ((1-p_unit[0])*(u100-u000) + p_unit[0]*(u101-u001)) + + + p_unit[1] + * + ((1-p_unit[0])*(u110-u010) + p_unit[0]*(u111-u011)) + )/dx[2]; + + return grad; + } + } + template InterpolatedTensorProductGridData:: InterpolatedTensorProductGridData(const std_cxx11::array,dim> &coordinate_values, @@ -2415,6 +2493,45 @@ namespace Functions + template + Tensor<1, dim> + InterpolatedTensorProductGridData::gradient(const Point &p, + const unsigned int component) const + { + (void)component; + Assert (component == 0, + ExcMessage ("This is a scalar function object, the component can only be zero.")); + + // find out where this data point lies + TableIndices ix; + for (unsigned int d=0; d 0) + --ix[d]; + } + + Point dx; + for (unsigned int d=0; d p_unit; + for (unsigned int d=0; d InterpolatedUniformGridData:: InterpolatedUniformGridData(const std_cxx11::array,dim> &interval_endpoints, diff --git a/tests/base/functions_10.cc b/tests/base/functions_10.cc index 9ffb6d3be6..7cb1ac0e47 100644 --- a/tests/base/functions_10.cc +++ b/tests/base/functions_10.cc @@ -70,7 +70,7 @@ void check () Functions::InterpolatedTensorProductGridData f(coordinates, data); // now choose a number of randomly chosen points inside the box and - // verify that the functions returned are correct + // verify that the function values and gradients returned are correct for (unsigned int i=0; i<10; ++i) { Point p; @@ -85,7 +85,18 @@ void check () AssertThrow (std::fabs (exact_value - f.value(p)) < 1e-12, ExcInternalError()); - } + + Tensor<1,dim> exact_gradient; + for (unsigned int d=0; d