diff --git a/applications/test/rawIOField/Make/files b/applications/test/rawIOField/Make/files
new file mode 100644
index 0000000000..a37b434282
--- /dev/null
+++ b/applications/test/rawIOField/Make/files
@@ -0,0 +1,3 @@
+Test-rawIOField.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-rawIOField
diff --git a/applications/test/rawIOField/Make/options b/applications/test/rawIOField/Make/options
new file mode 100644
index 0000000000..b6fa4bff16
--- /dev/null
+++ b/applications/test/rawIOField/Make/options
@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfileFormats \
+    -lsurfMesh \
+    -lmeshTools
diff --git a/applications/test/rawIOField/Test-rawIOField.C b/applications/test/rawIOField/Test-rawIOField.C
new file mode 100644
index 0000000000..43bd4979ce
--- /dev/null
+++ b/applications/test/rawIOField/Test-rawIOField.C
@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Application
+    Test-rawIOField
+
+Description
+    Reading rawIOField from disk
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "Switch.H"
+#include "primitiveFields.H"
+#include "pointField.H"
+#include "rawIOField.H"
+#include "exprTraits.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+using namespace Foam;
+
+#undef  USE_ROOT_CASE
+//#define USE_ROOT_CASE
+
+
+template<class Type>
+tmp<Field<Type>> readRawField
+(
+    const IOobject& io,
+    IOobjectOption::readOption withAverage
+)
+{
+    rawIOField<Type> raw(io, withAverage);
+
+    Info<< "File: " << io.objectPath() << nl
+        << "Read: " << raw.size()
+        << ' ' << pTraits<Type>::typeName << " entries" << nl
+        << "Average: " << Switch::name(raw.hasAverage())
+        << " = " << raw.average() << endl;
+
+    return tmp<Field<Type>>::New(std::move(static_cast<Field<Type>&>(raw)));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Test behaviour of rawIOField reading (writing?)"
+    );
+
+    argList::noCheckProcessorDirectories();
+
+    argList::addBoolOption("scalar", "Read scalar field");
+    argList::addBoolOption("vector", "Read vector field");
+    argList::addBoolOption("point", "Read point field");
+    argList::addBoolOption("average", "Require averaged value entry");
+    argList::addBoolOption("try-average", "Optional averaged value entry");
+
+    argList::addArgument("fileName");
+
+    #ifdef USE_ROOT_CASE
+        #include "setRootCase.H"
+        #include "createTime.H"
+    #else
+    // Without root case, or time
+    argList args(argc, argv);
+    #endif
+
+    fileName inputName = args.get<fileName>(1);
+
+    IOobjectOption::readOption withAverage = IOobjectOption::NO_READ;
+
+    if (args.found("average"))
+    {
+        withAverage = IOobjectOption::MUST_READ;
+    }
+    else if (args.found("try-average"))
+    {
+        withAverage = IOobjectOption::READ_IF_PRESENT;
+    }
+
+    Info<< "Using case: " << argList::envGlobalPath() << nl
+        << "Read file:  " << inputName << nl
+        << "with average: " << int(withAverage) << nl
+        << endl;
+
+
+    refPtr<Time> timePtr;
+
+    #ifdef USE_ROOT_CASE
+    timePtr.cref(runTime);
+    #endif
+
+    // Fallback (eg, no runTime)
+    if (!timePtr.good())
+    {
+        timePtr.reset(Time::New(argList::envGlobalPath()));
+    }
+
+
+    const auto& tm = timePtr();
+
+    fileName resolvedName(inputName);
+    resolvedName.toAbsolute();
+
+    IOobject io
+    (
+        resolvedName,   // absolute path
+        tm,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        IOobject::NO_REGISTER,
+        true            // is global object (currently not used)
+    );
+
+
+    if (args.found("scalar"))
+    {
+        auto tfield = readRawField<scalar>(io, withAverage);
+    }
+    else if (args.found("point"))
+    {
+        auto tfield = readRawField<point>(io, withAverage);
+    }
+    else if (args.found("vector"))
+    {
+        auto tfield = readRawField<vector>(io, withAverage);
+    }
+    else
+    {
+        Info<< "no data type specified!\n";
+    }
+
+    Info<< nl << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/db/IOobjects/GlobalIOField/GlobalIOField.C b/src/OpenFOAM/db/IOobjects/GlobalIOField/GlobalIOField.C
index 8f60905464..08853995bb 100644
--- a/src/OpenFOAM/db/IOobjects/GlobalIOField/GlobalIOField.C
+++ b/src/OpenFOAM/db/IOobjects/GlobalIOField/GlobalIOField.C
@@ -43,7 +43,7 @@ Foam::GlobalIOField<Type>::GlobalIOField(const IOobject& io)
 
 
 template<class Type>
-Foam::GlobalIOField<Type>::GlobalIOField(const IOobject& io, const label size)
+Foam::GlobalIOField<Type>::GlobalIOField(const IOobject& io, const label len)
 :
     regIOobject(io)
 {
@@ -52,7 +52,7 @@ Foam::GlobalIOField<Type>::GlobalIOField(const IOobject& io, const label size)
 
     if (!readHeaderOk(IOstreamOption::BINARY, typeName))
     {
-        Field<Type>::setSize(size);
+        Field<Type>::resize(len);
     }
 }
 
@@ -113,7 +113,7 @@ Foam::GlobalIOField<Type>::GlobalIOField
         Field<Type>::transfer(tfld.ref());
     }
 
-    if (!readHeaderOk(IOstreamOption::BINARY, typeName))
+    if (!readHeaderOk(IOstreamOption::BINARY, typeName) && !reuse)
     {
         Field<Type>::operator=(tfld());
     }
diff --git a/src/OpenFOAM/db/IOobjects/GlobalIOField/GlobalIOField.H b/src/OpenFOAM/db/IOobjects/GlobalIOField/GlobalIOField.H
index 165ea0f22a..95f7871ea0 100644
--- a/src/OpenFOAM/db/IOobjects/GlobalIOField/GlobalIOField.H
+++ b/src/OpenFOAM/db/IOobjects/GlobalIOField/GlobalIOField.H
@@ -35,11 +35,10 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef GlobalIOField_H
-#define GlobalIOField_H
+#ifndef Foam_GlobalIOField_H
+#define Foam_GlobalIOField_H
 
