diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C
index 6322c38f7c..cfdb76185d 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C
@@ -77,7 +77,7 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
 )
 :
     regIOobject(io),
-    Field<Type>(field),
+    DynamicField<Type>(field),
     mesh_(mesh),
     dimensions_(dims)
 {
@@ -95,7 +95,7 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
 )
 :
     regIOobject(io),
-    Field<Type>(std::move(field)),
+    DynamicField<Type>(std::move(field)),
     mesh_(mesh),
     dimensions_(dims)
 {
@@ -113,7 +113,7 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
 )
 :
     regIOobject(io),
-    Field<Type>(std::move(field)),
+    DynamicField<Type>(std::move(field)),
     mesh_(mesh),
     dimensions_(dims)
 {
@@ -131,7 +131,7 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
 )
 :
     regIOobject(io),
-    Field<Type>(tfield.constCast(), tfield.movable()),
+    DynamicField<Type>(tfield.constCast(), tfield.movable()),
     mesh_(mesh),
     dimensions_(dims)
 {
@@ -146,11 +146,21 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
     const IOobject& io,
     const Mesh& mesh,
     const dimensionSet& dims,
-    const bool checkIOFlags
+    const bool checkIOFlags,
+    const bool extraCapacity
 )
 :
     regIOobject(io),
-    Field<Type>(GeoMesh::size(mesh)),
+    DynamicField<Type>
+    (
+        // (size,capacity)
+        std::pair<label,label>
+        (
+            GeoMesh::size(mesh),
+            GeoMesh::size(mesh)
+          + (extraCapacity ? GeoMesh::boundary_size(mesh) : label(0))
+        )
+    ),
     mesh_(mesh),
     dimensions_(dims)
 {
@@ -168,11 +178,21 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
     const Mesh& mesh,
     const Type& value,
     const dimensionSet& dims,
-    const bool checkIOFlags
+    const bool checkIOFlags,
+    const bool extraCapacity
 )
 :
     regIOobject(io),
-    Field<Type>(GeoMesh::size(mesh)),
+    DynamicField<Type>
+    (
+        // (size,capacity)
+        std::pair<label,label>
+        (
+            GeoMesh::size(mesh),
+            GeoMesh::size(mesh)
+          + (extraCapacity ? GeoMesh::boundary_size(mesh) : label(0))
+        )
+    ),
     mesh_(mesh),
     dimensions_(dims)
 {
@@ -190,7 +210,8 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
     const IOobject& io,
     const Mesh& mesh,
     const dimensioned<Type>& dt,
-    const bool checkIOFlags
+    const bool checkIOFlags,
+    const bool extraCapacity
 )
 :
     DimensionedField<Type, GeoMesh>
@@ -199,7 +220,8 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
         mesh,
         dt.value(),
         dt.dimensions(),
-        checkIOFlags
+        checkIOFlags,
+        extraCapacity
     )
 {}
 
@@ -211,7 +233,7 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
 )
 :
     regIOobject(df),
-    Field<Type>(df),
+    DynamicField<Type>(df),
     mesh_(df.mesh_),
     dimensions_(df.dimensions_),
     oriented_(df.oriented_)
@@ -236,7 +258,7 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
 )
 :
     regIOobject(df, reuse),
-    Field<Type>(df, reuse),
+    DynamicField<Type>(df, reuse),
     mesh_(df.mesh_),
     dimensions_(df.dimensions_),
     oriented_(df.oriented_)
@@ -263,7 +285,7 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
 )
 :
     regIOobject(io),
-    Field<Type>(df),
+    DynamicField<Type>(df),
     mesh_(df.mesh_),
     dimensions_(df.dimensions_),
     oriented_(df.oriented_)
@@ -290,7 +312,7 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
 )
 :
     regIOobject(io, df),
-    Field<Type>(df, reuse),
+    DynamicField<Type>(df, reuse),
     mesh_(df.mesh_),
     dimensions_(df.dimensions_),
     oriented_(df.oriented_)
@@ -318,7 +340,7 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
 )
 :
     regIOobject(newName, df, newName != df.name()),
