From 78bd261ce98b915c6b0673f646518ed7e27b60de Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 17 Apr 2012 18:02:25 +0100
Subject: [PATCH] ENH: interpolationWeights: run-time selectable interpolation

---
 src/OpenFOAM/Make/files                       |   7 +
 .../GeometricField/uniformInterpolate.C       | 142 ++++++++++
 .../GeometricField/uniformInterpolate.H       |  83 ++++++
 .../interpolationWeights.C                    | 121 +++++++++
 .../interpolationWeights.H                    | 144 ++++++++++
 .../interpolationWeightsTemplates.C           | 174 ++++++++++++
 .../linearInterpolationWeights.C              | 249 ++++++++++++++++++
 .../linearInterpolationWeights.H              | 120 +++++++++
 .../splineInterpolationWeights.C              | 228 ++++++++++++++++
 .../splineInterpolationWeights.H              | 121 +++++++++
 .../functions/DataEntry/Table/TableBase.C     | 189 ++++++++-----
 .../functions/DataEntry/Table/TableBase.H     |  23 +-
 12 files changed, 1529 insertions(+), 72 deletions(-)
 create mode 100644 src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C
 create mode 100644 src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.H
 create mode 100644 src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C
 create mode 100644 src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.H
 create mode 100644 src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeightsTemplates.C
 create mode 100644 src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.C
 create mode 100644 src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.H
 create mode 100644 src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.C
 create mode 100644 src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.H

diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index c58cf1afc29..948ec4dbe7a 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -576,6 +576,13 @@ $(interpolations)/interpolationTable/tableReaders/tableReaders.C
 $(interpolations)/interpolationTable/tableReaders/openFoam/openFoamTableReaders.C
 $(interpolations)/interpolationTable/tableReaders/csv/csvTableReaders.C
 