-#include "regIOobject.H"
-#include "Field.H"
+#include "IOField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,9 +55,12 @@ class GlobalIOField
     public regIOobject,
     public Field<Type>
 {
-
 public:
 
+    //- The underlying content type
+    typedef Field<Type> content_type;
+
+    //- Runtime type information
     TypeName("Field");
 
 
@@ -70,17 +72,17 @@ public:
         //- Construct from IOobject
         explicit GlobalIOField(const IOobject& io);
 
-        //- Construct from IOobject and size (does not set values)
-        GlobalIOField(const IOobject& io, const label size);
+        //- Construct from IOobject and field size (if not read)
+        GlobalIOField(const IOobject& io, const label len);
 
-        //- Construct from IOobject and a List/Field content
+        //- Construct from IOobject and copy of List/Field content
         GlobalIOField(const IOobject& io, const UList<Type>& content);
 
         //- Construct by transferring the Field content
         GlobalIOField(const IOobject& io, Field<Type>&& content);
 
         //- Construct by copying/moving tmp content
-        GlobalIOField(const IOobject& io, const tmp<Field<Type>>& tf);
+        GlobalIOField(const IOobject& io, const tmp<Field<Type>>& tfld);
 
 
     //- Destructor
@@ -116,7 +118,6 @@ public:
 
         //- Copy or move assignment of entries
         using Field<Type>::operator=;
-
 };
 
 
diff --git a/src/OpenFOAM/db/IOobjects/GlobalIOList/GlobalIOList.H b/src/OpenFOAM/db/IOobjects/GlobalIOList/GlobalIOList.H
index a38a1097f4..666a02bb52 100644
--- a/src/OpenFOAM/db/IOobjects/GlobalIOList/GlobalIOList.H
+++ b/src/OpenFOAM/db/IOobjects/GlobalIOList/GlobalIOList.H
@@ -35,11 +35,10 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef GlobalIOList_H
-#define GlobalIOList_H
+#ifndef Foam_GlobalIOList_H
+#define Foam_GlobalIOList_H
 
