interpolation2DTable.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "IFstream.H"
27 #include "openFoamTableReader.H"
28 
29 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
30 
31 template<class Type>
33 {
34  fileName fName(fileName_);
35  fName.expand();
36 
37  // Read data from file
38  reader_()(fName, *this);
39 
40  if (this->empty())
41  {
43  (
44  "Foam::interpolation2DTable<Type>::readTable()"
45  ) << "table read from " << fName << " is empty" << nl
46  << exit(FatalError);
47  }
48 
49  // Check that the data are in ascending order
50  checkOrder();
51 }
52 
53 
54 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
55 
56 template<class Type>
58 :
59  List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >(),
60  boundsHandling_(interpolation2DTable::WARN),
61  fileName_("fileNameIsUndefined"),
62  reader_(NULL)
63 {}
64 
65 
66 template<class Type>
68 (
69  const List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >& values,
70  const boundsHandling bounds,
71  const fileName& fName
72 )
73 :
75  boundsHandling_(bounds),
76  fileName_(fName),
77  reader_(NULL)
78 {}
79 
80 
81 template<class Type>
83 :
84  List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >(),
85  boundsHandling_(interpolation2DTable::WARN),
86  fileName_(fName),
87  reader_(new openFoamTableReader<Type>(dictionary()))
88 {
89  readTable();
90 }
91 
92 
93 template<class Type>
95 :
96  List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >(),
97  boundsHandling_(wordToBoundsHandling(dict.lookup("outOfBounds"))),
98  fileName_(dict.lookup("fileName")),
99  reader_(tableReader<Type>::New(dict))
100 {
101  readTable();
102 }
103 
104 
105 template<class Type>
107 (
108  const interpolation2DTable& interpTable
109 )
110 :
112  boundsHandling_(interpTable.boundsHandling_),
113  fileName_(interpTable.fileName_),
114  reader_(interpTable.reader_) // note: steals reader. Used in write().
115 {}
116 
117 
118 
119 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
120 
121 template<class Type>
123 (
124  const List<Tuple2<scalar, Type> >& data,
125  const scalar lookupValue
126 ) const
127 {
128  label n = data.size();
129 
130  scalar minLimit = data.first().first();
131  scalar maxLimit = data.last().first();
132 
133  if (lookupValue < minLimit)
134  {
135  switch (boundsHandling_)
136  {
138  {
140  (
141  "Foam::interpolation2DTable<Type>::interpolateValue"
142  "("
143  "List<Tuple2<scalar, Type> >&, "
144  "const scalar"
145  ")"
146  ) << "value (" << lookupValue << ") less than lower "
147  << "bound (" << minLimit << ")" << nl
148  << exit(FatalError);
149  break;
150  }
152  {
153  WarningIn
154  (
155  "Foam::interpolation2DTable<Type>::interpolateValue"
156  "("
157  "List<Tuple2<scalar, Type> >&, "
158  "const scalar"
159  ")"
160  ) << "value (" << lookupValue << ") less than lower "
161  << "bound (" << minLimit << ")" << nl
162  << " Continuing with the first entry"
163  << endl;
164  // fall-through to 'CLAMP'
165  }
167  {
168  return data.first().second();
169  break;
170  }
171  }
172  }
173  else if (lookupValue >= maxLimit)
174  {
175  switch (boundsHandling_)
176  {
178  {
180  (
181  "Foam::interpolation2DTable<Type>::interpolateValue"
182  "("
183  "List<Tuple2<scalar, Type> >&, "
184  "const scalar"
185  ")"
186  ) << "value (" << lookupValue << ") greater than upper "
187  << "bound (" << maxLimit << ")" << nl
188  << exit(FatalError);
189  break;
190  }
192  {
193  WarningIn
194  (
195  "Foam::interpolation2DTable<Type>::interpolateValue"
196  "("
197  "List<Tuple2<scalar, Type> >&, "
198  "const scalar"
199  ")"
200  ) << "value (" << lookupValue << ") greater than upper "
201  << "bound (" << maxLimit << ")" << nl
202  << " Continuing with the last entry"
203  << endl;
204  // fall-through to 'CLAMP'
205  }
207  {
208  return data.last().second();
209  break;
210  }
211  }
212  }
213 
214  // look for the correct range in X
215  label lo = 0;
216  label hi = 0;
217 
218  for (label i = 0; i < n; ++i)
219  {
220  if (lookupValue >= data[i].first())
221  {
222  lo = hi = i;
223  }
224  else
225  {
226  hi = i;
227  break;
228  }
229  }
230 
231  if (lo == hi)
232  {
233  return data[lo].second();
234  }
235  else
236  {
237  Type m =
238  (data[hi].second() - data[lo].second())
239  /(data[hi].first() - data[lo].first());
240 
241  // normal interpolation
242  return data[lo].second() + m*(lookupValue - data[lo].first());
243  }
244 }
245 
246 
247 template<class Type>
248 template<class BinaryOp>
250 (
251  const BinaryOp& bop,
252  const scalar valueX,
253  const bool reverse
254 ) const
255 {
256  const table& t = *this;
257 
258  label limitI = 0;
259  if (reverse)
260  {
261  limitI = t.size() - 1;
262  }
263 
264  if (bop(valueX, t[limitI].first()))
265  {
266  switch (boundsHandling_)
267  {
269  {
271  (
272  "Foam::label Foam::interpolation2DTable<Type>::Xi"
273  "("
274  "const BinaryOp&, "
275  "const scalar, "
276  "const bool"
277  ") const"
278  ) << "value (" << valueX << ") out of bounds"
279  << exit(FatalError);
280  break;
281  }
283  {
284  WarningIn
285  (
286  "Foam::label Foam::interpolation2DTable<Type>::Xi"
287  "("
288  "const BinaryOp&, "
289  "const scalar, "
290  "const bool"
291  ) << "value (" << valueX << ") out of bounds"
292  << endl;
293  // fall-through to 'CLAMP'
294  }
296  {
297  return limitI;
298  }
299  default:
300  {
302  (
303  "Foam::label Foam::interpolation2DTable<Type>::Xi"
304  "("
305  "const BinaryOp&, "
306  "const scalar, "
307  "const bool"
308  ") const"
309  )
310  << "Un-handled enumeration " << boundsHandling_
311  << abort(FatalError);
312  }
313  }
314  }
315 
316  label i = 0;
317  if (reverse)
318  {
319  label nX = t.size();
320  i = 0;
321  while ((i < nX) && (valueX > t[i].first()))
322  {
323  i++;
324  }
325  }
326  else
327  {
328  i = t.size() - 1;
329  while ((i > 0) && (valueX < t[i].first()))
330  {
331  i--;
332  }
333  }
334 
335  return i;
336 }
337 
338 
339 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
340 
341 template<class Type>
343 (
344  const scalar valueX,
345  const scalar valueY
346 ) const
347 {
348  // Considers all of the list in Y being equal
349  label nX = this->size();
350 
351  const table& t = *this;
352 
353  if (nX == 0)
354  {
355  WarningIn
356  (
357  "Type Foam::interpolation2DTable<Type>::operator()"
358  "("
359  "const scalar, "
360  "const scalar"
361  ") const"
362  )
363  << "cannot interpolate a zero-sized table - returning zero" << endl;
364 
365  return pTraits<Type>::zero;
366  }
367  else if (nX == 1)
368  {
369  // only 1 column (in X) - interpolate to find Y value
370  return interpolateValue(t.first().second(), valueY);
371  }
372  else
373  {
374  // have 2-D data, interpolate
375 
376  // find low and high indices in the X range that bound valueX
377  label x0i = Xi(lessOp<scalar>(), valueX, false);
378  label x1i = Xi(greaterOp<scalar>(), valueX, true);
379 
380  if (x0i == x1i)
381  {
382  return interpolateValue(t[x0i].second(), valueY);
383  }
384  else
385  {
386  Type y0(interpolateValue(t[x0i].second(), valueY));
387  Type y1(interpolateValue(t[x1i].second(), valueY));
388 
389  // gradient in X
390  scalar x0 = t[x0i].first();
391  scalar x1 = t[x1i].first();
392  Type mX = (y1 - y0)/(x1 - x0);
393 
394  // interpolate
395  return y0 + mX*(valueX - x0);
396  }
397  }
398 }
399 
400 
401 template<class Type>
403 (
404  const boundsHandling& bound
405 ) const
406 {
407  word enumName("warn");
408 
409  switch (bound)
410  {
412  {
413  enumName = "error";
414  break;
415  }
417  {
418  enumName = "warn";
419  break;
420  }
422  {
423  enumName = "clamp";
424  break;
425  }
426  }
427 
428  return enumName;
429 }
430 
431 
432 template<class Type>
435 (
436  const word& bound
437 ) const
438 {
439  if (bound == "error")
440  {
442  }
443  else if (bound == "warn")
444  {
446  }
447  else if (bound == "clamp")
448  {
450  }
451  else
452  {
453  WarningIn
454  (
455  "Foam::interpolation2DTable<Type>::wordToBoundsHandling"
456  "("
457  " const word&"
458  ")"
459  ) << "bad outOfBounds specifier " << bound << " using 'warn'" << endl;
460 
462  }
463 }
464 
465 
466 template<class Type>
469 (
470  const boundsHandling& bound
471 )
472 {
473  boundsHandling prev = boundsHandling_;
474  boundsHandling_ = bound;
475  return prev;
476 }
477 
478 
479 template<class Type>
481 {
482  label n = this->size();
483  const table& t = *this;
484 
485  scalar prevValue = t[0].first();
486 
487  for (label i=1; i<n; ++i)
488  {
489  const scalar currValue = t[i].first();
490 
491  // avoid duplicate values (divide-by-zero error)
492  if (currValue <= prevValue)
493  {
495  (
496  "Foam::interpolation2DTable<Type>::checkOrder() const"
497  ) << "out-of-order value: "
498  << currValue << " at index " << i << nl
499  << exit(FatalError);
500  }
501  prevValue = currValue;
502  }
503 }
504 
505 
506 template<class Type>
508 {
509  os.writeKeyword("fileName")
510  << fileName_ << token::END_STATEMENT << nl;
511  os.writeKeyword("outOfBounds")
512  << boundsHandlingToWord(boundsHandling_) << token::END_STATEMENT << nl;
513 
514  *this >> os;
515 }
516 
517 
518 // ************************************************************************* //
T * last()
Return the last entry.
Definition: UILList.H:118
A 2-tuple for storing two objects of different types.
Definition: Tuple2.H:47
interpolation2DTable()
Construct null.
T * first()
Return the first entry.
Definition: UILList.H:106
string & expand(const bool allowEmpty=false)
Expand initial tildes and all occurences of environment variables.
Definition: string.C:98
dimensionedScalar y1(const dimensionedScalar &ds)
word boundsHandlingToWord(const boundsHandling &bound) const
Return the out-of-bounds handling as a word.
boundsHandling
Enumeration for handling out-of-bound values.
volScalarField & bound(volScalarField &, const dimensionedScalar &lowerBound)
Bound the given scalar field if it has gone unbounded.
Definition: bound.C:33
A class for handling words, derived from string.
Definition: word.H:59
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Database for solution data, solver performance and other reduced data.
Definition: data.H:52
Reads an interpolation table from a file - OpenFOAM-format.
friend Ostream & operator(Ostream &, const UList< T > &)
dimensionedScalar y0(const dimensionedScalar &ds)
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
label size() const
Return number of elements in list.
Definition: DLListBaseI.H:77
T & first()
Return the first element of the list.
Definition: UListI.H:117
autoPtr< BasicCompressibleTurbulenceModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleTurbulenceModel::transportModel &transport, const word &propertiesName)
void checkOrder() const
Check that list is monotonically increasing.
label n
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
boundsHandling outOfBounds(const boundsHandling &bound)
Set the out-of-bounds handling from enum, return previous setting.
dictionary dict
static const char nl
Definition: Ostream.H:260
2D table interpolation. The data must be in ascending order in both dimensions x and y...
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
stressControl lookup("compactNormalStress") >> compactNormalStress
#define WarningIn(functionName)
Report a warning using Foam::Warning.
Base class to read table data for the interpolationTable.
Definition: tableReader.H:57
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Ostream & writeKeyword(const keyType &)
Write the keyword followed by an appropriate indentation.
Definition: Ostream.C:59
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
Traits class for primitives.
Definition: pTraits.H:50
error FatalError
boundsHandling wordToBoundsHandling(const word &bound) const
Return the out-of-bounds handling as an enumeration.
A class for handling file names.
Definition: fileName.H:69
label size() const
Return the number of elements in the UList.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
void reverse(UList< T > &, const label n)
Definition: UListI.H:322
void write(Ostream &os) const
Write.