+interpolationWeights = $(interpolations)/interpolationWeights
+$(interpolationWeights)/interpolationWeights/interpolationWeights.C
+$(interpolationWeights)/linearInterpolationWeights/linearInterpolationWeights.C
+$(interpolationWeights)/splineInterpolationWeights/splineInterpolationWeights.C
+
+
+
 algorithms/indexedOctree/indexedOctreeName.C
 algorithms/indexedOctree/treeDataCell.C
 
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C
new file mode 100644
index 00000000000..968f5f4c0f0
--- /dev/null
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C
@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class GeoField>
+Foam::tmp<GeoField> Foam::uniformInterpolate
+(
+    const HashPtrTable<GeoField, label, Hash<label> >& fields,
+    const labelList& indices,
+    const scalarField& weights
+)
+{
+    const GeoField& field0 = *(*fields.begin());
+
+    // Interpolate
+    tmp<GeoField> tfld
+    (
+        new GeoField
+        (
+            IOobject
+            (
+                "uniformInterpolate(" + field0.name() + ')',
+                field0.time().timeName(),
+                field0.db(),
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            weights[0]*(*fields[indices[0]])
+        )
+    );
+    GeoField& fld = tfld();
+
+    for (label i = 1; i < indices.size(); ++i)
+    {
+        fld += weights[i]*(*fields[indices[i]]);
+    }
+
+    return tfld;
+}
+
+
+template<class GeoField>
+Foam::tmp<GeoField> Foam::uniformInterpolate
+(
+    const IOobject& fieldIO,
+    const word& fieldName,
+    const wordList& times,
+    const scalarField& weights,
+    const objectRegistry& fieldsCache
+)
+{
+    // Look up the first field
+    const objectRegistry& time0Fields = fieldsCache.lookupObject
+    <
+        const objectRegistry
+    >
+    (
+        times[0]
+    );
+    const GeoField& field0 = time0Fields.lookupObject
+    <
+        const GeoField
+    >
+    (
+        fieldName
+    );
+
+
+    // Interpolate
+    tmp<GeoField> tfld(new GeoField(fieldIO, weights[0]*field0));
+    GeoField& fld = tfld();
+
+    for (label i = 1; i < times.size(); ++i)
+    {
+        const objectRegistry& timeIFields = fieldsCache.lookupObject
+        <
+            const objectRegistry
+        >
+        (
+            times[i]
+        );
+        const GeoField& fieldI = timeIFields.lookupObject
+        <
+            const GeoField
+        >
+        (
+            fieldName
+        );
+
+        fld += weights[i]*fieldI;
+    }
+
+    return tfld;
+}
+
+
+template<class GeoField>
+Foam::tmp<GeoField> Foam::uniformInterpolate
+(
+    const IOobject& fieldIO,
+    const word& fieldName,
+    const wordList& times,
+    const scalarField& weights,
+    const word& registryName
+)
+{
+    return uniformInterpolate<GeoField>
+    (
+        fieldIO,
+        fieldName,
+        times,
+        weights,
+        fieldIO.db().subRegistry(registryName, true)
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.H
new file mode 100644
index 00000000000..f196f71c192
--- /dev/null
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.H
@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "GeometricField.H"
+#include "HashPtrTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Global functions  * * * * * * * * * * * * * //
+
+//- Interpolate selected fields (given by indices and corresponding weights)
+//  (boundary type becomes calculated). Fields stored per index. Field gets name
+//  "uniformInterpolate(" + fld.name() + ')'.
+template<class GeoField>
+tmp<GeoField> uniformInterpolate
+(
+    const HashPtrTable<GeoField, label, Hash<label> >& fields,
+    const labelList& indices,
+    const scalarField& weights
+);
+
+//- Interpolate fields. fieldsCache contains per timeName all loaded fields.
+//  Resulting field gets properties according to fieldIO
+template<class GeoField>
+tmp<GeoField> uniformInterpolate
+(
+    const IOobject& fieldIO,
+    const word& fieldName,
+    const wordList& times,
+    const scalarField& weights,
+    const objectRegistry& fieldsCache
+);
+
+//- Interpolate fields. fieldsCache contains per timeName all loaded fields.
+//  Resulting field gets properties according to fieldIO
+template<class GeoField>
+tmp<GeoField> uniformInterpolate
+(
+    const IOobject& fieldIO,
+    const word& fieldName,
+    const wordList& times,
+    const scalarField& weights,
+    const word& registryName = "fieldsCache"
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "uniformInterpolate.C"
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C
new file mode 100644
index 00000000000..30fc2fc0f69
--- /dev/null
+++ b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C
@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "interpolationWeights.H"
+#include "addToRunTimeSelectionTable.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(interpolationWeights, 0);
+defineRunTimeSelectionTable(interpolationWeights, word);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+interpolationWeights::interpolationWeights
+(
+    const scalarField& samples
+)
+:
+    samples_(samples)
+{}
+
+
+// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
+
+autoPtr<interpolationWeights> interpolationWeights::New
+(
+    const word& type,
+    const scalarField& samples
+)
+{
+    Info<< nl << "Selecting interpolationWeights "
+        << type << endl;
+
+    wordConstructorTable::iterator cstrIter =
+        wordConstructorTablePtr_->find(type);
+
+    if (cstrIter == wordConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "interpolationWeights::New(const word&, "
+            "const scalarField&)"
+        )   << "Unknown interpolationWeights type "
+            << type
+            << endl << endl
+            << "Valid interpolationWeights types are :" << endl
+            << wordConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<interpolationWeights>(cstrIter()(samples));
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+interpolationWeights::~interpolationWeights()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+//objectRegistry& interpolationWeights::registry
+//(
+//    const objectRegistry& obr,
+//    const word& name
+//)
+//{
+//    if (!obr.foundObject<objectRegistry>(name))
+//    {
+//        objectRegistry* fieldsCachePtr = new objectRegistry
+//        (
+//            IOobject
+//            (
+//                name,
+//                obr.time().constant(),
+//                obr,
+//                IOobject::NO_READ,
+//                IOobject::NO_WRITE
+//            )
+//        );
+//        fieldsCachePtr->store();
+//    }
+//    return const_cast<objectRegistry&>(obr.subRegistry(name));
+//}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.H b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.H
new file mode 100644
index 00000000000..889076a1ee2
--- /dev/null
+++ b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.H
@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::interpolationWeights
+
+Description
+    Abstract base class for interpolating in 1D
+
+SourceFiles
+    interpolationWeights.C
+    interpolationWeightsTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interpolationWeights_H
+#define interpolationWeights_H
+
+#include "scalarField.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "pointField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class fvMesh;
+class objectRegistry;
+
+/*---------------------------------------------------------------------------*\
+                       Class interpolationWeights Declaration
+\*---------------------------------------------------------------------------*/
+
+class interpolationWeights
+{
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        interpolationWeights(const interpolationWeights&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const interpolationWeights&);
+
+protected:
+
+        const scalarField& samples_;
+
+public:
+
+    //- Runtime type information
+    TypeName("interpolationWeights");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            interpolationWeights,
+            word,
+            (
+                const scalarField& samples
+            ),
+            (samples)
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        interpolationWeights(const scalarField& samples);
+
+
+    // Selectors
+
+        //- Return a reference to the selected interpolationWeights
+        static autoPtr<interpolationWeights> New
+        (
+            const word& type,
+            const scalarField& samples
+        );
+
+
+    //- Destructor
+    virtual ~interpolationWeights();
+
+
+    // Member Functions
+
+        //- Calculate weights and indices to calculate t from samples.
+        //  Returns true if indices changed.
+        virtual bool valueWeights
+        (
+            const scalar t,
+            labelList& indices,
+            scalarField& weights
+        ) const = 0;
+
+        //- Calculate weights and indices to calculate integrand of t1..t2
+        //  from samples. Returns true if indices changed.
+        virtual bool integrationWeights
+        (
+            const scalar t1,
+            const scalar t2,
+            labelList& indices,
+            scalarField& weights
+        ) const = 0;
+
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeightsTemplates.C b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeightsTemplates.C
new file mode 100644
index 00000000000..aac670138d9
--- /dev/null
+++ b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeightsTemplates.C
@@ -0,0 +1,174 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "interpolationWeights.H"
+#include "ListOps.H"
+#include "IOobject.H"
+#include "HashSet.H"
+#include "objectRegistry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+//template<class GeoField>
+//void interpolationWeights::readFields
+//(
+//    const word& fieldName,
+//    const typename GeoField::Mesh& mesh,
+//    const wordList& timeNames,
+//    objectRegistry& fieldsCache
+//)
+//{
+//    // Collect all times that are no longer used
+//    {
+//        HashSet<word> usedTimes(timeNames);
+//
+//        DynamicList<word> unusedTimes(fieldsCache.size());
+//
+//        forAllIter(objectRegistry, fieldsCache, timeIter)
+//        {
+//            const word& tm = timeIter.key();
+//            if (!usedTimes.found(tm))
+//            {
+//                unusedTimes.append(tm);
+//            }
+//        }
+//
+//        //Info<< "Unloading times " << unusedTimes << endl;
+//
+//        forAll(unusedTimes, i)
+//        {
+//            objectRegistry& timeCache = const_cast<objectRegistry&>
+//            (
+//                fieldsCache.lookupObject<objectRegistry>(unusedTimes[i])
+//            );
+//            fieldsCache.checkOut(timeCache);
+//        }
+//    }
+//
+//
+//    // Load any new fields
+//    forAll(timeNames, i)
+//    {
+//        const word& tm = timeNames[i];
+//
+//        // Create if not found
+//        if (!fieldsCache.found(tm))
+//        {
+//            //Info<< "Creating registry for time " << tm << endl;
+//
+//            // Create objectRegistry if not found
+//            objectRegistry* timeCachePtr = new objectRegistry
+//            (
+//                IOobject
+//                (
+//                    tm,
+//                    tm,
+//                    fieldsCache,
+//                    IOobject::NO_READ,
+//                    IOobject::NO_WRITE
+//                )
+//            );
+//            timeCachePtr->store();
+//        }
+//
+//        // Obtain cache for current time
+//        const objectRegistry& timeCache =
+//            fieldsCache.lookupObject<objectRegistry>
+//            (
+//                tm
+//            );
+//
+//        // Store field if not found
+//        if (!timeCache.found(fieldName))
+//        {
+//            //Info<< "Loading field " << fieldName
+//            //    << " for time " << tm << endl;
+//
+//            GeoField loadedFld
+//            (
+//                IOobject
+//                (
+//                    fieldName,
+//                    tm,
+//                    mesh.thisDb(),
+//                    IOobject::MUST_READ,
+//                    IOobject::NO_WRITE,
+//                    false
+//                ),
+//                mesh
+//            );
+//
+//            // Transfer to timeCache (new objectRegistry and store flag)
+//            GeoField* fldPtr = new GeoField
+//            (
+//                IOobject
+//                (
+//                    fieldName,
+//                    tm,
+//                    timeCache,
+//                    IOobject::NO_READ,
+//                    IOobject::NO_WRITE
+//                ),
+//                loadedFld
+//            );
+//            fldPtr->store();
+//        }
+//    }
+//}
+//
+//
+//template<class GeoField>
+//void interpolationWeights::readFields
+//(
+//    const word& fieldName,
+//    const typename GeoField::Mesh& mesh,
+//    const wordList& timeNames,
+//    const word& registryName
+//)
+//{
+//    readFields<GeoField>
+//    (
+//        fieldName,
+//        mesh,
+//        timeNames,
+//        //registry(mesh.thisDb(), registryName)
+//        const_cast<objectRegistry&>
+//        (
+//            mesh.thisDb().subRegistry(registryName, true)
+//        )
+//    );
+//}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.C b/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.C
new file mode 100644
index 00000000000..18233f1825e
--- /dev/null
+++ b/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.C
@@ -0,0 +1,249 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "linearInterpolationWeights.H"
+#include "addToRunTimeSelectionTable.H"
+#include "ListOps.H"
+#include "Pair.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(linearInterpolationWeights, 0);
+addToRunTimeSelectionTable
+(
+    interpolationWeights,
+    linearInterpolationWeights,
+    word
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::Pair<Foam::scalar> linearInterpolationWeights::integrationWeights
+(
+    const label i,
+    const scalar t
+) const
+{
+    // t is in range samples_[i] .. samples_[i+1]
+
+    scalar s = (t-samples_[i])/(samples_[i+1]-samples_[i]);
+
+    if (s < -SMALL || s > 1+SMALL)
+    {
+        FatalErrorIn("linearInterpolationWeights::integrationWeights(..)")
+            << "Value " << t << " outside range " << samples_[i]
+            << " .. " << samples_[i+1]
+            << exit(FatalError);
+    }
+
+    scalar d = samples_[i+1]-t;
+
+    return Pair<scalar>(d*0.5*(1-s), d*0.5*(1+s));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+linearInterpolationWeights::linearInterpolationWeights
+(
+    const scalarField& samples
+)
+:
+    interpolationWeights(samples),
+    index_(-1)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool linearInterpolationWeights::valueWeights
+(
+    const scalar t,
+    labelList& indices,
+    scalarField& weights
+) const
+{
+    bool indexChanged = false;
+
+    // Check if current timeIndex is still valid
+    if
+    (
+        index_ >= 0
+     && index_ < samples_.size()
+     && (
+            samples_[index_] <= t
+         && (index_ == samples_.size()-1 || t <= samples_[index_+1])
+        )
+    )
+    {
+        // index_ still at correct slot
+    }
+    else
+    {
+        // search for correct index
+        index_ = findLower(samples_, t);
+        indexChanged = true;
+    }
+
+
+    if (index_ == -1)
+    {
+        // Use first element only
+        indices.setSize(1);
+        weights.setSize(1);
+        indices[0] = 0;
+        weights[0] = 1.0;
+    }
+    else if (index_ == samples_.size()-1)
+    {
+        // Use last element only
+        indices.setSize(1);
+        weights.setSize(1);
+        indices[0] = samples_.size()-1;
+        weights[0] = 1.0;
+    }
+    else
+    {
+        // Interpolate
+        indices.setSize(2);
+        weights.setSize(2);
+
+        indices[0] = index_;
+        indices[1] = index_+1;
+
+        scalar t0 = samples_[indices[0]];
+        scalar t1 = samples_[indices[1]];
+        scalar deltaT = t1-t0;
+
+        weights[0] = (t1-t)/deltaT;
+        weights[1] = 1.0-weights[0];
+    }
+
+    return indexChanged;
+}
+
+
+bool linearInterpolationWeights::integrationWeights
+(
+    const scalar t1,
+    const scalar t2,
+    labelList& indices,
+    scalarField& weights
+) const
+{
+    if (t2 < t1-VSMALL)
+    {
+        FatalErrorIn("linearInterpolationWeights::integrationWeights(..)")
+            << "Integration should be in positive direction."
+            <<  " t1:" << t1 << " t2:" << t2
+            << exit(FatalError);
+    }
+
+    // Currently no fancy logic on cached index
+    label i1 = findLower(samples_, t1);
+    label i2 = findLower(samples_, t2);
+
+    // For now just fail if any outside table
+    if (i1 == -1 || i2 == samples_.size()-1)
+    {
+        FatalErrorIn("linearInterpolationWeights::integrationWeights(..)")
+            << "Integrating outside table " << samples_[0] << ".."
+            << samples_.last() << " not implemented."
+            << " t1:" << t1 << " t2:" << t2 << exit(FatalError);
+    }
+
+    label nIndices = i2-i1+2;
+
+
+    // Determine if indices already correct
+    bool anyChanged = false;
+
+    if (nIndices != indices.size())
+    {
+        anyChanged = true;
+    }
+    else
+    {
+        // Closer check
+
+        label index = i1;
+        forAll(indices, i)
+        {
+            if (indices[i] != index)
+            {
+                anyChanged = true;
+                break;
+            }
+            index++;
+        }
+    }
+
+    indices.setSize(nIndices);
+    weights.setSize(nIndices);
+    weights = 0.0;
+
+    // Sum from i1+1 to i2+1
+    for (label i = i1+1; i <= i2; i++)
+    {
+        scalar d = samples_[i+1]-samples_[i];
+        indices[i-i1] = i;
+        weights[i-i1] += 0.5*d;
+        indices[i+1-i1] = i+1;
+        weights[i+1-i1] += 0.5*d;
+    }
+
+    // Add from i1 to t1
+    {
+        Pair<scalar> i1Tot1 = integrationWeights(i1, t1);
+        indices[0] = i1;
+        weights[0] += i1Tot1.first();
+        indices[1] = i1+1;
+        weights[1] += i1Tot1.second();
+    }
+
+    // Subtract from t2 to i2+1
+    {
+        Pair<scalar> wghts = integrationWeights(i2, t2);
+        indices[i2-i1] = i2;
+        weights[i2-i1] += -wghts.first();
+        indices[i2-i1+1] = i2+1;
+        weights[i2-i1+1] += -wghts.second();
+    }
+
+    return anyChanged;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.H b/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.H
new file mode 100644
index 00000000000..871eb30599c
--- /dev/null
+++ b/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.H
@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::linearInterpolationWeights
+
+Description
+
+SourceFiles
+    linearInterpolationWeights.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef linearInterpolationWeights_H
+#define linearInterpolationWeights_H
+
+#include "interpolationWeights.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class linearInterpolationWeights Declaration
+\*---------------------------------------------------------------------------*/
+
+class linearInterpolationWeights
+:
+    public interpolationWeights
+{
+
+private:
+
+    // Private data
+
+        //- Cached index in samples from previous invocation
+        mutable label index_;
+
+    // Private Member Functions
+
+        //- Get weights of i and i+1 to calculate integration from t to
+        //  samples_[i+1]
+        Pair<scalar> integrationWeights
+        (
+            const label i,
+            const scalar t
+        ) const;
+
+public:
+
+    //- Runtime type information
+    TypeName("linear");
+
+    // Constructors
+
+        //- Construct from components
+        linearInterpolationWeights
+        (
+            const scalarField& samples
+        );
+
+
+    //- Destructor
+    virtual ~linearInterpolationWeights()
+    {}
+
+
+    // Member Functions
+
+        //- Calculate weights and indices to calculate t from samples.
+        //  Returns true if indices changed.
+        virtual bool valueWeights
+        (
+            const scalar t,
+            labelList& indices,
+            scalarField& weights
+        ) const;
+
+        //- Calculate weights and indices to calculate integrand of t1..t2
+        //  from samples. Returns true if indices changed.
+        virtual bool integrationWeights
+        (
+            const scalar t1,
+            const scalar t2,
+            labelList& indices,
+            scalarField& weights
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.C b/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.C
new file mode 100644
index 00000000000..854146cdff5
--- /dev/null
+++ b/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.C
@@ -0,0 +1,228 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "splineInterpolationWeights.H"
+#include "addToRunTimeSelectionTable.H"
+#include "ListOps.H"
+#include "linearInterpolationWeights.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(splineInterpolationWeights, 0);
+addToRunTimeSelectionTable
+(
+    interpolationWeights,
+    splineInterpolationWeights,
+    word
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+splineInterpolationWeights::splineInterpolationWeights
+(
+    const scalarField& samples,
+    const bool checkEqualDistance
+)
+:
+    interpolationWeights(samples),
+    index_(-1)
+{
+    if (checkEqualDistance && samples_.size() > 2)
+    {
+        const scalar interval = samples_[1]-samples[0];
+        for (label i = 2; i < samples_.size(); i++)
+        {
+            scalar d = samples_[i]-samples[i-1];
+
+            if (mag(d-interval) > SMALL)
+            {
+                WarningIn
+                (
+                    "splineInterpolationWeights::splineInterpolationWeights"
+                    "(const scalarField&)"
+                )   << "Spline interpolation only valid for constant intervals."
+                    << nl
+                    << "Interval 0-1 : " << interval << nl
+                    << "Interval " << i-1 << '-' << i << " : "
+                    << d << endl;
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool splineInterpolationWeights::valueWeights
+(
+    const scalar t,
+    labelList& indices,
+    scalarField& weights
+) const
+{
+    bool indexChanged = false;
+
+    // linear interpolation
+    if (samples_.size() <= 2)
+    {
+        return linearInterpolationWeights(samples_).valueWeights
+        (
+            t,
+            indices,
+            weights
+        );
+    }
+
+    // Check if current timeIndex is still valid
+    if
+    (
+        index_ >= 0
+     && index_ < samples_.size()
+     && (
+            samples_[index_] <= t
+         && (index_ == samples_.size()-1 || t <= samples_[index_+1])
+        )
+    )
+    {
+        // index_ still at correct slot
+    }
+    else
+    {
+        // search for correct index
+        index_ = findLower(samples_, t);
+        indexChanged = true;
+    }
+
+
+    // Clamp if outside table
+    if (index_ == -1)
+    {
+        indices.setSize(1);
+        weights.setSize(1);
+
+        indices[0] = 0;
+        weights[0] = 1;
+        return indexChanged;
+    }
+    else if (index_ == samples_.size()-1)
+    {
+        indices.setSize(1);
+        weights.setSize(1);
+
+        indices[0] = samples_.size()-1;
+        weights[0] = 1;
+        return indexChanged;
+    }
+
+
+
+    label lo = index_;
+    label hi = index_+1;
+
+    // weighting
+    scalar mu = (t - samples_[lo])/(samples_[hi] - samples_[lo]);
+
+    scalar w0 = 0.5*(mu*(-1+mu*(2-mu)));            // coeff of lo-1
+    scalar w1 = 0.5*(2+mu*(mu*(-5 + mu*(3))));      // coeff of lo
+    scalar w2 = 0.5*(mu*(1 + mu*(4 + mu*(-3))));    // coeff of hi
+    scalar w3 = 0.5*(mu*mu*(-1 + mu));              // coeff of hi+1
+
+    if (lo > 0)
+    {
+        if (hi < samples_.size()-1)
+        {
+            // Four points available
+            indices.setSize(4);
+            weights.setSize(4);
+
+            indices[0] = lo-1;
+            indices[1] = lo;
+            indices[2] = hi;
+            indices[3] = hi+1;
+
+            weights[0] = w0;
+            weights[1] = w1;
+            weights[2] = w2;
+            weights[3] = w3;
+        }
+        else
+        {
+            // No y3 available. Extrapolate: y3=3*y2-y1
+            indices.setSize(3);
+            weights.setSize(3);
+
+            indices[0] = lo-1;
+            indices[1] = lo;
+            indices[2] = hi;
+
+            weights[0] = w0;
+            weights[1] = w1 - w3;
+            weights[2] = w2 + 2*w3;
+        }
+    }
+    else
+    {
+        // No y0 available. Extrapolate: y0=2*y1-y2;
+        if (hi < samples_.size()-1)
+        {
+            indices.setSize(3);
+            weights.setSize(3);
+
+            indices[0] = lo;
+            indices[1] = hi;
+            indices[2] = hi+1;
+
+            weights[0] = w1 + 2*w0;
+            weights[1] = w2 - w0;
+            weights[2] = w3;
+        }
+        else
+        {
+            indices.setSize(2);
+            weights.setSize(2);
+
+            indices[0] = lo;
+            indices[1] = hi;
+
+            weights[0] = w1 + 2*w0 - w3;
+            weights[1] = w2 - w0 + 2*w3;
+        }
+    }
+
+    return indexChanged;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.H b/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.H
new file mode 100644
index 00000000000..15660ff90cf
--- /dev/null
+++ b/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.H
@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::splineInterpolationWeights
+
+Description
+    Catmull-Rom spline interpolation.
+
+SourceFiles
+    splineInterpolationWeights.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef splineInterpolationWeights_H
+#define splineInterpolationWeights_H
+
+#include "interpolationWeights.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class splineInterpolationWeights Declaration
+\*---------------------------------------------------------------------------*/
+
+class splineInterpolationWeights
+:
+    public interpolationWeights
+{
+
+private:
+
+    // Private data
+
+        //- Cached index in samples from previous invocation
+        mutable label index_;
+
+public:
+
+    //- Runtime type information
+    TypeName("spline");
+
+    // Constructors
+
+        //- Construct from components. By default make sure samples are
+        //  equidistant.
+        splineInterpolationWeights
+        (
+            const scalarField& samples,
+            const bool checkEqualDistance = true
+        );
+
+
+    //- Destructor
+    virtual ~splineInterpolationWeights()
+    {}
+
+
+    // Member Functions
+
+        //- Calculate weights and indices to calculate t from samples.
+        //  Returns true if indices changed.
+        virtual bool valueWeights
+        (
+            const scalar t,
+            labelList& indices,
+            scalarField& weights
+        ) const;
+
+        //- Calculate weights and indices to calculate integrand of t1..t2
+        //  from samples. Returns true if indices changed.
+        virtual bool integrationWeights
+        (
+            const scalar t1,
+            const scalar t2,
+            labelList& indices,
+            scalarField& weights
+        ) const
+        {
+            notImplemented
+            (
+                "splineInterpolationWeights::integrationWeights(..)"
+            );
+            return false;
+        }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C b/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C
index c977618e2a7..e59c5ec898d 100644
--- a/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C
+++ b/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C
@@ -26,6 +26,29 @@ License
 #include "TableBase.H"
 #include "Time.H"
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+const Foam::interpolationWeights& Foam::TableBase<Type>::interpolator() const
+{
+    if (interpolatorPtr_.empty())
+    {
+        // Re-work table into linear list
+        tableSamples_.setSize(table_.size());
+        forAll(table_, i)
+        {
+            tableSamples_[i] = table_[i].first();
+        }
+        interpolatorPtr_ = interpolationWeights::New
+        (
+            interpolationScheme_,
+            tableSamples_
+        );
+    }
+    return interpolatorPtr_();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type>
@@ -39,6 +62,10 @@ Foam::TableBase<Type>::TableBase(const word& name, const dictionary& dict)
             dict.lookupOrDefault<word>("outOfBounds", "clamp")
         )
     ),
+    interpolationScheme_
+    (
+        dict.lookupOrDefault<word>("interpolationScheme", "linear")
+    ),
     table_(),
     dimensions_(dimless)
 {}
@@ -49,8 +76,11 @@ Foam::TableBase<Type>::TableBase(const TableBase<Type>& tbl)
 :
     name_(tbl.name_),
     boundsHandling_(tbl.boundsHandling_),
+    interpolationScheme_(tbl.interpolationScheme_),
     table_(tbl.table_),
-    dimensions_(tbl.dimensions_)
+    dimensions_(tbl.dimensions_),
+    tableSamples_(tbl.tableSamples_),
+    interpolatorPtr_(tbl.interpolatorPtr_)
 {}
 
 
@@ -307,6 +337,9 @@ void Foam::TableBase<Type>::convertTimeBase(const Time& t)
         scalar value = table_[i].first();
         table_[i].first() = t.userTimeToTime(value);
     }
+
+    tableSamples_.clear();
+    interpolatorPtr_.clear();
 }
 
 
@@ -325,88 +358,104 @@ Type Foam::TableBase<Type>::value(const scalar x) const
         return table_.last().second();
     }
 
-    // Find i such that x(i) < xDash < x(i+1)
-    label i = 0;
-    while ((table_[i+1].first() < xDash) && (i+1 < table_.size()))
+    // Use interpolator
+    interpolator().valueWeights(x, currentIndices_, currentWeights_);
+
+    Type t = currentWeights_[0]*table_[currentIndices_[0]].second();
+    for (label i = 1; i < currentIndices_.size(); i++)
     {
-        i++;
+        t += currentWeights_[i]*table_[currentIndices_[i]].second();
     }
-
-Info <<
-     (xDash - table_[i].first())/(table_[i+1].first() - table_[i].first())
-      * (table_[i+1].second() - table_[i].second())
-      + table_[i].second() << endl;
-
-    // Linear interpolation to find value
-    return Type
-    (
-        (xDash - table_[i].first())/(table_[i+1].first() - table_[i].first())
-      * (table_[i+1].second() - table_[i].second())
-      + table_[i].second()
-    );
+    return t;
+
+    //// Find i such that x(i) < xDash < x(i+1)
+    //label i = 0;
+    //while ((table_[i+1].first() < xDash) && (i+1 < table_.size()))
+    //{
+    //    i++;
+    //}
+    //
+    //// Linear interpolation to find value
+    //return Type
+    //(
+    //    (xDash - table_[i].first())/(table_[i+1].first() - table_[i].first())
+    //  * (table_[i+1].second() - table_[i].second())
+    //  + table_[i].second()
+    //);
 }
 
 
 template<class Type>
 Type Foam::TableBase<Type>::integrate(const scalar x1, const scalar x2) const
 {
-    // Initialise return value
-    Type sum = pTraits<Type>::zero;
-
-    // Return zero if out of bounds
-    if ((x1 > table_.last().first()) || (x2 < table_[0].first()))
-    {
-        return sum;
-    }
-
-    // Find next index greater than x1
-    label id1 = 0;
-    while ((table_[id1].first() < x1) && (id1 < table_.size()))
-    {
-        id1++;
-    }
+    // Use interpolator
+    interpolator().integrationWeights(x1, x2, currentIndices_, currentWeights_);
 
-    // Find next index less than x2
-    label id2 = table_.size() - 1;
-    while ((table_[id2].first() > x2) && (id2 >= 1))
+    Type sum = currentWeights_[0]*table_[currentIndices_[0]].second();
+    for (label i = 1; i < currentIndices_.size(); i++)
     {
-        id2--;
+       sum += currentWeights_[i]*table_[currentIndices_[i]].second();
     }
+    return sum;
 
-    if ((id1 - id2) == 1)
-    {
-        // x1 and x2 lie within 1 interval
-        sum = 0.5*(value(x1) + value(x2))*(x2 - x1);
-    }
-    else
-    {
-        // x1 and x2 cross multiple intervals
-
-        // Integrate table body
-        for (label i=id1; i<id2; i++)
-        {
-            sum +=
-                (table_[i].second() + table_[i+1].second())
-              * (table_[i+1].first() - table_[i].first());
-        }
-        sum *= 0.5;
-
-        // Add table ends (partial segments)
-        if (id1 > 0)
-        {
-            sum += 0.5
-              * (value(x1) + table_[id1].second())
-              * (table_[id1].first() - x1);
-        }
-        if (id2 < table_.size() - 1)
-        {
-            sum += 0.5
-              * (table_[id2].second() + value(x2))
-              * (x2 - table_[id2].first());
-        }
-    }
 
-    return sum;
+    //// Initialise return value
+    //Type sum = pTraits<Type>::zero;
+    //
+    //// Return zero if out of bounds
+    //if ((x1 > table_.last().first()) || (x2 < table_[0].first()))
+    //{
+    //    return sum;
+    //}
+    //
+    //// Find next index greater than x1
+    //label id1 = 0;
+    //while ((table_[id1].first() < x1) && (id1 < table_.size()))
+    //{
+    //    id1++;
+    //}
+    //
+    //// Find next index less than x2
+    //label id2 = table_.size() - 1;
+    //while ((table_[id2].first() > x2) && (id2 >= 1))
+    //{
+    //    id2--;
+    //}
+    //
+    //if ((id1 - id2) == 1)
+    //{
+    //    // x1 and x2 lie within 1 interval
+    //    sum = 0.5*(value(x1) + value(x2))*(x2 - x1);
+    //}
+    //else
+    //{
+    //    // x1 and x2 cross multiple intervals
+    //
+    //    // Integrate table body
+    //    for (label i=id1; i<id2; i++)
+    //    {
+    //        sum +=
+    //            (table_[i].second() + table_[i+1].second())
+    //          * (table_[i+1].first() - table_[i].first());
+    //    }
+    //    sum *= 0.5;
+    //
+    //    // Add table ends (partial segments)
+    //    if (id1 > 0)
+    //    {
+    //        sum += 0.5
+    //          * (value(x1) + table_[id1].second())
+    //          * (table_[id1].first() - x1);
+    //    }
+    //    if (id2 < table_.size() - 1)
+    //    {
+    //        sum += 0.5
+    //          * (table_[id2].second() + value(x2))
+    //          * (x2 - table_[id2].first());
+    //    }
+    //}
+    //
+    //return sum;
 }
 
 
diff --git a/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.H b/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.H
index 15164b3d903..ce4c3d22879 100644
--- a/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.H
+++ b/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.H
@@ -38,6 +38,7 @@ SourceFiles
 #include "DataEntry.H"
 #include "Tuple2.H"
 #include "dimensionSet.H"
+#include "interpolationWeights.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,10 +81,13 @@ protected:
     // Protected data
 
         //- Table name
-        word name_;
+        const word name_;
 
         //- Enumeration for handling out-of-bound values
-        boundsHandling boundsHandling_;
+        const boundsHandling boundsHandling_;
+
+        //- Interpolation type
+        const word interpolationScheme_;
 
         //- Table data
         List<Tuple2<scalar, Type> > table_;
@@ -91,9 +95,24 @@ protected:
         //- The dimension set
         dimensionSet dimensions_;
 
+        //- Extracted values
+        mutable scalarField tableSamples_;
+
+        //- Interpolator method
+        mutable autoPtr<interpolationWeights> interpolatorPtr_;
+
+        //- Cached indices and weights
+        mutable labelList currentIndices_;
+
+        mutable scalarField currentWeights_;
+
 
     // Protected Member Functions
 
+
+        //- Return (demand driven) interpolator
+        const interpolationWeights& interpolator() const;
+
         //- Disallow default bitwise assignment
         void operator=(const TableBase<Type>&);
 
-- 
GitLab