-#include "List.H"
-#include "regIOobject.H"
+#include "IOList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,6 +57,10 @@ class GlobalIOList
 {
 public:
 
+    //- The underlying content type
+    typedef List<Type> content_type;
+
+    //- Runtime type information
     TypeName("List");
 
 
diff --git a/src/OpenFOAM/db/IOobjects/rawIOField/rawIOField.C b/src/OpenFOAM/db/IOobjects/rawIOField/rawIOField.C
index b084f6a985..35f6283f74 100644
--- a/src/OpenFOAM/db/IOobjects/rawIOField/rawIOField.C
+++ b/src/OpenFOAM/db/IOobjects/rawIOField/rawIOField.C
@@ -28,24 +28,54 @@ License
 #include "rawIOField.H"
 #include "IFstream.H"
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class Type>
-Foam::rawIOField<Type>::rawIOField(const IOobject& io, const bool readAverage)
-:
-    regIOobject(io),
-    hasAverage_(false),
-    average_(Zero)
+void Foam::rawIOField<Type>::readContents
+(
+    Istream& is,
+    IOobjectOption::readOption readAverage
+)
 {
-    // Check for MUST_READ_IF_MODIFIED
-    warnNoRereading<rawIOField<Type>>();
+    is >> static_cast<Field<Type>&>(*this);
 
-    if (io.isReadRequired() || io.isReadOptional())
+    if (IOobjectOption::isReadRequired(readAverage))
+    {
+        is >> average_;
+        hasAverage_ = true;
+    }
+    else if (IOobjectOption::isReadOptional(readAverage))
+    {
+        // Slightly heavy handed
+        const bool oldThrowingIOerr = FatalIOError.throwing(true);
+
+        try
+        {
+            is >> average_;
+            hasAverage_ = true;
+        }
+        catch (const Foam::IOerror& err)
+        {
+            average_ = Zero;
+            hasAverage_ = false;
+        }
+        FatalIOError.throwing(oldThrowingIOerr);
+    }
+}
+
+
+template<class Type>
+bool Foam::rawIOField<Type>::readContents
+(
+    IOobjectOption::readOption readAverage
+)
+{
+    if (isReadRequired() || isReadOptional())
     {
         bool haveFile = false;
         bool haveHeader = false;
 
-        // Replacement of regIOobject::headerok() since that one complains
+        // Replacement of regIOobject::headerOk() since that one complains
         // if there is no header. TBD - Move up to headerOk()/fileHandler.
         {
             const fileName fName(filePath());
@@ -66,7 +96,7 @@ Foam::rawIOField<Type>::rawIOField(const IOobject& io, const bool readAverage)
 
             if (debug)
             {
-                Pout<< "rawIOField : object:" << io.name()
+                Pout<< "rawIOField : object:" << name()
                     << " haveFile:" << haveFile
                     << " haveHeader:" << haveHeader << endl;
             }
@@ -81,35 +111,23 @@ Foam::rawIOField<Type>::rawIOField(const IOobject& io, const bool readAverage)
 
             if (is.good())
             {
-                is  >> static_cast<Field<Type>&>(*this);
-                if (readAverage)
-                {
-                    hasAverage_ = true;
-                    is >> average_;
-                }
+                readContents(is, readAverage);
                 close();
             }
         }
         else if (haveFile)
         {
             // Failed reading header - fall back to IFstream
-            autoPtr<ISstream> isPtr(fileHandler().NewIFstream(io.objectPath()));
+            autoPtr<ISstream> isPtr(fileHandler().NewIFstream(objectPath()));
 
             if (isPtr && isPtr->good())
             {
-                auto& is = *isPtr;
-
-                is  >> static_cast<Field<Type>&>(*this);
-                if (readAverage)
-                {
-                    hasAverage_ = true;
-                    is >> average_;
-                }
+                readContents(*isPtr, readAverage);
             }
             else
             {
                 // Error if required but missing
-                if (io.isReadRequired())
+                if (isReadRequired())
                 {
                     FatalIOErrorInFunction(*isPtr)
                         << "Trying to read raw field" << endl
@@ -120,23 +138,58 @@ Foam::rawIOField<Type>::rawIOField(const IOobject& io, const bool readAverage)
 
         if (debug)
         {
-            Pout<< "rawIOField : object:" << io.name()
+            Pout<< "rawIOField : object:" << name()
                 << " size:" << this->size() << endl;
         }
+
+        return true;
     }
+
+    return false;
 }
 
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::rawIOField<Type>::rawIOField
+(
+    const IOobject& io,
+    IOobjectOption::readOption readAverage
+)
+:
+    regIOobject(io),
+    hasAverage_(false),
+    average_(Zero)
+{
+    // Check for MUST_READ_IF_MODIFIED
+    warnNoRereading<rawIOField<Type>>();
+
+    readContents(readAverage);
+}
+
+
+template<class Type>
+Foam::rawIOField<Type>::rawIOField
+(
+    const IOobject& io,
+    const bool readAverage
+)
+:
+    rawIOField<Type>
+    (
+        io,
+        (
+            readAverage
+          ? IOobjectOption::readOption::MUST_READ
+          : IOobjectOption::readOption::NO_READ
+        )
+    )
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Type>
-void Foam::rawIOField<Type>::setAverage(const Type& val)
-{
-    hasAverage_ = true;
-    average_ = val;
-}
-
-
 template<class Type>
 bool Foam::rawIOField<Type>::writeData(Ostream& os) const
 {
diff --git a/src/OpenFOAM/db/IOobjects/rawIOField/rawIOField.H b/src/OpenFOAM/db/IOobjects/rawIOField/rawIOField.H
index df71f7ca65..3dfc3d0c6d 100644
--- a/src/OpenFOAM/db/IOobjects/rawIOField/rawIOField.H
+++ b/src/OpenFOAM/db/IOobjects/rawIOField/rawIOField.H
@@ -64,8 +64,21 @@ class rawIOField
         Type average_;
 
 
+    // Private Member Functions
+
+        //- Read contents and average
+        void readContents(Istream&, IOobjectOption::readOption readAverage);
+
+        //- Read if IOobject flags set. Return true if read.
+        bool readContents(IOobjectOption::readOption readAverage);
+
+
 public:
 
+    //- The underlying content type
+    typedef Field<Type> content_type;
+
+    //- Runtime type information
     TypeName("rawField");
 
 
@@ -74,8 +87,16 @@ public:
         //- Default copy construct
         rawIOField(const rawIOField&) = default;
 
-        //- Construct from IOobject
-        explicit rawIOField(const IOobject& io, const bool readAverage=false);
+        //- Construct from IOobject and specified request for reading average
+        //- (default: NO_READ)
+        explicit rawIOField
+        (
+            const IOobject& io,
+            IOobjectOption::readOption readAverage = IOobjectOption::NO_READ
+        );
+
+        //- Construct from IOobject and specified request for reading average
+        rawIOField(const IOobject& io, const bool readAverage);
 
 
     //- Destructor
@@ -84,10 +105,19 @@ public:
 
     // Member Functions
 
-        //- Has an average value
-        bool hasAverage() const noexcept
+        //- The writeData method for regIOobject write operation
+        virtual bool writeData(Ostream& os) const;
+
+        //- The field content
+        const Field<Type>& field() const noexcept
         {
-            return hasAverage_;
+            return static_cast<const Field<Type>&>(*this);
+        }
+
+        //- The field content, non-const access
+        Field<Type>& field() noexcept
+        {
+            return static_cast<Field<Type>&>(*this);
         }
 
         //- The average value (if any)
@@ -96,10 +126,25 @@ public:
             return average_;
         }
 
-        //- Set an average value
-        void setAverage(const Type& val);
+        //- Has an average value
+        bool hasAverage() const noexcept
+        {
+            return hasAverage_;
+        }
 
-        bool writeData(Ostream& os) const;
+        //- Set an average value
+        void setAverage(const Type& val)
+        {
+            hasAverage_ = true;
+            average_ = val;
+        }
+
+        //- Clear average value
+        void clearAverage()
+        {
+            hasAverage_ = false;
+            average_ = Zero;
+        }
 
 
     // Member Operators
diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C
index c98b89a3e2..c24a6b764e 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C
@@ -48,12 +48,9 @@ timeVaryingMappedFixedValuePointPatchField
     mapMethod_(),
     mapperPtr_(nullptr),
     sampleTimes_(),
-    begSampleIndex_(-1),
-    endSampleIndex_(-1),
-    begAverage_(Zero),
-    endAverage_(Zero),
-    begSampledValues_(),
-    endSampledValues_(),
+    sampleIndex_(-1, -1),
+    sampleAverage_(Zero, Zero),
+    sampleValues_(),
     offset_(nullptr)
 {}
 
@@ -75,12 +72,9 @@ timeVaryingMappedFixedValuePointPatchField
     mapMethod_(),
     mapperPtr_(nullptr),
     sampleTimes_(),
-    begSampleIndex_(-1),
-    endSampleIndex_(-1),
-    begAverage_(Zero),
-    endAverage_(Zero),
-    begSampledValues_(),
-    endSampledValues_(),
+    sampleIndex_(-1, -1),
+    sampleAverage_(Zero, Zero),
+    sampleValues_(),
     offset_
     (
         Function1<Type>::NewIfPresent("offset", dict, word::null, &this->db())
@@ -143,12 +137,9 @@ timeVaryingMappedFixedValuePointPatchField
     mapMethod_(ptf.mapMethod_),
     mapperPtr_(nullptr),
     sampleTimes_(),
-    begSampleIndex_(-1),
-    endSampleIndex_(-1),
-    begAverage_(Zero),
-    endAverage_(Zero),
-    begSampledValues_(),
-    endSampledValues_(),
+    sampleIndex_(-1, -1),
+    sampleAverage_(Zero, Zero),
+    sampleValues_(),
     offset_(ptf.offset_.clone())
 {}
 
@@ -168,12 +159,9 @@ timeVaryingMappedFixedValuePointPatchField
     mapMethod_(ptf.mapMethod_),
     mapperPtr_(ptf.mapperPtr_),
     sampleTimes_(ptf.sampleTimes_),
-    begSampleIndex_(ptf.begSampleIndex_),
-    endSampleIndex_(ptf.endSampleIndex_),
-    begAverage_(ptf.begAverage_),
-    endAverage_(ptf.endAverage_),
-    begSampledValues_(ptf.begSampledValues_),
-    endSampledValues_(ptf.endSampledValues_),
+    sampleIndex_(ptf.sampleIndex_),
+    sampleAverage_(ptf.sampleAverage_),
+    sampleValues_(ptf.sampleValues_),
     offset_(ptf.offset_.clone())
 {}
 
@@ -194,12 +182,9 @@ timeVaryingMappedFixedValuePointPatchField
     mapMethod_(ptf.mapMethod_),
     mapperPtr_(ptf.mapperPtr_),
     sampleTimes_(ptf.sampleTimes_),
-    begSampleIndex_(ptf.begSampleIndex_),
-    endSampleIndex_(ptf.endSampleIndex_),
-    begAverage_(ptf.begAverage_),
-    endAverage_(ptf.endAverage_),
-    begSampledValues_(ptf.begSampledValues_),
-    endSampledValues_(ptf.endSampledValues_),
+    sampleIndex_(ptf.sampleIndex_),
+    sampleAverage_(ptf.sampleAverage_),
+    sampleValues_(ptf.sampleValues_),
     offset_(ptf.offset_.clone())
 {}
 
@@ -214,20 +199,19 @@ void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::autoMap
 {
     fixedValuePointPatchField<Type>::autoMap(m);
 
-    if (begSampledValues_.size())
+    if (sampleValues_.first().size())
     {
-        begSampledValues_.autoMap(m);
+        sampleValues_.first().autoMap(m);
     }
 
-    if (endSampledValues_.size())
+    if (sampleValues_.second().size())
     {
-        endSampledValues_.autoMap(m);
+        sampleValues_.second().autoMap(m);
     }
 
     // Clear interpolator
     mapperPtr_.reset(nullptr);
-    begSampleIndex_ = -1;
-    endSampleIndex_ = -1;
+    sampleIndex_ = labelPair(-1, -1);
 }
 
 
@@ -243,87 +227,82 @@ void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::rmap
     const timeVaryingMappedFixedValuePointPatchField<Type>& tiptf =
         refCast<const timeVaryingMappedFixedValuePointPatchField<Type>>(ptf);
 
-    begSampledValues_.rmap(tiptf.begSampledValues_, addr);
-    endSampledValues_.rmap(tiptf.endSampledValues_, addr);
+    sampleValues_.first().rmap(tiptf.sampleValues_.first(), addr);
+    sampleValues_.second().rmap(tiptf.sampleValues_.second(), addr);
 
     // Clear interpolator
     mapperPtr_.reset(nullptr);
-    begSampleIndex_ = -1;
-    endSampleIndex_ = -1;
+    sampleIndex_ = labelPair(-1, -1);
 }
 
 
 template<class Type>
 void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::updateSampledValues
 (
-    const int whichEnd  // (0|1)
-)
+    const label sampleIndex,
+    Field<Type>& field,
+    Type& avg
+) const
 {
+    tmp<Field<Type>> tvalues;
+
     // Update sampled data fields
-    const Time& time = this->db().time();
-
-    const word& sampleTimeName =
-        sampleTimes_[(whichEnd ? endSampleIndex_ : begSampleIndex_)].name();
-
-    if (debug)
     {
-        Pout<< "checkTable : Reading values from "
-            <<
-            (
-                "boundaryData"
-              / this->patch().name()
-              / sampleTimeName
-              / fieldTableName_
-            ) << endl;
-    }
+        const Time& time = this->db().time();
 
-    // Reread values and interpolate
-    const fileName valsFile
-    (
-        time.caseConstant()
-        /"boundaryData"
-        /this->patch().name()
-        /sampleTimeName
-        /fieldTableName_
-    );
-
-    IOobject io
-    (
-        valsFile,   // absolute path
-        time,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE,
-        false,              // no need to register
-        true                // is global object (currently not used)
-    );
-
-    const rawIOField<Type> vals(io, setAverage_);
-
-    if (vals.size() != mapperPtr_().sourceSize())
-    {
-        FatalErrorInFunction
-            << "Number of values (" << vals.size()
-            << ") differs from the number of points ("
-            <<  mapperPtr_().sourceSize()
-            << ") in file " << valsFile << exit(FatalError);
-    }
-
-    if (whichEnd)
-    {
-        if (setAverage_)  // or vals.hasAverage()
+        if (debug)
         {
-            endAverage_ = vals.average();
+            Pout<< "checkTable : Reading values from "
+                <<
+                (
+                    "boundaryData"
+                  / this->patch().name()
+                  / sampleTimes_[sampleIndex].name()
+                  / fieldTableName_
+                ) << endl;
         }
-        endSampledValues_ = mapperPtr_().interpolate(vals);
-    }
-    else
-    {
-        if (setAverage_)  // or vals.hasAverage()
+
+        // Reread values and interpolate
+        const fileName valsFile
+        (
+            time.caseConstant()
+            /"boundaryData"
+            /this->patch().name()
+            /sampleTimes_[sampleIndex].name()
+            /fieldTableName_
+        );
+
+        IOobject io
+        (
+            valsFile,               // absolute path
+            time,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            IOobject::NO_REGISTER,
+            true                    // is global object (currently not used)
+        );
+
+        rawIOField<Type> vals(io, setAverage_);
+
+        if (vals.hasAverage())
         {
-            begAverage_ = vals.average();
+            avg = vals.average();
         }
-        begSampledValues_ = mapperPtr_().interpolate(vals);
+
+        if (vals.size() != mapperPtr_().sourceSize())
+        {
+            FatalErrorInFunction
+                << "Number of values (" << vals.size()
+                << ") differs from the number of points ("
+                <<  mapperPtr_().sourceSize()
+                << ") in file " << valsFile << exit(FatalError);
+        }
+
+        tvalues = tmp<Field<Type>>::New(std::move(vals.field()));
     }
+
+    // From input values to interpolated (sampled) positions
+    field = mapperPtr_().interpolate(tvalues);
 }
 
 
@@ -336,7 +315,7 @@ void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::checkTable
    const Time& time = this->db().time();
 
     // Initialise
-    if (begSampleIndex_ == -1 && endSampleIndex_ == -1)
+    if (sampleIndex_.first() == -1 && sampleIndex_.second() == -1)
     {
         const polyMesh& pMesh = this->patch().boundaryMesh().mesh()();
 
@@ -350,11 +329,9 @@ void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::checkTable
         else
         {
             // Load points from facesInstance
-            if (debug)
-            {
-                Info<< "Reloading points0 from " << pMesh.facesInstance()
-                    << endl;
-            }
+            DebugInFunction
+                << "Reloading points0 from " << pMesh.facesInstance()
+                << endl;
 
             pointIOField points0
             (
@@ -366,7 +343,7 @@ void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::checkTable
                     pMesh,
                     IOobject::MUST_READ,
                     IOobject::NO_WRITE,
-                    false
+                    IOobject::NO_REGISTER
                 )
             );
             meshPts = pointField(points0, this->patch().meshPoints());
@@ -387,12 +364,12 @@ void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::checkTable
             time,
             IOobject::MUST_READ,
             IOobject::NO_WRITE,
-            false,              // no need to register
+            IOobject::NO_REGISTER,
             true                // is global object (currently not used)
         );
 
         // Read data (no average value!)
-        const rawIOField<point> samplePoints(io, false);
+        const rawIOField<point> samplePoints(io);
 
         // tbd: run-time selection
         const bool nearestOnly =
@@ -417,9 +394,8 @@ void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::checkTable
         const fileName samplePointsDir = samplePointsFile.path();
         sampleTimes_ = Time::findTimes(samplePointsDir);
 
-        DebugInfo
-            << "timeVaryingMappedFixedValuePointPatchField : In directory "
-            << samplePointsDir << " found times "
+        DebugInFunction
+            << "Found times "
             << pointToPointPlanarInterpolation::timeNames(sampleTimes_)
             << endl;
     }
@@ -429,7 +405,7 @@ void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::checkTable
     (
         sampleTimes_,
         t,  // time.value(),
-        begSampleIndex_
+        sampleIndex_.first()
     );
 
     if (timeIndices.first() < 0)
@@ -449,48 +425,46 @@ void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::checkTable
 
     // Update sampled data fields.
 
-    if (begSampleIndex_ != timeIndices.first())
+    if (sampleIndex_.first() != timeIndices.first())
     {
-        begSampleIndex_ = timeIndices.first();
+        sampleIndex_.first() = timeIndices.first();
 
-        if (begSampleIndex_ == endSampleIndex_)
+        if (sampleIndex_.first() == sampleIndex_.second())
         {
-            // No need to reread since are end values
-            if (debug)
-            {
-                Pout<< "checkTable : Setting startValues to (already read) "
-                    <<   "boundaryData"
-                        /this->patch().name()
-                        /sampleTimes_[begSampleIndex_].name()
-                    << endl;
-            }
-            begAverage_ = endAverage_;
-            begSampledValues_ = endSampledValues_;
+            // Can reuse previous end values
+            sampleValues_.first() = sampleValues_.second();
+            sampleAverage_.first() = sampleAverage_.second();
         }
         else
         {
-            // Update begin values
-            this->updateSampledValues(0);
+            // Update first() values
+            this->updateSampledValues
+            (
+                sampleIndex_.first(),
+                sampleValues_.first(),
+                sampleAverage_.first()
+            );
         }
     }
 
-    if (endSampleIndex_ != timeIndices.second())
+    if (sampleIndex_.second() != timeIndices.second())
     {
-        endSampleIndex_ = timeIndices.second();
+        sampleIndex_.second() = timeIndices.second();
 
-        if (endSampleIndex_ == -1)
+        if (sampleIndex_.second() == -1)
         {
-            // endTime no longer valid. Might as well clear endValues.
-            if (debug)
-            {
-                Pout<< "checkTable : Clearing endValues" << endl;
-            }
-            endSampledValues_.clear();
+            // Index no longer valid - clear values accordingly
+            sampleValues_.second().clear();
         }
         else
         {
-            // Update end values
-            this->updateSampledValues(1);
+            // Update second() values
+            this->updateSampledValues
+            (
+                sampleIndex_.second(),
+                sampleValues_.second(),
+                sampleAverage_.second()
+            );
         }
     }
 }
@@ -513,36 +487,29 @@ void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::updateCoeffs()
     auto& fld = static_cast<Field<Type>&>(*this);
     Type wantedAverage;
 
-    if (endSampleIndex_ == -1)
+    if (sampleIndex_.second() == -1)
     {
         // Only start value
-        if (debug)
-        {
-            Pout<< "updateCoeffs : Sampled, non-interpolated values"
-                << " from start time:"
-                << sampleTimes_[begSampleIndex_].name() << nl;
-        }
-
-        fld = begSampledValues_;
-        wantedAverage = begAverage_;
+        fld = sampleValues_.first();
+        wantedAverage = sampleAverage_.first();
     }
     else
     {
-        const scalar beg = sampleTimes_[begSampleIndex_].value();
-        const scalar end = sampleTimes_[endSampleIndex_].value();
+        const scalar beg = sampleTimes_[sampleIndex_.first()].value();
+        const scalar end = sampleTimes_[sampleIndex_.second()].value();
         const scalar s = (x - beg)/(end - beg);
 
-        if (debug)
-        {
-            Pout<< "updateCoeffs : Sampled, interpolated values"
-                << " between start time:"
-                << sampleTimes_[begSampleIndex_].name()
-                << " and end time:" << sampleTimes_[endSampleIndex_].name()
-                << " with weight:" << s << endl;
-        }
+        fld = (1 - s)*sampleValues_.first() + s*sampleValues_.second();
 
-        fld = ((1 - s)*begSampledValues_ + s*endSampledValues_);
-        wantedAverage = (1 - s)*begAverage_ + s*endAverage_;
+        wantedAverage
+            = (1 - s)*sampleAverage_.first() + s*sampleAverage_.second();
+
+        DebugInFunction
+            << "Sampled, interpolated values"
+            << " between time:"
+            << sampleTimes_[sampleIndex_.first()].name()
+            << " and time:" << sampleTimes_[sampleIndex_.second()].name()
+            << " with weight:" << s << endl;
     }
 
     // Enforce average. Either by scaling (if scaling factor > 0.5) or by
diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H
index b2d146994a..2dedce3d91 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H
+++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H
@@ -53,9 +53,10 @@ SourceFiles
 #define Foam_timeVaryingMappedFixedValuePointPatchField_H
 
 #include "fixedValuePointPatchField.H"
+#include "Function1.H"
+#include "Pair.H"
 #include "instantList.H"
 #include "pointToPointPlanarInterpolation.H"
-#include "Function1.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -94,23 +95,14 @@ class timeVaryingMappedFixedValuePointPatchField
         //- List of boundaryData time directories
         instantList sampleTimes_;
 
-        //- Current start index in sampleTimes
-        label begSampleIndex_;
+        //- The first/second sample indices in sampleTimes
+        labelPair sampleIndex_;
 
-        //- Current end index in sampleTimes
-        label endSampleIndex_;
+        //- The average values at first/second sampleIndex_ (if setAverage)
+        Pair<Type> sampleAverage_;
 
-        //- Current average value at begSampleIndex (if setAverage)
-        Type begAverage_;
-
-        //- Current average value at endSampleIndex (if setAverage)
-        Type endAverage_;
-
-        //- Current interpolated values at begSampleIndex
-        Field<Type> begSampledValues_;
-
-        //- Current interpolated values at endSampleIndex
-        Field<Type> endSampledValues_;
+        //- The interpolated values at first/second sampleIndex_
+        Pair<Field<Type>> sampleValues_;
 
         //- Time varying offset values to interpolated data
         autoPtr<Function1<Type>> offset_;
@@ -118,12 +110,18 @@ class timeVaryingMappedFixedValuePointPatchField
 
     // Private Member Functions
 
-        //- Update start (0) or end (1) interpolated data
-        void updateSampledValues(const int whichEnd);
-
         //- Find boundary data between time 't' and interpolate
         void checkTable(const scalar t);
 
+        //- Update field and average value using interpolation at
+        //- given instant index
+        void updateSampledValues
+        (
+            const label sampleIndex,
+            Field<Type>& field,
+            Type& avg
+        ) const;
+
 
 public:
 
diff --git a/src/meshTools/PatchFunction1/MappedFile/MappedFile.C b/src/meshTools/PatchFunction1/MappedFile/MappedFile.C
index e3e8d4e681..00f118662e 100644
--- a/src/meshTools/PatchFunction1/MappedFile/MappedFile.C
+++ b/src/meshTools/PatchFunction1/MappedFile/MappedFile.C
@@ -49,12 +49,9 @@ Foam::PatchFunction1Types::MappedFile<Type>::MappedFile
     mapMethod_(),
     mapperPtr_(nullptr),
     sampleTimes_(),
-    begSampleIndex_(-1),
-    endSampleIndex_(-1),
-    begAverage_(Zero),
-    endAverage_(Zero),
-    begSampledValues_(),
-    endSampledValues_(),
+    sampleIndex_(-1, -1),
+    sampleAverage_(Zero, Zero),
+    sampleValues_(),
     offset_(Function1<Type>::NewIfPresent("offset", dict))
 {
     if
@@ -93,12 +90,9 @@ Foam::PatchFunction1Types::MappedFile<Type>::MappedFile
     mapMethod_(),
     mapperPtr_(nullptr),
     sampleTimes_(),
-    begSampleIndex_(-1),
-    endSampleIndex_(-1),
-    begAverage_(Zero),
-    endAverage_(Zero),
-    begSampledValues_(),
-    endSampledValues_(),
+    sampleIndex_(-1, -1),
+    sampleAverage_(Zero, Zero),
+    sampleValues_(),
     offset_(Function1<Type>::NewIfPresent("offset", dict))
 {
     if
@@ -144,12 +138,9 @@ Foam::PatchFunction1Types::MappedFile<Type>::MappedFile
     mapMethod_(rhs.mapMethod_),
     mapperPtr_(rhs.mapperPtr_.clone()),
     sampleTimes_(rhs.sampleTimes_),
-    begSampleIndex_(rhs.begSampleIndex_),
-    endSampleIndex_(rhs.endSampleIndex_),
-    begAverage_(rhs.begAverage_),
-    endAverage_(rhs.endAverage_),
-    begSampledValues_(rhs.begSampledValues_),
-    endSampledValues_(rhs.endSampledValues_),
+    sampleIndex_(rhs.sampleIndex_),
+    sampleAverage_(rhs.sampleAverage_),
+    sampleValues_(rhs.sampleValues_),
     offset_(rhs.offset_.clone())
 {}
 
@@ -164,19 +155,18 @@ void Foam::PatchFunction1Types::MappedFile<Type>::autoMap
 {
     PatchFunction1<Type>::autoMap(mapper);
 
-    if (begSampledValues_.size())
+    if (sampleValues_.first().size())
     {
-        begSampledValues_.autoMap(mapper);
+        sampleValues_.first().autoMap(mapper);
     }
-    if (endSampledValues_.size())
+    if (sampleValues_.second().size())
     {
-        endSampledValues_.autoMap(mapper);
+        sampleValues_.second().autoMap(mapper);
     }
 
     // Clear interpolator
     mapperPtr_.reset(nullptr);
-    begSampleIndex_ = -1;
-    endSampleIndex_ = -1;
+    sampleIndex_ = labelPair(-1, -1);
 }
 
 
@@ -189,103 +179,99 @@ void Foam::PatchFunction1Types::MappedFile<Type>::rmap
 {
     PatchFunction1<Type>::rmap(pf1, addr);
 
-    const PatchFunction1Types::MappedFile<Type>& tiptf =
+    const auto& tiptf =
         refCast<const PatchFunction1Types::MappedFile<Type>>(pf1);
 
-    if (tiptf.begSampledValues_.size())
+    if (tiptf.sampleValues_.first().size())
     {
-        begSampledValues_.resize(this->size());
-        begSampledValues_.rmap(tiptf.begSampledValues_, addr);
+        sampleValues_.first().resize(this->size());
+        sampleValues_.first().rmap(tiptf.sampleValues_.first(), addr);
     }
 
-    if (tiptf.endSampledValues_.size())
+    if (tiptf.sampleValues_.second().size())
     {
-        endSampledValues_.resize(this->size());
-        endSampledValues_.rmap(tiptf.endSampledValues_, addr);
+        sampleValues_.second().resize(this->size());
+        sampleValues_.second().rmap(tiptf.sampleValues_.second(), addr);
     }
 
     // Clear interpolator
     mapperPtr_.reset(nullptr);
-    begSampleIndex_ = -1;
-    endSampleIndex_ = -1;
+    sampleIndex_ = labelPair(-1, -1);
 }
 
 
 template<class Type>
 void Foam::PatchFunction1Types::MappedFile<Type>::updateSampledValues
 (
-    const int whichEnd  // (0|1)
+    const label sampleIndex,
+    Field<Type>& field,
+    Type& avg
 ) const
 {
+    tmp<Field<Type>> tvalues;
+
     // Update sampled data fields
-    const polyMesh& mesh = this->patch_.boundaryMesh().mesh();
-    const Time& time = mesh.time();
-
-    const word& sampleTimeName =
-        sampleTimes_[(whichEnd ? endSampleIndex_ : begSampleIndex_)].name();
-
-    if (debug)
     {
-        Pout<< "checkTable : Reading values from "
-            <<
-            (
-                "boundaryData"
-              / this->patch_.name()
-              / sampleTimeName
-              / fieldTableName_
-            ) << endl;
-    }
+        const polyMesh& mesh = this->patch_.boundaryMesh().mesh();
+        const Time& time = mesh.time();
 
-    // Reread values and interpolate
-    const fileName valsFile
-    (
-        time.globalPath()
-        /time.constant()
-        /mesh.dbDir()            // region
-        /"boundaryData"
-        /this->patch_.name()
-        /sampleTimeName
-        /fieldTableName_
-    );
-
-    IOobject io
-    (
-        valsFile,   // absolute path
-        time,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE,
-        false,              // no need to register
-        true                // is global object (currently not used)
-    );
-
-    const rawIOField<Type> vals(io, setAverage_);
-
-    if (vals.size() != mapperPtr_().sourceSize())
-    {
-        FatalErrorInFunction
-            << "Number of values (" << vals.size()
-            << ") differs from the number of points ("
-            <<  mapperPtr_().sourceSize()
-            << ") in file " << valsFile
-            << exit(FatalError);
-    }
-
-    if (whichEnd)
-    {
-        if (setAverage_)  // or vals.hasAverage()
+        if (debug)
         {
-            endAverage_ = vals.average();
+            Pout<< "checkTable : Reading values from "
+                <<
+                (
+                    "boundaryData"
+                  / this->patch_.name()
+                  / sampleTimes_[sampleIndex].name()
+                  / fieldTableName_
+                ) << endl;
         }
-        endSampledValues_ = mapperPtr_().interpolate(vals);
-    }
-    else
-    {
-        if (setAverage_)  // or vals.hasAverage()
+
+        // Reread values and interpolate
+        const fileName valsFile
+        (
+            time.globalPath()
+            /time.constant()
+            /mesh.dbDir()            // region
+            /"boundaryData"
+            /this->patch_.name()
+            /sampleTimes_[sampleIndex].name()
+            /fieldTableName_
+        );
+
+        IOobject io
+        (
+            valsFile,           // absolute path
+            time,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            IOobject::NO_REGISTER,
+            true                // is global object (currently not used)
+        );
+
+        rawIOField<Type> vals(io, setAverage_);
+
+        if (vals.hasAverage())
         {
-            begAverage_ = vals.average();
+            avg = vals.average();
         }
-        begSampledValues_ = mapperPtr_().interpolate(vals);
+
+        if (vals.size() != mapperPtr_().sourceSize())
+        {
+            FatalErrorInFunction
+                << "Number of values (" << vals.size()
+                << ") differs from the number of points ("
+                <<  mapperPtr_().sourceSize()
+                << ") in file " << valsFile
+                << exit(FatalError);
+        }
+
+        tvalues = tmp<Field<Type>>::New(std::move(vals.field()));
     }
+
+
+    // From input values to interpolated (sampled) positions
+    field = mapperPtr_().interpolate(tvalues);
 }
 
 
@@ -318,12 +304,12 @@ void Foam::PatchFunction1Types::MappedFile<Type>::checkTable
             time,
             IOobject::MUST_READ,
             IOobject::NO_WRITE,
-            false,              // no need to register
+            IOobject::NO_REGISTER,
             true                // is global object (currently not used)
         );
 
         // Read data (no average value!)
-        const rawIOField<point> samplePoints(io, false);
+        const rawIOField<point> samplePoints(io);
 
         // tbd: run-time selection
         const bool nearestOnly =
@@ -369,8 +355,7 @@ void Foam::PatchFunction1Types::MappedFile<Type>::checkTable
         sampleTimes_ = Time::findTimes(samplePointsDir);
 
         DebugInfo
-            << "In directory "
-            << samplePointsDir << " found times "
+            << "Found times "
             << pointToPointPlanarInterpolation::timeNames(sampleTimes_)
             << endl;
     }
@@ -381,7 +366,7 @@ void Foam::PatchFunction1Types::MappedFile<Type>::checkTable
     (
         sampleTimes_,
         t,  //mesh.time().value(),
-        begSampleIndex_
+        sampleIndex_.first()
     );
 
     if (timeIndices.first() < 0)
@@ -401,48 +386,46 @@ void Foam::PatchFunction1Types::MappedFile<Type>::checkTable
 
     // Update sampled data fields.
 
-    if (begSampleIndex_ != timeIndices.first())
+    if (sampleIndex_.first() != timeIndices.first())
     {
-        begSampleIndex_ = timeIndices.first();
+        sampleIndex_.first() = timeIndices.first();
 
-        if (begSampleIndex_ == endSampleIndex_)
+        if (sampleIndex_.first() == sampleIndex_.second())
         {
-            // No need to reread since are end values
-            if (debug)
-            {
-                Pout<< "checkTable : Setting startValues to (already read) "
-                    << "boundaryData"
-                      /this->patch_.name()
-                      /sampleTimes_[begSampleIndex_].name()
-                    << endl;
-            }
-            begAverage_ = endAverage_;
-            begSampledValues_ = endSampledValues_;
+            // Can reuse previous end values
+            sampleValues_.first() = sampleValues_.second();
+            sampleAverage_.first() = sampleAverage_.second();
         }
         else
         {
-            // Update begin values
-            this->updateSampledValues(0);
+            // Update first() values
+            this->updateSampledValues
+            (
+                sampleIndex_.first(),
+                sampleValues_.first(),
+                sampleAverage_.first()
+            );
         }
     }
 
-    if (endSampleIndex_ != timeIndices.second())
+    if (sampleIndex_.second() != timeIndices.second())
     {
-        endSampleIndex_ = timeIndices.second();
+        sampleIndex_.second() = timeIndices.second();
 
-        if (endSampleIndex_ == -1)
+        if (sampleIndex_.second() == -1)
         {
-            // endTime no longer valid. Might as well clear endValues.
-            if (debug)
-            {
-                Pout<< "checkTable : Clearing endValues" << endl;
-            }
-            endSampledValues_.clear();
+            // Index no longer valid - clear values accordingly
+            sampleValues_.second().clear();
         }
         else
         {
-            // Update end values
-            this->updateSampledValues(1);
+            // Update second() values
+            this->updateSampledValues
+            (
+                sampleIndex_.second(),
+                sampleValues_.second(),
+                sampleAverage_.second()
+            );
         }
     }
 }
