ENH: redistributePar - add support for pointFields (#2436)

applications/utilities/parallelProcessing/redistributePar/Make/files

@@ -1,6 +1,8 @@
 passivePositionParticleCloud.C
 parLagrangianRedistributor.C
+parPointFieldDistributor.C
 parFvFieldReconstructor.C
+
 loadOrCreateMesh.C
 redistributePar.C
 

applications/utilities/parallelProcessing/redistributePar/parPointFieldDistributor.C

@@ -0,0 +1,238 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "parPointFieldDistributor.H"
+#include "processorPointPatch.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+int Foam::parPointFieldDistributor::verbose_ = 0;
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::parPointFieldDistributor::parPointFieldDistributor
+(
+    const pointMesh& srcMesh,
+    const bool savePoints,
+    const bool isWriteProc
+)
+:
+    srcMesh_(srcMesh),
+    nOldPoints_(srcMesh.size()),
+    patchMeshPoints_(),
+    tgtMeshRef_(nullptr),
+    distMapRef_(nullptr),
+    patchPointMaps_(),
+    isWriteProc_(isWriteProc)
+{
+    if (savePoints)
+    {
+        saveMeshPoints();
+    }
+}
+
+
+Foam::parPointFieldDistributor::parPointFieldDistributor
+(
+    const pointMesh& srcMesh,
+    const pointMesh& tgtMesh,
+    const mapDistributePolyMesh& distMap,
+    const bool savePoints,
+    const bool isWriteProc
+)
+:
+    srcMesh_(srcMesh),
+    nOldPoints_(srcMesh.size()),
+    patchMeshPoints_(),
+    tgtMeshRef_(tgtMesh),
+    distMapRef_(distMap),
+    patchPointMaps_(),
+    isWriteProc_(isWriteProc)
+{
+    if (savePoints)
+    {
+        saveMeshPoints();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::parPointFieldDistributor::hasMeshPoints() const
+{
+    return !patchMeshPoints_.empty();
+}
+
+
+bool Foam::parPointFieldDistributor::hasPatchPointMaps() const
+{
+    return !patchPointMaps_.empty();
+}
+
+
+bool Foam::parPointFieldDistributor::hasTarget() const
+{
+    return (tgtMeshRef_ && distMapRef_);
+}
+
+
+void Foam::parPointFieldDistributor::clearMeshPoints()
+{
+    patchMeshPoints_.clear();
+}
+
+void Foam::parPointFieldDistributor::clearPatchPointMaps()
+{
+    patchPointMaps_.clear();
+}
+
+
+void Foam::parPointFieldDistributor::saveMeshPoints()
+{
+    const pointBoundaryMesh& patches = srcMesh_.boundary();
+
+    patchMeshPoints_.clear();
+    patchMeshPoints_.resize(patches.size());
+
+    forAll(patches, patchi)
+    {
+        if (!isA<processorPointPatch>(patches[patchi]))
+        {
+            // Copy meshPoints
+            patchMeshPoints_.set
+            (
+                patchi,
+                new labelList(patches[patchi].meshPoints())
+            );
+        }
+    }
+}
+
+
+void Foam::parPointFieldDistributor::createPatchPointMaps()
+{
+    if (!tgtMeshRef_ || !distMapRef_)
+    {
+        FatalErrorInFunction
+            << "Cannot create maps without target mesh and/or distribution!"
+            << abort(FatalError);
+    }
+
+    const auto& tgtMesh = tgtMeshRef_();
+    const auto& distMap = distMapRef_();
+
+    const auto& newPatches = tgtMesh.boundary();
+    const auto& oldPatches = srcMesh_.boundary();
+
+    patchPointMaps_.clear();
+    patchPointMaps_.resize(oldPatches.size());
+
+    // if (patchPointMaps_.size() != patchMeshPoints_.size())
+    // {
+    //     // Warn?
+    // }
+
+    forAll(oldPatches, patchi)
+    {
+        if (!isA<processorPointPatch>(oldPatches[patchi]))
+        {
+            // Create map for patch points only
+            labelList oldToNewSub;
+            labelList oldToNewConstruct;
+
+            // Copy point map
+            patchPointMaps_.set
+            (
+                patchi,
+                new mapDistributeBase(distMap.pointMap())
+            );
+
+            const labelList& oldMeshPoints =
+            (
+                patchMeshPoints_.test(patchi)
+              ? patchMeshPoints_[patchi]
+              : oldPatches[patchi].meshPoints()  // <- Questionable!
+            );
+
+            patchPointMaps_[patchi].compactData
+            (
+                oldMeshPoints,
+                newPatches[patchi].meshPoints(),
+                oldToNewSub,
+                oldToNewConstruct,
+                nOldPoints_,
+                UPstream::msgType()
+            );
+        }
+    }
+}
+
+
+void Foam::parPointFieldDistributor::resetTarget()
+{
+    tgtMeshRef_.reset(nullptr);
+    distMapRef_.reset(nullptr);
+
+    // Old maps are now invalid
+    clearPatchPointMaps();
+}
+
+
+void Foam::parPointFieldDistributor::resetTarget
+(
+    const pointMesh& tgtMesh,
+    const mapDistributePolyMesh& distMap
+)
+{
+    tgtMeshRef_.cref(tgtMesh);
+    distMapRef_.cref(distMap);
+
+    // Old maps are now invalid
+    clearPatchPointMaps();
+}
+
+
+Foam::label Foam::parPointFieldDistributor::distributeAllFields
+(
+    const IOobjectList& objects,
+    const wordRes& selected
+) const
+{
+    label nTotal = 0;
+
+    nTotal += distributePointFields<scalar>(objects, selected);
+    nTotal += distributePointFields<vector>(objects, selected);
+    nTotal += distributePointFields<symmTensor>(objects, selected);
+    nTotal += distributePointFields<sphericalTensor>(objects, selected);
+    nTotal += distributePointFields<tensor>(objects, selected);
+
+    return nTotal;
+}
+
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/redistributePar/parPointFieldDistributor.H

@@ -0,0 +1,255 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::parPointFieldDistributor
+
+Description
+    Distributor/redistributor for point fields,
+    uses a two (or three) stage construction.
+
+    The inconvenient multi-stage construction is needed since the
+    pointMesh is directly associated with a polyMesh, which will probably
+    have changed while creating the target mesh. This means that it is
+    necessary to save the size of the source mesh and all of its
+    patch meshPoints prior to making any changes (eg, creating the target
+    mesh).
+
+    -# Create with specified source mesh
+    -# Save the meshPoints (per boundary) for the source mesh
+    -# Attach a target mesh and mesh distribution
+    -# Map the point fields
+    .
+
+    Runs in parallel. Redistributes from srcMesh to tgtMesh.
+
+SourceFiles
+    parPointFieldDistributor.C
+    parPointFieldDistributorTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_parPointFieldDistributor_H
+#define Foam_parPointFieldDistributor_H
+
+#include "PtrList.H"
+#include "pointMesh.H"
+#include "pointFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class mapDistributePolyMesh;
+class mapDistributeBase;
+class IOobjectList;
+
+/*---------------------------------------------------------------------------*\
+                  Class parPointFieldDistributor Declaration
+\*---------------------------------------------------------------------------*/
+
+class parPointFieldDistributor
+{
+    // Private Data
+
+        //- The source mesh reference
+        const pointMesh& srcMesh_;
+
+        //- Number of points in the old (source) mesh
+        const label nOldPoints_;
+
+        //- The pointPatch mesh points
+        PtrList<labelList> patchMeshPoints_;
+
+        //- The target (destination) mesh reference
+        refPtr<pointMesh> tgtMeshRef_;
+
+        //- Distribution map reference
+        refPtr<mapDistributePolyMesh> distMapRef_;
+
+        //- Point patch mappers
+        PtrList<mapDistributeBase> patchPointMaps_;
+
+        //- Do I need to write (eg, master only for reconstruct)
+        bool isWriteProc_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        parPointFieldDistributor(const parPointFieldDistributor&) = delete;
+
+        //- No copy assignment
+        void operator=(const parPointFieldDistributor&) = delete;
+
+
+public:
+
+    //- Output verbosity when writing
+    static int verbose_;
+
+
+    // Constructors
+
+        //- Basic construction
+        //
+        //  \param srcMesh  The source pointMesh
+        //  \param savePoints  Call saveMeshPoints() immediately
+        //  \param isWriteProc  Tagged for output writing (on this proc)
+        explicit parPointFieldDistributor
+        (
+            const pointMesh& srcMesh,
+            const bool savePoints = false,
+            const bool isWriteProc = false
+        );
+
+        //- Full construction of source/target
+        //
+        //  \param srcMesh  The source pointMesh
+        //  \param tgtMesh  The target pointMesh
+        //  \param distMap  The distribution map
+        //  \param savePoints  Call saveMeshPoints() immediately
+        //  \param isWriteProc  Tagged for output writing (on this proc)
+        explicit parPointFieldDistributor
+        (
+            const pointMesh& srcMesh,
+            const pointMesh& tgtMesh,
+            const mapDistributePolyMesh& distMap,
+            const bool savePoints = false,
+            const bool isWriteProc = false
+        );
+
+
+    // Member Functions
+
+        //- True if meshPoints (per boundary) for the source mesh
+        //- have been saved
+        bool hasMeshPoints() const;
+
+        //- True if patch maps (per boundary) exist
+        bool hasPatchPointMaps() const;
+
+        //- True if a target mesh/distribution map has been attached
+        bool hasTarget() const;
+
+        //- Clear out meshPoints (per boundary) for the source mesh
+        void clearMeshPoints();
+
+        //- Clear out patch maps (per boundary)
+        void clearPatchPointMaps();
+
+        //- Create/recreate meshPoints (per boundary) for the source mesh
+        void saveMeshPoints();
+
+        //- Construct per-patch addressing
+        void createPatchPointMaps();
+
+        //- Clear target mesh / distribution map
+        void resetTarget();
+
+        //- Reset target mesh / distribution map
+        void resetTarget
+        (
+            const pointMesh& tgtMesh,
+            const mapDistributePolyMesh& distMap
+        );
+
+        //- Get status of write enabled (on this proc)
+        bool isWriteProc() const noexcept
+        {
+            return isWriteProc_;
+        }
+
+        //- Change status of write enabled (on this proc)
+        bool isWriteProc(const bool on) noexcept
+        {
+            bool old(isWriteProc_);
+            isWriteProc_ = on;
+            return old;
+        }
+
+
+    // Field Mapping
+
+        //- Read, distribute and write all/selected point field types
+        //- (scalar, vector, ... types)
+        label distributeAllFields
+        (
+            const IOobjectList& objects,
+            const wordRes& selectedFields = wordRes()
+        ) const;
+
+
+        //- Distribute point field
+        template<class Type>
+        tmp<GeometricField<Type, pointPatchField, pointMesh>>
+        distributeField
+        (
+            const GeometricField<Type, pointPatchField, pointMesh>& fld
+        ) const;
+
+        //- Read and distribute point field
+        template<class Type>
+        tmp<GeometricField<Type, pointPatchField, pointMesh>>
+        distributePointField
+        (
+            const IOobject& fieldObject
+        ) const;
+
+        //- Read, distribute and write all/selected point fields
+        template<class Type>
+        label distributePointFields
+        (
+            const IOobjectList& objects,
+            const wordRes& selectedFields = wordRes()
+        ) const;
+
+        //- Distributed each (unregistered!) point field
+        //- and store the result on its objectRegistry
+        template<class Type>
+        void distributeAndStore
+        (
+            const PtrList<GeometricField<Type, pointPatchField, pointMesh>>&
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "parPointFieldDistributorTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/redistributePar/parPointFieldDistributorTemplates.C

@@ -0,0 +1,239 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "Time.H"
+#include "emptyPointPatchField.H"
+#include "IOobjectList.H"
+#include "mapDistributePolyMesh.H"
+#include "distributedFieldMapper.H"
+#include "distributedPointPatchFieldMapper.H"
+#include "pointFields.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::pointPatchField, Foam::pointMesh>>
+Foam::parPointFieldDistributor::distributeField
+(
+    const GeometricField<Type, pointPatchField, pointMesh>& fld
+) const
+{
+    if (!tgtMeshRef_ || !distMapRef_)
+    {
+        FatalErrorInFunction
+            << "Cannot map field without target mesh and/or distribution!"
+            << abort(FatalError);
+    }
+    if (!hasPatchPointMaps())
+    {
+        const_cast<parPointFieldDistributor&>(*this).createPatchPointMaps();
+    }
+
+    const auto& tgtMesh = tgtMeshRef_();
+    const auto& distMap = distMapRef_();
+
+    // Create internalField by remote mapping
+    distributedFieldMapper mapper
+    (
+        labelUList::null(),
+        distMap.pointMap()
+    );
+
+    DimensionedField<Type, pointMesh> internalField
+    (
+        IOobject
+        (
+            fld.name(),
+            tgtMesh.time().timeName(),
+            fld.local(),
+            tgtMesh.thisDb(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        tgtMesh,
+        fld.dimensions(),
+        Field<Type>(fld.internalField(), mapper)
+    );
+
+    internalField.oriented() = fld.oriented();
+
+
+    // Create patchFields by remote mapping
+
+    PtrList<pointPatchField<Type>> newPatchFields(tgtMesh.boundary().size());
+
+    const auto& bfld = fld.boundaryField();
+
+    forAll(bfld, patchi)
+    {
+        if (patchPointMaps_.set(patchi))
+        {
+            // Clone local patch field
+
+            const distributedPointPatchFieldMapper mapper
+            (
+                labelUList::null(),
+                patchPointMaps_[patchi]
+            );
+
+            // Map into local copy
+            newPatchFields.set
+            (
+                patchi,
+                pointPatchField<Type>::New
+                (
+                    bfld[patchi],
+                    tgtMesh.boundary()[patchi],    // pointPatch
+                    DimensionedField<Type, pointMesh>::null(),
+                    mapper
+                )
+            );
+
+            // Note: below alternative, 'clone' method will not work since
+            // there is no clone to reset both internalField reference and
+            // patch reference. TBD.
+            //newPatchFields.set
+            //(
+            //    patchi,
+            //    bfld[patchi].clone
+            //    (
+            //        tgtMesh.boundary()[patchi],
+            //        DimensionedField<Type, pointMesh>::null(),
+            //        mapper
+            //    )
+            //);
+        }
+    }
+
+    // Add some empty patchFields on remaining patches (this also handles
+    // e.g. processorPatchFields or any other constraint type patches)
+    forAll(newPatchFields, patchi)
+    {
+        if (!newPatchFields.set(patchi))
+        {
+            newPatchFields.set
+            (
+                patchi,
+                pointPatchField<Type>::New
+                (
+                    emptyPointPatchField<Type>::typeName,
+                    tgtMesh.boundary()[patchi],
+                    DimensionedField<Type, pointMesh>::null()
+                )
+            );
+        }
+    }
+
+    return
+        tmp<GeometricField<Type, pointPatchField, pointMesh>>::New
+        (
+            std::move(internalField),
+            newPatchFields
+        );
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::pointPatchField, Foam::pointMesh>>
+Foam::parPointFieldDistributor::distributePointField
+(
+    const IOobject& fieldObject
+) const
+{
+    // Read field
+    GeometricField<Type, pointPatchField, pointMesh> fld
+    (
+        fieldObject,
+        srcMesh_
+    );
+
+    // Redistribute
+    return distributeField(fld);
+}
+
+
+template<class Type>
+Foam::label Foam::parPointFieldDistributor::distributePointFields
+(
+    const IOobjectList& objects,
+    const wordRes& selectedFields
+) const
+{
+    typedef GeometricField<Type, pointPatchField, pointMesh> fieldType;
+
+    UPtrList<const IOobject> fieldObjects
+    (
+        selectedFields.empty()
+      ? objects.sorted<fieldType>()
+      : objects.sorted<fieldType>(selectedFields)
+    );
+
+    label nFields = 0;
+    for (const IOobject& io : fieldObjects)
+    {
+        if (verbose_)
+        {
+            if (!nFields)
+            {
+                Info<< "    Reconstructing "
+                    << fieldType::typeName << "s\n" << nl;
+            }
+            Info<< "        " << io.name() << nl;
+        }
+        ++nFields;
+
+        tmp<fieldType> tfld(distributePointField<Type>(io));
+        if (isWriteProc_)
+        {
+            tfld().write();
+        }
+    }
+
+    if (nFields && verbose_) Info<< endl;
+    return nFields;
+}
+
+
+template<class Type>
+void Foam::parPointFieldDistributor::distributeAndStore
+(
+    const PtrList<GeometricField<Type, pointPatchField, pointMesh>>& fields
+) const
+{
+    for (const auto& fld : fields)
+    {
+        // Distribute and store
+        auto tfld = distributeField(fld);
+
+        tfld.ref().writeOpt(IOobject::AUTO_WRITE);
+
+        tfld().mesh().thisDb().store(tfld);
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/redistributePar/redistributePar.C

@@ -95,6 +95,7 @@ Usage
 #include "regionProperties.H"
 
 #include "parFvFieldReconstructor.H"
+#include "parPointFieldDistributor.H"
 #include "parLagrangianRedistributor.H"
 #include "unmappedPassivePositionParticleCloud.H"
 #include "hexRef8Data.H"
@@ -535,7 +536,7 @@ void determineDecomposition
 
 // Variant of GeometricField::correctBoundaryConditions that only
 // evaluates selected patch fields
-template<class GeoField, class CoupledPatchType>
+template<class CoupledPatchType, class GeoField>
 void correctCoupledBoundaryConditions(fvMesh& mesh)
 {
     HashTable<GeoField*> flds
@@ -576,6 +577,7 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
     //Info<< "Before distribution:" << endl;
     //printMeshData(mesh);
 
+    // Storage of fields
 
     PtrList<volScalarField> volScalarFields;
     PtrList<volVectorField> volVectorFields;
@@ -595,7 +597,24 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
     PtrList<DimensionedField<symmTensor, volMesh>> dimSymmTensorFields;
     PtrList<DimensionedField<tensor, volMesh>> dimTensorFields;
 
-    DynamicList<word> pointFieldNames;
+    PtrList<pointScalarField> pointScalarFields;
+    PtrList<pointVectorField> pointVectorFields;
+    PtrList<pointTensorField> pointTensorFields;
+    PtrList<pointSphericalTensorField> pointSphTensorFields;
+    PtrList<pointSymmTensorField> pointSymmTensorFields;
+
+    // Self-contained pointMesh for reading pointFields
+    const pointMesh oldPointMesh(mesh);
+
+    // Track how many (if any) pointFields are read/mapped
+    label nPointFields = 0;
+
+    parPointFieldDistributor pointDistributor
+    (
+        oldPointMesh,   // source mesh
+        false,          // savePoints=false (ie, delay until later)
+        false           // Do not write
+    );
 
 
     if (doReadFields)
@@ -656,190 +675,68 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
         }
 
 
-        // volFields
-
         if (Pstream::master() && decompose)
         {
             runTime.caseName() = baseRunTime.caseName();
             runTime.processorCase(false);
         }
 
-        fieldsDistributor::readFields
+        // Field reading
-        (
-            haveMesh,
-            mesh,
-            subsetterPtr,
-            objects,
-            volScalarFields
-        );
 
-        fieldsDistributor::readFields
+        #undef  doFieldReading
-        (
+        #define doFieldReading(Storage)                                       \
-            haveMesh,
+        {                                                                     \
-            mesh,
+            fieldsDistributor::readFields                                     \
-            subsetterPtr,
+            (                                                                 \
-            objects,
+                haveMesh, mesh, subsetterPtr, objects, Storage                \
-            volVectorFields
+            );                                                                \
-        );
+        }
 
-        fieldsDistributor::readFields
+        // volField
-        (
+        doFieldReading(volScalarFields);
-            haveMesh,
+        doFieldReading(volVectorFields);
-            mesh,
+        doFieldReading(volSphereTensorFields);
-            subsetterPtr,
+        doFieldReading(volSymmTensorFields);
-            objects,
+        doFieldReading(volTensorFields);
-            volSphereTensorFields
-        );
-
-        fieldsDistributor::readFields
-        (
-            haveMesh,
-            mesh,
-            subsetterPtr,
-            objects,
-            volSymmTensorFields
-        );
-
-        fieldsDistributor::readFields
-        (
-            haveMesh,
-            mesh,
-            subsetterPtr,
-            objects,
-            volTensorFields
-        );
-
-
-        // surfaceFields
-
-        fieldsDistributor::readFields
-        (
-            haveMesh,
-            mesh,
-            subsetterPtr,
-            objects,
-            surfScalarFields
-        );
-
-        fieldsDistributor::readFields
-        (
-            haveMesh,
-            mesh,
-            subsetterPtr,
-            objects,
-            surfVectorFields
-        );
-
-        fieldsDistributor::readFields
-        (
-            haveMesh,
-            mesh,
-            subsetterPtr,
-            objects,
-            surfSphereTensorFields
-        );
-
-        fieldsDistributor::readFields
-        (
-            haveMesh,
-            mesh,
-            subsetterPtr,
-            objects,
-            surfSymmTensorFields
-        );
-
-        fieldsDistributor::readFields
-        (
-            haveMesh,
-            mesh,
-            subsetterPtr,
-            objects,
-            surfTensorFields
-        );
 
+        // surfaceField
+        doFieldReading(surfScalarFields);
+        doFieldReading(surfVectorFields);
+        doFieldReading(surfSphereTensorFields);
+        doFieldReading(surfSymmTensorFields);
+        doFieldReading(surfTensorFields);
 
         // Dimensioned internal fields
-        fieldsDistributor::readFields
+        doFieldReading(dimScalarFields);
-        (
+        doFieldReading(dimVectorFields);
-            haveMesh,
+        doFieldReading(dimSphereTensorFields);
-            mesh,
+        doFieldReading(dimSymmTensorFields);
-            subsetterPtr,
+        doFieldReading(dimTensorFields);
-            objects,
-            dimScalarFields
-        );
 
-        fieldsDistributor::readFields
+        // pointFields
-        (
+        nPointFields = 0;
-            haveMesh,
-            mesh,
-            subsetterPtr,
-            objects,
-            dimVectorFields
-        );
 
-        fieldsDistributor::readFields
+        #undef  doFieldReading
-        (
+        #define doFieldReading(Storage)                                       \
-            haveMesh,
+        {                                                                     \
-            mesh,
+            fieldsDistributor::readFields                                     \
-            subsetterPtr,
+            (                                                                 \
-            objects,
+                haveMesh, oldPointMesh, subsetterPtr, objects, Storage,       \
-            dimSphereTensorFields
+                true  /* (deregister field) */                                \
-        );
+            );                                                                \
-
+            nPointFields += Storage.size();                                   \
-        fieldsDistributor::readFields
-        (
-            haveMesh,
-            mesh,
-            subsetterPtr,
-            objects,
-            dimSymmTensorFields
-        );
-
-        fieldsDistributor::readFields
-        (
-            haveMesh,
-            mesh,
-            subsetterPtr,
-            objects,
-            dimTensorFields
-        );
-
-
-        // pointFields currently not supported. Read their names so we
-        // can delete them.
-        {
-            // Get my objects of type
-            pointFieldNames.append
-            (
-                objects.lookupClass(pointScalarField::typeName).sortedNames()
-            );
-            pointFieldNames.append
-            (
-                objects.lookupClass(pointVectorField::typeName).sortedNames()
-            );
-            pointFieldNames.append
-            (
-                objects.lookupClass
-                (
-                    pointSphericalTensorField::typeName
-                ).sortedNames()
-            );
-            pointFieldNames.append
-            (
-                objects.lookupClass
-                (
-                    pointSymmTensorField::typeName
-                ).sortedNames()
-            );
-            pointFieldNames.append
-            (
-                objects.lookupClass(pointTensorField::typeName).sortedNames()
-            );
-
-            // Make sure all processors have the same set
-            Pstream::scatter(pointFieldNames);
         }
 
+        doFieldReading(pointScalarFields);
+        doFieldReading(pointVectorFields);
+        doFieldReading(pointSphTensorFields);
+        doFieldReading(pointSymmTensorFields);
+        doFieldReading(pointTensorFields);
+        #undef doFieldReading
+
+
+        // Done reading
+
         if (Pstream::master() && decompose)
         {
             runTime.caseName() = proc0CaseName;
@@ -847,6 +744,12 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
         }
     }
 
+    // Save pointMesh information before any topology changes occur!
+    if (nPointFields)
+    {
+        pointDistributor.saveMeshPoints();
+    }
+
 
     // Mesh distribution engine
     fvMeshDistribute distributor(mesh);
@@ -859,34 +762,40 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
     printMeshData(mesh);
 
     // Get other side of processor boundaries
-    correctCoupledBoundaryConditions
+    do
-    <
+    {
-        volScalarField,
+        #undef  doCorrectCoupled
-        processorFvPatch
+        #define doCorrectCoupled(FieldType)  \
-    >(mesh);
+        correctCoupledBoundaryConditions<processorFvPatch, FieldType>(mesh);
-    correctCoupledBoundaryConditions
+
-    <
+        doCorrectCoupled(volScalarField);
-        volVectorField,
+        doCorrectCoupled(volVectorField);
-        processorFvPatch
+        doCorrectCoupled(volSphericalTensorField);
-    >(mesh);
+        doCorrectCoupled(volSymmTensorField);
-    correctCoupledBoundaryConditions
+        doCorrectCoupled(volTensorField);
-    <
+        #undef doCorrectCoupled
-        volSphericalTensorField,
+    }
-        processorFvPatch
+    while (false);
-    >(mesh);
+
-    correctCoupledBoundaryConditions
-    <
-        volSymmTensorField,
-        processorFvPatch
-    >(mesh);
-    correctCoupledBoundaryConditions
-    <
-        volTensorField,
-        processorFvPatch
-    >(mesh);
     // No update surface fields
 
 
+    // Map pointFields
+    if (nPointFields)
+    {
+        // Construct new pointMesh from distributed mesh
+        const pointMesh& newPointMesh = pointMesh::New(mesh);
+
+        pointDistributor.resetTarget(newPointMesh, rawMap());
+
+        pointDistributor.distributeAndStore(pointScalarFields);
+        pointDistributor.distributeAndStore(pointVectorFields);
+        pointDistributor.distributeAndStore(pointSphTensorFields);
+        pointDistributor.distributeAndStore(pointSymmTensorFields);
+        pointDistributor.distributeAndStore(pointTensorFields);
+    }
+
+
     // Set the minimum write precision
     IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
 
@@ -925,24 +834,9 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
                 << " to write reconstructed mesh and fields." << endl;
             runTime.caseName() = baseRunTime.caseName();
             const bool oldProcCase(runTime.processorCase(false));
-            const bool oldParRun = Pstream::parRun(false);
 
             mesh.write();
             topoSet::removeFiles(mesh);
-            for (const word& fieldName : pointFieldNames)
-            {
-                IOobject io
-                (
-                    fieldName,
-                    runTime.timeName(),
-                    mesh
-                );
-
-                const fileName fieldFile(io.objectPath());
-                if (topoSet::debug) DebugVar(fieldFile);
-                rm(fieldFile);
-            }
-            Pstream::parRun(oldParRun);
 
             // Now we've written all. Reset caseName on master
             Info<< "Restoring caseName to " << proc0CaseName << endl;
@@ -965,19 +859,6 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
             writeHandler = fileHandler(std::move(defaultHandler));
         }
         topoSet::removeFiles(mesh);
-        for (const word& fieldName : pointFieldNames)
-        {
-            IOobject io
-            (
-                fieldName,
-                runTime.timeName(),
-                mesh
-            );
-
-            const fileName fieldFile(io.objectPath());
-            if (topoSet::debug) DebugVar(fieldFile);
-            rm(fieldFile);
-        }
     }
     Info<< "Written redistributed mesh to " << mesh.facesInstance() << nl
         << endl;
@@ -1724,6 +1605,7 @@ int main(int argc, char *argv[])
     #include "addOverwriteOption.H"
     argList::addBoolOption("decompose", "Decompose case");
     argList::addBoolOption("reconstruct", "Reconstruct case");
+    argList::addVerboseOption("Additional verbosity");
     argList::addDryRunOption
     (
         "Test without writing the decomposition. "
@@ -1794,6 +1676,12 @@ int main(int argc, char *argv[])
     const bool dryrun = args.dryRun();
     const bool newTimes = args.found("newTimes");
 
+    if (args.verbose())
+    {
+        // Report on output
+        parPointFieldDistributor::verbose_ = 1;
+    }
+
     bool decompose = args.found("decompose");
     bool overwrite = args.found("overwrite");
 
@@ -2240,17 +2128,28 @@ int main(int argc, char *argv[])
             distMap = fvMeshTools::readProcAddressing(mesh, baseMeshPtr);
 
             // Construct field mapper
-            autoPtr<parFvFieldReconstructor> fvReconstructorPtr
+            auto fvReconstructorPtr =
-            (
+                autoPtr<parFvFieldReconstructor>::New
-                new parFvFieldReconstructor
                 (
                     baseMeshPtr(),
                     mesh,
                     distMap(),
                     Pstream::master()       // do I need to write?
-                )
+                );
-            );
 
+            // Construct point field mapper
+            const auto& basePointMesh = pointMesh::New(baseMeshPtr());
+            const auto& procPointMesh = pointMesh::New(mesh);
+
+            auto pointFieldDistributorPtr =
+                autoPtr<parPointFieldDistributor>::New
+                (
+                    procPointMesh,   // source
+                    basePointMesh,   // target
+                    distMap(),
+                    false,           // delay
+                    UPstream::master()  // Write reconstructed on master
+                );
 
 
             // Since we start from Times[0] and not runTime.timeName() we
@@ -2324,7 +2223,8 @@ int main(int argc, char *argv[])
                     distMap =
                         fvMeshTools::readProcAddressing(mesh, baseMeshPtr);
 
-                    // Reset field mapper
+                    // Reset field mappers
+
                     fvReconstructorPtr.reset
                     (
                         new parFvFieldReconstructor
@@ -2335,6 +2235,23 @@ int main(int argc, char *argv[])
                             Pstream::master()
                         )
                     );
+
+                    // Construct point field mapper
+                    const auto& basePointMesh = pointMesh::New(baseMeshPtr());
+                    const auto& procPointMesh = pointMesh::New(mesh);
+
+                    pointFieldDistributorPtr.reset
+                    (
+                        new parPointFieldDistributor
+                        (
+                            procPointMesh,  // source
+                            basePointMesh,  // target
+                            distMap(),
+                            false,          // delay until later
+                            UPstream::master()  // Write reconstruct on master
+                        )
+                    );
+
                     lagrangianReconstructorPtr.clear();
                 }
 
@@ -2351,6 +2268,12 @@ int main(int argc, char *argv[])
                     selectedFields
                 );
 
+                // pointfields
+                // - distribute and write (verbose)
+                pointFieldDistributorPtr()
+                    .distributeAllFields(objects, selectedFields);
+
+
                 // Clouds (note: might not be present on all processors)
                 reconstructLagrangian
                 (