diff --git a/applications/test/FieldFields1/Make/files b/applications/test/FieldFields1/Make/files
new file mode 100644
index 0000000000..e90b498508
--- /dev/null
+++ b/applications/test/FieldFields1/Make/files
@@ -0,0 +1,3 @@
+Test-FieldFields1.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-FieldFields1
diff --git a/applications/test/tensorFieldFields1/Make/options b/applications/test/FieldFields1/Make/options
similarity index 100%
rename from applications/test/tensorFieldFields1/Make/options
rename to applications/test/FieldFields1/Make/options
diff --git a/applications/test/FieldFields1/Test-FieldFields1.C b/applications/test/FieldFields1/Test-FieldFields1.C
new file mode 100644
index 0000000000..53f3314c35
--- /dev/null
+++ b/applications/test/FieldFields1/Test-FieldFields1.C
@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Application
+    Test-FieldFields1
+
+\*---------------------------------------------------------------------------*/
+
+#include "symmTensorField.H"
+#include "tensorField.H"
+#include "FieldFields.H"
+#include "Random.H"
+
+using namespace Foam;
+
+
+template<class Cmpt>
+void printFieldField(const FieldField<Field, Cmpt>& ff)
+{
+    forAll(ff, i)
+    {
+        Info<< i << ": " << flatOutput(ff[i]) << nl;
+    }
+    Info<< nl;
+}
+
+
+template<class Type>
+tmp<Field<Type>> randomField(Random& rnd, label dim)
+{
+    auto tfld = tmp<Field<Type>>::New(dim);
+    auto& fld = tfld.ref();
+
+    for (Type& val : fld)
+    {
+        for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; ++cmpt)
+        {
+            setComponent(val, cmpt) = rnd.position<label>(0, 100);
+        }
+    }
+
+    return tfld;
+}
+
+
+template<class Type>
+tmp<Field<Type>> randomField(label dim)
+{
+    Random rnd;
+    return randomField<Type>(rnd, dim);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    // scalarField
+    {
+        Info<< nl << "scalarFieldField" << nl;
+
+        Random rnd;
+
+        FieldField<Field, scalar> sff1(6);
+        forAll(sff1, i)
+        {
+            sff1.set(i, randomField<scalar>(rnd, 8));
+        }
+
+        printFieldField(sff1);
+
+        Info<< nl << "indexing:" << nl;
+
+        {
+            labelPair index;
+            const label range1 = sff1.size()-1;
+            const label range2 = sff1[0].size()-1;
+
+            for (label iter = 0; iter < 10; ++iter)
+            {
+                index.first() = rnd.position<label>(0, range1);
+                index.second() = rnd.position<label>(0, range2);
+
+                Info<< index << " => " << sff1[index] << nl;
+            }
+        }
+    }
+
+    Info<< nl << "End\n" << nl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/FieldFields2/Make/files b/applications/test/FieldFields2/Make/files
new file mode 100644
index 0000000000..d7f632aaf2
--- /dev/null
+++ b/applications/test/FieldFields2/Make/files
@@ -0,0 +1,3 @@
+Test-FieldFields2.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-FieldFields2
diff --git a/applications/test/FieldFields2/Make/options b/applications/test/FieldFields2/Make/options
new file mode 100644
index 0000000000..18e6fe47af
--- /dev/null
+++ b/applications/test/FieldFields2/Make/options
@@ -0,0 +1,2 @@
+/* EXE_INC = */
+/* EXE_LIBS = */
diff --git a/applications/test/tensorFieldFields1/Test-tensorFieldFields1.C b/applications/test/FieldFields2/Test-FieldFields2.C
similarity index 91%
rename from applications/test/tensorFieldFields1/Test-tensorFieldFields1.C
rename to applications/test/FieldFields2/Test-FieldFields2.C
index c7be55e686..4cf802a4f3 100644
--- a/applications/test/tensorFieldFields1/Test-tensorFieldFields1.C
+++ b/applications/test/FieldFields2/Test-FieldFields2.C
@@ -24,7 +24,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    Test-tensorFieldFields1
+    Test-FieldFields2
 
 \*---------------------------------------------------------------------------*/
 
@@ -73,21 +73,29 @@ void allocComponents
 }
 
 