@@ -457,40 +440,35 @@ Foam::PatchFunction1Types::MappedFile<Type>::value
 {
     checkTable(x);
 
-    auto tfld = tmp<Field<Type>>::New(begSampledValues_.size());
+    // Interpolate between the sampled data
+    auto tfld = tmp<Field<Type>>::New();
     auto& fld = tfld.ref();
     Type wantedAverage;
 
-    if (endSampleIndex_ == -1)
+    if (sampleIndex_.second() == -1)
     {
         // Only start value
-        if (debug)
-        {
-            Pout<< "MappedFile<Type>::value : Sampled, non-interpolated values"
-                << " from start time:"
-                << sampleTimes_[begSampleIndex_].name() << nl;
-        }
-
-        fld = begSampledValues_;
-        wantedAverage = begAverage_;
+        fld = sampleValues_.first();
+        wantedAverage = sampleAverage_.first();
     }
     else
     {
-        const scalar beg = sampleTimes_[begSampleIndex_].value();
-        const scalar end = sampleTimes_[endSampleIndex_].value();
+        const scalar beg = sampleTimes_[sampleIndex_.first()].value();
+        const scalar end = sampleTimes_[sampleIndex_.second()].value();
         const scalar s = (x - beg)/(end - beg);
 
-        if (debug)
-        {
-            Pout<< "MappedFile<Type>::value : Sampled, interpolated values"
-                << " between start time:"
-                << sampleTimes_[begSampleIndex_].name()
-                << " and end time:" << sampleTimes_[endSampleIndex_].name()
-                << " with weight:" << s << endl;
-        }
+        fld = (1 - s)*sampleValues_.first() + s*sampleValues_.second();
 
-        fld = ((1 - s)*begSampledValues_ + s*endSampledValues_);
-        wantedAverage = (1 - s)*begAverage_ + s*endAverage_;
+        wantedAverage
+            = (1 - s)*sampleAverage_.first() + s*sampleAverage_.second();
+
+        DebugInfo
+            << "MappedFile<Type>::value : "
+            << "Sampled, interpolated values"
+            << " between time:"
+            << sampleTimes_[sampleIndex_.first()].name()
+            << " and time:" << sampleTimes_[sampleIndex_.second()].name()
+            << " with weight:" << s << endl;
     }
 
     // Enforce average. Either by scaling (if scaling factor > 0.5) or by
@@ -498,6 +476,7 @@ Foam::PatchFunction1Types::MappedFile<Type>::value
     if (setAverage_)
     {
         Type averagePsi;
+
         if (this->faceValues())
         {
             const scalarField magSf(mag(this->patch_.faceAreas()));
diff --git a/src/meshTools/PatchFunction1/MappedFile/MappedFile.H b/src/meshTools/PatchFunction1/MappedFile/MappedFile.H
index 1948af097a..b45a7cca8e 100644
--- a/src/meshTools/PatchFunction1/MappedFile/MappedFile.H
+++ b/src/meshTools/PatchFunction1/MappedFile/MappedFile.H
@@ -64,6 +64,8 @@ SourceFiles
 
 #include "PatchFunction1.H"
 #include "Function1.H"
+#include "Pair.H"
+#include "instantList.H"
 #include "pointToPointPlanarInterpolation.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -108,23 +110,14 @@ class MappedFile
         //- List of boundaryData time directories
         mutable instantList sampleTimes_;
 
-        //- Current start index in sampleTimes
-        mutable label begSampleIndex_;
+        //- The first/second sample indices in sampleTimes
+        mutable labelPair sampleIndex_;
 
-        //- Current end index in sampleTimes
-        mutable label endSampleIndex_;
+        //- The average values at first/second sampleIndex_ (if setAverage)
+        mutable Pair<Type> sampleAverage_;
 
-        //- Current average value at begSampleIndex (if setAverage)
-        mutable Type begAverage_;
-
-        //- Current average value at endSampleIndex (if setAverage)
-        mutable Type endAverage_;
-
-        //- Current interpolated values at begSampleIndex
-        mutable Field<Type> begSampledValues_;
-
-        //- Current interpolated values at endSampleIndex
-        mutable Field<Type> endSampledValues_;
+        //- The interpolated values at first/second sampleIndex_
+        mutable Pair<Field<Type>> sampleValues_;
 
         //- Time varying offset values to interpolated data
         autoPtr<Function1<Type>> offset_;
@@ -132,12 +125,18 @@ class MappedFile
 
     // Private Member Functions
 
-        //- Update start (0) or end (1) interpolated data
-        void updateSampledValues(const int whichEnd) const;
-
         //- Find boundary data between time 't' and interpolate
         void checkTable(const scalar t) const;
 
+        //- Update field and average value using interpolation at
+        //- given instant index
+        void updateSampledValues
+        (
+            const label sampleIndex,
+            Field<Type>& field,
+            Type& avg
+        ) const;
+
 public:
 
     //- Runtime type information
diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H
index 80842319d2..93a297d0c6 100644
--- a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H
+++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H
@@ -196,6 +196,10 @@ public:
         //- Interpolate from field on source points to dest points
         template<class Type>
         tmp<Field<Type>> interpolate(const Field<Type>& sourceFld) const;
+
+        //- Interpolate from field on source points to dest points
+        template<class Type>
+        tmp<Field<Type>> interpolate(const tmp<Field<Type>>& tsource) const;
 };
 
 
diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolationTemplates.C b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolationTemplates.C
index f571becd8e..aae7dcfa3d 100644
--- a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolationTemplates.C
+++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolationTemplates.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -81,4 +82,18 @@ Foam::pointToPointPlanarInterpolation::interpolate
 }
 
 
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::pointToPointPlanarInterpolation::interpolate
+(
+    const tmp<Field<Type>>& tsource
+) const
+{
+    tmp<Field<Type>> tresult = this->interpolate(tsource());
+    tsource.clear();
+
+    return tresult;
+}
+
+
 // ************************************************************************* //
diff --git a/src/surfMesh/writers/boundary/boundaryDataSurfaceWriter.C b/src/surfMesh/writers/boundary/boundaryDataSurfaceWriter.C
index 0746958e26..ff76eccbcb 100644
--- a/src/surfMesh/writers/boundary/boundaryDataSurfaceWriter.C
+++ b/src/surfMesh/writers/boundary/boundaryDataSurfaceWriter.C
@@ -171,7 +171,7 @@ void Foam::surfaceWriters::boundaryDataWriter::serialWriteGeometry
                 iopts.db(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
-                false  // Do not register
+                IOobject::NO_REGISTER
             )
         );
         iofld.note() = "face data";
