www.mooseframework.org
PiecewiseMultilinear.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 #include "PiecewiseMultilinear.h"
16 #include "GriddedData.h"
17 
18 template <>
21 {
23  params.addParam<FileName>(
24  "data_file",
25  "File holding data for use with PiecewiseMultilinear. Format: any empty line and any line "
26  "beginning with # are ignored, all other lines are assumed to contain relevant information. "
27  "The file must begin with specification of the grid. This is done through lines containing "
28  "the keywords: AXIS X; AXIS Y; AXIS Z; or AXIS T. Immediately following the keyword line "
29  "must be a space-separated line of real numbers which define the grid along the specified "
30  "axis. These data must be monotonically increasing. After all the axes and their grids "
31  "have been specified, there must be a line that is DATA. Following that line, function "
32  "values are given in the correct order (they may be on indivicual lines, or be "
33  "space-separated on a number of lines). When the function is evaluated, f[i,j,k,l] "
34  "corresponds to the i + j*Ni + k*Ni*Nj + l*Ni*Nj*Nk data value. Here i>=0 corresponding to "
35  "the index along the first AXIS, j>=0 corresponding to the index along the second AXIS, etc, "
36  "and Ni = number of grid points along the first AXIS, etc.");
37  params.addClassDescription("PiecewiseMultilinear performs interpolation on 1D, 2D, 3D or 4D "
38  "data. The data_file specifies the axes directions and the function "
39  "values. If a point lies outside the data range, the appropriate end "
40  "value is used.");
41  return params;
42 }
43 
45  : Function(parameters),
46  _gridded_data(libmesh_make_unique<GriddedData>(getParam<FileName>("data_file"))),
47  _dim(_gridded_data->getDim())
48 {
49  _gridded_data->getAxes(_axes);
50  _gridded_data->getGrid(_grid);
51 
52  // GriddedData does not require monotonicity of axes, but we do
53  for (unsigned int i = 0; i < _dim; ++i)
54  for (unsigned int j = 1; j < _grid[i].size(); ++j)
55  if (_grid[i][j - 1] >= _grid[i][j])
56  mooseError("PiecewiseMultilinear needs monotonically-increasing axis data. Axis ",
57  i,
58  " contains non-monotonicity at value ",
59  _grid[i][j]);
60 
61  // GriddedData does not demand that each axis is independent, but we do
62  std::set<int> s(_axes.begin(), _axes.end());
63  if (s.size() != _dim)
64  mooseError("PiecewiseMultilinear needs the AXES to be independent. Check the AXIS lines in "
65  "your data file.");
66 }
67 
69 
70 Real
71 PiecewiseMultilinear::value(Real t, const Point & p)
72 {
73  // convert the inputs to an input to the sample function using _axes
74  std::vector<Real> pt_in_grid(_dim);
75  for (unsigned int i = 0; i < _dim; ++i)
76  {
77  if (_axes[i] < 3)
78  pt_in_grid[i] = p(_axes[i]);
79  else if (_axes[i] == 3) // the time direction
80  pt_in_grid[i] = t;
81  }
82  return sample(pt_in_grid);
83 }
84 
85 Real
86 PiecewiseMultilinear::sample(const std::vector<Real> & pt)
87 {
88  /*
89  * left contains the indices of the point to the 'left', 'down', etc, of pt
90  * right contains the indices of the point to the 'right', 'up', etc, of pt
91  * Hence, left and right define the vertices of the hypercube containing pt
92  */
93  std::vector<unsigned int> left(_dim);
94  std::vector<unsigned int> right(_dim);
95  for (unsigned int i = 0; i < _dim; ++i)
96  {
97  getNeighborIndices(_grid[i], pt[i], left[i], right[i]);
98  }
99 
100  /*
101  * The following just loops through all the vertices of the
102  * hypercube containing pt, evaluating the function at all
103  * those vertices, and weighting the contributions to the
104  * final result depending on the distance of pt from the vertex
105  */
106  Real f = 0;
107  Real weight;
108  std::vector<unsigned int> arg(_dim);
109  for (unsigned int i = 0; i < std::pow(2.0, int(_dim));
110  ++i) // number of points in hypercube = 2^_dim
111  {
112  weight = 1;
113  for (unsigned int j = 0; j < _dim; ++j)
114  if ((i >> j) % 2 ==
115  0) // shift i j-bits to the right and see if the result has a 0 as its right-most bit
116  {
117  arg[j] = left[j];
118  if (left[j] != right[j])
119  weight *= std::abs(pt[j] - _grid[j][right[j]]);
120  else // unusual "end condition" case. weight by 0.5 because we will encounter this twice
121  weight *= 0.5;
122  }
123  else
124  {
125  arg[j] = right[j];
126  if (left[j] != right[j])
127  weight *= std::abs(pt[j] - _grid[j][left[j]]);
128  else // unusual "end condition" case. weight by 0.5 because we will encounter this twice
129  weight *= 0.5;
130  }
131  f += _gridded_data->evaluateFcn(arg) * weight;
132  }
133 
134  /*
135  * finally divide by the volume of the hypercube
136  */
137  weight = 1;
138  for (unsigned int dim = 0; dim < pt.size(); ++dim)
139  if (left[dim] != right[dim])
140  weight *= _grid[dim][right[dim]] - _grid[dim][left[dim]];
141  else // unusual "end condition" case. weight by 1 to cancel the two 0.5 encountered previously
142  weight *= 1;
143 
144  return f / weight;
145 }
146 
147 void
149  Real x,
150  unsigned int & lower_x,
151  unsigned int & upper_x)
152 {
153  int N = in_arr.size();
154  if (x <= in_arr[0])
155  {
156  lower_x = 0;
157  upper_x = 0;
158  }
159  else if (x >= in_arr[N - 1])
160  {
161  lower_x = N - 1;
162  upper_x = N - 1;
163  }
164  else
165  {
166  // returns up which points at the first element in inArr that is not less than x
167  std::vector<double>::iterator up = std::lower_bound(in_arr.begin(), in_arr.end(), x);
168 
169  // std::distance returns std::difference_type, which can be negative in theory, but
170  // in this context will always be >=0. Therefore the explicit cast is just to shut
171  // the compiler up.
172  upper_x = static_cast<unsigned int>(std::distance(in_arr.begin(), up));
173  if (in_arr[upper_x] == x)
174  lower_x = upper_x;
175  else
176  lower_x = upper_x - 1;
177  }
178 }
Container for holding a function defined on a grid of arbitrary dimension.
Definition: GriddedData.h:39
Base class for function objects.
Definition: Function.h:46
PetscInt N
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
static PetscErrorCode Vec x
PiecewiseMultilinear(const InputParameters &parameters)
Create new PiecewiseMultilinear object.
void getNeighborIndices(std::vector< Real > in_arr, Real x, unsigned int &lower_x, unsigned int &upper_x)
Operates on monotonically increasing in_arr.
InputParameters validParams< PiecewiseMultilinear >()
std::unique_ptr< GriddedData > _gridded_data
object to provide function evaluations at points on the grid
std::vector< int > _axes
_axes specifies how to embed the grid into the MOOSE coordinate frame if _axes[i] = 0 then the i_th a...
virtual Real value(Real t, const Point &pt) override
Given t and p, return the interpolated value.
void addClassDescription(const std::string &doc_string)
This method adds a description of the class that will be displayed in the input file syntax dump...
unsigned int _dim
dimension of the grid
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an option parameter and a documentation string to the InputParameters object...
void mooseError(Args &&...args) const
Definition: MooseObject.h:80
InputParameters validParams< Function >()
Definition: Function.C:19
std::vector< std::vector< Real > > _grid
the grid
Real sample(const std::vector< Real > &pt)
This does the core work.