From 9335b641f1d38bfa56a22c844ad64810b0bb1a8a Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Mon, 13 May 2024 12:31:57 +0200
Subject: [PATCH] ENH: GeometricBoundaryField evaluate with specified
 communication type

- evaluate()
  the communication type is exposed as a parameter to allow for more
  tuning, but default parameter remains defaultCommsType so there is
  no change in behaviour for existing code

- evaluate_if()
  supports a general selection predicate

- evaluateSelected()
  now does initEvaluate() for all patches, waits and then calls
  evaluate(). This avoids potential deadlocks when multiple patches
  are inter-communicating.

ENH: align DimensionedField reading with GeometricField treatment

- use localIOdictionary to obtain the dictionary contents

ENH: update GeometricField code

- change GeometricField writeData() as primary output method
  (not operator<<) for better clarity of purpose

- use unique_ptr for GeometricField demand-driven data
---
 .../DimensionedField/DimensionedField.H       |  10 +-
 .../DimensionedField/DimensionedFieldIO.C     |  58 +++--
 .../GeometricField/GeometricBoundaryField.C   | 218 +++++++++++-------
 .../GeometricField/GeometricBoundaryField.H   |  58 +++--
 .../GeometricField/GeometricField.C           | 159 ++++---------
 .../GeometricField/GeometricField.H           |  21 +-
 6 files changed, 282 insertions(+), 242 deletions(-)

diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H
index e7b4d21c6f..d37b651781 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H
@@ -115,8 +115,12 @@ private:
         //- Assert that non-zero field size == mesh size
         void checkFieldSize() const;
 
+        //- Read from file if it is present. Checks the readOption flag
         void readIfPresent(const word& fieldDictEntry = "value");
 
+        //- Read the field - create the field dictionary on-the-fly
+        void readField(const word& fieldDictEntry = "value");
+
         //- Implementation for 'New' with specified registerObject preference.
         //  For LEGACY_REGISTER, registration is determined by
         //  objectRegistry::is_cacheTemporaryObject().
@@ -580,8 +584,8 @@ public:
         bool writeData(Ostream& os, const word& fieldDictEntry) const;
 
         //- The writeData function (required by regIOobject),
-        //- call writeData with dictionary entry name = "value"
-        bool writeData(Ostream& os) const;
+        //- calls writeData with dictionary entry name = "value"
+        bool writeData(Ostream& os) const { return writeData(os, "value"); }
 
 
     // Member Operators
@@ -613,12 +617,14 @@ public:
 
     // Ostream Operators
 
+        //- Calls DimensionedField::writeData()
         friend Ostream& operator<< <Type, GeoMesh>
         (
             Ostream& os,
             const DimensionedField<Type, GeoMesh>& df
         );
 
+        //- Calls DimensionedField::writeData()
         friend Ostream& operator<< <Type, GeoMesh>
         (
             Ostream& os,
diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.C b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.C
index 57a0da9e66..7ad8544616 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.C
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.C
@@ -28,6 +28,7 @@ License
 
 #include "DimensionedField.H"
 #include "IOstreams.H"
+#include "localIOdictionary.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -58,15 +59,49 @@ void Foam::DimensionedField<Type, GeoMesh>::readField
 }
 
 
+template<class Type, class GeoMesh>
+void Foam::DimensionedField<Type, GeoMesh>::readField
+(
+    const word& fieldDictEntry
+)
+{
+    dictionary dict
+    (
+        localIOdictionary::readContents
+        (
+            IOobject
+            (
+                this->name(),
+                this->instance(),
+                this->local(),
+                this->db(),
+                IOobjectOption::MUST_READ,
+                IOobjectOption::NO_WRITE,
+                IOobjectOption::NO_REGISTER
+            ),
+            typeName
+        )
+    );
+
+    this->close();
+
+    readField(dict, fieldDictEntry);
+}
+
+
 template<class Type, class GeoMesh>
 void Foam::DimensionedField<Type, GeoMesh>::readIfPresent
 (
     const word& fieldDictEntry
 )
 {
-    if (this->isReadRequired() || (this->isReadOptional() && this->headerOk()))
+    if
+    (
+        this->isReadRequired()
+     || (this->isReadOptional() && this->headerOk())
+    )
     {
-        readField(dictionary(readStream(typeName)), fieldDictEntry);
+        readField(fieldDictEntry);
     }
 }
 