-    Field<Type>(df),
+    DynamicField<Type>(df),
     mesh_(df.mesh_),
     dimensions_(df.dimensions_),
     oriented_(df.oriented_)
@@ -345,7 +367,7 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
 )
 :
     regIOobject(newName, df, true),
-    Field<Type>(df, reuse),
+    DynamicField<Type>(df, reuse),
     mesh_(df.mesh_),
     dimensions_(df.dimensions_),
     oriented_(df.oriented_)
@@ -421,7 +443,7 @@ void Foam::DimensionedField<Type, GeoMesh>::replace
     >& df
 )
 {
-    Field<Type>::replace(d, df);
+    DynamicField<Type>::replace(d, df);
 }
 
 
@@ -523,7 +545,7 @@ void Foam::DimensionedField<Type, GeoMesh>::operator=
 
     dimensions_ = df.dimensions();
     oriented_ = df.oriented();
-    Field<Type>::operator=(df);
+    DynamicField<Type>::operator=(df);
 }
 
 
@@ -556,14 +578,14 @@ void Foam::DimensionedField<Type, GeoMesh>::operator=
 )
 {
     dimensions_ = dt.dimensions();
-    Field<Type>::operator=(dt.value());
+    DynamicField<Type>::operator=(dt.value());
 }
 
 
 template<class Type, class GeoMesh>
 void Foam::DimensionedField<Type, GeoMesh>::operator=(Foam::zero)
 {
-    Field<Type>::operator=(Foam::zero{});
+    DynamicField<Type>::operator=(Foam::zero{});
 }
 
 
@@ -579,7 +601,7 @@ void Foam::DimensionedField<Type, GeoMesh>::operator op                        \
                                                                                \
     dimensions_ op df.dimensions();                                            \
     oriented_ op df.oriented();                                                \
-    Field<Type>::operator op(df);                                              \
+    DynamicField<Type>::operator op(df);                                       \
 }                                                                              \
                                                                                \
 template<class Type, class GeoMesh>                                            \
@@ -599,7 +621,7 @@ void Foam::DimensionedField<Type, GeoMesh>::operator op                        \
 )                                                                              \
 {                                                                              \
     dimensions_ op dt.dimensions();                                            \
-    Field<Type>::operator op(dt.value());                                      \
+    DynamicField<Type>::operator op(dt.value());                               \
 }
 
 COMPUTED_ASSIGNMENT(Type, +=)
diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H
index a2741b89e9..3112500947 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H
@@ -42,7 +42,7 @@ SourceFiles
 #define Foam_DimensionedField_H
 
 #include "regIOobject.H"
-#include "Field.H"
+#include "DynamicField.H"
 #include "dimensionedType.H"
 #include "orientedType.H"
 
@@ -77,7 +77,7 @@ template<class Type, class GeoMesh>
 class DimensionedField
 :
     public regIOobject,
-    public Field<Type>
+    public DynamicField<Type>
 {
 public:
 
@@ -93,6 +93,7 @@ public:
         typedef DimensionedField<Type, GeoMesh> Internal;
 
         //- Type of the field from which this DimensionedField is derived
+        //  Use Field (not DynamicField) here.
         typedef Field<Type> FieldType;
 
         //- Component type of the field elements
@@ -197,7 +198,10 @@ public:
             const IOobject& io,
             const Mesh& mesh,
             const dimensionSet& dims,
-            const bool checkIOFlags = true
+            //! Handle mandatory/optional reading?
+            const bool checkIOFlags = true,
+            //! Additional space for geometric boundary size?
+            const bool extraCapacity = false
         );
 
         //- Construct from components, setting dimensions and initial
@@ -208,7 +212,10 @@ public:
             const Mesh& mesh,
             const Type& value,
             const dimensionSet& dims,
-            const bool checkIOFlags = true
+            //! Handle mandatory/optional reading?
+            const bool checkIOFlags = true,
+            //! Additional space for geometric boundary size?
+            const bool extraCapacity = false
         );
 
         //- Construct from components, setting dimensions and initial
@@ -219,15 +226,20 @@ public:
             const IOobject& io,
             const Mesh& mesh,
             const dimensioned<Type>& dt,