-vectorField randomVectorField(label size)
+template<class Type>
+tmp<Field<Type>> randomField(Random& rnd, label dim)
 {
-    Random rnd;
+    auto tfld = tmp<Field<Type>>::New(dim);
+    auto& fld = tfld.ref();
 
-    vectorField vf(size);
-
-    forAll(vf, i)
+    for (Type& val : fld)
     {
-        for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
+        for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; ++cmpt)
         {
-            vf[i][cmpt] = rnd.position<label>(0, 100);
+            setComponent(val, cmpt) = rnd.position<label>(0, 100);
         }
     }
 
-    return vf;
+    return tfld;
+}
+
+
+template<class Type>
+tmp<Field<Type>> randomField(label dim)
+{
+    Random rnd;
+    return randomField<Type>(rnd, dim);
 }
 
 
@@ -191,7 +199,7 @@ int main(int argc, char *argv[])
         printFieldField(sf1);
 
         Info<< nl;
-        for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
+        for (direction cmpt = 0; cmpt < pTraits<vector>::nComponents; ++cmpt)
         {
             unzipRow(sf1, cmpt, slice[0]);
 
@@ -253,8 +261,9 @@ int main(int argc, char *argv[])
     {
         Info<< nl << "vectorField" << nl;
 
+        Random rnd;
         FieldField<Field, vector> vf1(1);
-        vf1.set(0, new vectorField(randomVectorField(4)));
+        vf1.set(0, randomField<vector>(rnd, 4));
 
         FixedList<FieldField<Field, scalar>, 3> cmpts;
         allocComponents(cmpts, 4);
diff --git a/applications/test/tensorFieldFields1/Make/files b/applications/test/tensorFieldFields1/Make/files
deleted file mode 100644
index 0bc78a8cea..0000000000
--- a/applications/test/tensorFieldFields1/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-Test-tensorFieldFields1.C
-
-EXE = $(FOAM_USER_APPBIN)/Test-tensorFieldFields1
diff --git a/applications/test/tensorFields1/Test-tensorFields1.C b/applications/test/tensorFields1/Test-tensorFields1.C
index 60dadb2629..e92689422c 100644
--- a/applications/test/tensorFields1/Test-tensorFields1.C
+++ b/applications/test/tensorFields1/Test-tensorFields1.C
@@ -34,21 +34,29 @@ Application
 
 using namespace Foam;
 
-vectorField randomVectorField(label size)
+template<class Type>
+tmp<Field<Type>> randomField(Random& rnd, label dim)
 {
-    Random rnd;
+    auto tfld = tmp<Field<Type>>::New(dim);
+    auto& fld = tfld.ref();
 
-    vectorField vf(size);
-
-    forAll(vf, i)
+    for (Type& val : fld)
     {
-        for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
+        for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; ++cmpt)
         {
-            vf[i][cmpt] = rnd.position<label>(0, 100);
+            setComponent(val, cmpt) = rnd.position<label>(0, 100);
         }
     }
 
-    return vf;
+    return tfld;
+}
+
+
+template<class Type>
+tmp<Field<Type>> randomField(label dim)
+{
+    Random rnd;
+    return randomField<Type>(rnd, dim);
 }
 
 
@@ -251,7 +259,7 @@ int main(int argc, char *argv[])
     {
         Info<< nl << "vectorField" << nl;
 
-        vectorField vf1(randomVectorField(4));
+        vectorField vf1(randomField<vector>(4));
         FixedList<scalarField, 3> cmpts(scalarField(vf1.size()));
 
         Info<< nl
diff --git a/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.C b/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.C
index 4905f9c03a..16d16264e4 100644
--- a/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.C
+++ b/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.C
@@ -297,6 +297,20 @@ tmp<FieldField<Field, Type>> FieldField<Field, Type>::T() const
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
+template<template<class> class Field, class Type>
+const Type& FieldField<Field, Type>::operator[](const labelPair& index) const
+{
+    return this->operator[](index.first())[index.second()];
+}
+
+
+template<template<class> class Field, class Type>
+Type& FieldField<Field, Type>::operator[](const labelPair& index)
+{
+    return this->operator[](index.first())[index.second()];
+}
+
+
 template<template<class> class Field, class Type>
 void FieldField<Field, Type>::operator=(const FieldField<Field, Type>& ff)
 {
@@ -356,6 +370,16 @@ void FieldField<Field, Type>::operator=(const Type& val)
 }
 
 
+template<template<class> class Field, class Type>
+void FieldField<Field, Type>::operator=(const Foam::zero)
+{
+    forAll(*this, i)
+    {
+        this->operator[](i) = Foam::zero{};
+    }
+}
+
+
 #define COMPUTED_ASSIGNMENT(TYPE, op)                                          \
                                                                                \
 template<template<class> class Field, class Type>                              \
diff --git a/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.H b/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.H
index 7b408ceca5..2227f295fa 100644
--- a/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.H
+++ b/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.H
@@ -39,6 +39,7 @@ SourceFiles
 #define Foam_FieldField_H
 
 #include "tmp.H"
+#include "Pair.H"
 #include "PtrList.H"
 #include "scalar.H"
 #include "direction.H"
@@ -79,7 +80,6 @@ class FieldField
     public refCount,
     public PtrList<Field<Type>>
 {
-
 public:
 
     //- Component type
@@ -153,7 +153,16 @@ public:
         tmp<FieldField<Field, Type>> T() const;
 
 
-    // Member operators
+    // Member Operators
+
+        //- Const or non-const access to a field
+        using PtrList<Field<Type>>::operator[];
+
+        //- Const access to a single field element via (fieldi, elemi)
+        inline const Type& operator[](const labelPair& index) const;
+
+        //- Non-const access to a single field element via (fieldi, elemi)
+        inline Type& operator[](const labelPair& index);
 
         //- Copy assignment
         void operator=(const FieldField<Field, Type>&);
@@ -167,6 +176,9 @@ public:
         //- Assign uniform value
         void operator=(const Type& val);
 
+        //- Assign uniform zero
+        void operator=(const Foam::zero);
+
         void operator+=(const FieldField<Field, Type>&);
         void operator+=(const tmp<FieldField<Field, Type>>&);
 
diff --git a/src/OpenFOAM/fields/FieldFields/oneFieldField/oneFieldField.H b/src/OpenFOAM/fields/FieldFields/oneFieldField/oneFieldField.H
index c6d31e6aca..c68c1a0ede 100644
--- a/src/OpenFOAM/fields/FieldFields/oneFieldField/oneFieldField.H
+++ b/src/OpenFOAM/fields/FieldFields/oneFieldField/oneFieldField.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,10 +37,11 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef oneFieldField_H
-#define oneFieldField_H
+#ifndef Foam_oneFieldField_H
+#define Foam_oneFieldField_H
 
 #include "oneField.H"
+#include "Pair.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,10 +66,17 @@ public:
 
     // Member Operators
 
+        //- Return the field
         oneField operator[](const label) const noexcept
         {
             return oneField{};
         }
+
+        //- Return the field element
+        one operator[](const labelPair&) const noexcept
+        {
+            return one{};
+        }
 };
 
 
diff --git a/src/OpenFOAM/fields/FieldFields/zeroFieldField/zeroFieldField.H b/src/OpenFOAM/fields/FieldFields/zeroFieldField/zeroFieldField.H
index 1a2aa0f5a5..4834de4e12 100644
--- a/src/OpenFOAM/fields/FieldFields/zeroFieldField/zeroFieldField.H
+++ b/src/OpenFOAM/fields/FieldFields/zeroFieldField/zeroFieldField.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,10 +37,11 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef zeroFieldField_H
-#define zeroFieldField_H
+#ifndef Foam_zeroFieldField_H
+#define Foam_zeroFieldField_H
 
 #include "zeroField.H"
+#include "Pair.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,10 +66,17 @@ public:
 
     // Member Operators
 
+        //- Return the field
         zeroField operator[](const label) const noexcept
         {
             return zeroField{};
         }
+
+        //- Return the field element
+        zero operator[](const labelPair&) const noexcept
+        {
+            return zero{};
+        }
 };
 
 
diff --git a/src/OpenFOAM/fields/Fields/oneField/oneField.H b/src/OpenFOAM/fields/Fields/oneField/oneField.H
index 523437756e..6358ec3f38 100644
--- a/src/OpenFOAM/fields/Fields/oneField/oneField.H
+++ b/src/OpenFOAM/fields/Fields/oneField/oneField.H
@@ -36,11 +36,11 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef oneField_H
-#define oneField_H
+#ifndef Foam_oneField_H
+#define Foam_oneField_H
 
 #include "one.H"
-#include "scalar.H"
+#include "scalarFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/fields/Fields/zeroField/zeroField.H b/src/OpenFOAM/fields/Fields/zeroField/zeroField.H
index 2969ec1393..70d8ab0b3f 100644
--- a/src/OpenFOAM/fields/Fields/zeroField/zeroField.H
+++ b/src/OpenFOAM/fields/Fields/zeroField/zeroField.H
@@ -36,11 +36,11 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef zeroField_H
-#define zeroField_H
+#ifndef Foam_zeroField_H
+#define Foam_zeroField_H
 
 #include "zero.H"
-#include "scalar.H"
+#include "scalarFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/primitives/tuples/Pair.H b/src/OpenFOAM/primitives/tuples/Pair.H
index 9a9458bc22..01307ba7c7 100644
--- a/src/OpenFOAM/primitives/tuples/Pair.H
+++ b/src/OpenFOAM/primitives/tuples/Pair.H
@@ -39,8 +39,8 @@ See also
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Pair_H
-#define Pair_H
+#ifndef Foam_Pair_H
+#define Foam_Pair_H
 
 #include "FixedList.H"
 #include "Istream.H"
diff --git a/src/OpenFOAM/primitives/tuples/Tuple2.H b/src/OpenFOAM/primitives/tuples/Tuple2.H
index 874f34945d..bd63390228 100644
--- a/src/OpenFOAM/primitives/tuples/Tuple2.H
+++ b/src/OpenFOAM/primitives/tuples/Tuple2.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,13 +37,12 @@ See also
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Tuple2_H
-#define Tuple2_H
+#ifndef Foam_Tuple2_H
+#define Foam_Tuple2_H
 
 #include "Istream.H"
 #include "Ostream.H"
 #include "Pair.H"
-#include <utility>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -79,80 +78,61 @@ public:
         Tuple2() = default;
 
         //- Copy construct from components
-        inline Tuple2(const T1& f, const T2& s)
+        Tuple2(const T1& f, const T2& s)
         :
             f_(f),
             s_(s)
         {}
 
         //- Move construct from components
-        inline Tuple2(T1&& f, T2&& s)
+        Tuple2(T1&& f, T2&& s)
         :
             f_(std::move(f)),
             s_(std::move(s))
         {}
 
         //- Copy construct from std::pair
-        inline Tuple2(const std::pair<T1,T2>& vals)
+        Tuple2(const std::pair<T1,T2>& vals)
         :
             f_(vals.first),
             s_(vals.second)
         {}
 
         //- Move construct from std::pair
-        inline Tuple2(std::pair<T1,T2>&& vals)
+        Tuple2(std::pair<T1,T2>&& vals)
         :
             f_(std::move(vals.first)),
             s_(std::move(vals.second))
         {}
 
         //- Construct from Istream
-        inline explicit Tuple2(Istream& is)
+        explicit Tuple2(Istream& is)
         {
-            is >> *this;
+            is.readBegin("Tuple2");
+            is >> f_ >> s_;
+            is.readEnd("Tuple2");
+            is.check(FUNCTION_NAME);
         }
 
 
     // Member Functions
 
         //- Return first
-        inline const T1& first() const noexcept
-        {
-            return f_;
-        }
+        const T1& first() const noexcept { return f_; }
 
         //- Return first
-        inline T1& first() noexcept
-        {
-            return f_;
-        }
+        T1& first() noexcept { return f_; }
 
         //- Return second
-        inline const T2& second() const noexcept
-        {
-            return s_;
-        }
+        const T2& second() const noexcept { return s_; }
 
         //- Return second
-        inline T2& second() noexcept
-        {
-            return s_;
-        }
+        T2& second() noexcept { return s_; }
 };
 
 
 // * * * * * * * * * * * * * * * * * Traits  * * * * * * * * * * * * * * * * //
 
