www.mooseframework.org
LinearInterpolation.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 "LinearInterpolation.h"
11 
12 #include "DualRealOps.h"
13 #include "ChainedReal.h"
14 
15 #include <cassert>
16 #include <fstream>
17 #include <stdexcept>
18 
19 LinearInterpolation::LinearInterpolation(const std::vector<Real> & x,
20  const std::vector<Real> & y,
21  const bool extrap)
22  : _x(x), _y(y), _extrap(extrap)
23 {
24  errorCheck();
25 }
26 
27 void
29 {
30  if (_x.size() != _y.size())
31  throw std::domain_error("Vectors are not the same length");
32 
33  for (unsigned int i = 0; i + 1 < _x.size(); ++i)
34  if (_x[i] >= _x[i + 1])
35  {
36  std::ostringstream oss;
37  oss << "x-values are not strictly increasing: x[" << i << "]: " << _x[i] << " x[" << i + 1
38  << "]: " << _x[i + 1];
39  throw std::domain_error(oss.str());
40  }
41 }
42 
43 template <typename T>
44 T
45 LinearInterpolation::sample(const T & x) const
46 {
47  // this is a hard error
48  if (_x.empty())
49  mooseError("Trying to evaluate an empty LinearInterpolation");
50 
51  // special case for single function nodes
52  if (_x.size() == 1)
53  return _y[0];
54 
55  // endpoint cases
56  if (_extrap)
57  {
58  if (x < _x[0])
59  return _y[0] + (x - _x[0]) / (_x[1] - _x[0]) * (_y[1] - _y[0]);
60  if (x >= _x.back())
61  return _y.back() +
62  (x - _x.back()) / (_x[_x.size() - 2] - _x.back()) * (_y[_y.size() - 2] - _y.back());
63  }
64  else
65  {
66  if (x < _x[0])
67  return _y[0];
68  if (x >= _x.back())
69  return _y.back();
70  }
71 
72  auto upper = std::upper_bound(_x.begin(), _x.end(), x);
73  int i = std::distance(_x.begin(), upper) - 1;
74  return _y[i] + (_y[i + 1] - _y[i]) * (x - _x[i]) / (_x[i + 1] - _x[i]);
75 
76  // If this point is reached, x must be a NaN.
77  mooseException("Sample point in LinearInterpolation is a NaN.");
78 }
79 
80 template Real LinearInterpolation::sample<Real>(const Real &) const;
81 template ADReal LinearInterpolation::sample<ADReal>(const ADReal &) const;
82 template ChainedReal LinearInterpolation::sample<ChainedReal>(const ChainedReal &) const;
83 
84 template <typename T>
85 T
87 {
88  // endpoint cases
89  if (_extrap)
90  {
91  if (x <= _x[0])
92  return (_y[1] - _y[0]) / (_x[1] - _x[0]);
93  if (x >= _x.back())
94  return (_y[_y.size() - 2] - _y.back()) / (_x[_x.size() - 2] - _x.back());
95  }
96  else
97  {
98  if (x < _x[0])
99  return 0.0;
100  if (x >= _x.back())
101  return 0.0;
102  }
103 
104  auto upper = std::upper_bound(_x.begin(), _x.end(), x);
105  int i = std::distance(_x.begin(), upper) - 1;
106  return (_y[i + 1] - _y[i]) / (_x[i + 1] - _x[i]);
107 
108  // If this point is reached, x must be a NaN.
109  mooseException("Sample point in LinearInterpolation is a NaN.");
110 }
111 
112 template Real LinearInterpolation::sampleDerivative<Real>(const Real &) const;
113 template ADReal LinearInterpolation::sampleDerivative<ADReal>(const ADReal &) const;
114 template ChainedReal LinearInterpolation::sampleDerivative<ChainedReal>(const ChainedReal &) const;
115 
116 Real
118 {
119  Real answer = 0;
120  for (unsigned int i = 1; i < _x.size(); ++i)
121  answer += 0.5 * (_y[i] + _y[i - 1]) * (_x[i] - _x[i - 1]);
122 
123  return answer;
124 }
125 
126 Real
128 {
129  return _x[i];
130 }
131 
132 Real
134 {
135  return _y[i];
136 }
137 
138 unsigned int
140 {
141  return _x.size();
142 }
std::vector< Real > _x
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:284
DualNumber< Real, Real > ChainedReal
Definition: ChainedReal.h:30
Real integrate()
This function returns the integral of the function.
T sample(const T &x) const
This function will take an independent variable input and will return the dependent variable based on...
std::vector< Real > _y
unsigned int getSampleSize() const
This function returns the size of the array holding the points, i.e.
Real range(int i) const
Real domain(int i) const
T sampleDerivative(const T &x) const
This function will take an independent variable input and will return the derivative of the dependent...
DualReal ADReal
Definition: ADRealForward.h:14
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real