-            const bool checkIOFlags = true
+            //! Handle mandatory/optional reading?
+            const bool checkIOFlags = true,
+            //! Additional space for geometric boundary size?
+            const bool extraCapacity = false
         );
 
-        //- Construct from Istream.
+        //- Construct from Istream (uses an intermediate dictionary)
         DimensionedField
         (
             const IOobject& io,
             const Mesh& mesh,
-            const word& fieldDictEntry = "value"
+            const word& fieldDictEntry = "value",
+            //! Additional space for geometric boundary size?
+            const bool extraCapacity = false
         );
 
         //- Construct from dictionary
@@ -236,7 +248,9 @@ public:
             const IOobject& io,
             const Mesh& mesh,
             const dictionary& fieldDict,
-            const word& fieldDictEntry = "value"
+            const word& fieldDictEntry = "value",
+            //! Additional space for geometric boundary size?
+            const bool extraCapacity = false
         );
 
         //- Copy construct
@@ -572,6 +586,15 @@ public:
         //- Return reference to the primitive field values
         Field<Type>& field() noexcept { return *this; }
 
+        //- Return const-reference to the primitive field storage.
+        //- Use sparingly.
+        const DynamicField<Type>& dyn_field() const noexcept { return *this; }
+
+        //- Return reference to the primitive field storage.
+        //- Use sparingly.
+        DynamicField<Type>& dyn_field() noexcept { return *this; }
+
+
         //- Return a component field of the field
         tmp<DimensionedField<cmptType, GeoMesh>> component
         (
diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.C b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.C
index b72557a381..8d68de252e 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.C
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.C
@@ -50,12 +50,19 @@ void Foam::DimensionedField<Type, GeoMesh>::readField
         oriented_.read(fieldDict);  // The "oriented" entry (if present)
     }
 
+    const label meshSize = GeoMesh::size(mesh_);
 
-    // The primitive field
-    auto& fld = static_cast<Field<Type>&>(*this);
+    // The primitive field : dyn_field()
+    auto& fld = static_cast<DynamicField<Type>&>(*this);
 
-    fld.resize_nocopy(GeoMesh::size(mesh_));
-    fld.assign(fieldDictEntry, fieldDict, fld.size());  // <- MUST_READ
+    // Resize primitive field to make space for the internal field
+    // - avoid size doubling
+    // - without copying any existing content
+
+    fld.clear();  // ie, ignore any existing content
+    fld.reserve_exact(meshSize);
+    fld.resize_nocopy(meshSize);
+    fld.assign(fieldDictEntry, fieldDict, meshSize);  // <- MUST_READ
 }
 
 
@@ -116,12 +123,20 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
 (
     const IOobject& io,
     const Mesh& mesh,
-    const word& fieldDictEntry
+    const word& fieldDictEntry,
+    const bool extraCapacity
 )
 :
     regIOobject(io),
     mesh_(mesh)
 {
+    if (extraCapacity)
+    {
+        DynamicField<Type>::reserve_exact
+        (
+            GeoMesh::size(mesh) + GeoMesh::boundary_size(mesh)
+        );
+    }
     readField(fieldDictEntry);
 }
 
@@ -132,12 +147,20 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
     const IOobject& io,
     const Mesh& mesh,
     const dictionary& fieldDict,
-    const word& fieldDictEntry
+    const word& fieldDictEntry,
+    const bool extraCapacity
 )
 :
     regIOobject(io),
     mesh_(mesh)
 {
+    if (extraCapacity)
+    {
+        DynamicField<Type>::reserve_exact
+        (
+            GeoMesh::size(mesh) + GeoMesh::boundary_size(mesh)
+        );
+    }
     readField(fieldDict, fieldDictEntry);
 }
 
@@ -159,7 +182,7 @@ bool Foam::DimensionedField<Type, GeoMesh>::writeData
         os << nl;
     }
 
-    Field<Type>::writeEntry(fieldDictEntry, os);
+    DynamicField<Type>::writeEntry(fieldDictEntry, os);
 
     os.check(FUNCTION_NAME);
     return os.good();
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
index e6fe9c335e..69d9dc30f3 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
@@ -182,7 +182,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const word& patchFieldType
 )
 :