@@ -203,7 +203,7 @@ Foam::fileName Foam::surfaceWriters::boundaryDataWriter::write()
     fileName surfaceDir = outputPath_;
 
     // Dummy Time to use as objectRegistry
-    autoPtr<Time> dummyTimePtr(Time::New(argList::envGlobalPath()));
+    refPtr<Time> timePtr(Time::New(argList::envGlobalPath()));
 
     // const meshedSurf& surf = surface();
     const meshedSurfRef& surf = adjustSurface();
@@ -221,10 +221,10 @@ Foam::fileName Foam::surfaceWriters::boundaryDataWriter::write()
             IOobject
             (
                 surfaceDir/"points",
-                *dummyTimePtr,
+                *timePtr,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
-                false  // Do not register
+                IOobject::NO_REGISTER
             )
         );
         iopts.note() = (this->isPointData() ? "point data" : "face data");
@@ -265,7 +265,7 @@ Foam::fileName Foam::surfaceWriters::boundaryDataWriter::writeTemplate
 
 
     // Dummy Time to use as objectRegistry
-    autoPtr<Time> dummyTimePtr(Time::New(argList::envGlobalPath()));
+    refPtr<Time> timePtr(Time::New(argList::envGlobalPath()));
 
     // const meshedSurf& surf = surface();
     const meshedSurfRef& surf = adjustSurface();
