www.mooseframework.org
LinearInterpolation.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 "LinearInterpolation.h"
16 
17 #include <cassert>
18 #include <fstream>
19 #include <stdexcept>
20 
22 
23 LinearInterpolation::LinearInterpolation(const std::vector<Real> & x, const std::vector<Real> & y)
24  : _x(x), _y(y)
25 {
26  errorCheck();
27 }
28 
29 void
31 {
32  if (_x.size() != _y.size())
33  throw std::domain_error("Vectors are not the same length");
34 
35  for (unsigned int i = 0; i + 1 < _x.size(); ++i)
36  if (_x[i] >= _x[i + 1])
37  {
38  std::ostringstream oss;
39  oss << "x-values are not strictly increasing: x[" << i << "]: " << _x[i] << " x[" << i + 1
40  << "]: " << _x[i + 1];
41  throw std::domain_error(oss.str());
42  }
43 }
44 
45 Real
47 {
48  // sanity check (empty LinearInterpolations get constructed in many places
49  // so we cannot put this into the errorCheck)
50  assert(_x.size() > 0);
51 
52  // endpoint cases
53  if (x <= _x[0])
54  return _y[0];
55  if (x >= _x.back())
56  return _y.back();
57 
58  for (unsigned int i = 0; i + 1 < _x.size(); ++i)
59  if (x >= _x[i] && x < _x[i + 1])
60  return _y[i] + (_y[i + 1] - _y[i]) * (x - _x[i]) / (_x[i + 1] - _x[i]);
61 
62  throw std::out_of_range("Unreachable");
63  return 0;
64 }
65 
66 Real
68 {
69  // endpoint cases
70  if (x < _x[0])
71  return 0.0;
72  if (x >= _x[_x.size() - 1])
73  return 0.0;
74 
75  for (unsigned int i = 0; i + 1 < _x.size(); ++i)
76  if (x >= _x[i] && x < _x[i + 1])
77  return (_y[i + 1] - _y[i]) / (_x[i + 1] - _x[i]);
78 
79  throw std::out_of_range("Unreachable");
80  return 0;
81 }
82 
83 Real
85 {
86  Real answer = 0;
87  for (unsigned int i = 1; i < _x.size(); ++i)
88  answer += 0.5 * (_y[i] + _y[i - 1]) * (_x[i] - _x[i - 1]);
89 
90  return answer;
91 }
92 
93 Real
95 {
96  return _x[i];
97 }
98 
99 Real
101 {
102  return _y[i];
103 }
104 
105 void
107  std::string x_label,
108  std::string y_label,
109  Real xmin,
110  Real xmax,
111  Real ymin,
112  Real ymax)
113 {
114  std::stringstream filename, filename_pts;
115  const unsigned char fill_character = '0';
116  const unsigned int field_width = 4;
117 
118  filename.fill(fill_character);
119  filename << base_name;
120  filename.width(field_width);
121  filename << _file_number << ".plt";
122 
123  filename_pts.fill(fill_character);
124  filename_pts << base_name << "_pts";
125  filename_pts.width(field_width);
126  filename_pts << _file_number << ".dat";
127 
128  /* First dump the GNUPLOT file with the Piecewise Linear Equations */
129  std::ofstream out(filename.str().c_str());
130  out.precision(15);
131  out.fill(fill_character);
132 
133  out << "set terminal postscript color enhanced\n"
134  << "set output \"" << base_name;
135  out.width(field_width);
136  out << _file_number << ".eps\"\n"
137  << "set xlabel \"" << x_label << "\"\n"
138  << "set ylabel \"" << y_label << "\"\n";
139  if (xmin != 0 && xmax != 0)
140  out << "set xrange [" << xmin << ":" << xmax << "]\n";
141  if (ymin != 0 && ymax != 0)
142  out << "set yrange [" << ymin << ":" << ymax << "]\n";
143  out << "set key left top\n"
144  << "f(x)=";
145 
146  for (unsigned int i = 1; i < _x.size(); ++i)
147  {
148  Real m = (_y[i] - _y[i - 1]) / (_x[i] - _x[i - 1]);
149  Real b = (_y[i] - m * _x[i]);
150 
151  out << _x[i - 1] << "<=x && x<" << _x[i] << " ? " << m << "*x+(" << b << ") : ";
152  }
153  out << " 1/0\n";
154 
155  out << "\nplot f(x) with lines, '" << filename_pts.str() << "' using 1:2 title \"Points\"\n";
156  out.close();
157 
158  assert(_x.size() == _y.size());
159 
160  out.open(filename_pts.str().c_str());
161  /* Next dump the data points into a separate file */
162  for (unsigned int i = 0; i < _x.size(); ++i)
163  out << _x[i] << " " << _y[i] << "\n";
164  out << std::endl;
165 
166  ++_file_number;
167  out.close();
168 }
169 
170 unsigned int
172 {
173  return _x.size();
174 }
std::vector< Real > _x
unsigned int getSampleSize()
This function returns the size of the array holding the points, i.e.
Real integrate()
This function returns the integral of the function.
std::vector< Real > _y
Real sample(Real x) const
This function will take an independent variable input and will return the dependent variable based on...
static PetscErrorCode Vec x
Real sampleDerivative(Real x) const
This function will take an independent variable input and will return the derivative of the dependent...
PetscInt m
Real range(int i) const
void dumpSampleFile(std::string base_name, std::string x_label="X", std::string y_label="Y", Real xmin=0, Real xmax=0, Real ymin=0, Real ymax=0)
This function will dump GNUPLOT input files that can be run to show the data points and function fits...
Real domain(int i) const