-    Internal(io, mesh, dims, false),
+    Internal(io, mesh, dims, false, FieldBase::unifiedGeometricField),
     timeIndex_(this->time().timeIndex()),
     boundaryField_(mesh.boundary(), *this, patchFieldType)
 {
@@ -203,7 +203,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const wordList& actualPatchTypes
 )
 :
-    Internal(io, mesh, dims, false),
+    Internal(io, mesh, dims, false, FieldBase::unifiedGeometricField),
     timeIndex_(this->time().timeIndex()),
     boundaryField_(mesh.boundary(), *this, patchFieldTypes, actualPatchTypes)
 {
@@ -224,7 +224,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const word& patchFieldType
 )
 :
-    Internal(io, mesh, value, dims, false),
+    Internal(io, mesh, value, dims, false, FieldBase::unifiedGeometricField),
     timeIndex_(this->time().timeIndex()),
     boundaryField_(mesh.boundary(), *this, patchFieldType)
 {
@@ -248,7 +248,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const wordList& actualPatchTypes
 )
 :
-    Internal(io, mesh, value, dims, false),
+    Internal(io, mesh, value, dims, false, FieldBase::unifiedGeometricField),
     timeIndex_(this->time().timeIndex()),
     boundaryField_(mesh.boundary(), *this, patchFieldTypes, actualPatchTypes)
 {
@@ -530,7 +530,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const bool readOldTime
 )
 :
-    Internal(io, mesh, dimless, false),
+    Internal(io, mesh, dimless, false, FieldBase::unifiedGeometricField),
     timeIndex_(this->time().timeIndex()),
     boundaryField_(mesh.boundary())
 {
@@ -567,7 +567,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const dictionary& dict
 )
 :
-    Internal(io, mesh, dimless, false),
+    Internal(io, mesh, dimless, false, FieldBase::unifiedGeometricField),
     timeIndex_(this->time().timeIndex()),
     boundaryField_(mesh.boundary())
 {
@@ -1087,6 +1087,46 @@ correctLocalBoundaryConditions()
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+template<class Cop>
+Foam::label Foam::GeometricField<Type, PatchField, GeoMesh>::
+boundaryEvaluate(const Cop& cop)
+{
+    const label meshSize = GeoMesh::size(this->mesh());
+    const label totalSize = GeoMesh::boundary_size(this->mesh()) + meshSize;
+
+    if (FOAM_UNLIKELY(meshSize != this->size() && totalSize != this->size()))
+    {
+        FatalErrorInFunction
+            << "Problem : field:" << this->name()
+            << " size:" << this->size()
+            << " capacity:" << this->capacity()
+            << " is not mesh size:" << meshSize
+            << " or total size:" << totalSize
+            << exit(FatalError);
+    }
+
+    auto& fld = static_cast<DynamicField<Type>&>(*this);
+
+    // Resize primitive field to allow for internal+boundary fields
+    // - avoid size doubling
+    // - only copying the internal field without boundaries
+
+    fld.reserve_exact(totalSize);
+    fld.resize_copy(meshSize, totalSize);
+
+    // Populate the extra space with the flattened boundary values:
+    for (const auto& pfld : this->boundaryField())
+    {
+        const label start = (meshSize + pfld.patch().offset());
+        SubList<Type> slice(fld, pfld.size(), start);
+        cop(pfld, slice);
+    }
+
+    return meshSize;
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 bool Foam::GeometricField<Type, PatchField, GeoMesh>::needReference() const
 {
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
index dc4b0caaf2..9f42d3f79f 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
@@ -923,6 +923,14 @@ public:
         //  Is dummy for processor boundary conditions etc
         void correctLocalBoundaryConditions();
 
+        //- Evaluate boundary functions using the field storage:
+        //  - extends the internal field storage
+        //  - applies passed-in operator to fill in the slots
+        //  .
+        //  \returns The nominal geometric field size
+        template<class Cop>
+        label boundaryEvaluate(const Cop& cop);
+
         //- Does the field need a reference level for solution
         bool needReference() const;