diff --git a/applications/utilities/preProcessing/smoothSurfaceData/Make/files b/applications/utilities/preProcessing/smoothSurfaceData/Make/files
new file mode 100644
index 0000000000..cb031ed9d6
--- /dev/null
+++ b/applications/utilities/preProcessing/smoothSurfaceData/Make/files
@@ -0,0 +1,3 @@
+smoothSurfaceData.C
+
+EXE = $(FOAM_USER_APPBIN)/smoothSurfaceData
diff --git a/applications/utilities/preProcessing/smoothSurfaceData/Make/options b/applications/utilities/preProcessing/smoothSurfaceData/Make/options
new file mode 100644
index 0000000000..128a319f08
--- /dev/null
+++ b/applications/utilities/preProcessing/smoothSurfaceData/Make/options
@@ -0,0 +1,11 @@
+/* Can compile with -fopenmp etc */
+
+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/utilities/preProcessing/smoothSurfaceData/smoothSurfaceData.C b/applications/utilities/preProcessing/smoothSurfaceData/smoothSurfaceData.C
new file mode 100644
index 0000000000..7ca090ed43
--- /dev/null
+++ b/applications/utilities/preProcessing/smoothSurfaceData/smoothSurfaceData.C
@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    smoothSurfaceData
+
+Description
+    Pre-processing, filtering of surface field data.
+    Currently this is just used to test filter settings
+    (for MappedFile) and only generates vtk output, which can be
+    easily loaded in ParaView.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "clockTime.H"
+#include "primitiveFields.H"
+#include "surfaceReader.H"
+#include "surfaceWriter.H"
+#include "foamVtkSurfaceWriter.H"
+#include "MappedFileFilterField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+using namespace Foam;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Testing, pre-processing, filtering of surface field data"
+    );
+
+    argList::noCheckProcessorDirectories();
+    argList::addVerboseOption("Additional verbosity");
+
+    argList::addArgument("input", "The input surface file");
+
+    argList::addOption("radius", "m", "Specify filter radius [metres]");
+    argList::addOption("sweeps", "N", "Number of median filter stages");
+    argList::addOption
+    (
+        "field",
+        "name",
+        "Field <scalar> to process (default: T)"
+    );
+
+    argList::addOption
+    (
+        "read-format",
+        "type",
+        "Input format (default: ensight)"
+    );
+
+    argList args(argc, argv);
+    // #include "setRootCase.H"
+
+    const int optVerbose = args.verbose();
+
+    const word readFileType
+    (
+        args.getOrDefault<word>("read-format", "ensight")
+    );
+
+    // Constant radius searching
+    label filterSweeps_(1);
+    scalar filterRadius_(0);
+    args.readIfPresent("sweeps", filterSweeps_);
+    args.readIfPresent("radius", filterRadius_);
+
+    Info<< nl
+        << "Filter: radius=" << filterRadius_
+        << " sweeps=" << filterSweeps_ << endl;
+
+    // Simple sanity check
+    if ((filterSweeps_ < 1) || (filterRadius_ <= VSMALL))
+    {
+        Info<< nl << "Nothing to do. Exiting..." << nl << endl;
+        return 0;
+    }
+
+    word fieldName("T");
+    args.readIfPresent("field", fieldName);
+
+
+    const fileName importName = args.get<fileName>(1);
+
+    auto readerPtr_ = surfaceReader::New(readFileType, importName);
+
+    auto& reader = readerPtr_();
+
+    const label fieldIndex = reader.fieldNames(0).find(fieldName);
+    if (fieldIndex == -1)
+    {
+        FatalErrorInFunction
+            << "Unable to find field name: " << fieldName
+            << " in list of available fields: " << reader.fieldNames(0)
+            << exit(FatalError);
+    }
+
+    clockTime timing;
+
+    const meshedSurface& geom = reader.geometry(0);
+
+    Info<< nl << "Read " << geom.nFaces() << " faces and "
+        << geom.nPoints() << " points in "
+        << timing.timeIncrement() << "s" << endl;
+
+    PatchFunction1Types::FilterField fieldFilter;
+
+    PatchFunction1Types::FilterField::debug = optVerbose;
+
+    fieldFilter.reset(geom, filterRadius_);
+
+    Info<< nl << "Built weights/addressing "
+        << timing.timeIncrement() << "s" << endl;
+
+    Info<< nl << "Processing " << reader.times().size() << " times" << nl;
+
+    const instantList times(reader.times());
+
+    forAll(times, timeIndex)
+    {
+        tmp<scalarField> tfield
+        (
+            reader.field(timeIndex, fieldIndex, pTraits<scalar>::zero)
+        );
+
+        tfield = fieldFilter.evaluate(tfield, filterSweeps_);
+
+        vtk::surfaceWriter writer
+        (
+            geom.points(),
+            geom,
+            word::printf("filtered_%06d", timeIndex)
+        );
+
+        writer.piece(geom.points(), geom);
+
+        writer.writeGeometry();
+        writer.beginCellData();
+        writer.writeCellData(fieldName, tfield());
+    }
+
+    Info<< nl << "Smoothing/writing "
+            << timing.timeIncrement() << "s" << endl;
+
+    Info<< nl << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index 116463e51c..5f3a5d7d25 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -327,6 +327,8 @@ PatchFunction1/makePatchFunction1s.C
 PatchFunction1/coordinateScaling/coordinateScalings.C
 PatchFunction1/CodedField/makeCodedFields.C
 
+PatchFunction1/MappedFile/MappedFileFilterField.C
+
 meshStructure/meshStructure.C
 
 coupling/externalFileCoupler.C
diff --git a/src/meshTools/PatchFunction1/MappedFile/MappedFile.C b/src/meshTools/PatchFunction1/MappedFile/MappedFile.C
index 00f118662e..270ceaeaf1 100644
--- a/src/meshTools/PatchFunction1/MappedFile/MappedFile.C
+++ b/src/meshTools/PatchFunction1/MappedFile/MappedFile.C
@@ -47,6 +47,12 @@ Foam::PatchFunction1Types::MappedFile<Type>::MappedFile
     fieldTableName_(dict.getOrDefault<word>("fieldTable", entryName)),
     pointsName_(dict.getOrDefault<word>("points", "points")),
     mapMethod_(),
