Commit 87330972 authored by Mark OLESEN's avatar Mark OLESEN Committed by Andrew Heather
Browse files

ENH: interpolationTable improvements

- reduce code duplication, support returning multiple interpolations
  as a Field
parent 784d3ad5
......@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2016 OpenCFD Ltd.
Copyright (C) 2016-2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -47,7 +47,7 @@ void Foam::interpolation2DTable<Type>::readTable()
}
// Check that the data are in ascending order
checkOrder();
check();
}
......@@ -56,7 +56,7 @@ void Foam::interpolation2DTable<Type>::readTable()
template<class Type>
Foam::interpolation2DTable<Type>::interpolation2DTable()
:
List<Tuple2<scalar, List<Tuple2<scalar, Type>>>>(),
List<value_type>(),
bounding_(bounds::normalBounding::WARN),
fileName_("fileNameIsUndefined"),
reader_(nullptr)
......@@ -71,7 +71,7 @@ Foam::interpolation2DTable<Type>::interpolation2DTable
const fileName& fName
)
:
List<Tuple2<scalar, List<Tuple2<scalar, Type>>>>(values),
List<value_type>(values),
bounding_(bounding),
fileName_(fName),
reader_(nullptr)
......@@ -81,7 +81,7 @@ Foam::interpolation2DTable<Type>::interpolation2DTable
template<class Type>
Foam::interpolation2DTable<Type>::interpolation2DTable(const fileName& fName)
:
List<Tuple2<scalar, List<Tuple2<scalar, Type>>>>(),
List<value_type>(),
bounding_(bounds::normalBounding::WARN),
fileName_(fName),
reader_(new openFoamTableReader<Type>(dictionary()))
......@@ -93,7 +93,7 @@ Foam::interpolation2DTable<Type>::interpolation2DTable(const fileName& fName)
template<class Type>
Foam::interpolation2DTable<Type>::interpolation2DTable(const dictionary& dict)
:
List<Tuple2<scalar, List<Tuple2<scalar, Type>>>>(),
List<value_type>(),
bounding_
(
bounds::normalBoundingNames.lookupOrDefault
......@@ -104,7 +104,7 @@ Foam::interpolation2DTable<Type>::interpolation2DTable(const dictionary& dict)
true // Failsafe behaviour
)
),
fileName_(dict.lookup("file")),
fileName_(dict.get<fileName>("file")),
reader_(tableReader<Type>::New(dict))
{
readTable();
......@@ -117,121 +117,28 @@ Foam::interpolation2DTable<Type>::interpolation2DTable
const interpolation2DTable& interpTable
)
:
List<Tuple2<scalar, List<Tuple2<scalar, Type>>>>(interpTable),
List<value_type>(interpTable),
bounding_(interpTable.bounding_),
fileName_(interpTable.fileName_),
reader_(interpTable.reader_) // note: steals reader. Used in write().
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Type Foam::interpolation2DTable<Type>::interpolateValue
(
const List<Tuple2<scalar, Type>>& data,
const scalar lookupValue
const List<Tuple2<scalar, Type>>& list,
scalar lookupValue
) const
{
const label n = data.size();
const scalar minLimit = data.first().first();
const scalar maxLimit = data.last().first();
if (lookupValue < minLimit)
{
switch (bounding_)
{
case bounds::normalBounding::ERROR:
{
FatalErrorInFunction
<< "value (" << lookupValue << ") less than lower "
<< "bound (" << minLimit << ")" << nl
<< exit(FatalError);
break;
}
case bounds::normalBounding::WARN:
{
WarningInFunction
<< "value (" << lookupValue << ") less than lower "
<< "bound (" << minLimit << ")" << nl
<< " Continuing with the first entry"
<< endl;
// Behaviour as per CLAMP
return data.first().second();
break;
}
case bounds::normalBounding::CLAMP:
{
return data.first().second();
break;
}
}
}
else if (lookupValue >= maxLimit)
{
switch (bounding_)
{
case bounds::normalBounding::ERROR:
{
FatalErrorInFunction
<< "value (" << lookupValue << ") greater than upper "
<< "bound (" << maxLimit << ")" << nl
<< exit(FatalError);
break;
}
case bounds::normalBounding::WARN:
{
WarningInFunction
<< "value (" << lookupValue << ") greater than upper "
<< "bound (" << maxLimit << ")" << nl
<< " Continuing with the last entry"
<< endl;
// Behaviour as per CLAMP
return data.last().second();
break;
}
case bounds::normalBounding::CLAMP:
{
return data.last().second();
break;
}
}
}
// look for the correct range in X
label lo = 0;
label hi = 0;
for (label i = 0; i < n; ++i)
{
if (lookupValue >= data[i].first())
{
lo = hi = i;
}
else
{
hi = i;
break;
}
}
if (lo == hi)
{
return data[lo].second();
}
else
{
Type m =
(data[hi].second() - data[lo].second())
/(data[hi].first() - data[lo].first());
// normal interpolation
return data[lo].second() + m*(lookupValue - data[lo].first());
}
return interpolationTable<Type>::interpolateValue
(
list,
lookupValue,
bounds::repeatableBounding(bounding_)
);
}
......@@ -244,7 +151,7 @@ Foam::label Foam::interpolation2DTable<Type>::Xi
const bool reverse
) const
{
const table& t = *this;
const List<value_type>& t = *this;
label limitI = 0;
if (reverse)
......@@ -259,15 +166,14 @@ Foam::label Foam::interpolation2DTable<Type>::Xi
case bounds::normalBounding::ERROR:
{
FatalErrorInFunction
<< "value (" << valueX << ") out of bounds"
<< "value (" << valueX << ") out of bounds" << nl
<< exit(FatalError);
break;
}
case bounds::normalBounding::WARN:
{
WarningInFunction
<< "value (" << valueX << ") out of bounds"
<< endl;
<< "value (" << valueX << ") out of bounds" << nl;
// Behaviour as per CLAMP
return limitI;
......@@ -290,11 +196,11 @@ Foam::label Foam::interpolation2DTable<Type>::Xi
label i = 0;
if (reverse)
{
label nX = t.size();
const label nX = t.size();
i = 0;
while ((i < nX) && (valueX > t[i].first()))
{
i++;
++i;
}
}
else
......@@ -302,7 +208,7 @@ Foam::label Foam::interpolation2DTable<Type>::Xi
i = t.size() - 1;
while ((i > 0) && (valueX < t[i].first()))
{
i--;
--i;
}
}
......@@ -319,66 +225,61 @@ Type Foam::interpolation2DTable<Type>::operator()
const scalar valueY
) const
{
// Considers all of the list in Y being equal
const label nX = this->size();
const List<value_type>& t = *this;
const table& t = *this;
// Assumes all of the list in Y being equal length
const label nX = t.size();
if (nX == 0)
{
WarningInFunction
<< "cannot interpolate a zero-sized table - returning zero" << endl;
<< "Cannot interpolate zero-sized table - returning zero" << nl;
return Zero;
}
else if (nX == 1)
{
// only 1 column (in X) - interpolate to find Y value
// Only 1 column (in X) - simply interpolate to find Y value
return interpolateValue(t.first().second(), valueY);
}
else
// Find low and high indices in the X range that bound valueX
const label lo = Xi(lessOp<scalar>(), valueX, false);
const label hi = Xi(greaterOp<scalar>(), valueX, true);
if (lo == hi)
{
// have 2-D data, interpolate
return interpolateValue(t[lo].second(), valueY);
}
// find low and high indices in the X range that bound valueX
const label x0i = Xi(lessOp<scalar>(), valueX, false);
const label x1i = Xi(greaterOp<scalar>(), valueX, true);
if (x0i == x1i)
{
return interpolateValue(t[x0i].second(), valueY);
}
else
{
Type y0(interpolateValue(t[x0i].second(), valueY));
Type y1(interpolateValue(t[x1i].second(), valueY));
// Normal interpolation
// gradient in X
const scalar x0 = t[x0i].first();
const scalar x1 = t[x1i].first();
Type mX = (y1 - y0)/(x1 - x0);
const Type y0(interpolateValue(t[lo].second(), valueY));
const Type y1(interpolateValue(t[hi].second(), valueY));
// interpolate
return y0 + mX*(valueX - x0);
}
}
const scalar& x0 = t[lo].first();
const scalar& x1 = t[hi].first();
return (y0 + (y1 - y0)*(valueX - x0)/(x1 - x0));
}
template<class Type>
void Foam::interpolation2DTable<Type>::checkOrder() const
void Foam::interpolation2DTable<Type>::check() const
{
const label n = this->size();
const table& t = *this;
const List<value_type>& list = *this;
scalar prevValue = t[0].first();
scalar prevValue(0);
for (label i=1; i<n; ++i)
label i = 0;
for (const auto& item : list)
{
const scalar currValue = t[i].first();
const scalar& currValue = item.first();
// avoid duplicate values (divide-by-zero error)
if (currValue <= prevValue)
// Avoid duplicate values (divide-by-zero error)
if (i && currValue <= prevValue)
{
FatalErrorInFunction
<< "out-of-order value: "
......@@ -386,6 +287,7 @@ void Foam::interpolation2DTable<Type>::checkOrder() const
<< exit(FatalError);
}
prevValue = currValue;
++i;
}
}
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -27,8 +28,8 @@ Class
Foam::interpolation2DTable
Description
2D table interpolation. The data must be in ascending order in both
dimensions x and y.
2D table interpolation.
The data must be in ascending order in both dimensions x and y.
SourceFiles
interpolation2DTable.C
......@@ -38,10 +39,7 @@ SourceFiles
#ifndef interpolation2DTable_H
#define interpolation2DTable_H
#include "List.H"
#include "Tuple2.H"
#include "tableBounds.H"
#include "tableReader.H"
#include "interpolationTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -57,17 +55,7 @@ class interpolation2DTable
:
public List<Tuple2<scalar, List<Tuple2<scalar, Type>>>>
{
public:
// Public data types
//- Convenience typedef
typedef List<Tuple2<scalar, List<Tuple2<scalar, Type>>>> table;
private:
// Private data
// Private Data
//- Handling for out-of-bound values
bounds::normalBounding bounding_;
......@@ -84,11 +72,11 @@ private:
//- Read the table of data from file
void readTable();
//- Return interpolated value in List
//- Interpolated value within 1-D list
Type interpolateValue
(
const List<Tuple2<scalar, Type>>& data,
const scalar
const List<Tuple2<scalar, Type>>& list,
scalar lookupValue
) const;
//- Return an X index from the matrix
......@@ -103,6 +91,15 @@ private:
public:
// Public Data Types
//- The element data type
typedef Tuple2<scalar, List<Tuple2<scalar, Type>>> value_type;
//- Convenience typedef
typedef List<Tuple2<scalar, List<Tuple2<scalar, Type>>>> table;
// Constructors
//- Construct null
......@@ -117,12 +114,12 @@ public:
);
//- Construct given the name of the file containing the table of data
interpolation2DTable(const fileName& fName);
explicit interpolation2DTable(const fileName& fName);
//- Construct by reading file name and outOfBounds from dictionary
interpolation2DTable(const dictionary& dict);
explicit interpolation2DTable(const dictionary& dict);
//- Construct copy
//- Copy construct
interpolation2DTable(const interpolation2DTable& interpTable);
......@@ -130,7 +127,7 @@ public:
//- Check that list is monotonically increasing
// Exit with a FatalError if there is a problem
void checkOrder() const;
void check() const;
//- Write
void write(Ostream& os) const;
......@@ -138,11 +135,15 @@ public:
// Member Operators
//- Return an element of constant List<Tuple2<scalar, Type>>
const List<Tuple2<scalar, Type>>& operator[](const label) const;
//- Return an interpolated value
Type operator()(const scalar valueX, const scalar valueY) const;
// Housekeeping
//- Deprecated(2019-08) check list is monotonically increasing
// \deprecated(2019-08) - older name for check() method
void checkOrder() const { check(); }
};
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -36,7 +37,6 @@ Description
If \c repeat is chosen for the out-of-bounds handling, the final time
value is treated as being equivalent to time=0 for the following periods.
The construct from dictionary reads a filename from a dictionary and
has an optional readerType. Default is to read OpenFOAM format. The only
other format is csv (comma separated values):
......@@ -65,10 +65,10 @@ SourceFiles
#include "List.H"
#include "Tuple2.H"
#include "tableBounds.H"
#include "tableReader.H"
#include "autoPtr.H"
#include "Field.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -84,7 +84,7 @@ class interpolationTable
:
public List<Tuple2<scalar, Type>>
{
// Private data
// Private Data
//- Handling for out-of-bound values
bounds::repeatableBounding bounding_;
......@@ -104,6 +104,15 @@ class interpolationTable
public:
// Public Data Types
//- The element data type
typedef Tuple2<scalar, Type> value_type;
//- The mapped data type
typedef Type mapped_type;
// Constructors
//- Construct null
......@@ -118,19 +127,28 @@ public:
);
//- Construct given the name of the file containing the table of data
interpolationTable(const fileName& fName);
explicit interpolationTable(const fileName& fName);
//- Construct by reading file name and outOfBounds from dictionary
// and read the table from that file.
//- and read the table from that file.
// This is a specialised constructor used by patchFields
interpolationTable(const dictionary& dict);
explicit interpolationTable(const dictionary& dict);
//- Construct copy
//- Copy construct
interpolationTable(const interpolationTable& interpTable);
// Member Functions
//- Return an interpolated value in List
static Type interpolateValue
(
const List<Tuple2<scalar, Type>>& list,
scalar lookupValue,
bounds::repeatableBounding = bounds::repeatableBounding::CLAMP
);
//- Check that list is monotonically increasing
// Exit with a FatalError if there is a problem
void check() const;
......@@ -139,17 +157,26 @@ public:
void write(Ostream& os) const;
//- Return the rate of change at the interpolation location
// for the give value
Type rateOfChange(const scalar) const;
//- for the given lookup value
Type rateOfChange(scalar lookupValue) const;
//- Return an interpolated value
Type interpolateValue(scalar lookupValue) const;
//- Return multiple interpolated values
tmp<Field<Type>> interpolateValues
(
const UList<scalar>& vals
) const;
// Member Operators
//- Return an element of constant Tuple2<scalar, Type>
const Tuple2<scalar, Type>& operator[](const label) const;
const Tuple2<scalar, Type>& operator[](label idx) const;
//- Return an interpolated value
Type operator()(const scalar) const;
Type operator()(scalar lookupValue) const;
};
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -38,7 +39,7 @@ Foam::Function1Types::Table<Type>::Table
:
TableBase<Type>(entryName, dict)
{
Istream& is(dict.lookup(entryName));
Istream& is = dict.lookup(entryName);
word entryType(is);
is >> this->table_;
TableBase<Type>::check();
......@@ -52,11 +53,4 @@ Foam::Function1Types::Table<Type>::Table(const Table<Type>& tbl)
{}