@@ -87,7 +122,7 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
     dimensions_(dimless),
     oriented_()
 {
-    readField(dictionary(readStream(typeName)), fieldDictEntry);
+    readField(fieldDictEntry);
 }
 
 
@@ -134,23 +169,16 @@ bool Foam::DimensionedField<Type, GeoMesh>::writeData
 }
 
 
-template<class Type, class GeoMesh>
-bool Foam::DimensionedField<Type, GeoMesh>::writeData(Ostream& os) const
-{
-    return writeData(os, "value");
-}
-
-
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 template<class Type, class GeoMesh>
 Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
-    const DimensionedField<Type, GeoMesh>& df
+    const DimensionedField<Type, GeoMesh>& fld
 )
 {
-    df.writeData(os);
+    fld.writeData(os);
 
     return os;
 }
@@ -160,11 +188,11 @@ template<class Type, class GeoMesh>
 Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
-    const tmp<DimensionedField<Type, GeoMesh>>& tdf
+    const tmp<DimensionedField<Type, GeoMesh>>& tfld
 )
 {
-    tdf().writeData(os);
-    tdf.clear();
+    tfld().writeData(os);
+    tfld.clear();
 
     return os;
 }
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C
index 70d5ce2556..aff9204eef 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C
@@ -91,7 +91,7 @@ bool Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::checkConsistency
             }
         }
 
-        // Wait for outstanding requests
+        // Wait for outstanding requests (non-blocking)
         UPstream::waitRequests(startOfRequests);
 
         for (auto& pfld : bfld)