@@ -284,10 +284,10 @@ Foam::fileName Foam::surfaceWriters::boundaryDataWriter::writeTemplate
                 IOobject
                 (
                     surfaceDir/"points",
-                    *dummyTimePtr,
+                    *timePtr,
                     IOobject::NO_READ,
                     IOobject::NO_WRITE,
-                    false  // Do not register
+                    IOobject::NO_REGISTER
                 )
             );
             iopts.note() = (this->isPointData() ? "point data" : "face data");
@@ -302,10 +302,10 @@ Foam::fileName Foam::surfaceWriters::boundaryDataWriter::writeTemplate
                 IOobject
                 (
                     outputFile,
-                    *dummyTimePtr,
+                    *timePtr,
                     IOobject::NO_READ,
                     IOobject::NO_WRITE,
-                    false  // Do not register
+                    IOobject::NO_REGISTER
                 )
             );
             iofld.note() = (this->isPointData() ? "point data" : "face data");
diff --git a/src/surfMesh/writers/boundary/boundaryDataSurfaceWriter.H b/src/surfMesh/writers/boundary/boundaryDataSurfaceWriter.H
index f48d188aa6..8e42d89ba7 100644
--- a/src/surfMesh/writers/boundary/boundaryDataSurfaceWriter.H
+++ b/src/surfMesh/writers/boundary/boundaryDataSurfaceWriter.H
@@ -36,9 +36,10 @@ Description
     {
         boundaryData
         {
-            header      false;
+            header      true;
             format      ascii;
             compression false;
+            normal      false;
         }
     }
     \endverbatim