diff --git a/applications/test/surface-sampling/Make/files b/applications/test/surface-sampling/Make/files
new file mode 100644
index 0000000000..e6db680175
--- /dev/null
+++ b/applications/test/surface-sampling/Make/files
@@ -0,0 +1,4 @@
+mydebugSurfaceWriter.C
+Test-surface-sampling.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-surface-sampling
diff --git a/applications/test/surface-sampling/Make/options b/applications/test/surface-sampling/Make/options
new file mode 100644
index 0000000000..4a8ebb6560
--- /dev/null
+++ b/applications/test/surface-sampling/Make/options
@@ -0,0 +1,22 @@
+EXE_INC = \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/fileFormats/lnInclude \
+ -I$(LIB_SRC)/surfMesh/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/sampling/lnInclude \
+ -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+ -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+ -I$(LIB_SRC)/transportModels \
+ -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel
+
+EXE_LIBS = \
+ -lfiniteVolume \
+ -lfvOptions \
+ -lfileFormats \
+ -lsurfMesh \
+ -lmeshTools \
+ -lsampling \
+ -lturbulenceModels \
+ -lincompressibleTurbulenceModels \
+ -lincompressibleTransportModels \
+ -latmosphericModels
diff --git a/applications/test/surface-sampling/Test-surface-sampling.C b/applications/test/surface-sampling/Test-surface-sampling.C
new file mode 100644
index 0000000000..bec2a2a71c
--- /dev/null
+++ b/applications/test/surface-sampling/Test-surface-sampling.C
@@ -0,0 +1,275 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023 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 .
+
+Application
+ Test-surface-sampling
+
+Description
+ Simple test of surface sampling, including timings
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "profiling.H"
+#include "clockTime.H"
+
+#include "fileName.H"
+#include "sampledSurfaces.H"
+#include "IOstreams.H"
+#include "OSspecific.H"
+#include "profiling.H"
+#include "ReadFields.H"
+#include "volFields.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+ argList::addNote
+ (
+ "Loads mesh and fields from latest time and writes multiple times"
+ );
+
+ argList::noFunctionObjects(); // Disallow function objects
+ argList::addVerboseOption("additional verbosity");
+ #include "addProfilingOption.H"
+
+ argList::addOption
+ (
+ "sample",
+ "file",
+ "Name of surface sampling to use (default: test-sample)"
+ );
+ argList::addOption("output", "Begin output iteration (default: 10000)");
+ argList::addOption("count", "Number of writes (default: 1)");
+
+ #include "setRootCase.H"
+ #include "createTime.H"
+
+ const int verbose = args.verbose();
+ const label firstOutput = args.getOrDefault("output", 10000);
+ const label nOutput = args.getOrDefault("count", 1);
+
+ // Select latestTime, including 0 and constant
+ {
+ const auto& times = runTime.times();
+ const label timeIndex = (times.size()-1);
+
+ if (timeIndex < 0)
+ {
+ FatalErrorInFunction
+ << "No times!"
+ << exit(FatalError);
+ }
+
+ runTime.setTime(times[timeIndex], timeIndex);
+ }
+
+ // #include "createMesh.H"
+
+ Info << "Create mesh time = " << runTime.timeName() << nl;
+
+ fvMesh mesh
+ (
+ IOobject
+ (
+ polyMesh::defaultRegion,
+ runTime.timeName(),
+ runTime,
+ IOobject::MUST_READ
+ ),
+ false
+ );
+ mesh.init(true); // initialise all (lower levels and current)
+ Info<< endl;
+
+ // Like "setSystemMeshDictionaryIO.H"
+ IOobject dictIO = IOobject::selectIO
+ (
+ IOobject
+ (
+ "test-sample",
+ runTime.system(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE,
+ IOobject::NO_REGISTER
+ ),
+ args.getOrDefault("sample", "")
+ );
+
+ dictionary dictContents = IOdictionary(dictIO);
+ const dictionary* sampleDict = nullptr;
+
+ if (!dictContents.empty())
+ {
+ // Either have 'regular form' (from functionObjects)
+ // in which the sample surfaces are buried one layer deep
+ // or a flattened dictionary
+
+ if (dictContents.front()->dictPtr())
+ {
+ // Appears to be a dictionary of contents
+ // - get the first sub-dictionary with the correct "type"
+
+ for (const entry& e : dictContents)
+ {
+ const dictionary* dptr = e.dictPtr();
+
+ if
+ (
+ dptr
+ &&
+ (
+ sampledSurfaces::typeName
+ == dptr->getOrDefault("type", word::null)
+ )
+ )
+ {
+ sampleDict = dptr;
+ break;
+ }
+ }
+ }
+ else
+ {
+ // Probably a flattened dictionary,
+ // just check directly
+
+ if
+ (
+ sampledSurfaces::typeName
+ == dictContents.getOrDefault("type", word::null)
+ )
+ {
+ sampleDict = &dictContents;
+ }
+ }
+ }
+
+ if (!sampleDict)
+ {
+ FatalErrorInFunction
+ << "Dictionary does not appear to contain type:"
+ << sampledSurfaces::typeName << nl
+ << " " << dictIO.objectRelPath() << nl
+ << exit(FatalError);
+ }
+
+ // Construct from Time and dictionary, without loadFromFiles
+ sampledSurfaces sampling("test-sample", runTime, *sampleDict);
+
+ Info<< "Loaded " << sampling.size() << " surface samplers" << nl;
+
+ if (sampling.empty())
+ {
+ FatalErrorInFunction
+ << "No surface samplers loaded" << nl
+ << " " << dictIO.objectRelPath() << nl
+ << exit(FatalError);
+ }
+
+
+ // Manually read and load files
+
+ // Read objects in time directory
+ IOobjectList objects(mesh, runTime.timeName());
+
+ // Read GeometricFields
+ Info<< nl << "Load fields" << nl;
+
+ #if (OPENFOAM <= 2306)
+ // List of stored objects to clear after (as required)
+ wordHashSet allFields(objects.names());
+ LIFOStack storedObjects;
+
+ #undef ReadFields
+ #define ReadFields(FieldType) \
+ readFields(mesh, objects, allFields, storedObjects);
+ #else
+ // List of stored objects to clear after (as required)
+ DynamicList storedObjects;
+
+ #undef ReadFields
+ #define ReadFields(FieldType) \
+ readFields(mesh, objects, predicates::always{}, storedObjects);
+ #endif
+
+ // Read volFields
+ ReadFields(volScalarField);
+ ReadFields(volVectorField);
+ ReadFields(volSphericalTensorField);
+ ReadFields(volSymmTensorField);
+ ReadFields(volTensorField);
+
+ // Set fields to AUTO_WRITE (not really necessary for sampling...)
+ for (regIOobject* io : storedObjects)
+ {
+ io->writeOpt(IOobjectOption::AUTO_WRITE);
+ }
+
+ Info<< nl
+ << "Start " << nOutput << " times starting at "
+ << firstOutput << nl;
+
+ clockTime timing;
+
+ if (verbose) Info<< "Time:";
+
+ for
+ (
+ label timeIndex = firstOutput, count = 0;
+ count < nOutput;
+ ++timeIndex, ++count
+ )
+ {
+ runTime.setTime(timeIndex, timeIndex);
+ if (verbose) Info<< ' ' << runTime.timeName() << flush;
+ sampling.write();
+ }
+
+ if (verbose) Info<< nl;
+ Info<< nl << "Writing took "
+ << timing.timeIncrement() << "s" << endl;
+
+ //TBD profiling::writeNow();
+
+ // Info<< nl
+ // << "Cleanup newly generated files with" << nl << nl
+ // << " foamListTimes -rm -time "
+ // << firstOutput << ":" << nl
+ // << " foamListTimes -processor -rm -time "
+ // << firstOutput << ":" << nl;
+
+
+ Info<< "\nEnd\n" << endl;
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/surface-sampling/mydebugSurfaceWriter.C b/applications/test/surface-sampling/mydebugSurfaceWriter.C
new file mode 100644
index 0000000000..3a740aa697
--- /dev/null
+++ b/applications/test/surface-sampling/mydebugSurfaceWriter.C
@@ -0,0 +1,473 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2022-2023 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 .
+
+\*---------------------------------------------------------------------------*/
+
+#include "mydebugSurfaceWriter.H"
+#include "globalIndex.H"
+#include "argList.H"
+#include "OFstream.H"
+#include "OSspecific.H"
+#include "IOmanip.H"
+#include "Time.H"
+#include "pointIOField.H"
+#include "primitivePatch.H"
+#include "profiling.H"
+#include "surfaceWriterMethods.H"
+#include "PrecisionAdaptor.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceWriters
+{
+ defineTypeName(mydebugWriter);
+ addToRunTimeSelectionTable(surfaceWriter, mydebugWriter, word);
+ addToRunTimeSelectionTable(surfaceWriter, mydebugWriter, wordDict);
+}
+}
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Type narrowing - base implementation is pass-through
+template struct narrowType
+{
+ typedef Type type;
+};
+
+template<> struct narrowType
+{
+ typedef float type;
+};
+
+template<> struct narrowType>
+{
+ typedef Vector type;
+};
+
+template<> struct narrowType>
+{
+ typedef SphericalTensor type;
+};
+
+template<> struct narrowType>
+{
+ typedef SymmTensor type;
+};
+
+// FIXME: Not sure why this one seems to be broken...
+//
+// template<> struct narrowType>
+// {
+// typedef Tensor type;
+// };
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+template
+Foam::tmp>
+Foam::surfaceWriters::mydebugWriter::mergeField
+(
+ const Field& fld
+) const
+{
+ addProfiling(merge, "debugWriter::merge-field");
+
+ // This is largely identical to surfaceWriter::mergeField()
+ // but with narrowing for communication
+ if (narrowTransfer_ && parallel_ && UPstream::parRun())
+ {
+ // The narrowed type
+ typedef typename narrowType::type narrowedType;
+
+ // Ensure geometry is also merged
+ merge();
+
+ // Gather all values
+ auto tfield = tmp>::New();
+ auto& allFld = tfield.ref();
+
+ // Update any expired global index (as required)
+
+ const globalIndex& globIndex =
+ (
+ this->isPointData()
+ ? mergedSurf_.pointGlobalIndex()
+ : mergedSurf_.faceGlobalIndex()
+ );
+
+ ConstPrecisionAdaptor input(fld);
+ PrecisionAdaptor output(allFld);
+
+ globIndex.gather
+ (
+ input.cref(), // fld,
+ output.ref(), // allFld,
+ UPstream::msgType(),
+ commType_,
+ UPstream::worldComm
+ );
+
+ // Commit adapted content changes
+ input.commit();
+ output.commit();
+
+ // Discard adaptors
+ input.clear();
+ output.clear();
+
+ // Renumber (point data) to correspond to merged points
+ if
+ (
+ UPstream::master()
+ && this->isPointData()
+ && mergedSurf_.pointsMap().size()
+ )
+ {
+ inplaceReorder(mergedSurf_.pointsMap(), allFld);
+ allFld.resize(mergedSurf_.points().size());
+ }
+
+ return tfield;
+ }
+
+ return surfaceWriter::mergeField(fld);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::surfaceWriters::mydebugWriter::mydebugWriter()
+:
+ surfaceWriter(),
+ enableMerge_(true),
+ enableWrite_(false),
+ header_(true),
+ narrowTransfer_(false),
+ streamOpt_(IOstreamOption::BINARY)
+{}
+
+
+Foam::surfaceWriters::mydebugWriter::mydebugWriter
+(
+ const dictionary& options
+)
+:
+ surfaceWriter(options),
+ enableMerge_(options.getOrDefault("merge", true)),
+ enableWrite_(options.getOrDefault("write", false)),
+ header_(true),
+ narrowTransfer_(options.getOrDefault("narrow", false)),
+ streamOpt_(IOstreamOption::BINARY)
+{
+ Info<< "Using debug surface writer ("
+ << (this->isPointData() ? "point" : "face") << " data):"
+ << " commsType=" << UPstream::commsTypeNames[commType_]
+ << " merge=" << Switch::name(enableMerge_)
+ << " write=" << Switch::name(enableWrite_)
+ << " narrow=" << Switch::name(narrowTransfer_)
+ << endl;
+}
+
+
+Foam::surfaceWriters::mydebugWriter::mydebugWriter
+(
+ const meshedSurf& surf,
+ const fileName& outputPath,
+ bool parallel,
+ const dictionary& options
+)
+:
+ mydebugWriter(options)
+{
+ open(surf, outputPath, parallel);
+}
+
+
+Foam::surfaceWriters::mydebugWriter::mydebugWriter
+(
+ const pointField& points,
+ const faceList& faces,
+ const fileName& outputPath,
+ bool parallel,
+ const dictionary& options
+)
+:
+ mydebugWriter(options)
+{
+ open(points, faces, outputPath, parallel);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::surfaceWriters::mydebugWriter::serialWriteGeometry
+(
+ const regIOobject& iopts,
+ const meshedSurf& surf
+)
+{
+ const pointField& points = surf.points();
+ const faceList& faces = surf.faces();
+
+ if (verbose_)
+ {
+ if (this->isPointData())
+ {
+ Info<< "Writing points: " << iopts.objectPath() << endl;
+ }
+ else
+ {
+ Info<< "Writing face centres: " << iopts.objectPath() << endl;
+ }
+ }
+
+ // Like regIOobject::writeObject without instance() adaptation
+ // since this would write to e.g. 0/ instead of postProcessing/
+
+ autoPtr ppPtr;
+
+ {
+ OFstream os(iopts.objectPath(), streamOpt_);
+
+ if (header_)
+ {
+ iopts.writeHeader(os);
+ }
+
+ if (this->isPointData())
+ {
+ // Just like writeData, but without copying beforehand
+ os << points;
+ }
+ else
+ {
+ ppPtr.reset(new primitivePatch(SubList(faces), points));
+
+ // Just like writeData, but without copying beforehand
+ os << ppPtr().faceCentres();
+ }
+
+ if (header_)
+ {
+ IOobject::writeEndDivider(os);
+ }
+ }
+}
+
+
+Foam::fileName Foam::surfaceWriters::mydebugWriter::write()
+{
+ checkOpen();
+
+ // Geometry: rootdir/surfaceName/"points"
+ // Field: rootdir/surfaceName/