From a7ef33da6b93708e2891d5ba674583ae59978530 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Fri, 23 Sep 2022 19:20:57 +0200
Subject: [PATCH] ENH: add finite-area support to setFields (#2591)

- for example,

    defaultFieldValues
    (
        areaScalarFieldValue h 0.00014
    );

    regions
    (
        clipPlaneToFace
        {
            point  (0 0 0);
            normal (1 0 0);

            fieldValues
            (
                areaScalarFieldValue h 0.00015
            );
        }
    );

ENH: additional clipPlaneTo{Cell,Face,Point} topo sets

- less cumbersome than defining a semi-infinite bounding box
---
 .../preProcessing/setFields/Make/options      |   2 +
 .../preProcessing/setFields/setFields.C       | 734 +++++++++++++-----
 src/meshTools/Make/files                      |   3 +
 .../clipPlaneToCell/clipPlaneToCell.C         | 157 ++++
 .../clipPlaneToCell/clipPlaneToCell.H         | 172 ++++
 .../clipPlaneToFace/clipPlaneToFace.C         | 157 ++++
 .../clipPlaneToFace/clipPlaneToFace.H         | 172 ++++
 .../clipPlaneToPoint/clipPlaneToPoint.C       | 157 ++++
 .../clipPlaneToPoint/clipPlaneToPoint.H       | 172 ++++
 .../finiteArea/liquidFilmFoam/cylinder/Allrun |   2 +
 .../liquidFilmFoam/cylinder/Allrun-parallel   |   2 +
 .../cylinder/system/setFieldsDict             |  37 +
 .../laminar/fluidisedBed/system/setFieldsDict |   1 +
 13 files changed, 1571 insertions(+), 197 deletions(-)
 create mode 100644 src/meshTools/topoSet/cellSources/clipPlaneToCell/clipPlaneToCell.C
 create mode 100644 src/meshTools/topoSet/cellSources/clipPlaneToCell/clipPlaneToCell.H
 create mode 100644 src/meshTools/topoSet/faceSources/clipPlaneToFace/clipPlaneToFace.C
 create mode 100644 src/meshTools/topoSet/faceSources/clipPlaneToFace/clipPlaneToFace.H
 create mode 100644 src/meshTools/topoSet/pointSources/clipPlaneToPoint/clipPlaneToPoint.C
 create mode 100644 src/meshTools/topoSet/pointSources/clipPlaneToPoint/clipPlaneToPoint.H
 create mode 100644 tutorials/finiteArea/liquidFilmFoam/cylinder/system/setFieldsDict

diff --git a/applications/utilities/preProcessing/setFields/Make/options b/applications/utilities/preProcessing/setFields/Make/options
index 969020c4af..b5844533a1 100644
--- a/applications/utilities/preProcessing/setFields/Make/options
+++ b/applications/utilities/preProcessing/setFields/Make/options
@@ -1,8 +1,10 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfiniteArea \
     -lmeshTools \
     -lgenericPatchFields
diff --git a/applications/utilities/preProcessing/setFields/setFields.C b/applications/utilities/preProcessing/setFields/setFields.C
index 0af8c05ec0..a5c984cfad 100644
--- a/applications/utilities/preProcessing/setFields/setFields.C
+++ b/applications/utilities/preProcessing/setFields/setFields.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,214 +38,311 @@ Description
 #include "argList.H"
 #include "Time.H"
 #include "fvMesh.H"
+#include "faMesh.H"
 #include "topoSetSource.H"
 #include "cellSet.H"
 #include "faceSet.H"
 #include "volFields.H"
+#include "areaFields.H"
 
 using namespace Foam;
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Simple tuple of field type and name (read from stream)
+class fieldDescription
+{
+    word type_;
+    word name_;
+
+public:
+
+    const word& type() const noexcept { return type_; }
+    const word& name() const noexcept { return name_; }
+
+    explicit fieldDescription(Istream& is)
+    {
+        is >> type_;
+        is >> name_;
+
+        // Eg, read as "volScalarFieldValue", but change to "volScalarField"
+        if (type_.ends_with("Value"))
+        {
+            type_.erase(type_.size()-5);
+        }
+    }
+};
+
+
+// Consume unused field information
+template<class Type>
+bool consumeUnusedType(const fieldDescription& fieldDesc, Istream& is)
+{
+    typedef GeometricField<Type, faPatchField, areaMesh> fieldType1;
+    typedef GeometricField<Type, fvPatchField, volMesh>  fieldType2;
+    //? typedef GeometricField<Type, faePatchField, areaMesh> fieldType3;
+    //? typedef GeometricField<Type, fvsPatchField, volMesh>  fieldType4;
+
+    if
+    (
+        fieldDesc.type() == fieldType1::typeName
+     || fieldDesc.type() == fieldType2::typeName
+    )
+    {
+        (void) pTraits<Type>(is);
+        return true;
+    }
+
+    return false;
+}
+
+
+// Consume unused field information
+static bool consumeUnused(const fieldDescription& fieldDesc, Istream& is)
+{
+    return
+    (
+        consumeUnusedType<scalar>(fieldDesc, is)
+     || consumeUnusedType<vector>(fieldDesc, is)
+     || consumeUnusedType<sphericalTensor>(fieldDesc, is)
+     || consumeUnusedType<symmTensor>(fieldDesc, is)
+     || consumeUnusedType<tensor>(fieldDesc, is)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Setting volume fields
 template<class Type>
 bool setCellFieldType
 (
-    const word& fieldTypeDesc,
+    const fieldDescription& fieldDesc,
     const fvMesh& mesh,
     const labelList& selectedCells,
-    Istream& fieldValueStream
+    Istream& is
 )
 {
     typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
 
-    if (fieldTypeDesc != fieldType::typeName + "Value")
+    if (fieldDesc.type() != fieldType::typeName)
     {
         return false;
     }
 
-    word fieldName(fieldValueStream);
+    // Get value from stream
+    const Type fieldValue = pTraits<Type>(is);
+
 
     // Check the current time directory
     IOobject fieldHeader
     (
-        fieldName,
-        mesh.time().timeName(),
-        mesh,
+        fieldDesc.name(),
+        mesh.thisDb().time().timeName(),
+        mesh.thisDb(),
         IOobject::MUST_READ
     );
 
-    // Check the "constant" directory
-    if (!fieldHeader.typeHeaderOk<fieldType>(true))
+    bool found = fieldHeader.typeHeaderOk<fieldType>(true);
+
+    if (!found)
     {
+        // Fallback to "constant" directory
         fieldHeader = IOobject
         (
-            fieldName,
-            mesh.time().constant(),
-            mesh,
+            fieldDesc.name(),
+            mesh.thisDb().time().constant(),
+            mesh.thisDb(),
             IOobject::MUST_READ
         );
+        found = fieldHeader.typeHeaderOk<fieldType>(true);
     }
 
-    // Check field exists
-    if (fieldHeader.typeHeaderOk<fieldType>(true))
+    // Field exists
+    if (found)
     {
-        Info<< "    Setting internal values of "
+        Info<< "    - set internal values of "
             << fieldHeader.headerClassName()
-            << " " << fieldName << endl;
+            << ": " << fieldDesc.name()
+            << " = " << fieldValue << endl;
 
         fieldType field(fieldHeader, mesh, false);
 
-        const Type value = pTraits<Type>(fieldValueStream);
-
-        if (selectedCells.size() == field.size())
+        if (isNull(selectedCells) || selectedCells.size() == field.size())
         {
-            field.primitiveFieldRef() = value;
+            field.primitiveFieldRef() = fieldValue;
         }
         else
         {
-            forAll(selectedCells, celli)
+            for (const label celli : selectedCells)
             {
-                field[selectedCells[celli]] = value;
+                field[celli] = fieldValue;
             }
         }
 
-        typename GeometricField<Type, fvPatchField, volMesh>::
-            Boundary& fieldBf = field.boundaryFieldRef();
-
-        forAll(field.boundaryField(), patchi)
+        // Make boundary fields consistent - treat like zeroGradient
+        for (auto& pfld : field.boundaryFieldRef())
         {
-            fieldBf[patchi] = fieldBf[patchi].patchInternalField();
+            pfld = pfld.patchInternalField();
         }
 
         if (!field.write())
         {
             FatalErrorInFunction
-              << "Failed writing field " << fieldName << endl;
+                << "Failed writing field " << field.name() << endl;
         }
     }
     else
     {
-        WarningInFunction
-          << "Field " << fieldName << " not found" << endl;
-
-        // Consume value
-        (void)pTraits<Type>(fieldValueStream);
+        Warning
+            << "Field " << fieldDesc.name() << " not found" << endl;
     }
 
     return true;
 }
 
 
-class setCellField
-{
-
-public:
-
-    setCellField()
-    {}
-
-    autoPtr<setCellField> clone() const
-    {
-        return autoPtr<setCellField>::New();
-    }
-
-    class iNew
-    {
-        const fvMesh& mesh_;
-        const labelList& selectedCells_;
-
-    public:
-
-        iNew(const fvMesh& mesh, const labelList& selectedCells)
-        :
-            mesh_(mesh),
-            selectedCells_(selectedCells)
-        {}
-
-        iNew(const fvMesh& mesh, labelList&& selectedCells)
-        :
-            mesh_(mesh),
-            selectedCells_(std::move(selectedCells))
-        {}
-
-        autoPtr<setCellField> operator()(Istream& fieldValues) const
-        {
-            word fieldType(fieldValues);
-
-            if
-            (
-               !(
-                    setCellFieldType<scalar>
-                        (fieldType, mesh_, selectedCells_, fieldValues)
-                 || setCellFieldType<vector>
-                        (fieldType, mesh_, selectedCells_, fieldValues)
-                 || setCellFieldType<sphericalTensor>
-                        (fieldType, mesh_, selectedCells_, fieldValues)
-                 || setCellFieldType<symmTensor>
-                        (fieldType, mesh_, selectedCells_, fieldValues)
-                 || setCellFieldType<tensor>
-                        (fieldType, mesh_, selectedCells_, fieldValues)
-                )
-            )
-            {
-                WarningInFunction
-                    << "field type " << fieldType << " not currently supported"
-                    << endl;
-            }
-
-            return autoPtr<setCellField>::New();
-        }
-    };
-};
-
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Setting finite-area face fields
 template<class Type>
-bool setFaceFieldType
+bool setAreaFieldType
 (
-    const word& fieldTypeDesc,
-    const fvMesh& mesh,
+    const fieldDescription& fieldDesc,
+    const faMesh& mesh,
     const labelList& selectedFaces,
-    Istream& fieldValueStream
+    Istream& is
 )
 {
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+    typedef GeometricField<Type, faPatchField, areaMesh> fieldType;
 
-    if (fieldTypeDesc != fieldType::typeName + "Value")
+    if (fieldDesc.type() != fieldType::typeName)
     {
         return false;
     }
 
-    word fieldName(fieldValueStream);
+    // Get value from stream
+    const Type fieldValue = pTraits<Type>(is);
+
 
     // Check the current time directory
     IOobject fieldHeader
     (
-        fieldName,
-        mesh.time().timeName(),
-        mesh,
+        fieldDesc.name(),
+        mesh.thisDb().time().timeName(),
+        mesh.thisDb(),
         IOobject::MUST_READ
     );
 
-    // Check the "constant" directory
-    if (!fieldHeader.typeHeaderOk<fieldType>(true))
+    bool found = fieldHeader.typeHeaderOk<fieldType>(true);
+
+    if (!found)
     {
+        // Fallback to "constant" directory
         fieldHeader = IOobject
         (
-            fieldName,
-            mesh.time().constant(),
-            mesh,
+            fieldDesc.name(),
+            mesh.thisDb().time().constant(),
+            mesh.thisDb(),
             IOobject::MUST_READ
         );
+        found = fieldHeader.typeHeaderOk<fieldType>(true);
     }
 
-    // Check field exists
-    if (fieldHeader.typeHeaderOk<fieldType>(true))
+    // Field exists
+    if (found)
     {
-        Info<< "    Setting patchField values of "
+        Info<< "    - set internal values of "
             << fieldHeader.headerClassName()
-            << " " << fieldName << endl;
+            << ": " << fieldDesc.name()
+            << " = " << fieldValue << endl;
 
         fieldType field(fieldHeader, mesh);
 
-        const Type value = pTraits<Type>(fieldValueStream);
+        if (isNull(selectedFaces) || selectedFaces.size() == field.size())
+        {
+            field.primitiveFieldRef() = fieldValue;
+        }
+        else
+        {
+            for (const label facei : selectedFaces)
+            {
+                field[facei] = fieldValue;
+            }
+        }
+
+        if (!field.write())
+        {
+            FatalErrorInFunction
+                << "Failed writing field " << field.name() << endl;
+        }
+    }
+    else
+    {
+        Warning
+            << "Field " << fieldDesc.name() << " not found" << endl;
+    }
+
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Setting volume boundary fields
+template<class Type>
+bool setFaceFieldType
+(
+    const fieldDescription& fieldDesc,
+    const fvMesh& mesh,
+    const labelList& selectedFaces,
+    Istream& is
+)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if (fieldDesc.type() != fieldType::typeName)
+    {
+        return false;
+    }
+
+    // Get value from stream
+    const Type fieldValue = pTraits<Type>(is);
+
+
+    // Check the current time directory
+    IOobject fieldHeader
+    (
+        fieldDesc.name(),
+        mesh.thisDb().time().timeName(),
+        mesh.thisDb(),
+        IOobject::MUST_READ
+    );
+
+    bool found = fieldHeader.typeHeaderOk<fieldType>(true);
+
+    if (!found)
+    {
+        // Fallback to "constant" directory
+        fieldHeader = IOobject
+        (
+            fieldDesc.name(),
+            mesh.thisDb().time().constant(),
+            mesh.thisDb(),
+            IOobject::MUST_READ
+        );
+        found = fieldHeader.typeHeaderOk<fieldType>(true);
+    }
+
+    // Field exists
+    if (found)
+    {
+        Info<< "    - set boundary values of "
+            << fieldHeader.headerClassName()
+            << ": " << fieldDesc.name()
+            << " = " << fieldValue << endl;
+
+        fieldType field(fieldHeader, mesh);
 
         // Create flat list of selected faces and their value.
         Field<Type> allBoundaryValues(mesh.nBoundaryFaces());
@@ -260,34 +358,47 @@ bool setFaceFieldType
         }
 
         // Override
-        bool hasWarned = false;
+        unsigned hasWarned = 0;
         labelList nChanged
         (
             returnReduce(field.boundaryField().size(), maxOp<label>()),
-            0
+            Zero
         );
-        forAll(selectedFaces, i)
+
+        for (const label facei : selectedFaces)
         {
-            label facei = selectedFaces[i];
-            if (mesh.isInternalFace(facei))
+            const label bFacei = facei-mesh.nInternalFaces();
+
+            if (bFacei < 0)
             {
-                if (!hasWarned)
+                if (!(hasWarned & 1))
                 {
-                    hasWarned = true;
+                    hasWarned |= 1;
                     WarningInFunction
                         << "Ignoring internal face " << facei
                         << ". Suppressing further warnings." << endl;
                 }
             }
+            else if (bFacei >= mesh.nBoundaryFaces())
+            {
+                if (!(hasWarned & 2))
+                {
+                    hasWarned |= 2;
+                    WarningInFunction
+                        << "Ignoring out-of-range face " << facei
+                        << ". Suppressing further warnings." << endl;
+                }
+            }
             else
             {
-                label bFacei = facei-mesh.nInternalFaces();
-                allBoundaryValues[bFacei] = value;
-                nChanged[mesh.boundaryMesh().patchID()[bFacei]]++;
+                label patchi = mesh.boundaryMesh().patchID()[bFacei];
+
+                allBoundaryValues[bFacei] = fieldValue;
+                ++nChanged[patchi];
             }
         }
 
-        Pstream::listCombineAllGather(nChanged, plusEqOp<label>());
+        Pstream::listCombineReduce(nChanged, plusEqOp<label>());
 
         auto& fieldBf = field.boundaryFieldRef();
 
@@ -299,6 +410,7 @@ bool setFaceFieldType
                 Info<< "    On patch "
                     << field.boundaryField()[patchi].patch().name()
                     << " set " << nChanged[patchi] << " values" << endl;
+
                 fieldBf[patchi] == SubField<Type>
                 (
                     allBoundaryValues,
@@ -312,85 +424,222 @@ bool setFaceFieldType
         if (!field.write())
         {
             FatalErrorInFunction
-                << "Failed writing field " << field.name() << exit(FatalError);
+                << "Failed writing field " << field.name() << endl;
         }
     }
     else
     {
-        WarningInFunction
-          << "Field " << fieldName << " not found" << endl;
-
-        // Consume value
-        (void)pTraits<Type>(fieldValueStream);
+        Warning
+            << "Field " << fieldDesc.name() << " not found" << endl;
     }
 
     return true;
 }
 
 
-class setFaceField
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Dispatcher for setting volume fields
+struct setCellField
 {
+    autoPtr<setCellField> clone() const { return nullptr; }  // placeholder
 
-public:
-
-    setFaceField()
-    {}
-
-    autoPtr<setFaceField> clone() const
+    static bool apply
+    (
+        const fieldDescription& fieldDesc,
+        const fvMesh& m,
+        const labelList& selectedCells,
+        Istream& is
+    )
     {
-        return autoPtr<setFaceField>::New();
+        return
+        (
+            setCellFieldType<scalar>(fieldDesc, m, selectedCells, is)
+         || setCellFieldType<vector>(fieldDesc, m, selectedCells, is)
+         || setCellFieldType<sphericalTensor>(fieldDesc, m, selectedCells, is)
+         || setCellFieldType<symmTensor>(fieldDesc, m, selectedCells, is)
+         || setCellFieldType<tensor>(fieldDesc, m, selectedCells, is)
+        );
     }
 
     class iNew
     {
         const fvMesh& mesh_;
-        const labelList& selectedFaces_;
+        const labelList& selected_;
+
+    public:
+
+        iNew(const fvMesh& mesh, const labelList& selectedCells)
+        :
+            mesh_(mesh),
+            selected_(selectedCells)
+        {}
+
+        autoPtr<setCellField> operator()(Istream& is) const
+        {
+            const fieldDescription fieldDesc(is);
+
+            bool ok = setCellField::apply(fieldDesc, mesh_, selected_, is);
+
+            if (!ok)
+            {
+                ok = consumeUnused(fieldDesc, is);
+
+                if (ok)
+                {
+                    // Not meant for us
+                    Info<< "Skip " << fieldDesc.type()
+                        << " for finite-volume" << nl;
+                }
+                else
+                {
+                    WarningInFunction
+                        << "Unsupported field type: "
+                        << fieldDesc.type() << endl;
+                }
+            }
+
+            return nullptr;  // Irrelevant return value
+        }
+    };
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Dispatcher for setting volume boundary fields
+struct setFaceField
+{
+    autoPtr<setFaceField> clone() const { return nullptr; }  // placeholder
+
+    static bool apply
+    (
+        const fieldDescription& fieldDesc,
+        const fvMesh& m,
+        const labelList& selectedFaces,
+        Istream& is
+    )
+    {
+        return
+        (
+            setFaceFieldType<scalar>(fieldDesc, m, selectedFaces, is)
+         || setFaceFieldType<vector>(fieldDesc, m, selectedFaces, is)
+         || setFaceFieldType<sphericalTensor>(fieldDesc, m, selectedFaces, is)
+         || setFaceFieldType<symmTensor>(fieldDesc, m, selectedFaces, is)
+         || setFaceFieldType<tensor>(fieldDesc, m, selectedFaces, is)
+        );
+    }
+
+    class iNew
+    {
+        const fvMesh& mesh_;
+        const labelList& selected_;
 
     public:
 
         iNew(const fvMesh& mesh, const labelList& selectedFaces)
         :
             mesh_(mesh),
-            selectedFaces_(selectedFaces)
+            selected_(selectedFaces)
         {}
 
-        iNew(const fvMesh& mesh, labelList&& selectedFaces)
-        :
-            mesh_(mesh),
-            selectedFaces_(std::move(selectedFaces))
-        {}
-
-        autoPtr<setFaceField> operator()(Istream& fieldValues) const
+        autoPtr<setFaceField> operator()(Istream& is) const
         {
-            word fieldType(fieldValues);
+            const fieldDescription fieldDesc(is);
 
-            if
-            (
-               !(
-                    setFaceFieldType<scalar>
-                        (fieldType, mesh_, selectedFaces_, fieldValues)
-                 || setFaceFieldType<vector>
-                        (fieldType, mesh_, selectedFaces_, fieldValues)
-                 || setFaceFieldType<sphericalTensor>
-                        (fieldType, mesh_, selectedFaces_, fieldValues)
-                 || setFaceFieldType<symmTensor>
-                        (fieldType, mesh_, selectedFaces_, fieldValues)
-                 || setFaceFieldType<tensor>
-                        (fieldType, mesh_, selectedFaces_, fieldValues)
-                )
-            )
+            bool ok = setFaceField::apply(fieldDesc, mesh_, selected_, is);
+
+            if (!ok)
             {
-                WarningInFunction
-                    << "field type " << fieldType << " not currently supported"
-                    << endl;
+                ok = consumeUnused(fieldDesc, is);
+
+                if (ok)
+                {
+                    // Not meant for us
+                    Info<< "Skip " << fieldDesc.type()
+                        << " for finite-volume" << nl;
+                }
+                else
+                {
+                    WarningInFunction
+                        << "Unsupported field type: "
+                        << fieldDesc.type() << endl;
+                }
             }
 
-            return autoPtr<setFaceField>::New();
+            return nullptr;  // Irrelevant return value
         }
     };
 };
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Dispatcher for setting area face fields
+struct setAreaField
+{
+    autoPtr<setAreaField> clone() const { return nullptr; }  // placeholder
+
+    static bool apply
+    (
+        const fieldDescription& fieldDesc,
+        const faMesh& m,
+        const labelList& selectedFaces,
+        Istream& is
+    )
+    {
+        return
+        (
+            setAreaFieldType<scalar>(fieldDesc, m, selectedFaces, is)
+         || setAreaFieldType<vector>(fieldDesc, m, selectedFaces, is)
+         || setAreaFieldType<sphericalTensor>(fieldDesc, m, selectedFaces, is)
+         || setAreaFieldType<symmTensor>(fieldDesc, m, selectedFaces, is)
+         || setAreaFieldType<tensor>(fieldDesc, m, selectedFaces, is)
+        );
+    }
+
+    class iNew
+    {
+        const faMesh& mesh_;
+        const labelList& selected_;
+
+    public:
+
+        iNew(const faMesh& mesh, const labelList& selectedFaces)
+        :
+            mesh_(mesh),
+            selected_(selectedFaces)
+        {}
+
+        autoPtr<setAreaField> operator()(Istream& is) const
+        {
+            const fieldDescription fieldDesc(is);
+
+            bool ok = setAreaField::apply(fieldDesc, mesh_, selected_, is);
+
+            if (!ok)
+            {
+                ok = consumeUnused(fieldDesc, is);
+
+                if (ok)
+                {
+                    // Not meant for us
+                    Info<< "Skip " << fieldDesc.type()
+                        << " for finite-volume" << nl;
+                }
+                else
+                {
+                    WarningInFunction
+                        << "Unsupported field type: "
+                        << fieldDesc.type() << endl;
+                }
+            }
+
+            return nullptr;  // Irrelevant return value
+        }
+    };
+};
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -403,11 +652,28 @@ int main(int argc, char *argv[])
 
     argList::addOption("dict", "file", "Alternative setFieldsDict");
 
+    argList::addBoolOption
+    (
+        "no-finite-area",
+        "Suppress handling of finite-area mesh/fields"
+    );
+
     #include "addRegionOption.H"
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createNamedMesh.H"
 
+    autoPtr<faMesh> faMeshPtr;
+
+    if (!args.found("no-finite-area"))
+    {
+        faMeshPtr = faMesh::TryNew(mesh);
+    }
+    if (faMeshPtr)
+    {
+        Info<< "Detected finite-area mesh" << nl;
+    }
+
     const word dictName("setFieldsDict");
     #include "setSystemMeshDictionaryIO.H"
 
@@ -415,70 +681,144 @@ int main(int argc, char *argv[])
 
     IOdictionary setFieldsDict(dictIO);
 
-    if (setFieldsDict.found("defaultFieldValues"))
+
+    // Note.
+    // The PtrList + iNew mechanism is used to trigger the callbacks
+    // and perform the desired actions. The contents of the PtrList
+    // itself are actually irrelevant.
+
+
+    // Default field values
     {
-        Info<< "Setting field default values" << endl;
-        PtrList<setCellField> defaultFieldValues
-        (
-            setFieldsDict.lookup("defaultFieldValues"),
-            setCellField::iNew(mesh, labelList(mesh.nCells()))
-        );
-        Info<< endl;
+        const entry* eptr =
+            setFieldsDict.findEntry("defaultFieldValues", keyType::LITERAL);
+
+        if (eptr)
+        {
+            ITstream& is = eptr->stream();
+
+            Info<< "Setting volume field default values" << endl;
+
+            PtrList<setCellField> defaultFieldValues
+            (
+                is,
+                setCellField::iNew(mesh, labelList::null())
+            );
+
+            if (faMeshPtr)
+            {
+                const faMesh& areaMesh = faMeshPtr();
+                is.rewind();
+
+                Info<< "Setting area field default values" << endl;
+
+                PtrList<setAreaField> defaultFieldValues
+                (
+                    is,
+                    setAreaField::iNew(areaMesh, labelList::null())
+                );
+            }
+            Info<< endl;
+        }
     }
 
 
-    Info<< "Setting field region values" << endl;
+    Info<< "Setting field region values" << nl << endl;
 
     PtrList<entry> regions(setFieldsDict.lookup("regions"));
 
-    forAll(regions, regionI)
+    for (const entry& region : regions)
     {
-        const entry& region = regions[regionI];
-
         autoPtr<topoSetSource> source =
             topoSetSource::New(region.keyword(), mesh, region.dict());
 
         if (source().setType() == topoSetSource::CELLSET_SOURCE)
         {
-            cellSet selectedCellSet
+            cellSet subset
             (
                 mesh,
                 "cellSet",
                 mesh.nCells()/10+1  // Reasonable size estimate.
             );
 
-            source->applyToSet
-            (
-                topoSetSource::NEW,
-                selectedCellSet
-            );
+            source->applyToSet(topoSetSource::NEW, subset);
+
+            labelList selectedCells(subset.sortedToc());
+
+            Info<< "    Selected "
+                << returnReduce(selectedCells.size(), sumOp<label>())
+                << '/'
+                << returnReduce(mesh.nCells(), sumOp<label>())
+                << " cells" << nl;
+
+            ITstream& is = region.dict().lookup("fieldValues");
 
             PtrList<setCellField> fieldValues
             (
-                region.dict().lookup("fieldValues"),
-                setCellField::iNew(mesh, selectedCellSet.sortedToc())
+                is,
+                setCellField::iNew(mesh, selectedCells)
             );
         }
         else if (source().setType() == topoSetSource::FACESET_SOURCE)
         {
-            faceSet selectedFaceSet
+            faceSet subset
             (
                 mesh,
                 "faceSet",
                 mesh.nBoundaryFaces()/10+1
             );
 
-            source->applyToSet
-            (
-                topoSetSource::NEW,
-                selectedFaceSet
-            );
+            source->applyToSet(topoSetSource::NEW, subset);
+
+            labelList selectedFaces(subset.sortedToc());
+
+            Info<< "    Selected " << selectedFaces.size()
+                << " faces" << nl;
+
+            ITstream& is = region.dict().lookup("fieldValues");
 
             PtrList<setFaceField> fieldValues
             (
-                region.dict().lookup("fieldValues"),
-                setFaceField::iNew(mesh, selectedFaceSet.sortedToc())
+                is,
+                setFaceField::iNew(mesh, selectedFaces)
             );
+
+            if (faMeshPtr)
+            {
+                const faMesh& areaMesh = faMeshPtr();
+
+                const labelUList& faceLabels = areaMesh.faceLabels();
+
+                // Transcribe from mesh faces to finite-area addressing
+                labelList areaFaces(faceLabels.size());
+
+                label nUsed = 0;
+                forAll(faceLabels, facei)
+                {
+                    const label meshFacei = faceLabels[facei];
+
+                    if (subset.test(meshFacei))
+                    {
+                        areaFaces[nUsed] = facei;
+                        ++nUsed;
+                    }
+                }
+                areaFaces.resize(nUsed);
+
+                Info<< "    Selected "
+                    << returnReduce(areaFaces.size(), sumOp<label>())
+                    << '/'
+                    << returnReduce(faceLabels.size(), sumOp<label>())
+                    << " area faces" << nl;
+
+                is.rewind();
+
+                PtrList<setAreaField> fieldValues
+                (
+                    is,
+                    setAreaField::iNew(areaMesh, areaFaces)
+                );
+            }
         }
     }
 
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index 7713661d99..4194866ded 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -163,6 +163,7 @@ $(cellSources)/topoSetCellSource/topoSetCellSource.C
 $(cellSources)/boundaryToCell/boundaryToCell.C
 $(cellSources)/boxToCell/boxToCell.C
 $(cellSources)/cellToCell/cellToCell.C
+$(cellSources)/clipPlaneToCell/clipPlaneToCell.C
 $(cellSources)/cylinderAnnulusToCell/cylinderAnnulusToCell.C
 $(cellSources)/cylinderToCell/cylinderToCell.C
 $(cellSources)/faceToCell/faceToCell.C
@@ -188,6 +189,7 @@ $(faceSources)/topoSetFaceSource/topoSetFaceSource.C
 $(faceSources)/boundaryToFace/boundaryToFace.C
 $(faceSources)/boxToFace/boxToFace.C
 $(faceSources)/cellToFace/cellToFace.C
+$(faceSources)/clipPlaneToFace/clipPlaneToFace.C
 $(faceSources)/cylinderAnnulusToFace/cylinderAnnulusToFace.C
 $(faceSources)/cylinderToFace/cylinderToFace.C
 $(faceSources)/faceToFace/faceToFace.C
@@ -205,6 +207,7 @@ pointSources = topoSet/pointSources
 $(pointSources)/topoSetPointSource/topoSetPointSource.C
 $(pointSources)/boxToPoint/boxToPoint.C
 $(pointSources)/cellToPoint/cellToPoint.C
+$(pointSources)/clipPlaneToPoint/clipPlaneToPoint.C
 $(pointSources)/cylinderToPoint/cylinderToPoint.C
 $(pointSources)/faceToPoint/faceToPoint.C
 $(pointSources)/labelToPoint/labelToPoint.C
diff --git a/src/meshTools/topoSet/cellSources/clipPlaneToCell/clipPlaneToCell.C b/src/meshTools/topoSet/cellSources/clipPlaneToCell/clipPlaneToCell.C
new file mode 100644
index 0000000000..521ba39942
--- /dev/null
+++ b/src/meshTools/topoSet/cellSources/clipPlaneToCell/clipPlaneToCell.C
@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "clipPlaneToCell.H"
+#include "polyMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(clipPlaneToCell, 0);
+    addToRunTimeSelectionTable(topoSetSource, clipPlaneToCell, word);
+    addToRunTimeSelectionTable(topoSetSource, clipPlaneToCell, istream);
+    addToRunTimeSelectionTable(topoSetCellSource, clipPlaneToCell, word);
+    addToRunTimeSelectionTable(topoSetCellSource, clipPlaneToCell, istream);
+    addNamedToRunTimeSelectionTable
+    (
+        topoSetCellSource,
+        clipPlaneToCell,
+        word,
+        clipPlane
+    );
+    addNamedToRunTimeSelectionTable
+    (
+        topoSetCellSource,
+        clipPlaneToCell,
+        istream,
+        clipPlane
+    );
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::clipPlaneToCell::usage_
+(
+    clipPlaneToCell::typeName,
+    "\n    Usage: clipPlaneToCell (px py pz) (nx ny nz)\n\n"
+    "    Select cells with centres above the plane\n\n"
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::clipPlaneToCell::combine(topoSet& set, const bool add) const
+{
+    // Cell centres above the plane
+
+    const pointField& ctrs = mesh_.cellCentres();
+
+    forAll(ctrs, elemi)
+    {
+        if (((ctrs[elemi] - point_) & normal_) > 0)
+        {
+            addOrDelete(set, elemi, add);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::clipPlaneToCell::clipPlaneToCell
+(
+    const polyMesh& mesh,
+    const point& basePoint,
+    const vector& normal
+)
+:
+    topoSetCellSource(mesh),
+    point_(basePoint),
+    normal_(normal)
+{}
+
+
+Foam::clipPlaneToCell::clipPlaneToCell
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    clipPlaneToCell
+    (
+        mesh,
+        dict.get<vector>("point"),
+        dict.get<vector>("normal")
+    )
+{}
+
+
+Foam::clipPlaneToCell::clipPlaneToCell
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetCellSource(mesh),
+    point_(checkIs(is)),
+    normal_(checkIs(is))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::clipPlaneToCell::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (action == topoSetSource::ADD || action == topoSetSource::NEW)
+    {
+        if (verbose_)
+        {
+            Info<< "    Adding cells with centres above plane at "
+                << point_ << " with normal " << normal_ << endl;
+        }
+
+        combine(set, true);
+    }
+    else if (action == topoSetSource::SUBTRACT)
+    {
+        if (verbose_)
+        {
+            Info<< "    Removing cells with centres above plane at "
+                << point_ << " with normal " << normal_ << endl;
+        }
+
+        combine(set, false);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/topoSet/cellSources/clipPlaneToCell/clipPlaneToCell.H b/src/meshTools/topoSet/cellSources/clipPlaneToCell/clipPlaneToCell.H
new file mode 100644
index 0000000000..544a82e2c0
--- /dev/null
+++ b/src/meshTools/topoSet/cellSources/clipPlaneToCell/clipPlaneToCell.H
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::clipPlaneToCell
+
+Description
+    A \c topoSetCellSource to select all cells
+    whose cell centre is above the specified plane
+
+    Operands:
+    \table
+      Operand   | Type    | Location
+      output    | cellSet | $FOAM_CASE/constant/polyMesh/sets/\<set\>
+    \endtable
+
+Usage
+    Minimal example by using \c system/topoSetDict.actions:
+    \verbatim
+    {
+        // Mandatory (inherited) entries
+        name        <name>;
+        type        cellSet;
+        action      <action>;
+
+        // Mandatory entries
+        source      clipPlaneToCell;
+
+        point  (x y z);
+        normal (x y z);
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property   | Description                         | Type | Req'd | Dflt
+      name       | Name of cellSet                     | word |  yes  | -
+      type       | Type name: cellSet                  | word |  yes  | -
+      action     | Action applied on cells - see below | word |  yes  | -
+      source     | Source name: clipPlaneToCell        | word |  yes  | -
+    \endtable
+
+    Options for the \c action entry:
+    \verbatim
+      new      | Create a new cellSet from selected cells
+      add      | Add selected cells into this cellSet
+      subtract | Remove selected cells from this cellSet
+    \endverbatim
+
+    Options for the conditional mandatory entries:
+    \verbatim
+      Entry    | Description                        | Type   | Req'd  | Dflt
+      point    | The reference point for the plane  | vector | yes    | -
+      normal   | The normal of the plane            | vector | yes    | -
+    \endverbatim
+
+See also
+    - Foam::topoSetSource
+    - Foam::topoSetCellSource
+
+SourceFiles
+    clipPlaneToCell.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_clipPlaneToCell_H
+#define Foam_clipPlaneToCell_H
+
+#include "topoSetCellSource.H"
+#include "point.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class clipPlaneToCell Declaration
+\*---------------------------------------------------------------------------*/
+
+class clipPlaneToCell
+:
+    public topoSetCellSource
+{
+    // Private Data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Point on the plane
+        const vector point_;
+
+        //- Normal to the plane
+        const vector normal_;
+
+
+    // Private Member Functions
+
+        void combine(topoSet& set, const bool add) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("clipPlaneToCell");
+
+
+    // Constructors
+
+        //- No default construct
+        clipPlaneToCell() = delete;
+
+        //- Construct from components
+        clipPlaneToCell
+        (
+            const polyMesh& mesh,
+            const point& basePoint,
+            const vector& normal
+        );
+
+        //- Construct from dictionary
+        clipPlaneToCell(const polyMesh& mesh, const dictionary& dict);
+
+        //- Construct from stream - a single box.
+        clipPlaneToCell(const polyMesh& mesh, Istream& is);
+
+
+    //- Destructor
+    virtual ~clipPlaneToCell() = default;
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet& set
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/topoSet/faceSources/clipPlaneToFace/clipPlaneToFace.C b/src/meshTools/topoSet/faceSources/clipPlaneToFace/clipPlaneToFace.C
new file mode 100644
index 0000000000..004a4d859a
--- /dev/null
+++ b/src/meshTools/topoSet/faceSources/clipPlaneToFace/clipPlaneToFace.C
@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "clipPlaneToFace.H"
+#include "polyMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(clipPlaneToFace, 0);
+    addToRunTimeSelectionTable(topoSetSource, clipPlaneToFace, word);
+    addToRunTimeSelectionTable(topoSetSource, clipPlaneToFace, istream);
+    addToRunTimeSelectionTable(topoSetFaceSource, clipPlaneToFace, word);
+    addToRunTimeSelectionTable(topoSetFaceSource, clipPlaneToFace, istream);
+    addNamedToRunTimeSelectionTable
+    (
+        topoSetFaceSource,
+        clipPlaneToFace,
+        word,
+        clipPlane
+    );
+    addNamedToRunTimeSelectionTable
+    (
+        topoSetFaceSource,
+        clipPlaneToFace,
+        istream,
+        clipPlane
+    );
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::clipPlaneToFace::usage_
+(
+    clipPlaneToFace::typeName,
+    "\n    Usage: clipPlaneToFace (px py pz) (nx ny nz)\n\n"
+    "    Select faces with centres above the plane\n\n"
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::clipPlaneToFace::combine(topoSet& set, const bool add) const
+{
+    // Face centres above the plane
+
+    const pointField& ctrs = mesh_.faceCentres();
+
+    forAll(ctrs, elemi)
+    {
+        if (((ctrs[elemi] - point_) & normal_) > 0)
+        {
+            addOrDelete(set, elemi, add);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::clipPlaneToFace::clipPlaneToFace
+(
+    const polyMesh& mesh,
+    const point& basePoint,
+    const vector& normal
+)
+:
+    topoSetFaceSource(mesh),
+    point_(basePoint),
+    normal_(normal)
+{}
+
+
+Foam::clipPlaneToFace::clipPlaneToFace
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    clipPlaneToFace
+    (
+        mesh,
+        dict.get<vector>("point"),
+        dict.get<vector>("normal")
+    )
+{}
+
+
+Foam::clipPlaneToFace::clipPlaneToFace
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetFaceSource(mesh),
+    point_(checkIs(is)),
+    normal_(checkIs(is))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::clipPlaneToFace::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (action == topoSetSource::ADD || action == topoSetSource::NEW)
+    {
+        if (verbose_)
+        {
+            Info<< "    Adding faces with centres above plane at "
+                << point_ << " with normal " << normal_ << endl;
+        }
+
+        combine(set, true);
+    }
+    else if (action == topoSetSource::SUBTRACT)
+    {
+        if (verbose_)
+        {
+            Info<< "    Removing faces with centres above plane at "
+                << point_ << " with normal " << normal_ << endl;
+        }
+
+        combine(set, false);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/topoSet/faceSources/clipPlaneToFace/clipPlaneToFace.H b/src/meshTools/topoSet/faceSources/clipPlaneToFace/clipPlaneToFace.H
new file mode 100644
index 0000000000..1842fce072
--- /dev/null
+++ b/src/meshTools/topoSet/faceSources/clipPlaneToFace/clipPlaneToFace.H
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::clipPlaneToFace
+
+Description
+    A \c topoSetFaceSource to select all faces
+    whose face centre is above the specified plane
+
+    Operands:
+    \table
+      Operand   | Type    | Location
+      output    | faceSet | $FOAM_CASE/constant/polyMesh/sets/\<set\>
+    \endtable
+
+Usage
+    Minimal example by using \c system/topoSetDict.actions:
+    \verbatim
+    {
+        // Mandatory (inherited) entries
+        name        <name>;
+        type        faceSet;
+        action      <action>;
+
+        // Mandatory entries
+        source      clipPlaneToFace;
+
+        point  (x y z);
+        normal (x y z);
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property   | Description                         | Type | Req'd | Dflt
+      name       | Name of faceSet                     | word |  yes  | -
+      type       | Type name: faceSet                  | word |  yes  | -
+      action     | Action applied on faces - see below | word |  yes  | -
+      source     | Source name: clipPlaneToFace        | word |  yes  | -
+    \endtable
+
+    Options for the \c action entry:
+    \verbatim
+      new      | Create a new faceSet from selected faces
+      add      | Add selected faces into this faceSet
+      subtract | Remove selected faces from this faceSet
+    \endverbatim
+
+    Options for the conditional mandatory entries:
+    \verbatim
+      Entry    | Description                        | Type   | Req'd  | Dflt
+      point    | The reference point for the plane  | vector | yes    | -
+      normal   | The normal of the plane            | vector | yes    | -
+    \endverbatim
+
+See also
+    - Foam::topoSetSource
+    - Foam::topoSetFaceSource
+
+SourceFiles
+    clipPlaneToFace.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_clipPlaneToFace_H
+#define Foam_clipPlaneToFace_H
+
+#include "topoSetFaceSource.H"
+#include "point.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class clipPlaneToFace Declaration
+\*---------------------------------------------------------------------------*/
+
+class clipPlaneToFace
+:
+    public topoSetFaceSource
+{
+    // Private Data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Point on the plane
+        const vector point_;
+
+        //- Normal to the plane
+        const vector normal_;
+
+
+    // Private Member Functions
+
+        void combine(topoSet& set, const bool add) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("clipPlaneToFace");
+
+
+    // Constructors
+
+        //- No default construct
+        clipPlaneToFace() = delete;
+
+        //- Construct from components
+        clipPlaneToFace
+        (
+            const polyMesh& mesh,
+            const point& basePoint,
+            const vector& normal
+        );
+
+        //- Construct from dictionary
+        clipPlaneToFace(const polyMesh& mesh, const dictionary& dict);
+
+        //- Construct from stream - a single box.
+        clipPlaneToFace(const polyMesh& mesh, Istream& is);
+
+
+    //- Destructor
+    virtual ~clipPlaneToFace() = default;
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet& set
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/topoSet/pointSources/clipPlaneToPoint/clipPlaneToPoint.C b/src/meshTools/topoSet/pointSources/clipPlaneToPoint/clipPlaneToPoint.C
new file mode 100644
index 0000000000..336983da34
--- /dev/null
+++ b/src/meshTools/topoSet/pointSources/clipPlaneToPoint/clipPlaneToPoint.C
@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "clipPlaneToPoint.H"
+#include "polyMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(clipPlaneToPoint, 0);
+    addToRunTimeSelectionTable(topoSetSource, clipPlaneToPoint, word);
+    addToRunTimeSelectionTable(topoSetSource, clipPlaneToPoint, istream);
+    addToRunTimeSelectionTable(topoSetPointSource, clipPlaneToPoint, word);
+    addToRunTimeSelectionTable(topoSetPointSource, clipPlaneToPoint, istream);
+    addNamedToRunTimeSelectionTable
+    (
+        topoSetPointSource,
+        clipPlaneToPoint,
+        word,
+        clipPlane
+    );
+    addNamedToRunTimeSelectionTable
+    (
+        topoSetPointSource,
+        clipPlaneToPoint,
+        istream,
+        clipPlane
+    );
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::clipPlaneToPoint::usage_
+(
+    clipPlaneToPoint::typeName,
+    "\n    Usage: clipPlaneToPoint (px py pz) (nx ny nz)\n\n"
+    "    Select points above the plane\n\n"
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::clipPlaneToPoint::combine(topoSet& set, const bool add) const
+{
+    // Mesh points above plane
+
+    const pointField& ctrs = mesh_.points();
+
+    forAll(ctrs, elemi)
+    {
+        if (((ctrs[elemi] - point_) & normal_) > 0)
+        {
+            addOrDelete(set, elemi, add);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::clipPlaneToPoint::clipPlaneToPoint
+(
+    const polyMesh& mesh,
+    const point& basePoint,
+    const vector& normal
+)
+:
+    topoSetPointSource(mesh),
+    point_(basePoint),
+    normal_(normal)
+{}
+
+
+Foam::clipPlaneToPoint::clipPlaneToPoint
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    clipPlaneToPoint
+    (
+        mesh,
+        dict.get<vector>("point"),
+        dict.get<vector>("normal")
+    )
+{}
+
+
+Foam::clipPlaneToPoint::clipPlaneToPoint
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetPointSource(mesh),
+    point_(checkIs(is)),
+    normal_(checkIs(is))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::clipPlaneToPoint::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (action == topoSetSource::ADD || action == topoSetSource::NEW)
+    {
+        if (verbose_)
+        {
+            Info<< "    Adding points above plane at "
+                << point_ << " with normal " << normal_ << endl;
+        }
+
+        combine(set, true);
+    }
+    else if (action == topoSetSource::SUBTRACT)
+    {
+        if (verbose_)
+        {
+            Info<< "    Removing points above plane at "
+                << point_ << " with normal " << normal_ << endl;
+        }
+
+        combine(set, false);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/topoSet/pointSources/clipPlaneToPoint/clipPlaneToPoint.H b/src/meshTools/topoSet/pointSources/clipPlaneToPoint/clipPlaneToPoint.H
new file mode 100644
index 0000000000..707f533339
--- /dev/null
+++ b/src/meshTools/topoSet/pointSources/clipPlaneToPoint/clipPlaneToPoint.H
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::clipPlaneToPoint
+
+Description
+    A \c topoSetPointSource to select all points
+    above the specified plane
+
+    Operands:
+    \table
+      Operand   | Type     | Location
+      output    | pointSet | $FOAM_CASE/constant/polyMesh/sets/\<set\>
+    \endtable
+
+Usage
+    Minimal example by using \c system/topoSetDict.actions:
+    \verbatim
+    {
+        // Mandatory (inherited) entries
+        name        <name>;
+        type        pointSet;
+        action      <action>;
+
+        // Mandatory entries
+        source      clipPlaneToPoint;
+
+        point  (x y z);
+        normal (x y z);
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property   | Description                         | Type | Req'd | Dflt
+      name       | Name of pointSet                    | word |  yes  | -
+      type       | Type name: pointSet                 | word |  yes  | -
+      action     | Action applied on points - see below| word |  yes  | -
+      source     | Source name: clipPlaneToPoint        | word |  yes  | -
+    \endtable
+
+    Options for the \c action entry:
+    \verbatim
+      new      | Create a new pointSet from selected points
+      add      | Add selected points into this pointSet
+      subtract | Remove selected cells from this pointSet
+    \endverbatim
+
+    Options for the conditional mandatory entries:
+    \verbatim
+      Entry    | Description                        | Type   | Req'd  | Dflt
+      point    | The reference point for the plane  | vector | yes    | -
+      normal   | The unit normal of the plane       | vector | yes    | -
+    \endverbatim
+
+See also
+    - Foam::topoSetSource
+    - Foam::topoSetPointSource
+
+SourceFiles
+    clipPlaneToPoint.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_clipPlaneToPoint_H
+#define Foam_clipPlaneToPoint_H
+
+#include "topoSetPointSource.H"
+#include "point.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class clipPlaneToPoint Declaration
+\*---------------------------------------------------------------------------*/
+
+class clipPlaneToPoint
+:
+    public topoSetPointSource
+{
+    // Private Data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Point on the plane
+        const vector point_;
+
+        //- Normal to the plane
+        const vector normal_;
+
+
+    // Private Member Functions
+
+        void combine(topoSet& set, const bool add) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("clipPlaneToPoint");
+
+
+    // Constructors
+
+        //- No default construct
+        clipPlaneToPoint() = delete;
+
+        //- Construct from components
+        clipPlaneToPoint
+        (
+            const polyMesh& mesh,
+            const point& basePoint,
+            const vector& normal
+        );
+
+        //- Construct from dictionary
+        clipPlaneToPoint(const polyMesh& mesh, const dictionary& dict);
+
+        //- Construct from stream - a single box.
+        clipPlaneToPoint(const polyMesh& mesh, Istream& is);
+
+
+    //- Destructor
+    virtual ~clipPlaneToPoint() = default;
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet& set
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/finiteArea/liquidFilmFoam/cylinder/Allrun b/tutorials/finiteArea/liquidFilmFoam/cylinder/Allrun
index 1da825daf3..6c2d0c8c0e 100755
--- a/tutorials/finiteArea/liquidFilmFoam/cylinder/Allrun
+++ b/tutorials/finiteArea/liquidFilmFoam/cylinder/Allrun
@@ -9,6 +9,8 @@ runApplication blockMesh
 
 runApplication makeFaMesh
 
+runApplication setFields
+
 runApplication $(getApplication)
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/finiteArea/liquidFilmFoam/cylinder/Allrun-parallel b/tutorials/finiteArea/liquidFilmFoam/cylinder/Allrun-parallel
index 0876a09fcf..b5691f584b 100755
--- a/tutorials/finiteArea/liquidFilmFoam/cylinder/Allrun-parallel
+++ b/tutorials/finiteArea/liquidFilmFoam/cylinder/Allrun-parallel
@@ -11,6 +11,8 @@ runApplication decomposePar -no-finite-area
 
 runParallel makeFaMesh
 
+runParallel setFields
+
 runParallel $(getApplication)
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/finiteArea/liquidFilmFoam/cylinder/system/setFieldsDict b/tutorials/finiteArea/liquidFilmFoam/cylinder/system/setFieldsDict
new file mode 100644
index 0000000000..5a85687929
--- /dev/null
+++ b/tutorials/finiteArea/liquidFilmFoam/cylinder/system/setFieldsDict
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2212                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    areaScalarFieldValue h 0.00014
+);
+
+regions
+(
+    clipPlaneToFace
+    {
+        point  (0 0 0);
+        normal (1 0 0);
+
+        fieldValues
+        (
+            areaScalarFieldValue h 0.00015
+        );
+    }
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/fluidisedBed/system/setFieldsDict b/tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/fluidisedBed/system/setFieldsDict
index d61ff2eb8f..724a5d1789 100644
--- a/tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/fluidisedBed/system/setFieldsDict
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/fluidisedBed/system/setFieldsDict
@@ -18,6 +18,7 @@ defaultFieldValues
 (
     volScalarFieldValue alpha.air 1
     volScalarFieldValue alpha.particles 0
+    // areaScalarFieldValue h 1e-3
 );
 
 regions