value (const Point<dim> &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<dim> &p,
+ const unsigned int component = 0) const;
+
private:
/**
* The set of coordinate values in each of the coordinate directions.
+
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 <int dim>
InterpolatedTensorProductGridData<dim>::
InterpolatedTensorProductGridData(const std_cxx11::array<std::vector<double>,dim> &coordinate_values,
+ template <int dim>
+ Tensor<1, dim>
+ InterpolatedTensorProductGridData<dim>::gradient(const Point<dim> &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<dim> ix;
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ ix[d] = (std::lower_bound (coordinate_values[d].begin(),
+ coordinate_values[d].end(),
+ p[d])
+ - coordinate_values[d].begin());
+
+ if (ix[d] == coordinate_values[d].size())
+ ix[d] = coordinate_values[d].size()-2;
+ else if (ix[d] > 0)
+ --ix[d];
+ }
+
+ Point<dim> dx;
+ for (unsigned int d=0; d<dim; ++d)
+ dx[d] = coordinate_values[d][ix[d]+1]-coordinate_values[d][ix[d]];
+
+ 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]]) / dx[d],
+ 1.),
+ 0.0);
+
+ return gradient_interpolate(data_values, ix, p_unit, dx);
+ }
+
+
+
template <int dim>
InterpolatedUniformGridData<dim>::
InterpolatedUniformGridData(const std_cxx11::array<std::pair<double,double>,dim> &interval_endpoints,
Functions::InterpolatedTensorProductGridData<dim> 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<dim> p;
AssertThrow (std::fabs (exact_value - f.value(p)) < 1e-12,
ExcInternalError());
- }
+
+ Tensor<1,dim> exact_gradient;
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ exact_gradient[d] = 1.0;
+ for (unsigned int k=0; k<dim; ++k)
+ exact_gradient[d] *= (k == d) ? 1.0 : p[k];
+ }
+
+ AssertThrow ((exact_gradient - f.gradient(p)).norm() < 1e-12,
+ ExcInternalError());
+ }
// now also verify that it computes values outside the box correctly, as
// documented