-//- Hashing for Tuple2 data
-template<class T1, class T2>
-struct Hash<Tuple2<T1, T2>>
-{
-    unsigned operator()(const Tuple2<T1, T2>& obj, unsigned seed=0) const
-    {
-        return Hash<T2>()(obj.second(), Hash<T1>()(obj.first(), seed));
-    }
-};
-
 //- Hashing for std::pair data
 template<class T1, class T2>
 struct Hash<std::pair<T1, T2>>
@@ -164,6 +144,17 @@ struct Hash<std::pair<T1, T2>>
 };
 
 
+//- Hashing for Tuple2 data
+template<class T1, class T2>
+struct Hash<Tuple2<T1, T2>>
+{
+    unsigned operator()(const Tuple2<T1, T2>& obj, unsigned seed=0) const
+    {
+        return Hash<T2>()(obj.second(), Hash<T1>()(obj.first(), seed));
+    }
+};
+
+
 // * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * * //
 
 //- Return reverse of a Tuple2
@@ -300,7 +291,20 @@ struct maxFirstEqOp
 
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
-//- Read Tuple2 from Istream, discarding contents of existing Tuple2.
+//- Read std::pair from Istream
+template<class T1, class T2>
+inline Istream& operator>>(Istream& is, std::pair<T1,T2>& t)
+{
+    is.readBegin("pair");
+    is >> t.first >> t.second;
+    is.readEnd("pair");
+
+    is.check(FUNCTION_NAME);
+    return is;
+}
+
+
+//- Read Tuple2 from Istream
 template<class T1, class T2>
 inline Istream& operator>>(Istream& is, Tuple2<T1,T2>& t)
 {
@@ -313,16 +317,14 @@ inline Istream& operator>>(Istream& is, Tuple2<T1,T2>& t)
 }
 
 
-//- Read std::pair from Istream
+//- Write std::pair to Ostream.
 template<class T1, class T2>
-inline Istream& operator>>(Istream& is, std::pair<T1,T2>& t)
+inline Ostream& operator<<(Ostream& os, const std::pair<T1,T2>& t)
 {
-    is.readBegin("std::pair");
-    is >> t.first >> t.second;
-    is.readEnd("std::pair");
-
-    is.check(FUNCTION_NAME);
-    return is;
+    os  << token::BEGIN_LIST
+        << t.first << token::SPACE << t.second
+        << token::END_LIST;
+    return os;
 }
 
 
@@ -333,19 +335,6 @@ inline Ostream& operator<<(Ostream& os, const Tuple2<T1,T2>& t)
     os  << token::BEGIN_LIST
         << t.first() << token::SPACE << t.second()
         << token::END_LIST;
-
-    return os;
-}
-
-
-//- Write std::pair to Ostream.
-template<class T1, class T2>
-inline Ostream& operator<<(Ostream& os, const std::pair<T1,T2>& t)
-{
-    os  << token::BEGIN_LIST
-        << t.first << token::SPACE << t.second
-        << token::END_LIST;
-
     return os;
 }