@@ -193,10 +193,7 @@ void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::readField
     const dictionary& dict
 )
 {
-    ///if (GeometricField<Type, PatchField, GeoMesh::debug)
-    ///{
-    ///    InfoInFunction << nl;
-    ///}
+    // DebugInFunction << nl;
 
     // Clear the boundary field if already initialised
     this->clear();
@@ -363,10 +360,7 @@ Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::GeometricBoundaryField
     FieldField<PatchField, Type>(bmesh.size()),
     bmesh_(bmesh)
 {
-    ///if (GeometricField<Type, PatchField, GeoMesh::debug)
-    ///{
-    ///    InfoInFunction << nl;
-    ///}
+    // DebugInFunction << nl;
 
     forAll(bmesh_, patchi)
     {
@@ -396,10 +390,7 @@ Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::GeometricBoundaryField
     FieldField<PatchField, Type>(bmesh.size()),
     bmesh_(bmesh)
 {
-    ///if (GeometricField<Type, PatchField, GeoMesh::debug)
-    ///{
-    ///    InfoInFunction << nl;
-    ///}
+    // DebugInFunction << nl;
 
     if
     (
@@ -462,10 +453,7 @@ Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::GeometricBoundaryField
     FieldField<PatchField, Type>(bmesh.size()),
     bmesh_(bmesh)
 {
-    ///if (GeometricField<Type, PatchField, GeoMesh::debug)
-    ///{
-    ///    InfoInFunction << nl;
-    ///}
+    // DebugInFunction << nl;
 
     forAll(bmesh_, patchi)
     {
@@ -484,10 +472,7 @@ Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::GeometricBoundaryField
     FieldField<PatchField, Type>(btf.size()),
     bmesh_(btf.bmesh_)
 {
-    ///if (GeometricField<Type, PatchField, GeoMesh::debug)
-    ///{
-    ///    InfoInFunction << nl;
-    ///}
+    // DebugInFunction << nl;
 
     forAll(bmesh_, patchi)
     {
@@ -508,10 +493,7 @@ Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::GeometricBoundaryField
     FieldField<PatchField, Type>(btf.size()),
     bmesh_(btf.bmesh_)
 {
-    ///if (GeometricField<Type, PatchField, GeoMesh::debug)
-    ///{
-    ///    InfoInFunction << nl;
-    ///}
+    // DebugInFunction << nl;
 
     for (const label patchi : patchIDs)
     {
@@ -545,12 +527,7 @@ Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::GeometricBoundaryField
 :
     FieldField<PatchField, Type>(btf),
     bmesh_(btf.bmesh_)
-{
-    ///if (GeometricField<Type, PatchField, GeoMesh::debug)
-    ///{
-    ///    InfoInFunction << nl;
-    ///}
-}
+{}
 
 
 template<class Type, template<class> class PatchField, class GeoMesh>
@@ -573,10 +550,7 @@ Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::GeometricBoundaryField
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::updateCoeffs()
 {
-    ///if (GeometricField<Type, PatchField, GeoMesh::debug)
-    ///{
-    ///    InfoInFunction << nl;
-    ///}
+    // DebugInFunction << nl;
 
     for (auto& pfld : *this)
     {
@@ -586,32 +560,26 @@ void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::updateCoeffs()
 
 
 template<class Type, template<class> class PatchField, class GeoMesh>
-void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::evaluate()
+void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::evaluate
+(
+    const UPstream::commsTypes commsType
+)
 {
-    ///if (GeometricField<Type, PatchField, GeoMesh::debug)
-    ///{
-    ///    InfoInFunction << nl;
-    ///}
-
-    const UPstream::commsTypes commsType = UPstream::defaultCommsType;
-    const label startOfRequests = UPstream::nRequests();
-
     if
     (
         commsType == UPstream::commsTypes::buffered
      || commsType == UPstream::commsTypes::nonBlocking
     )
     {
+        const label startOfRequests = UPstream::nRequests();
+
         for (auto& pfld : *this)
         {
             pfld.initEvaluate(commsType);
         }
 
-        // Wait for outstanding requests
-        if (commsType == Pstream::commsTypes::nonBlocking)
-        {
-            UPstream::waitRequests(startOfRequests);
-        }
+        // Wait for outstanding requests (non-blocking)
+        UPstream::waitRequests(startOfRequests);
 
         for (auto& pfld : *this)
         {
@@ -641,8 +609,74 @@ void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::evaluate()
     else
     {
         FatalErrorInFunction
-            << "Unsupported communications type "
-            << UPstream::commsTypeNames[commsType]
+            << "Unsupported communications type " << int(commsType) << nl
+            << exit(FatalError);
+    }
+}
+
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+template<class UnaryPredicate>
+void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::evaluate_if
+(
+    const UnaryPredicate& pred,
+    const UPstream::commsTypes commsType
+)
+{
+    if
+    (
+        commsType == UPstream::commsTypes::buffered
+     || commsType == UPstream::commsTypes::nonBlocking
+    )
+    {
+        const label startOfRequests = UPstream::nRequests();
+
+        for (auto& pfld : *this)
+        {
+            if (pred(pfld))
+            {
+                pfld.initEvaluate(commsType);
+            }
+        }
+
+        // Wait for outstanding requests (non-blocking)
+        UPstream::waitRequests(startOfRequests);
+
+        for (auto& pfld : *this)
+        {
+            if (pred(pfld))
+            {
+                pfld.evaluate(commsType);
+            }
+        }
+    }
+    else if (commsType == UPstream::commsTypes::scheduled)
+    {
+        const lduSchedule& patchSchedule =
+            bmesh_.mesh().globalData().patchSchedule();
+
+        for (const auto& schedEval : patchSchedule)
+        {
+            const label patchi = schedEval.patch;
+            auto& pfld = (*this)[patchi];
+
+            if (pred(pfld))
+            {
+                if (schedEval.init)
+                {
+                    pfld.initEvaluate(commsType);
+                }
+                else
+                {
+                    pfld.evaluate(commsType);
+                }
+            }
+        }
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Unsupported communications type " << int(commsType) << nl
             << exit(FatalError);
     }
 }
@@ -651,58 +685,60 @@ void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::evaluate()
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::evaluateSelected
 (
-    const UList<label>& patchIDs
+    const labelUList& patchIDs
 )
 {
+    const label startOfRequests = UPstream::nRequests();
+
     for (const label patchi : patchIDs)
     {
-        auto& pf = (*this)[patchi];
+        auto& pfld = (*this)[patchi];
 
-        DebugInfo<< "Updating " << pf.patch().name() << endl;
+        DebugInfo<< "Updating " << pfld.patch().name() << endl;
 
-        const label startOfRequests = UPstream::nRequests();
+        pfld.initEvaluate(UPstream::commsTypes::nonBlocking);
+    }
 
-        pf.initEvaluate(UPstream::commsTypes::nonBlocking);
+    // Wait for outstanding requests (non-blocking)
+    UPstream::waitRequests(startOfRequests);
 
-        UPstream::waitRequests(startOfRequests);
+    for (const label patchi : patchIDs)
+    {
+        auto& pfld = (*this)[patchi];
 
-        pf.evaluate(UPstream::commsTypes::nonBlocking);
+        pfld.evaluate(UPstream::commsTypes::nonBlocking);
     }
 }
 
 
 template<class Type, template<class> class PatchField, class GeoMesh>
-void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::evaluateLocal()
+void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::evaluateLocal
+(
+    const UPstream::commsTypes commsType
+)
 {
-    ///if (GeometricField<Type, PatchField, GeoMesh::debug)
-    ///{
-    ///    InfoInFunction << nl;
-    ///}
+    // DebugInFunction << nl;
 
     if (!localConsistency)
     {
         return;
     }
 
-    const UPstream::commsTypes commsType = UPstream::defaultCommsType;
-    const label startOfRequests = UPstream::nRequests();
-
     if
     (
         commsType == UPstream::commsTypes::buffered
      || commsType == UPstream::commsTypes::nonBlocking
     )
     {
+        const label startOfRequests = UPstream::nRequests();
+
         for (auto& pfld : *this)
         {
             pfld.initEvaluateLocal(commsType);
         }
 
-        // Wait for outstanding requests
-        if (commsType == Pstream::commsTypes::nonBlocking)
-        {
-            UPstream::waitRequests(startOfRequests);
-        }
+        // Wait for outstanding requests (non-blocking)
+        UPstream::waitRequests(startOfRequests);
 
         for (auto& pfld : *this)
         {
@@ -732,8 +768,7 @@ void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::evaluateLocal()
     else
     {
         FatalErrorInFunction
-            << "Unsupported communications type "
-            << UPstream::commsTypeNames[commsType]
+            << "Unsupported communications type " << int(commsType) << nl
             << exit(FatalError);
     }
 }
@@ -741,17 +776,24 @@ void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::evaluateLocal()
 
 template<class Type, template<class> class PatchField, class GeoMesh>
 template<class CoupledPatchType>
-void
-Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>
-::evaluateCoupled()
+void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::evaluateCoupled
+(
+    const UPstream::commsTypes commsType
+)
 {
-    ///if (GeometricField<Type, PatchField, GeoMesh::debug)
-    ///{
-    ///    InfoInFunction << nl;
-    ///}
+    // Alternative (C++14)
+    //
+    // this->evaluate_if
+    // (
+    //     [](const auto& pfld) -> bool
+    //     {
+    //         const auto* cpp = isA<CoupledPatchType>(pfld.patch());
+    //         return (cpp && cpp->coupled());
+    //     },
+    //     commsType
+    // );
 
-    const UPstream::commsTypes commsType = UPstream::defaultCommsType;
-    const label startOfRequests = UPstream::nRequests();
+    // DebugInFunction << nl;
 
     if
     (
@@ -759,6 +801,8 @@ Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>
      || commsType == UPstream::commsTypes::nonBlocking
     )
     {
+        const label startOfRequests = UPstream::nRequests();
+
         for (auto& pfld : *this)
         {
             const auto* cpp = isA<CoupledPatchType>(pfld.patch());
@@ -769,11 +813,8 @@ Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>
             }
         }
 
-        // Wait for outstanding requests
-        if (commsType == UPstream::commsTypes::nonBlocking)
-        {
-            UPstream::waitRequests(startOfRequests);
-        }
+        // Wait for outstanding requests (non-blocking)
+        UPstream::waitRequests(startOfRequests);
 
         for (auto& pfld : *this)
         {
@@ -813,8 +854,7 @@ Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>
     else
     {
         FatalErrorInFunction
-            << "Unsupported communications type "
-            << UPstream::commsTypeNames[commsType]
+            << "Unsupported communications type " << int(commsType) << nl
             << exit(FatalError);
     }
 }
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.H
index 3e32718409..28834c39ce 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017,2022 OpenFOAM Foundation
-    Copyright (C) 2015-2023 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -49,12 +49,6 @@ SourceFiles
 namespace Foam
 {
 
-// Forward Declarations
-class dictionary;
-
-template<class Type, template<class> class PatchField, class GeoMesh>
-class GeometricField;
-
 /*---------------------------------------------------------------------------*\
                    Class GeometricBoundaryField Declaration
 \*---------------------------------------------------------------------------*/
@@ -74,6 +68,9 @@ public:
         //- The internal field type associated with the boundary fields
         typedef DimensionedField<Type, GeoMesh> Internal;
 
+        //- This boundary field type
+        typedef GeometricBoundaryField<Type, PatchField, GeoMesh> Boundary;
+
         //- The patch field type for the boundary fields
         typedef PatchField<Type> Patch;
 
@@ -182,18 +179,40 @@ public:
         //- Update the boundary condition coefficients
         void updateCoeffs();
 
-        //- Evaluate boundary conditions
-        void evaluate();
+        //- Evaluate boundary conditions for each patch field.
+        //- Uses specified or default comms.
+        void evaluate
+        (
+            const UPstream::commsTypes commsType = UPstream::defaultCommsType
+        );
 
-        //- Evaluate boundary conditions after change in local values
-        void evaluateLocal();
+        //- Evaluate boundary conditions for patch fields matching the
+        //- given predicate. Uses specified or default comms.
+        template<class UnaryPredicate>
+        void evaluate_if
+        (
+            const UnaryPredicate& pred,
+            const UPstream::commsTypes commsType = UPstream::defaultCommsType
+        );
 
-        //- Evaluate boundary conditions on a subset of coupled patches
+        //- Evaluate boundary conditions after change in local values.
+        //- Uses specified or default comms.
+        void evaluateLocal
+        (
+            const UPstream::commsTypes commsType = UPstream::defaultCommsType
+        );
+
+        //- Evaluate boundary conditions on coupled patches of given type.
+        //- Uses specified or default comms.
         template<class CoupledPatchType>
-        void evaluateCoupled();
+        void evaluateCoupled
+        (
+            const UPstream::commsTypes commsType = UPstream::defaultCommsType
+        );
 
-        //- Evaluate boundary conditions for selected patches
-        void evaluateSelected(const UList<label>& patchIDs);
+        //- Evaluate boundary conditions for selected patches.
+        //- Uses non-blocking comms.
+        void evaluateSelected(const labelUList& patchIDs);
 
         //- Return a list of the patch types
         wordList types() const;
@@ -209,15 +228,18 @@ public:
         //- pointing to interfaces being set
         lduInterfaceFieldPtrsList scalarInterfaces() const;
 
+        //- Helper: check if field has been evaluated. See instantiations.
+        bool check() const;
+
+
+    // Write
+
         //- Write boundary field as dictionary entry
         void writeEntry(const word& keyword, Ostream& os) const;
 
         //- Write dictionary entries of the individual boundary fields.
         void writeEntries(Ostream& os) const;
 
-        //- Helper: check if field has been evaluated. See instantiations.
-        bool check() const;
-
 
     // Member Operators
 
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
index 12603e72ef..2ec59ff2d7 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2023 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,6 @@ License
 
 #include "GeometricField.H"
 #include "Time.H"
-#include "demandDrivenData.H"
 #include "dictionary.H"
 #include "localIOdictionary.H"
 #include "meshState.H"
@@ -77,19 +76,22 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::readFields
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::readFields()
 {
-    const localIOdictionary dict
+    dictionary dict
     (
-        IOobject
+        localIOdictionary::readContents
         (
-            this->name(),
-            this->instance(),
-            this->local(),
-            this->db(),
-            IOobjectOption::MUST_READ,
-            IOobjectOption::NO_WRITE,
-            IOobjectOption::NO_REGISTER
-        ),
-        typeName
+            IOobject
+            (
+                this->name(),
+                this->instance(),
+                this->local(),
+                this->db(),
+                IOobjectOption::MUST_READ,
+                IOobjectOption::NO_WRITE,
+                IOobjectOption::NO_REGISTER
+            ),
+            typeName
+        )
     );
 
     this->close();
@@ -111,10 +113,7 @@ bool Foam::GeometricField<Type, PatchField, GeoMesh>::readIfPresent()
     else if
     (
         this->isReadOptional()
-     && this->template typeHeaderOk<GeometricField<Type, PatchField, GeoMesh>>
-        (
-            true
-        )
+     && this->template typeHeaderOk<this_type>(true)  // checkType = true
     )
     {
         readFields();
@@ -143,20 +142,13 @@ bool Foam::GeometricField<Type, PatchField, GeoMesh>::readOldTimeIfPresent()
 
     if
     (
-        field0.template typeHeaderOk<GeometricField<Type, PatchField, GeoMesh>>
-        (
-            true
-        )
+        field0.template typeHeaderOk<this_type>(true)  // checkType = true
     )
     {
         DebugInFunction
             << "Reading old time level for field" << nl << this->info() << endl;
 
-        field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
-        (
-            field0,
-            this->mesh()
-        );
+        field0Ptr_ = std::make_unique<this_type>(field0, this->mesh());
 
         // Ensure the old time field oriented flag is set to the parent's state
         // Note: required for backwards compatibility in case of restarting from
@@ -190,8 +182,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, mesh, dims, false),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(mesh.boundary(), *this, patchFieldType)
 {
     DebugInFunction
@@ -213,8 +203,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, mesh, dims, false),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(mesh.boundary(), *this, patchFieldTypes, actualPatchTypes)
 {
     DebugInFunction
@@ -236,8 +224,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, mesh, value, dims, false),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(mesh.boundary(), *this, patchFieldType)
 {
     DebugInFunction
@@ -262,8 +248,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, mesh, value, dims, false),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(mesh.boundary(), *this, patchFieldTypes, actualPatchTypes)
 {
     DebugInFunction
@@ -327,8 +311,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, diField),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(this->mesh().boundary(), *this, ptfl)
 {
     DebugInFunction
@@ -348,8 +330,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, std::move(diField)),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(this->mesh().boundary(), *this, ptfl)
 {
     DebugInFunction
@@ -369,8 +349,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, tfield),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(this->mesh().boundary(), *this, ptfl)
 {
     DebugInFunction
@@ -389,8 +367,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(diField),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(this->mesh().boundary(), *this, ptfl)
 {
     DebugInFunction
@@ -409,8 +385,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(std::move(diField)),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(this->mesh().boundary(), *this, ptfl)
 {
     DebugInFunction
@@ -432,8 +406,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, mesh, dims, iField),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(mesh.boundary(), *this, patchFieldType)
 {
     DebugInFunction
@@ -455,8 +427,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, mesh, dims, std::move(iField)),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(mesh.boundary(), *this, patchFieldType)
 {
     DebugInFunction
@@ -478,8 +448,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, mesh, dims, iField),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(mesh.boundary(), *this, ptfl)
 {
     DebugInFunction
@@ -501,8 +469,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, mesh, dims, std::move(iField)),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(mesh.boundary(), *this, ptfl)
 {
     DebugInFunction
@@ -524,8 +490,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, mesh, dims, tfield),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(mesh.boundary(), *this, ptfl)
 {
     DebugInFunction
@@ -545,8 +509,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, mesh, dimless, false),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(mesh.boundary())
 {
     DebugInFunction
@@ -584,8 +546,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, mesh, dimless, false),
     timeIndex_(this->time().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(mesh.boundary())
 {
     readFields(dict);
@@ -603,8 +563,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(gf),
     timeIndex_(gf.timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(*this, gf.boundaryField_)
 {
     DebugInFunction
@@ -612,10 +570,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 
     if (gf.field0Ptr_)
     {
-        field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
-        (
-            *gf.field0Ptr_
-        );
+        field0Ptr_ = std::make_unique<this_type>(*gf.field0Ptr_);
     }
 
     this->writeOpt(IOobjectOption::NO_WRITE);
@@ -630,8 +585,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(tgf.constCast(), tgf.movable()),
     timeIndex_(tgf().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(*this, tgf().boundaryField_)
 {
     DebugInFunction
@@ -652,8 +605,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, gf),
     timeIndex_(gf.timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(*this, gf.boundaryField_)
 {
     DebugInFunction
@@ -662,7 +613,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 
     if (!readIfPresent() && gf.field0Ptr_)
     {
-        field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
+        field0Ptr_ = std::make_unique<this_type>
         (
             io.name() + "_0",
             *gf.field0Ptr_
@@ -680,8 +631,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, tgf.constCast(), tgf.movable()),
     timeIndex_(tgf().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(*this, tgf().boundaryField_)
 {
     DebugInFunction
@@ -703,8 +652,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(newName, gf),
     timeIndex_(gf.timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(*this, gf.boundaryField_)
 {
     DebugInFunction
@@ -713,7 +660,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 
     if (!readIfPresent() && gf.field0Ptr_)
     {
-        field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
+        field0Ptr_ = std::make_unique<this_type>
         (
             newName + "_0",
             *gf.field0Ptr_
@@ -731,8 +678,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(newName, tgf.constCast(), tgf.movable()),
     timeIndex_(tgf().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(*this, tgf().boundaryField_)
 {
     DebugInFunction
@@ -753,8 +698,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, gf),
     timeIndex_(gf.timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(this->mesh().boundary(), *this, patchFieldType)
 {
     DebugInFunction
@@ -765,7 +708,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 
     if (!readIfPresent() && gf.field0Ptr_)
     {
-        field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
+        field0Ptr_ = std::make_unique<this_type>
         (
             io.name() + "_0",
             *gf.field0Ptr_
@@ -785,8 +728,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, gf),
     timeIndex_(gf.timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_
     (
         this->mesh().boundary(),
@@ -803,7 +744,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 
     if (!readIfPresent() && gf.field0Ptr_)
     {
-        field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
+        field0Ptr_ = std::make_unique<this_type>
         (
             io.name() + "_0",
             *gf.field0Ptr_
@@ -823,8 +764,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, gf),
     timeIndex_(gf.timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_(*this, gf.boundaryField_, patchIDs, patchFieldType)
 {
     DebugInFunction
@@ -834,7 +773,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 
     if (!readIfPresent() && gf.field0Ptr_)
     {
-        field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
+        field0Ptr_ = std::make_unique<this_type>
         (
             io.name() + "_0",
             *gf.field0Ptr_
@@ -854,8 +793,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 :
     Internal(io, tgf.constCast(), tgf.movable()),
     timeIndex_(tgf().timeIndex()),
-    field0Ptr_(nullptr),
-    fieldPrevIterPtr_(nullptr),
     boundaryField_
     (
         this->mesh().boundary(),
@@ -897,13 +834,8 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::~GeometricField()
     }
     */
 
-    deleteDemandDrivenData(field0Ptr_);
-    deleteDemandDrivenData(fieldPrevIterPtr_);
-
     // FUTURE: register cache field info
     // // this->db().cacheTemporaryObject(*this);
-
-    clearOldTimes();
 }
 
 
@@ -1014,7 +946,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::oldTime() const
 {
     if (!field0Ptr_)
     {
-        field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
+        field0Ptr_ = std::make_unique<this_type>
         (
             IOobject
             (
@@ -1068,7 +1000,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::storePrevIter() const
             << "Allocating previous iteration field" << nl
             << this->info() << endl;
 
-        fieldPrevIterPtr_ = new GeometricField<Type, PatchField, GeoMesh>
+        fieldPrevIterPtr_ = std::make_unique<this_type>
         (
             this->name() + "PrevIter",
             *this
@@ -1101,8 +1033,8 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::prevIter() const
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::clearOldTimes()
 {
-    deleteDemandDrivenData(field0Ptr_);
-    deleteDemandDrivenData(fieldPrevIterPtr_);
+    field0Ptr_.reset(nullptr);
+    fieldPrevIterPtr_.reset(nullptr);
 }
 
 
@@ -1110,8 +1042,11 @@ template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::
 correctBoundaryConditions()
 {
-    this->setUpToDate();
-    storeOldTimes();
+    // updateAccessTime
+    {
+        this->setUpToDate();
+        storeOldTimes();
+    }
     boundaryField_.evaluate();
 }
 
@@ -1120,8 +1055,11 @@ template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::
 correctLocalBoundaryConditions()
 {
-    this->setUpToDate();
-    storeOldTimes();
+    // updateAccessTime
+    {
+        this->setUpToDate();
+        storeOldTimes();
+    }
     boundaryField_.evaluateLocal();
 }
 
@@ -1208,7 +1146,11 @@ template<class Type, template<class> class PatchField, class GeoMesh>
 bool Foam::GeometricField<Type, PatchField, GeoMesh>::
 writeData(Ostream& os) const
 {
-    os << *this;
+    this->internalField().writeData(os, "internalField");
+    os  << nl;
+    this->boundaryField().writeEntry("boundaryField", os);
+
+    os.check(FUNCTION_NAME);
     return os.good();
 }
 
@@ -1540,14 +1482,10 @@ template<class Type, template<class> class PatchField, class GeoMesh>
 Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
-    const GeometricField<Type, PatchField, GeoMesh>& gf
+    const GeometricField<Type, PatchField, GeoMesh>& fld
 )
 {
-    gf.internalField().writeData(os, "internalField");
-    os  << nl;
-    gf.boundaryField().writeEntry("boundaryField", os);
-
-    os.check(FUNCTION_NAME);
+    fld.writeData(os);
     return os;
 }
 
@@ -1556,11 +1494,12 @@ template<class Type, template<class> class PatchField, class GeoMesh>
 Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
-    const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
+    const tmp<GeometricField<Type, PatchField, GeoMesh>>& tfld
 )
 {
-    os << tgf();
-    tgf.clear();
+    tfld().writeData(os);
+    tfld.clear();
+
     return os;
 }
 
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
index 427009b245..8cb500f580 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
@@ -41,7 +41,6 @@ SourceFiles
 #ifndef Foam_GeometricField_H
 #define Foam_GeometricField_H
 
-#include "regIOobject.H"
 #include "GeometricBoundaryField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -50,7 +49,6 @@ namespace Foam
 {
 
 // Forward Declarations
-class dictionary;
 
 template<class Type, template<class> class PatchField, class GeoMesh>
 class GeometricField;
@@ -104,6 +102,12 @@ public:
 
 private:
 
+    // Private Typedefs
+
+        //- The GeometricField type
+        typedef GeometricField<Type, PatchField, GeoMesh> this_type;
+
+
     // Private Data
 
         //- Current time index.
@@ -111,10 +115,10 @@ private:
         mutable label timeIndex_;
 
         //- Pointer to old time field
-        mutable GeometricField<Type, PatchField, GeoMesh>* field0Ptr_;
+        mutable std::unique_ptr<this_type> field0Ptr_;
 
         //- Pointer to previous iteration (used for under-relaxation)
-        mutable GeometricField<Type, PatchField, GeoMesh>* fieldPrevIterPtr_;
+        mutable std::unique_ptr<this_type> fieldPrevIterPtr_;
 
         //- Boundary field containing boundary field values
         Boundary boundaryField_;
@@ -871,8 +875,8 @@ public:
         //- Correct boundary field
         void correctBoundaryConditions();
 
-        //- Correct boundary conditions after a purely local operation. Is
-        //  dummy for e.g. processor boundary conditions
+        //- Correct boundary conditions after a purely local operation.
+        //  Is dummy for processor boundary conditions etc
         void correctLocalBoundaryConditions();
 
         //- Does the field need a reference level for solution
@@ -1008,19 +1012,20 @@ public:
 
     // Write
 
-        //- The writeData function (required by regIOobject),
-        //- calls operator<<
+        //- The writeData function (required by regIOobject)
         bool writeData(Ostream& os) const;
 
 
     // Ostream Operators
 
+        //- Calls GeometricField::writeData()
         friend Ostream& operator<< <Type, PatchField, GeoMesh>
         (
             Ostream&,
             const GeometricField<Type, PatchField, GeoMesh>&
         );
 
+        //- Calls GeometricField::writeData()
         friend Ostream& operator<< <Type, PatchField, GeoMesh>
         (
             Ostream&,