interpolationTable.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2019-2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "interpolationTable.H"
30#include "openFoamTableReader.H"
31
32// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
33
34template<class Type>
36{
37 // preserve the original (unexpanded) fileName to avoid absolute paths
38 // appearing subsequently in the write() method
39 fileName fName(fileName_);
40
41 fName.expand();
42
43 // Read data from file
44 reader_()(fName, *this);
45
46 if (this->empty())
47 {
49 << "table read from " << fName << " is empty" << nl
50 << exit(FatalError);
51 }
52
53 // Check that the data are okay
54 check();
55}
56
57
58// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
59
60template<class Type>
62:
64 bounding_(bounds::repeatableBounding::WARN),
65 fileName_("fileNameIsUndefined"),
66 reader_(nullptr)
67{}
68
69
70template<class Type>
72(
73 const List<Tuple2<scalar, Type>>& values,
74 const bounds::repeatableBounding bounding,
75 const fileName& fName
76)
77:
78 List<value_type>(values),
79 bounding_(bounding),
80 fileName_(fName),
81 reader_(nullptr)
82{}
83
84
85template<class Type>
87:
89 bounding_(bounds::repeatableBounding::WARN),
90 fileName_(fName),
91 reader_(new openFoamTableReader<Type>())
92{
93 readTable();
94}
95
96
97template<class Type>
99:
100 List<value_type>(),
101 bounding_
102 (
103 bounds::repeatableBoundingNames.getOrDefault
104 (
105 "outOfBounds",
106 dict,
107 bounds::repeatableBounding::WARN,
108 true // Failsafe behaviour
109 )
110 ),
111 fileName_(dict.get<fileName>("file")),
112 reader_(tableReader<Type>::New(dict))
113{
114 readTable();
115}
116
118template<class Type>
121 const interpolationTable& tbl
122)
123:
124 List<value_type>(tbl),
125 bounding_(tbl.bounding_),
126 fileName_(tbl.fileName_),
127 reader_(tbl.reader_.clone())
129
130
131// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
132
133template<class Type>
135{
136 const List<value_type>& list = *this;
137
138 scalar prevValue(0);
139
140 label i = 0;
141 for (const auto& item : list)
143 const scalar& currValue = item.first();
144
145 // Avoid duplicate values (divide-by-zero error)
146 if (i && currValue <= prevValue)
147 {
149 << "out-of-order value: "
150 << currValue << " at index " << i << nl
151 << exit(FatalError);
153 prevValue = currValue;
154 ++i;
156}
157
158
159template<class Type>
161{
162 os.writeEntry("file", fileName_);
163 os.writeEntry("outOfBounds", bounds::repeatableBoundingNames[bounding_]);
164 if (reader_)
166 reader_->write(os);
167 }
168}
169
170
171template<class Type>
173{
174 const List<value_type>& list = *this;
175
176 const label n = list.size();
178 if (n <= 1)
179 {
180 // Not enough entries for a rate of change
181 return Zero;
182 }
183
184 const scalar minLimit = list.first().first();
185 const scalar maxLimit = list.last().first();
186
187 if (lookupValue < minLimit)
188 {
189 switch (bounding_)
190 {
192 {
194 << "value (" << lookupValue << ") less than lower "
195 << "bound (" << minLimit << ")\n"
196 << exit(FatalError);
197 break;
198 }
200 {
202 << "value (" << lookupValue << ") less than lower "
203 << "bound (" << minLimit << ")\n"
204 << " Zero rate of change." << endl;
205
206 // Behaviour as per CLAMP
207 return Zero;
208 break;
209 }
211 {
212 return Zero;
213 break;
214 }
216 {
217 // Adjust lookupValue to >= minLimit
218 scalar span = maxLimit-minLimit;
219 lookupValue = fmod(lookupValue - minLimit, span) + minLimit;
220 break;
221 }
222 }
223 }
224 else if (lookupValue >= maxLimit)
225 {
226 switch (bounding_)
227 {
229 {
231 << "value (" << lookupValue << ") greater than upper "
232 << "bound (" << maxLimit << ")\n"
233 << exit(FatalError);
234 break;
235 }
237 {
239 << "value (" << lookupValue << ") greater than upper "
240 << "bound (" << maxLimit << ")\n"
241 << " Zero rate of change." << endl;
242
243 // Behaviour as per CLAMP
244 return Zero;
245 break;
246 }
248 {
249 return Zero;
250 break;
251 }
253 {
254 // Adjust lookupValue <= maxLimit
255 scalar span = maxLimit-minLimit;
256 lookupValue = fmod(lookupValue - minLimit, span) + minLimit;
257 break;
258 }
259 }
260 }
261
262 label lo = 0;
263 label hi = 0;
264
265 // Look for the correct range
266 for (label i = 0; i < n; ++i)
267 {
268 if (lookupValue >= list[i].first())
269 {
270 lo = hi = i;
271 }
272 else
273 {
274 hi = i;
275 break;
276 }
277 }
278
279 if (lo == hi)
280 {
281 return Zero;
282 }
283 else if (hi == 0)
284 {
285 // This treatment should only occur under these conditions:
286 // -> the 'REPEAT' treatment
287 // -> (0 <= value <= minLimit)
288 // -> minLimit > 0
289 // Use the value at maxLimit as the value for value=0
290 lo = n - 1;
291
292 return
293 (
294 (list[hi].second() - list[lo].second())
295 / (list[hi].first() + minLimit - list[lo].first())
296 );
297 }
298
299
300 // Normal rate of change
301 return
302 (
303 (list[hi].second() - list[lo].second())
304 / (list[hi].first() - list[lo].first())
305 );
306}
307
308
309template<class Type>
311(
312 const List<Tuple2<scalar, Type>>& list,
313 scalar lookupValue,
315)
316{
317 const label n = list.size();
318
319 if (n <= 1)
320 {
321 #ifdef FULLDEBUG
322 if (!n)
323 {
325 << "Cannot interpolate from zero-sized table" << nl
326 << exit(FatalError);
327 }
328 #endif
329
330 return list.first().second();
331 }
332
333 const scalar minLimit = list.first().first();
334 const scalar maxLimit = list.last().first();
335
336 if (lookupValue < minLimit)
337 {
338 switch (bounding)
339 {
341 {
343 << "value (" << lookupValue << ") less than lower "
344 << "bound (" << minLimit << ")\n"
345 << exit(FatalError);
346 break;
347 }
349 {
351 << "value (" << lookupValue << ") less than lower "
352 << "bound (" << minLimit << ")\n"
353 << " Continuing with the first entry" << endl;
354
355 // Behaviour as per CLAMP
356 return list.first().second();
357 break;
358 }
360 {
361 return list.first().second();
362 break;
363 }
365 {
366 // adjust lookupValue to >= minLimit
367 const scalar span = maxLimit-minLimit;
368 lookupValue = fmod(lookupValue - minLimit, span) + minLimit;
369 break;
370 }
371 }
372 }
373 else if (lookupValue >= maxLimit)
374 {
375 switch (bounding)
376 {
378 {
380 << "value (" << lookupValue << ") greater than upper "
381 << "bound (" << maxLimit << ")\n"
382 << exit(FatalError);
383 break;
384 }
386 {
388 << "value (" << lookupValue << ") greater than upper "
389 << "bound (" << maxLimit << ")\n"
390 << " Continuing with the last entry" << endl;
391
392 // Behaviour as per 'CLAMP'
393 return list.last().second();
394 break;
395 }
397 {
398 return list.last().second();
399 break;
400 }
402 {
403 // Adjust lookupValue <= maxLimit
404 const scalar span = maxLimit-minLimit;
405 lookupValue = fmod(lookupValue - minLimit, span) + minLimit;
406 break;
407 }
408 }
409 }
410
411
412 label lo = 0;
413 label hi = 0;
414
415 // Look for the correct range
416 for (label i = 0; i < n; ++i)
417 {
418 if (lookupValue >= list[i].first())
419 {
420 lo = hi = i;
421 }
422 else
423 {
424 hi = i;
425 break;
426 }
427 }
428
429 if (lo == hi)
430 {
431 return list[hi].second();
432 }
433 else if (hi == 0)
434 {
435 // This treatment should only occur under these conditions:
436 // -> the 'REPEAT' treatment
437 // -> (0 <= value <= minLimit)
438 // -> minLimit > 0
439 // Use the value at maxLimit as the value for value=0
440 lo = n - 1;
441
442 return
443 (
444 list[lo].second()
445 + (list[hi].second() - list[lo].second())
446 * (lookupValue / minLimit)
447 );
448 }
449
450
451 // Normal interpolation
452 return
453 (
454 list[lo].second()
455 + (list[hi].second() - list[lo].second())
456 * (lookupValue - list[lo].first())
457 / (list[hi].first() - list[lo].first())
458 );
459}
460
461
462template<class Type>
464(
465 scalar lookupValue
466) const
467{
468 return interpolateValue(*this, lookupValue, bounding_);
469}
470
471
472template<class Type>
475(
476 const UList<scalar>& vals
477) const
478{
479 auto tfld = tmp<Field<Type>>::New(vals.size());
480 auto& fld = tfld.ref();
481
482 forAll(fld, i)
483 {
484 fld[i] = interpolateValue(vals[i]);
485 }
486
487 return tfld;
488}
489
490
491// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
492
493template<class Type>
495(
496 const interpolationTable<Type>& rhs
497)
498{
499 if (this == &rhs)
500 {
501 return;
502 }
503
504 static_cast<List<value_type>&>(*this) = rhs;
505 bounding_ = rhs.bounding_;
506 fileName_ = rhs.fileName_;
507 reader_.reset(rhs.reader_.clone());
508}
509
510
511template<class Type>
514{
515 const List<value_type>& list = *this;
516 const label n = list.size();
517
518 if (n <= 1)
519 {
520 idx = 0;
521
522 #ifdef FULLDEBUG
523 if (!n)
524 {
526 << "Cannot interpolate from zero-sized table" << nl
527 << exit(FatalError);
528 }
529 #endif
530 }
531 else if (idx < 0)
532 {
533 switch (bounding_)
534 {
536 {
538 << "index (" << idx << ") underflow" << nl
539 << exit(FatalError);
540 break;
541 }
543 {
545 << "index (" << idx << ") underflow" << nl
546 << " Continuing with the first entry" << nl;
547
548 // Behaviour as per 'CLAMP'
549 idx = 0;
550 break;
551 }
553 {
554 idx = 0;
555 break;
556 }
558 {
559 while (idx < 0)
560 {
561 idx += n;
562 }
563 break;
564 }
565 }
566 }
567 else if (idx >= n)
568 {
569 switch (bounding_)
570 {
572 {
574 << "index (" << idx << ") overflow" << nl
575 << exit(FatalError);
576 break;
577 }
579 {
581 << "index (" << idx << ") overflow" << nl
582 << " Continuing with the last entry" << nl;
583
584 // Behaviour as per 'CLAMP'
585 idx = n - 1;
586 break;
587 }
589 {
590 idx = n - 1;
591 break;
592 }
594 {
595 while (idx >= n)
596 {
597 idx -= n;
598 }
599 break;
600 }
601 }
602 }
603
604 return list[idx];
605}
606
607
608template<class Type>
609Type Foam::interpolationTable<Type>::operator()(scalar lookupValue) const
610{
611 return interpolateValue(*this, lookupValue, bounding_);
612}
613
614
615// ************************************************************************* //
label n
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
Ostream & writeEntry(const keyType &key, const T &value)
Write a keyword/value entry.
Definition: Ostream.H:239
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: UList.H:94
T & first()
Return the first element of the list.
Definition: UListI.H:202
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
T & last()
Return the last element of the list.
Definition: UListI.H:216
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
Ostream & operator()() const
Output stream (master only).
Definition: ensightCaseI.H:74
A class for handling file names.
Definition: fileName.H:76
virtual bool write()
Write the output fields.
An interpolation/look-up table of scalar vs <Type> values. The reference scalar values must be monoto...
void check() const
Check that list is monotonically increasing.
interpolationTable()
Default construct.
tmp< Field< Type > > interpolateValues(const UList< scalar > &vals) const
Return multiple interpolated values.
static Type interpolateValue(const List< Tuple2< scalar, Type > > &list, scalar lookupValue, bounds::repeatableBounding=bounds::repeatableBounding::CLAMP)
Return an interpolated value in List.
Type rateOfChange(scalar lookupValue) const
const Tuple2< scalar, Type > & operator[](label idx) const
Return an element of constant Tuple2<scalar, Type>
Reads an interpolation table from a file - OpenFOAM-format.
Base class to read table data for the interpolationTable.
Definition: tableReader.H:61
A class for managing temporary objects.
Definition: tmp.H:65
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
OBJstream os(runTime.globalPath()/outputName)
#define WarningInFunction
Report a warning using Foam::Warning.
repeatableBounding
Enumeration for handling out-of-bound values that are repeatable.
Definition: tableBounds.H:62
@ WARN
Issue warning and clamp value (this is a good default)
@ REPEAT
Treat as a repeating list.
@ ERROR
Exit with a FatalError.
@ CLAMP
Clamp value to the start/end value.
const Foam::Enum< repeatableBounding > repeatableBoundingNames
Strings corresponding to the repeatableBounding.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh > > &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
error FatalError
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333