www.mooseframework.org
PiecewiseMultilinear.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "PiecewiseMultilinear.h"
11 #include "GriddedData.h"
12 
14 
17 {
19  params.addClassDescription(
20  "PiecewiseMultilinear performs linear interpolation on 1D, 2D, 3D or 4D "
21  "data. The data_file specifies the axes directions and the function "
22  "values. If a point lies outside the data range, the appropriate end "
23  "value is used.");
24  params.addParam<Real>(
25  "epsilon", 1e-12, "Finite differencing parameter for gradient and time derivative");
26  return params;
27 }
28 
30  : PiecewiseMultiInterpolation(parameters), _epsilon(getParam<Real>("epsilon"))
31 {
32 }
33 
34 Real
36 {
37  return sampleInternal<false>(pt);
38 }
39 
40 ADReal
42 {
43  return sampleInternal<true>(pt);
44 }
45 
46 template <bool is_ad>
49 {
50  /*
51  * left contains the indices of the point to the 'left', 'down', etc, of pt
52  * right contains the indices of the point to the 'right', 'up', etc, of pt
53  * Hence, left and right define the vertices of the hypercube containing pt
54  */
55  GridIndex left(_dim);
56  GridIndex right(_dim);
57  for (unsigned int i = 0; i < _dim; ++i)
58  getNeighborIndices(_grid[i], MetaPhysicL::raw_value(pt[i]), left[i], right[i]);
59 
60  /*
61  * The following just loops through all the vertices of the
62  * hypercube containing pt, evaluating the function at all
63  * those vertices, and weighting the contributions to the
64  * final result depending on the distance of pt from the vertex
65  */
68  GridIndex arg(_dim);
69  // number of points in hypercube = 2^_dim
70  for (unsigned int i = 0; i < (1u << _dim); ++i)
71  {
72  weight = 1;
73  for (unsigned int j = 0; j < _dim; ++j)
74  // shift i j-bits to the right and see if the result has a 0 as its right-most bit
75  if ((i >> j) % 2 == 0)
76  {
77  arg[j] = left[j];
78  if (left[j] != right[j])
79  weight *= std::abs(pt[j] - _grid[j][right[j]]);
80  else
81  // unusual "end condition" case. weight by 0.5 because we will encounter this twice
82  weight *= 0.5;
83  }
84  else
85  {
86  arg[j] = right[j];
87  if (left[j] != right[j])
88  weight *= std::abs(pt[j] - _grid[j][left[j]]);
89  else
90  // unusual "end condition" case. weight by 0.5 because we will encounter this twice
91  weight *= 0.5;
92  }
93  f += _gridded_data->evaluateFcn(arg) * weight;
94  }
95 
96  /*
97  * finally divide by the volume of the hypercube
98  */
99  weight = 1;
100  for (unsigned int dim = 0; dim < pt.size(); ++dim)
101  if (left[dim] != right[dim])
102  weight *= _grid[dim][right[dim]] - _grid[dim][left[dim]];
103  else
104  // unusual "end condition" case. weight by 1 to cancel the two 0.5 encountered previously
105  weight *= 1;
106 
107  return f / weight;
108 }
109 
111 PiecewiseMultilinear::gradient(Real t, const Point & p) const
112 {
113  RealGradient grad;
114 
115  // sample center point
116  auto s1 = sample(pointInGrid<false>(t, p));
117 
118  // sample epsilon steps in all directions
119  for (const auto dir : _axes)
120  if (dir < 3)
121  {
122  Point pp = p;
123  pp(dir) += _epsilon;
124  grad(dir) = (sample(pointInGrid<false>(t, pp)) - s1) / _epsilon;
125  }
126 
127  return grad;
128 }
129 
130 Real
131 PiecewiseMultilinear::timeDerivative(Real t, const Point & p) const
132 {
133  return (sample(pointInGrid<false>(t + _epsilon, p)) - sample(pointInGrid<false>(t, p))) /
134  _epsilon;
135 }
Uses GriddedData to define data on a grid, and does linear interpolation on that data to provide func...
static InputParameters validParams()
static InputParameters validParams()
Create new PiecewiseMultiInterpolation object.
std::unique_ptr< GriddedData > _gridded_data
object to provide function evaluations at points on the grid
std::vector< std::vector< Real > > _grid
the grid
Uses GriddedData to define data on a grid, and does linear interpolation on that data to provide func...
auto raw_value(const Eigen::Map< T > &in)
Definition: ADReal.h:73
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
Definition: Moose.h:148
Utility class template for a semidynamic vector with a maximum size N and a chosen dynamic size...
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
void getNeighborIndices(std::vector< Real > in_arr, Real x, unsigned int &lower_x, unsigned int &upper_x) const
Operates on monotonically increasing in_arr.
virtual RealGradient gradient(Real t, const Point &p) const override
Function objects can optionally provide a gradient at a point.
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
PiecewiseMultilinear(const InputParameters &parameters)
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...
MooseADWrapper< Real, is_ad > sampleInternal(const MooseADWrapper< GridPoint, is_ad > pt) const
unsigned int _dim
dimension of the grid
DualReal ADReal
Definition: ADRealForward.h:14
virtual Real timeDerivative(Real t, const Point &p) const override
Get the time derivative of the function.
registerMooseObject("MooseApp", PiecewiseMultilinear)
virtual Real sample(const GridPoint &pt) const override
This does the core work.
typename MooseADWrapperStruct< T, is_ad >::type MooseADWrapper
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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...
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...