+    filterRadius_(dict.getOrDefault<scalar>("filterRadius", 0)),
+    filterSweeps_(dict.getOrDefault<label>("filterSweeps", 0)),
+    filterFieldPtr_(nullptr),
+    readerFormat_(),
+    readerFile_(),
+    readerPtr_(nullptr),
     mapperPtr_(nullptr),
     sampleTimes_(),
     sampleIndex_(-1, -1),
@@ -54,6 +60,40 @@ Foam::PatchFunction1Types::MappedFile<Type>::MappedFile
     sampleValues_(),
     offset_(Function1<Type>::NewIfPresent("offset", dict))
 {
+    // Simple sanity check
+    if ((filterSweeps_ < 1) || (filterRadius_ <= VSMALL))
+    {
+        filterRadius_ = 0;
+        filterSweeps_ = 0;
+    }
+
+    if (dict.readIfPresent("sampleFormat", readerFormat_))
+    {
+        dict.readEntry("sampleFile", readerFile_);
+
+        fileName fName(readerFile_);
+        fName.expand();
+
+        readerPtr_ = surfaceReader::New(readerFormat_, fName);
+    }
+
+    if (debug)
+    {
+        Info<< "mappedFile:" << nl;
+        if (readerFormat_.empty())
+        {
+            Info<< "    boundary format" << nl;
+        }
+        else
+        {
+            Info<< "    format:" << readerFormat_
+                << " file:" << readerFile_ << nl;
+        }
+
+        Info<< "    filter radius=" << filterRadius_
+            << " sweeps=" << filterSweeps_ << endl;
+    }
+
     if
     (
         dict.readIfPresent("mapMethod", mapMethod_)
@@ -88,6 +128,12 @@ Foam::PatchFunction1Types::MappedFile<Type>::MappedFile
     fieldTableName_(fieldTableName),
     pointsName_(dict.getOrDefault<word>("points", "points")),
     mapMethod_(),
+    filterRadius_(dict.getOrDefault<scalar>("filterRadius", 0)),
+    filterSweeps_(dict.getOrDefault<label>("filterSweeps", 0)),
+    filterFieldPtr_(nullptr),
+    readerFormat_(),
+    readerFile_(),
+    readerPtr_(nullptr),
     mapperPtr_(nullptr),
     sampleTimes_(),
     sampleIndex_(-1, -1),
@@ -95,6 +141,40 @@ Foam::PatchFunction1Types::MappedFile<Type>::MappedFile
     sampleValues_(),
     offset_(Function1<Type>::NewIfPresent("offset", dict))
 {
+    // Simple sanity check
+    if ((filterSweeps_ < 1) || (filterRadius_ <= VSMALL))
+    {
+        filterRadius_ = 0;
+        filterSweeps_ = 0;
+    }
+
+    if (dict.readIfPresent("sampleFormat", readerFormat_))
+    {
+        dict.readEntry("sampleFile", readerFile_);
+
+        fileName fName(readerFile_);
+        fName.expand();
+
+        readerPtr_ = surfaceReader::New(readerFormat_, fName);
+    }
+
+    if (debug)
+    {
+        Info<< "mappedFile:" << nl;
+        if (readerFormat_.empty())
+        {
+            Info<< "    boundary format" << nl;
+        }
+        else
+        {
+            Info<< "    format:" << readerFormat_
+                << " file:" << readerFile_ << nl;
+        }
+
+        Info<< "    filter radius=" << filterRadius_
+            << " sweeps=" << filterSweeps_ << endl;
+    }
+
     if
     (
         dict.readIfPresent("mapMethod", mapMethod_)
@@ -136,13 +216,27 @@ Foam::PatchFunction1Types::MappedFile<Type>::MappedFile
     fieldTableName_(rhs.fieldTableName_),
     pointsName_(rhs.pointsName_),
     mapMethod_(rhs.mapMethod_),
+    filterRadius_(rhs.filterRadius_),
+    filterSweeps_(rhs.filterSweeps_),
+    filterFieldPtr_(nullptr),
+    readerFormat_(rhs.readerFormat_),
+    readerFile_(rhs.readerFile_),
+    readerPtr_(nullptr),
     mapperPtr_(rhs.mapperPtr_.clone()),
     sampleTimes_(rhs.sampleTimes_),
     sampleIndex_(rhs.sampleIndex_),
     sampleAverage_(rhs.sampleAverage_),
     sampleValues_(rhs.sampleValues_),
     offset_(rhs.offset_.clone())
-{}
+{
+    if (!readerFormat_.empty() && !readerFile_.empty())
+    {
+        fileName fName(readerFile_);
+        fName.expand();
+
+        readerPtr_ = surfaceReader::New(readerFormat_, fName);
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -165,6 +259,7 @@ void Foam::PatchFunction1Types::MappedFile<Type>::autoMap
     }
 
     // Clear interpolator
+    filterFieldPtr_.reset(nullptr);
     mapperPtr_.reset(nullptr);
     sampleIndex_ = labelPair(-1, -1);
 }
@@ -195,6 +290,7 @@ void Foam::PatchFunction1Types::MappedFile<Type>::rmap
     }
 
     // Clear interpolator
+    filterFieldPtr_.reset(nullptr);
     mapperPtr_.reset(nullptr);
     sampleIndex_ = labelPair(-1, -1);
 }
@@ -211,13 +307,41 @@ void Foam::PatchFunction1Types::MappedFile<Type>::updateSampledValues
     tmp<Field<Type>> tvalues;
 
     // Update sampled data fields
+    if (readerPtr_)
+    {
+        wordList fieldNames = readerPtr_->fieldNames(sampleIndex);
+
+        label fieldIndex = fieldNames.find(fieldTableName_);
+
+        DebugInfo
+            << "checkTable : Update index=" << sampleIndex
+            << " field=" << fieldNames[fieldIndex] << endl;
+
+        tvalues = readerPtr_->field
+        (
+            sampleIndex,
+            fieldIndex,
+            pTraits<Type>::zero
+        );
+
+        if (tvalues().size() != mapperPtr_().sourceSize())
+        {
+            FatalErrorInFunction
+                << "Number of values (" << tvalues().size()
+                << ") differs from the number of points ("
+                <<  mapperPtr_().sourceSize() << ")"
+                << exit(FatalError);
+        }
+    }
+    else
     {
         const polyMesh& mesh = this->patch_.boundaryMesh().mesh();
         const Time& time = mesh.time();
 
         if (debug)
         {
-            Pout<< "checkTable : Reading values from "
+            Pout<< "checkTable : Update index=" << sampleIndex
+                << " Reading values from "
                 <<
                 (
                     "boundaryData"
@@ -269,6 +393,13 @@ void Foam::PatchFunction1Types::MappedFile<Type>::updateSampledValues
         tvalues = tmp<Field<Type>>::New(std::move(vals.field()));
     }
 
+    if (filterFieldPtr_)
+    {
+        DebugInfo
+            << "apply " << filterSweeps_ << " filter sweeps" << endl;
+
+        tvalues = filterFieldPtr_().evaluate(tvalues, filterSweeps_);
+    }
 
     // From input values to interpolated (sampled) positions
     field = mapperPtr_().interpolate(tvalues);
@@ -285,7 +416,74 @@ void Foam::PatchFunction1Types::MappedFile<Type>::checkTable
     const Time& time = mesh.time();
 
     // Initialise
-    if (!mapperPtr_)
+    if (!mapperPtr_ && readerPtr_)
+    {
+        auto& reader = readerPtr_();
+
+        const meshedSurface& geom = reader.geometry(0);
+
+        sampleTimes_ = reader.times();
+
+        // We may have been passed in geometry without any faces
+        // eg, boundaryData
+
+        const pointField& samplePoints =
+        (
+            geom.nFaces() ? geom.faceCentres() : geom.points()
+        );
+
+        DebugInfo
+            << "Read " << samplePoints.size() << " sample points from "
+            << readerFile_ << endl
+            << "Found times "
+            << pointToPointPlanarInterpolation::timeNames(sampleTimes_) << nl;
+
+        // tbd: run-time selection
+        const bool nearestOnly =
+        (
+            !mapMethod_.empty() && !mapMethod_.starts_with("planar")
+        );
+
+
+        // Allocate the interpolator
+        if (this->faceValues())
+        {
+            mapperPtr_.reset
+            (
+                new pointToPointPlanarInterpolation
+                (
+                    samplePoints,
+                    this->localPosition(this->patch_.faceCentres()),
+                    perturb_,
+                    nearestOnly
+                )
+            );
+        }
+        else
+        {
+            mapperPtr_.reset
+            (
+                new pointToPointPlanarInterpolation
+                (
+                    samplePoints,
+                    this->localPosition(this->patch_.localPoints()),
+                    perturb_,
+                    nearestOnly
+                )
+            );
+        }
+
+        // Setup median filter (if any)
+        if (filterSweeps_ > 0)
+        {
+            filterFieldPtr_.reset(new FilterField(geom, filterRadius_));
+        }
+        else
+        {
+            filterFieldPtr_.reset(nullptr);
+        }
+    }
+    else if (!mapperPtr_)
     {
         // Reread values and interpolate
         const fileName samplePointsFile
@@ -349,7 +547,6 @@ void Foam::PatchFunction1Types::MappedFile<Type>::checkTable
             );
         }
 
-
         // Read the times for which data is available
         const fileName samplePointsDir = samplePointsFile.path();
         sampleTimes_ = Time::findTimes(samplePointsDir);
@@ -358,11 +555,21 @@ void Foam::PatchFunction1Types::MappedFile<Type>::checkTable
             << "Found times "
             << pointToPointPlanarInterpolation::timeNames(sampleTimes_)
             << endl;
+
+        // Setup median filter (if any)
+        if (filterSweeps_ > 0)
+        {
+            filterFieldPtr_.reset(new FilterField(samplePoints, filterRadius_));
+        }
+        else
+        {
+            filterFieldPtr_.reset(nullptr);
+        }
     }
 
 
     // Find range of current time indices in sampleTimes
-    Pair<label> timeIndices = instant::findRange
+    labelPair timeIndices = instant::findRange
     (
         sampleTimes_,
         t,  //mesh.time().value(),
@@ -555,6 +762,12 @@ void Foam::PatchFunction1Types::MappedFile<Type>::writeEntries
     Ostream& os
 ) const
 {
+    if (!readerFormat_.empty() && !readerFile_.empty())
+    {
+        os.writeEntry("readerFormat", readerFormat_);
+        os.writeEntry("readerFile", readerFile_);
+    }
+
     os.writeEntryIfDifferent
     (
         "fieldTable",
@@ -579,6 +792,12 @@ void Foam::PatchFunction1Types::MappedFile<Type>::writeEntries
 
     os.writeEntryIfDifferent<scalar>("perturb", 1e-5, perturb_);
 
+    if (filterSweeps_ >= 1)
+    {
+        os.writeEntry("filterRadius", filterRadius_);
+        os.writeEntry("filterSweeps", filterSweeps_);
+    }
+
     if (offset_)
     {
         offset_->writeData(os);
diff --git a/src/meshTools/PatchFunction1/MappedFile/MappedFile.H b/src/meshTools/PatchFunction1/MappedFile/MappedFile.H
index b45a7cca8e..b8ef66b9d8 100644
--- a/src/meshTools/PatchFunction1/MappedFile/MappedFile.H
+++ b/src/meshTools/PatchFunction1/MappedFile/MappedFile.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -49,10 +49,26 @@ Description
     \endtable
 
     Options for the \c mapMethod entry:
-    \verbatim
+    \table
+      Property    | Description
       nearest     | Use nearest points only (avoids triangulation)
       planar      | Interpolation using 2D Delaunay triangulation
-    \endverbatim
+    \endtable
+
+    Options for reading and filtering:
+    \table
+      Property     | Description
+      sampleFormat | The surfaceReader format (eg, ensight)
+      sampleFile   | <case>/foo/bar/window.case
+      filterRadius | Search radius [m] for median filter neighbours
+      filterSweeps | Filter sweeps for median filter
+    \endtable
+
+Note
+    The MappedFile handling currently has two different reading types:
+    a builtin boundaryData reader (default) and a surfaceReader
+    (specified by sampleFormat and sampleFile keywords). The generic
+    surfaceReader handling does \em not support field averaging.
 
 SourceFiles
     MappedFile.C
@@ -66,6 +82,8 @@ SourceFiles
 #include "Function1.H"
 #include "Pair.H"
 #include "instantList.H"
+#include "surfaceReader.H"
+#include "MappedFileFilterField.H"
 #include "pointToPointPlanarInterpolation.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -104,6 +122,24 @@ class MappedFile
         //- Interpolation scheme to use (default is empty == "planar")
         word mapMethod_;
 
+        //- Radius for filter
+        scalar filterRadius_;
+
+        //- Number of median filter sweeps
+        label filterSweeps_;
+
+        //- Median filtering (for input values)
+        mutable autoPtr<FilterField> filterFieldPtr_;
+
+        //- Format name for surfaceReader
+        word readerFormat_;
+
+        //- File name associated with surfaceReader
+        fileName readerFile_;
+
+        //- Meshed surface with fields (for input values)
+        mutable autoPtr<surfaceReader> readerPtr_;
+
         //- 2D interpolation (for 'planar' mapMethod)
         mutable autoPtr<pointToPointPlanarInterpolation> mapperPtr_;
 
diff --git a/src/meshTools/PatchFunction1/MappedFile/MappedFileFilterField.C b/src/meshTools/PatchFunction1/MappedFile/MappedFileFilterField.C
new file mode 100644
index 0000000000..59cde08c01
--- /dev/null
+++ b/src/meshTools/PatchFunction1/MappedFile/MappedFileFilterField.C
@@ -0,0 +1,381 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "MappedFileFilterField.H"
+#include "oneField.H"
+#include "MeshedSurfaces.H"
+#include "MinMax.H"
+#include "indexedOctree.H"
+#include "treeDataPoint.H"
+#include <algorithm>
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace PatchFunction1Types
+{
+    defineTypeNameAndDebug(FilterField, 0);
+}
+}
+
+
+const Foam::Enum
+<
+    Foam::PatchFunction1Types::FilterField::RBF_type
+>
+Foam::PatchFunction1Types::FilterField::RBF_typeNames_
+{
+    { RBF_type::RBF_linear, "linear" },
+    { RBF_type::RBF_quadratic, "quadratic" },
+    { RBF_type::RBF_linear, "default" },
+};
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+//- Linear basis function: 1 - (r/r0)
+static inline scalar linearWeight
+(
+    const point& p,
+    const point& p0,
+    const scalar radiusSqr
+)
+{
+    return (1 - ::sqrt(p.distSqr(p0) / radiusSqr));
+}
+
+
+//- Quadratic basis function: 1 - (r/r0)^2
+static inline scalar quadraticWeight
+(
+    const point& p,
+    const point& p0,
+    const scalar radiusSqr
+)
+{
+    return (1 - p.distSqr(p0) / radiusSqr);
+}
+
+
+//- Construct search tree for points
+static autoPtr<indexedOctree<treeDataPoint>>
+createTree(const pointField& points)
+{
+    // Avoid bounding box error in case of 2D cases
+    treeBoundBox bb(points);
+    bb.grow(1e-4);
+
+    const int debug = PatchFunction1Types::FilterField::debug;
+
+    const int oldDebug = indexedOctreeBase::debug;
+
+    if (debug & 2)
+    {
+        indexedOctreeBase::debug = 1;
+    }
+
+    auto treePtr = autoPtr<indexedOctree<treeDataPoint>>::New
+    (
+        treeDataPoint(points),
+        bb,             // overall search domain
+        points.size(),  // maxLevel
+        16,             // leafsize
+        1.0             // duplicity
+    );
+
+    indexedOctreeBase::debug = oldDebug;
+
+    if (debug & 2)
+    {
+        const auto& tree = treePtr();
+
+        OFstream os("indexedOctree.obj");
+        tree.writeOBJ(os);
+        Info<< "=================" << endl;
+        Info<< "have " << tree.nodes().size() << " nodes" << nl;
+        Info<< "=================" << endl;
+    }
+
+    return treePtr;
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template
+<
+    class TreeType,
+    class RadiusSqrOp,
+    class BasisFunctionOp,
+    class PointWeightFieldType
+>
+void Foam::PatchFunction1Types::FilterField::buildWeightsImpl
+(
+    const TreeType& tree,
+    const RadiusSqrOp& radiusSqrOp,
+    const BasisFunctionOp& basisFuncOp,
+    const PointWeightFieldType& pointWeightFld
+)
+{
+    tmp<pointField> tpoints = tree.shapes().centres();
+
+    const auto& points = tpoints();
+
+    const label nPoints = points.size();
+
+    weights_.clear();
+    addressing_.clear();
+
+    weights_.resize(nPoints);
+    addressing_.resize(nPoints);
+
+    labelHashSet neighbours(2*128);
+
+    // Catch degenerate case where weight/addressing not actually needed
+    bool usesNeighbours = false;
+
+    for (label pointi = 0; pointi < nPoints; ++pointi)
+    {
+        const point& p0 = points[pointi];
+        auto& currAddr = addressing_[pointi];
+        auto& currWeight = weights_[pointi];
+
+        // Search with local sphere
+        const scalar radiusSqr = radiusSqrOp(pointi);
+
+        tree.findSphere(p0, radiusSqr, neighbours);
+
+        // Paranoid, enforce identity weighting as a minimum
+        if (neighbours.size() < 2)
+        {
+            currAddr.resize(1, pointi);
+            currWeight.resize(1, scalar(1));
+            continue;
+        }
+
+        usesNeighbours = true;
+
+        currAddr = neighbours.sortedToc();
+        currWeight.resize_nocopy(currAddr.size());
+
+        scalar sumWeight = 0;
+
+        forAll(currAddr, i)
+        {
+            const point& p = points[currAddr[i]];
+
+            currWeight[i] = basisFuncOp(p, p0, radiusSqr);
+
+            // Eg, face area weighting.
+            // - cast required for oneField()
+            currWeight[i] *= static_cast<scalar>(pointWeightFld[i]);
+            sumWeight += currWeight[i];
+        }
+
+        for (auto& w : currWeight)
+        {
+            w /= sumWeight;
+        }
+    }
+
+    if (!usesNeighbours)
+    {
+        // No addressing/weighting required
+        reset();
+    }
+
+    if (debug && !addressing_.empty())
+    {
+        labelMinMax limits(addressing_[0].size());
+
+        label total = 0;
+
+        for (const labelList& addr : addressing_)
+        {
+            const label n = addr.size();
+
+            limits.add(n);
+            total += n;
+        }
+
+        Info<< "Weight neighbours: min=" << limits.min()
+            << " avg=" << (total / addressing_.size())
+            << " max=" << limits.max() << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::PatchFunction1Types::FilterField::FilterField
+(
+    const pointField& points,
+    const scalar radius,
+    const RBF_type interp
+)
+{
+    reset(points, radius, interp);
+}
+
+
+Foam::PatchFunction1Types::FilterField::FilterField
+(
+    const meshedSurface& geom,
+    const scalar radius,
+    const bool relative,
+    const RBF_type interp
+)
+{
+    reset(geom, radius, relative, interp);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::PatchFunction1Types::FilterField::reset()
+{
+    addressing_.clear();
+    weights_.clear();
+}
+
+
+void Foam::PatchFunction1Types::FilterField::reset
+(
+    const pointField& points,
+    const scalar radius,
+    const RBF_type interp
+)
+{
+    scalarField dummyWeights;
+
+    if (debug)
+    {
+        Info<< "Apply " << RBF_typeNames_[interp] << " filter,"
+            << " radius=" << radius << nl
+            << "Create tree..." << endl;
+    }
+
+    autoPtr<indexedOctree<treeDataPoint>> treePtr
+        = createTree(points);
+
+    const scalar radiusSqr = sqr(radius);
+
+    auto weightFunc = linearWeight;
+
+    if (interp == RBF_type::RBF_quadratic)
+    {
+        weightFunc = quadraticWeight;
+    }
+
+
+    // Use (constant) absolute size
+    buildWeightsImpl
+    (
+        treePtr(),
+        [=](label index) -> scalar { return radiusSqr; },
+        weightFunc,
+        oneField()
+    );
+}
+
+
+void Foam::PatchFunction1Types::FilterField::reset
+(
+    const meshedSurface& geom,
+    const scalar radius,
+    const bool relative,
+    const RBF_type interp
+)
+{
+    if (!geom.nFaces())
+    {
+        // Received geometry without any faces eg, boundaryData
+        reset(geom.points(), radius, interp);
+        return;
+    }
+
+    if (debug)
+    {
+        Info<< "Apply " << RBF_typeNames_[interp] << " filter,";
+
+        if (relative)
+        {
+            Info<< " relative";
+        }
+
+        Info<< " radius=" << radius << endl;
+        Info<< "Create tree..." << endl;
+    }
+
+    autoPtr<indexedOctree<treeDataPoint>> treePtr
+        = createTree(geom.faceCentres());
+
+    const scalar radiusSqr = sqr(radius);
+
+    auto weightFunc = linearWeight;
+
+    if (interp == RBF_type::RBF_quadratic)
+    {
+        weightFunc = quadraticWeight;
+    }
+
+
+    if (relative)
+    {
+        // Use relative size
+        buildWeightsImpl
+        (
+            treePtr(),
+            [&](label index) -> scalar
+            {
+                return (radiusSqr * geom.sphere(index));
+            },
+            weightFunc,
+            geom.magSf()
+        );
+    }
+    else
+    {
+        // Use (constant) absolute size
+        buildWeightsImpl
+        (
+            treePtr(),
+            [=](label index) -> scalar { return radiusSqr; },
+            weightFunc,
+            geom.magSf()
+        );
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/PatchFunction1/MappedFile/MappedFileFilterField.H b/src/meshTools/PatchFunction1/MappedFile/MappedFileFilterField.H
new file mode 100644
index 0000000000..2cc99df6ac
--- /dev/null
+++ b/src/meshTools/PatchFunction1/MappedFile/MappedFileFilterField.H
@@ -0,0 +1,213 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::PatchFunction1Types::FilterField
+
+Description
+    The FilterField helper class provides a multi-sweep median filter
+    for a Field of data associated with a geometric point cloud.
+
+    The points can be freestanding or the faceCentres (or points)
+    of a meshedSurface, for example.
+
+    Using an initial specified search radius, the nearest point
+    neighbours are gathered and addressing/weights are built for them.
+    This currently uses an area-weighted, linear RBF interpolator
+    with provision for quadratic RBF interpolator etc.
+
+    After the weights and addressing are established,
+    the evaluate() method can be called to apply a median filter
+    to data fields, with a specified number of sweeps.
+
+Note
+    When handling large search radii and/or an extensive number of
+    filter sweeps, compiling with openmp can yield some speedup.
+
+SourceFiles
+    MappedFileFilterField.C
+    MappedFileFilterFieldTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_PatchFunction1Types_MappedFileFilterField_H
+#define Foam_PatchFunction1Types_MappedFileFilterField_H
+
+#include "primitiveFields.H"
+#include "pointField.H"
+#include "Enum.H"
+#include "className.H"
+#include "MeshedSurfacesFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace PatchFunction1Types
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class FilterField Declaration
+\*---------------------------------------------------------------------------*/
+
+class FilterField
+{
+public:
+
+     // Data Types
+
+        //- Basis type
+        enum RBF_type
+        {
+            RBF_linear,
+            RBF_quadratic
+        };
+
+
+private:
+
+    // Private Data
+
+        //- Names for RBF basis types
+        static const Enum<RBF_type> RBF_typeNames_;
+
+        //- Addressing
+        List<labelList> addressing_;
+
+        //- Weights
+        List<scalarField> weights_;
+
+
+    // Private Member Functions
+
+        //- Construct neighbour point weights and addressing based
+        //- on specified search radius
+        //
+        // Search with local sphere
+        //
+        // RadiusSqrOp (index) -> radiusSqr
+        // BasisFunctionOp(point, point, radius^2) -> RBF
+        // PointWeightFieldType [index] -> weighting
+        template
+        <
+            class TreeType,
+            class RadiusSqrOp,
+            class BasisFunctionOp,
+            class PointWeightFieldType
+        >
+        void buildWeightsImpl
+        (
+            const TreeType& tree,
+            const RadiusSqrOp& radiusSqrOp,
+            const BasisFunctionOp& basisFuncOp,
+            const PointWeightFieldType& pointWeightFld
+        );
+
+
+public:
+
+    //- Runtime type information
+    ClassName("filterField");
+
+
+    // Constructors
+
+        //- Default construct
+        FilterField() noexcept = default;
+
+        //- Construct with weights for a field of points
+        //- (constant search radius)
+        FilterField
+        (
+            const pointField& points,
+            const scalar radius,
+            const RBF_type interp = RBF_type::RBF_linear
+        );
+
+        //- Construct with weights for faceCentres of a meshedSurface
+        //- using a constant search radius search
+        //- (or relative radius multiplier for the face bounding spheres).
+        FilterField
+        (
+            const meshedSurface& geom,
+            const scalar radius,
+            const bool relative = false,
+            const RBF_type interp = RBF_type::RBF_linear
+        );
+
+
+    // Member Functions
+
+        //- Reset to unweighted (pass-through)
+        void reset();
+
+        //- Create weights for field of points (constant search radius)
+        void reset
+        (
+            const pointField& points,
+            const scalar radius,
+            const RBF_type interp = RBF_type::RBF_linear
+        );
+
+        //- Create weights for meshedSurface using a constant search radius
+        //- or optionally with a search radius multiplier for the
+        //- face bounding spheres.
+        void reset
+        (
+            const meshedSurface& geom,
+            const scalar radius,
+            const bool relative = false,
+            const RBF_type interp = RBF_type::RBF_linear
+        );
+
+
+    // Evaluation
+
+        //- Return the median smoothed field
+        template<class Type>
+        tmp<Field<Type>> evaluate
+        (
+            const tmp<Field<Type>>& tinput,
+            const label nSweeps
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace PatchFunction1Types
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "MappedFileFilterFieldTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/PatchFunction1/MappedFile/MappedFileFilterFieldTemplates.C b/src/meshTools/PatchFunction1/MappedFile/MappedFileFilterFieldTemplates.C
new file mode 100644
index 0000000000..92180ade66
--- /dev/null
+++ b/src/meshTools/PatchFunction1/MappedFile/MappedFileFilterFieldTemplates.C
@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "indexedOctree.H"
+#include "treeDataPoint.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::PatchFunction1Types::FilterField::evaluate
+(
+    const tmp<Field<Type>>& tinput,
+    const label nSweeps
+) const
+{
+    if (nSweeps < 1 || !tinput.good())
+    {
+        // Nothing to do
+        return tinput;
+    }
+
+    label nPoints = tinput().size();
+
+    const label nAddr = addressing_.size();
+
+    if (!nPoints || !nAddr)
+    {
+        // No input or using an identity mapping
+        return tinput;
+    }
+
+    // Output
+    auto toutput = tmp<Field<Type>>::New(nPoints);
+
+    if (nAddr < nPoints)
+    {
+        WarningInFunction
+            << "Addressing/weights shorter than input field"
+            << endl;
+
+        // Restrict addressing space within loop
+        nPoints = nAddr;
+
+        // Fill trailing portion with the input values
+        toutput.ref().slice(nAddr) = tinput().slice(nAddr);
+    }
+
+    // Intermediate buffer
+    tmp<Field<Type>> tbuffer;
+
+    if (nSweeps == 1)
+    {
+        // Simple reference ie enough
+        tbuffer.cref(tinput);
+    }
+    else
+    {
+        // Need buffer for swap - copy or steal contents
+        tbuffer.reset(tinput.ptr());
+    }
+    tinput.clear();
+
+    // If there are any 'senseless' sweeps
+    // (ie, with identity mapping of values), they will have
+    // been detected during addressing construction, can ignore here
+
+    for (label sweep = 0; sweep < nSweeps; ++sweep)
+    {
+        if (sweep > 0)
+        {
+            // Reuse previous output for subsequent sweeps
+            tbuffer.swap(toutput);
+        }
+
+        const auto& input = tbuffer();
+        auto& output = toutput.ref();
+
+        #pragma omp parallel for if (nPoints > 1000)
+        for (label pointi = 0; pointi < nPoints; ++pointi)
+        {
+            const auto& addr = addressing_[pointi];
+            const auto& weight = weights_[pointi];
+
+            auto& interp = output[pointi];
+
+            if (addr.empty())
+            {
+                // Could happen if the search radius was really small
+                interp = input[pointi];
+            }
+            else
+            {
+                interp = Zero;
+
+                forAll(addr, i)
+                {
+                    interp += (weight[i] * input[addr[i]]);
+                }
+            }
+        }
+    }
+
+    return toutput;
+}
+
+
+// ************************************************************************* //
diff --git a/src/surfMesh/Make/files b/src/surfMesh/Make/files
index 9d786fb75a..1216f4f491 100644
--- a/src/surfMesh/Make/files
+++ b/src/surfMesh/Make/files
@@ -60,6 +60,7 @@ triSurface/patches/surfacePatch.C
 readers = readers
 
 $(readers)/common/surfaceReader.C
+$(readers)/boundary/boundaryDataSurfaceReader.C
 $(readers)/ensight/ensightSurfaceReader.C
 
 writers = writers
diff --git a/src/surfMesh/readers/boundary/boundaryDataSurfaceReader.C b/src/surfMesh/readers/boundary/boundaryDataSurfaceReader.C
new file mode 100644
index 0000000000..2f51ee0383
--- /dev/null
+++ b/src/surfMesh/readers/boundary/boundaryDataSurfaceReader.C
@@ -0,0 +1,281 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "boundaryDataSurfaceReader.H"
+#include "rawIOField.H"
+#include "argList.H"
+#include "Time.H"
+#include "addToRunTimeSelectionTable.H"
+#include <algorithm>
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(boundaryDataSurfaceReader, 0);
+    addToRunTimeSelectionTable
+    (
+        surfaceReader,
+        boundaryDataSurfaceReader,
+        fileName
+    );
+}
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::pointField Foam::boundaryDataSurfaceReader::readPoints
+(
+    const Time& runTime,
+    const fileName& baseDir,
+    const word& pointsName
+)
+{
+    fileName pointsFile
+    (
+        baseDir / (pointsName.empty() ? "points" : pointsName)
+    );
+    pointsFile.toAbsolute();
+
+    IOobject io
+    (
+        pointsFile,             // absolute path
+        runTime,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        IOobject::NO_REGISTER,
+        true                    // global object (currently not used)
+    );
+
+    DebugInfo
+        << "File: " << io.objectPath() << endl;
+
+    // Read data (no average value!)
+    rawIOField<point> rawData(io);
+
+    pointField points(std::move(rawData.field()));
+
+    DebugInfo
+        << "File: " << io.objectPath()
+        << " " << points.size() << " points" << endl;
+
+    return points;
+}
+
+
+Foam::pointField Foam::boundaryDataSurfaceReader::readPoints
+(
+    const fileName& dirName,
+    const word& pointsName
+)
+{
+    refPtr<Time> timePtr(Time::New(argList::envGlobalPath()));
+
+    return readPoints(*timePtr, dirName, pointsName);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::boundaryDataSurfaceReader::readCase()
+{
+    DebugInFunction << endl;
+
+    timeValues_ = TimePaths::findTimes(baseDir_);
+}
+
+
+void Foam::boundaryDataSurfaceReader::readGeometry
+(
+    meshedSurface& surf,
+    const label timeIndex
+)
+{
+    surf.clear();
+
+    pointField points(this->readPoints(baseDir_, pointsName_));
+
+    surf = meshedSurface(std::move(points), faceList());
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::boundaryDataSurfaceReader::readFieldTemplate
+(
+    const label timeIndex,
+    const label fieldIndex
+) const
+{
+    Type dummyAvg;
+
+    return readField<Type>
+    (
+        baseDir_,
+        timeValues_[timeIndex],
+        fieldNames_[fieldIndex],
+        dummyAvg
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::boundaryDataSurfaceReader::boundaryDataSurfaceReader
+(
+    const fileName& fName,
+    const word& pointsName
+)
+:
+    surfaceReader(fName),
+    baseDir_(fName.path()),
+    pointsName_(pointsName),
+    timeValues_(),
+    fieldNames_(),
+    surfPtr_(nullptr)
+{
+    baseDir_.toAbsolute();
+    debug = 1;
+    DebugInFunction << endl;
+    Info<< "create with " << baseDir_ << endl;
+
+    readCase();
+}
+
+
+// * * * * * * * * * * * * * Public Member Functions   * * * * * * * * * * * //
+
+const Foam::meshedSurface& Foam::boundaryDataSurfaceReader::geometry
+(
+    const label timeIndex
+)
+{
+    DebugInFunction << endl;
+
+    if (!surfPtr_)
+    {
+        surfPtr_.reset(new meshedSurface);
+        readGeometry(*surfPtr_, timeIndex);
+    }
+
+    return *surfPtr_;
+}
+
+
+Foam::instantList Foam::boundaryDataSurfaceReader::times() const
+{
+    return timeValues_;
+}
+
+
+Foam::wordList Foam::boundaryDataSurfaceReader::fieldNames
+(
+    const label timeIndex
+) const
+{
+    if (timeValues_.empty() || timeIndex >= timeValues_.size())
+    {
+        fieldNames_.clear();
+        return wordList();
+    }
+
+    fileNameList items =
+        fileHandler().readDir(baseDir_/timeValues_[timeIndex].name());
+
+    fieldNames_.resize_nocopy(items.size());
+
+    std::transform
+    (
+        items.begin(),
+        items.end(),
+        fieldNames_.begin(),
+        [](const fileName& f) { return word(f); }
+    );
+
+    Foam::sort(fieldNames_);
+
+    return fieldNames_;
+}
+
+
+Foam::tmp<Foam::Field<Foam::scalar>> Foam::boundaryDataSurfaceReader::field
+(
+    const label timeIndex,
+    const label fieldIndex,
+    const scalar& refValue
+) const
+{
+    return readFieldTemplate<scalar>(timeIndex, fieldIndex);
+}
+
+
+Foam::tmp<Foam::Field<Foam::vector>> Foam::boundaryDataSurfaceReader::field
+(
+    const label timeIndex,
+    const label fieldIndex,
+    const vector& refValue
+) const
+{
+    return readFieldTemplate<vector>(timeIndex, fieldIndex);
+}
+
+
+Foam::tmp<Foam::Field<Foam::sphericalTensor>>
+Foam::boundaryDataSurfaceReader::field
+(
+    const label timeIndex,
+    const label fieldIndex,
+    const sphericalTensor& refValue
+) const
+{
+    return readFieldTemplate<sphericalTensor>(timeIndex, fieldIndex);
+}
+
+
+Foam::tmp<Foam::Field<Foam::symmTensor>> Foam::boundaryDataSurfaceReader::field
+(
+    const label timeIndex,
+    const label fieldIndex,
+    const symmTensor& refValue
+) const
+{
+    return readFieldTemplate<symmTensor>(timeIndex, fieldIndex);
+}
+
+
+Foam::tmp<Foam::Field<Foam::tensor>> Foam::boundaryDataSurfaceReader::field
+(
+    const label timeIndex,
+    const label fieldIndex,
+    const tensor& refValue
+) const
+{
+    return readFieldTemplate<tensor>(timeIndex, fieldIndex);
+}
+
+
+// ************************************************************************* //
diff --git a/src/surfMesh/readers/boundary/boundaryDataSurfaceReader.H b/src/surfMesh/readers/boundary/boundaryDataSurfaceReader.H
new file mode 100644
index 0000000000..76821c298d
--- /dev/null
+++ b/src/surfMesh/readers/boundary/boundaryDataSurfaceReader.H
@@ -0,0 +1,227 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::boundaryDataSurfaceReader
+
+Description
+    A boundaryData format surface reader.
+    However, the "surface" represented by boundaryData is actually only
+    point data!
+
+    // Points
+    <case>/constant/region0/"boundaryData"/patchName/points
+
+    // Values
+    <case>/constant/region0/"boundaryData"/patchName/TIME/field
+
+SourceFiles
+    boundaryDataSurfaceReader.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_boundaryDataSurfaceReader_H
+#define Foam_boundaryDataSurfaceReader_H
+
+#include "surfaceReader.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class Time;
+
+/*---------------------------------------------------------------------------*\
+                  Class boundaryDataSurfaceReader Declaration
+\*---------------------------------------------------------------------------*/
+
+class boundaryDataSurfaceReader
+:
+    public surfaceReader
+{
+    // Private Data
+
+        //- Base directory
+        fileName baseDir_;
+
+        //- Name of the 'points' file
+        word pointsName_;
+
+        //- Times
+        instantList timeValues_;
+
+        //- Field names
+        mutable List<word> fieldNames_;
+
+        //- Pointer to the surface
+        autoPtr<meshedSurface> surfPtr_;
+
+
+    // Private Member Functions
+
+        //- Populate time instances
+        void readCase();
+
+        //- Read geometry
+        void readGeometry(meshedSurface& surf, const label timeIndex);
+
+        //- Read and return given field
+        template<class Type>
+        tmp<Field<Type>> readFieldTemplate
+        (
+            const label timeIndex,
+            const label fieldIndex
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("boundaryData");
+
+
+    // Constructors
+
+        //- Construct from fileName
+        explicit boundaryDataSurfaceReader
+        (
+            const fileName& fName,
+            const word& pointsName = "points"
+        );
+
+
+    //- Destructor
+    virtual ~boundaryDataSurfaceReader() = default;
+
+
+    // Static Member Functions
+
+        //- Read points file
+        static pointField readPoints
+        (
+            const Time& runTime,
+            const fileName& baseDir,
+            const word& pointsName = "points"
+        );
+
+        //- Read points file
+        static pointField readPoints
+        (
+            const fileName& baseDir,
+            const word& pointsName = "points"
+        );
+
+        //- Read and return given field
+        template<class Type>
+        static tmp<Field<Type>> readField
+        (
+            const Time& runTime,
+            const fileName& baseDir,
+            const instant& timeDir,
+            const word& fieldName,
+            Type& avg
+        );
+
+        //- Read and return given field
+        template<class Type>
+        static tmp<Field<Type>> readField
+        (
+            const fileName& baseDir,
+            const instant& timeDir,
+            const word& fieldName,
+            Type& avg
+        );
+
+
+    // Member Functions
+
+        //- Return a reference to the surface geometry
+        virtual const meshedSurface& geometry(const label timeIndex);
+
+        //- Return a list of the available times
+        virtual instantList times() const;
+
+        //- Return a list of the available fields at a given time
+        virtual wordList fieldNames(const label timeIndex) const;
+
+        //- Return a scalar field at a given time
+        virtual tmp<Field<scalar>> field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const scalar& refValue = pTraits<scalar>::zero
+        ) const;
+
+        //- Return a vector field at a given time
+        virtual tmp<Field<vector>> field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const vector& refValue = pTraits<vector>::zero
+        ) const;
+
+        //- Return a sphericalTensor field at a given time
+        virtual tmp<Field<sphericalTensor>> field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const sphericalTensor& refValue = pTraits<sphericalTensor>::zero
+        ) const;
+
+        //- Return a symmTensor field at a given time
+        virtual tmp<Field<symmTensor>> field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const symmTensor& refValue = pTraits<symmTensor>::zero
+        ) const;
+
+        //- Return a tensor field at a given time
+        virtual tmp<Field<tensor>> field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const tensor& refValue = pTraits<tensor>::zero
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "boundaryDataSurfaceReaderTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfMesh/readers/boundary/boundaryDataSurfaceReaderTemplates.C b/src/surfMesh/readers/boundary/boundaryDataSurfaceReaderTemplates.C
new file mode 100644
index 0000000000..8494099738
--- /dev/null
+++ b/src/surfMesh/readers/boundary/boundaryDataSurfaceReaderTemplates.C
@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "boundaryDataSurfaceReader.H"
+#include "rawIOField.H"
+#include "argList.H"
+#include "Time.H"
+#include <algorithm>
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::boundaryDataSurfaceReader::readField
+(
+    const Time& runTime,
+    const fileName& baseDir,
+    const instant& timeDir,
+    const word& fieldName,
+    Type& avg
+)
+{
+    // Reread values and interpolate
+    fileName valuesFile(baseDir/timeDir.name()/fieldName);
+    valuesFile.toAbsolute();
+
+    IOobject io
+    (
+        valuesFile,   // absolute path
+        runTime,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        IOobject::NO_REGISTER,
+        true                // global object (currently not used)
+    );
+
+    DebugInfo
+        << "File: " << io.objectPath() << endl;
+
+    // Read data (TDB: setAverage)
+    rawIOField<Type> rawData(io, IOobjectOption::READ_IF_PRESENT);
+
+    if (rawData.hasAverage())
+    {
+        avg = rawData.average();
+    }
+
+    DebugInfo
+        << "File: " << io.objectPath()
+        << " " << rawData.size() << " values" << endl;
+
+    return tmp<Field<Type>>::New(std::move(rawData.field()));
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::boundaryDataSurfaceReader::readField
+(
+    const fileName& baseDir,
+    const instant& timeDir,
+    const word& fieldName,
+    Type& avg
+)
+{
+    refPtr<Time> timePtr(Time::New(argList::envGlobalPath()));
+
+    return readField<Type>(baseDir, timeDir, fieldName, avg);
+}
+
+
+// ************************************************************************* //