diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldReuseFunctions.H b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldReuseFunctions.H
index a4647cb1c82d52e11e2c2c76500457d2848ed1c9..ce01f1b6b3a81738fe6d8f2063b43a697e978f7a 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldReuseFunctions.H
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldReuseFunctions.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,6 +40,15 @@ namespace Foam
 template<class TypeR, class Type1, class GeoMesh>
 struct reuseTmpDimensionedField
 {
+    //- Dissimilar types: non-reusable
+    static bool reusable
+    (
+        const tmp<DimensionedField<Type1, GeoMesh>>& tf1
+    )
+    {
+        return false;
+    }
+
     //- Pass-through to New DimensionedField
     static tmp<DimensionedField<TypeR, GeoMesh>> New
     (
@@ -79,6 +88,15 @@ struct reuseTmpDimensionedField
 template<class TypeR, class GeoMesh>
 struct reuseTmpDimensionedField<TypeR, TypeR, GeoMesh>
 {
+    //- Identical types: possibly reusable
+    static bool reusable
+    (
+        const tmp<DimensionedField<TypeR, GeoMesh>>& tf1
+    )
+    {
+        return tf1.movable();
+    }
+
     //- Allow optional copy assignment of the initial content
     //- for identical input and output types
     static tmp<DimensionedField<TypeR, GeoMesh>> New
@@ -137,9 +155,22 @@ tmp<DimensionedField<TypeR, GeoMesh>> New
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Two-parameter versions
+
 template<class TypeR, class Type1, class Type12, class Type2, class GeoMesh>
 struct reuseTmpTmpDimensionedField
 {
+    //- Dissimilar types: non-reusable
+    static bool reusable
+    (
+        const tmp<DimensionedField<Type1, GeoMesh>>& tf1,
+        const tmp<DimensionedField<Type2, GeoMesh>>& tf2
+    )
+    {
+        return false;
+    }
+
     static tmp<DimensionedField<TypeR, GeoMesh>> New
     (
         const tmp<DimensionedField<Type1, GeoMesh>>& tf1,
@@ -163,6 +194,17 @@ struct reuseTmpTmpDimensionedField
 template<class TypeR, class Type1, class Type12, class GeoMesh>
 struct reuseTmpTmpDimensionedField<TypeR, Type1, Type12, TypeR, GeoMesh>
 {
+    //- Second input has return type: possibly reusable
+    static bool reusable
+    (
+        const tmp<DimensionedField<Type1, GeoMesh>>& tf1,
+        const tmp<DimensionedField<TypeR, GeoMesh>>& tf2
+    )
+    {
+        return tf2.movable();
+    }
+
+    //- Second input has return type
     static tmp<DimensionedField<TypeR, GeoMesh>> New
     (
         const tmp<DimensionedField<Type1, GeoMesh>>& tf1,
@@ -195,6 +237,17 @@ struct reuseTmpTmpDimensionedField<TypeR, Type1, Type12, TypeR, GeoMesh>
 template<class TypeR, class Type2, class GeoMesh>
 struct reuseTmpTmpDimensionedField<TypeR, TypeR, TypeR, Type2, GeoMesh>
 {
+    //- First input has return type: possibly reusable
+    static bool reusable
+    (
+        const tmp<DimensionedField<TypeR, GeoMesh>>& tf1,
+        const tmp<DimensionedField<Type2, GeoMesh>>& tf2
+    )
+    {
+        return tf1.movable();
+    }
+
+    //- First input has return type
     static tmp<DimensionedField<TypeR, GeoMesh>> New
     (
         const tmp<DimensionedField<TypeR, GeoMesh>>& tf1,
@@ -227,6 +280,17 @@ struct reuseTmpTmpDimensionedField<TypeR, TypeR, TypeR, Type2, GeoMesh>
 template<class TypeR, class GeoMesh>
 struct reuseTmpTmpDimensionedField<TypeR, TypeR, TypeR, TypeR, GeoMesh>
 {
+    //- Both inputs have return type: possibly reusable
+    static bool reusable
+    (
+        const tmp<DimensionedField<TypeR, GeoMesh>>& tf1,
+        const tmp<DimensionedField<TypeR, GeoMesh>>& tf2
+    )
+    {
+        return (tf1.movable() || tf2.movable());
+    }
+
+    //- Both inputs have return type
     static tmp<DimensionedField<TypeR, GeoMesh>> New
     (
         const tmp<DimensionedField<TypeR, GeoMesh>>& tf1,
diff --git a/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctions.C b/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctions.C
index 54a489deac158c79b665e8ca553560c3338defb9..2b0fa83aa54488188d22a2130d1017305b6f5dad 100644
--- a/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctions.C
+++ b/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctions.C
@@ -47,9 +47,9 @@ void component
     const direction d
 )
 {
-    const label loopLen = (sf).size();
+    const label loop_len = (sf).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         component(sf[i], f[i], d);
     }
@@ -59,9 +59,9 @@ void component
 template<template<class> class Field, class Type>
 void T(FieldField<Field, Type>& f1, const FieldField<Field, Type>& f2)
 {
-    const label loopLen = (f1).size();
+    const label loop_len = (f1).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         T(f1[i], f2[i]);
     }
@@ -75,9 +75,9 @@ void pow
     const FieldField<Field, Type>& vf
 )
 {
-    const label loopLen = (f).size();
+    const label loop_len = (f).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         pow(f[i], vf[i]);
     }
@@ -122,9 +122,9 @@ void sqr
     const FieldField<Field, Type>& vf
 )
 {
-    const label loopLen = (f).size();
+    const label loop_len = (f).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         sqr(f[i], vf[i]);
     }
@@ -163,9 +163,9 @@ void magSqr
     const FieldField<Field, Type>& f
 )
 {
-    const label loopLen = (sf).size();
+    const label loop_len = (sf).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         magSqr(sf[i], f[i]);
     }
@@ -204,9 +204,9 @@ void mag
     const FieldField<Field, Type>& f
 )
 {
-    const label loopLen = (sf).size();
+    const label loop_len = (sf).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         mag(sf[i], f[i]);
     }
@@ -245,9 +245,9 @@ void cmptMax
     const FieldField<Field, Type>& f
 )
 {
-    const label loopLen = (cf).size();
+    const label loop_len = (cf).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         cmptMax(cf[i], f[i]);
     }
@@ -290,9 +290,9 @@ void cmptMin
     const FieldField<Field, Type>& f
 )
 {
-    const label loopLen = (cf).size();
+    const label loop_len = (cf).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         cmptMin(cf[i], f[i]);
     }
@@ -335,9 +335,9 @@ void cmptAv
     const FieldField<Field, Type>& f
 )
 {
-    const label loopLen = (cf).size();
+    const label loop_len = (cf).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         cmptAv(cf[i], f[i]);
     }
@@ -380,9 +380,9 @@ void cmptMag
     const FieldField<Field, Type>& f
 )
 {
-    const label loopLen = (cf).size();
+    const label loop_len = (cf).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         cmptMag(cf[i], f[i]);
     }
@@ -428,9 +428,9 @@ Type max(const FieldField<Field, Type>& f)
 {
     Type result = pTraits<Type>::min;
 
-    const label loopLen = (f).size();
+    const label loop_len = (f).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         if (f[i].size())
         {
@@ -450,9 +450,9 @@ Type min(const FieldField<Field, Type>& f)
 {
     Type result = pTraits<Type>::max;
 
-    const label loopLen = (f).size();
+    const label loop_len = (f).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         if (f[i].size())
         {
@@ -471,9 +471,9 @@ Type sum(const FieldField<Field, Type>& f)
 {
     Type result = Zero;
 
-    const label loopLen = (f).size();
+    const label loop_len = (f).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         result += sum(f[i]);
     }
@@ -490,9 +490,9 @@ typename typeOfMag<Type>::type sumMag(const FieldField<Field, Type>& f)
 
     resultType result = Zero;
 
-    const label loopLen = (f).size();
+    const label loop_len = (f).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         result += sumMag(f[i]);
     }
@@ -507,9 +507,9 @@ Type average(const FieldField<Field, Type>& f)
 {
     label n = 0;
 
-    const label loopLen = (f).size();
+    const label loop_len = (f).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         n += f[i].size();
     }
@@ -535,9 +535,9 @@ MinMax<Type> minMax(const FieldField<Field, Type>& f)
 {
     MinMax<Type> result;
 
-    const label loopLen = (f).size();
+    const label loop_len = (f).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         result += minMax(f[i]);
     }
@@ -552,9 +552,9 @@ scalarMinMax minMaxMag(const FieldField<Field, Type>& f)
 {
     scalarMinMax result;
 
-    const label loopLen = (f).size();
+    const label loop_len = (f).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         result += minMaxMag(f[i]);
     }
@@ -593,9 +593,9 @@ Type gAverage(const FieldField<Field, Type>& f)
 {
     label n = 0;
 
-    const label loopLen = (f).size();
+    const label loop_len = (f).size();
 
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         n += f[i].size();
     }
@@ -672,9 +672,9 @@ void OpFunc                                                                    \
     const FieldField<Field2, Type2>& f2                                        \
 )                                                                              \
 {                                                                              \
-    const label loopLen = (f).size();                                          \
+    const label loop_len = (f).size();                                         \
                                                                                \
-    for (label i = 0; i < loopLen; ++i)                                        \
+    for (label i = 0; i < loop_len; ++i)                                       \
     {                                                                          \
         OpFunc(f[i], f1[i], f2[i]);                                            \
     }                                                                          \
@@ -798,9 +798,9 @@ void OpFunc                                                                    \
     const VectorSpace<Form,Cmpt,nCmpt>& vs                                     \
 )                                                                              \
 {                                                                              \
-    const label loopLen = (result).size();                                     \
+    const label loop_len = (result).size();                                    \
                                                                                \
-    for (label i = 0; i < loopLen; ++i)                                        \
+    for (label i = 0; i < loop_len; ++i)                                       \
     {                                                                          \
         OpFunc(result[i], f1[i], vs);                                          \
     }                                                                          \
@@ -864,9 +864,9 @@ void OpFunc                                                                    \
     const FieldField<Field, Type>& f1                                          \
 )                                                                              \
 {                                                                              \
-    const label loopLen = (result).size();                                     \
+    const label loop_len = (result).size();                                    \
                                                                                \
-    for (label i = 0; i < loopLen; ++i)                                        \
+    for (label i = 0; i < loop_len; ++i)                                       \
     {                                                                          \
         OpFunc(result[i], vs, f1[i]);                                          \
     }                                                                          \
diff --git a/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctionsM.C b/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctionsM.C
index 5becaa47539331891d4e5981d306a6d3cdfda988..b3d2ee3e6a8e44a72e3ffe6deaa9f7b8bc74481e 100644
--- a/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctionsM.C
+++ b/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctionsM.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2022-2023 OpenCFD Ltd.
+    Copyright (C) 2022-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,9 +40,9 @@ void Func                                                                      \
     const FieldField<Field, Type1>& f1                                         \
 )                                                                              \
 {                                                                              \
-    const label loopLen = (result).size();                                     \
+    const label loop_len = (result).size();                                    \
                                                                                \
-    for (label i = 0; i < loopLen; ++i)                                        \
+    for (label i = 0; i < loop_len; ++i)                                       \
     {                                                                          \
         Func(result[i], f1[i]);                                                \
     }                                                                          \
@@ -83,9 +83,9 @@ void OpFunc                                                                    \
     const FieldField<Field, Type1>& f1                                         \
 )                                                                              \
 {                                                                              \
-    const label loopLen = (result).size();                                     \
+    const label loop_len = (result).size();                                    \
                                                                                \
-    for (label i = 0; i < loopLen; ++i)                                        \
+    for (label i = 0; i < loop_len; ++i)                                       \
     {                                                                          \
         OpFunc(result[i], f1[i]);                                              \
     }                                                                          \
@@ -127,9 +127,9 @@ void Func                                                                      \
     const FieldField<Field, Type2>& f2                                         \
 )                                                                              \
 {                                                                              \
-    const label loopLen = (result).size();                                     \
+    const label loop_len = (result).size();                                    \
                                                                                \
-    for (label i = 0; i < loopLen; ++i)                                        \
+    for (label i = 0; i < loop_len; ++i)                                       \
     {                                                                          \
         Func(result[i], f1[i], f2[i]);                                         \
     }                                                                          \
@@ -204,9 +204,9 @@ void Func                                                                      \
     const FieldField<Field, Type2>& f2                                         \
 )                                                                              \
 {                                                                              \
-    const label loopLen = (result).size();                                     \
+    const label loop_len = (result).size();                                    \
                                                                                \
-    for (label i = 0; i < loopLen; ++i)                                        \
+    for (label i = 0; i < loop_len; ++i)                                       \
     {                                                                          \
         Func(result[i], s1, f2[i]);                                            \
     }                                                                          \
@@ -248,9 +248,9 @@ void Func                                                                      \
     const Type2& s2                                                            \
 )                                                                              \
 {                                                                              \
-    const label loopLen = (result).size();                                     \
+    const label loop_len = (result).size();                                    \
                                                                                \
-    for (label i = 0; i < loopLen; ++i)                                        \
+    for (label i = 0; i < loop_len; ++i)                                       \
     {                                                                          \
         Func(result[i], f1[i], s2);                                            \
     }                                                                          \
@@ -299,9 +299,9 @@ void OpFunc                                                                    \
     const FieldField<Field, Type2>& f2                                         \
 )                                                                              \
 {                                                                              \
-    const label loopLen = (result).size();                                     \
+    const label loop_len = (result).size();                                    \
                                                                                \
-    for (label i = 0; i < loopLen; ++i)                                        \
+    for (label i = 0; i < loop_len; ++i)                                       \
     {                                                                          \
         OpFunc(result[i], f1[i], f2[i]);                                       \
     }                                                                          \
@@ -376,9 +376,9 @@ void OpFunc                                                                    \
     const FieldField<Field, Type2>& f2                                         \
 )                                                                              \
 {                                                                              \
-    const label loopLen = (result).size();                                     \
+    const label loop_len = (result).size();                                    \
                                                                                \
-    for (label i = 0; i < loopLen; ++i)                                        \
+    for (label i = 0; i < loop_len; ++i)                                       \
     {                                                                          \
         OpFunc(result[i], s1, f2[i]);                                          \
     }                                                                          \
@@ -420,9 +420,9 @@ void OpFunc                                                                    \
     const Type2& s2                                                            \
 )                                                                              \
 {                                                                              \
-    const label loopLen = (result).size();                                     \
+    const label loop_len = (result).size();                                    \
                                                                                \
-    for (label i = 0; i < loopLen; ++i)                                        \
+    for (label i = 0; i < loop_len; ++i)                                       \
     {                                                                          \
         OpFunc(result[i], f1[i], s2);                                          \
     }                                                                          \
@@ -472,9 +472,9 @@ void Func                                                                      \
     const FieldField<Field, Type3>& f3                                         \
 )                                                                              \
 {                                                                              \
-    const label loopLen = (result).size();                                     \
+    const label loop_len = (result).size();                                    \
                                                                                \
-    for (label i = 0; i < loopLen; ++i)                                        \
+    for (label i = 0; i < loop_len; ++i)                                       \
     {                                                                          \
         Func(result[i], f1[i], f2[i], f3[i]);                                  \
     }                                                                          \
@@ -625,9 +625,9 @@ void Func                                                                      \
     const Type3& s3                                                            \
 )                                                                              \
 {                                                                              \
-    const label loopLen = (result).size();                                     \
+    const label loop_len = (result).size();                                    \
                                                                                \
-    for (label i = 0; i < loopLen; ++i)                                        \
+    for (label i = 0; i < loop_len; ++i)                                       \
     {                                                                          \
         Func(result[i], f1[i], f2[i], s3);                                     \
     }                                                                          \
diff --git a/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldReuseFunctions.H b/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldReuseFunctions.H
index 7e51d03daf9d5b04c0d116bc8ac1a80729ea246e..c81e1069aa5456a35d9262ffc80b50b423638bdb 100644
--- a/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldReuseFunctions.H
+++ b/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldReuseFunctions.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,6 +40,15 @@ namespace Foam
 template<template<class> class Field, class TypeR, class Type1>
 struct reuseTmpFieldField
 {
+    //- Dissimilar types: non-reusable
+    static bool reusable
+    (
+        const tmp<FieldField<Field, Type1>>& tf1
+    )
+    {
+        return false;
+    }
+
     //- Pass-through to NewCalculatedType
     static tmp<FieldField<Field, TypeR>> New
     (
@@ -63,6 +72,15 @@ struct reuseTmpFieldField
 template<template<class> class Field, class TypeR>
 struct reuseTmpFieldField<Field, TypeR, TypeR>
 {
+    //- Identical types: possibly reusable
+    static bool reusable
+    (
+        const tmp<FieldField<Field, TypeR>>& tf1
+    )
+    {
+        return tf1.movable();
+    }
+
     //- Identical input and return types:
     //- allow optional copy assignment of the initial content
     static tmp<FieldField<Field, TypeR>> New
@@ -76,11 +94,13 @@ struct reuseTmpFieldField<Field, TypeR, TypeR>
             return tf1;
         }
 
-        auto tresult = FieldField<Field, TypeR>::NewCalculatedType(tf1());
+        const auto& f1 = tf1();
+
+        auto tresult = FieldField<Field, TypeR>::NewCalculatedType(f1);
 
         if (initCopy)
         {
-            tresult.ref() = tf1();
+            tresult.ref() = f1;
         }
 
         return tresult;
@@ -113,6 +133,16 @@ template
 >
 struct reuseTmpTmpFieldField
 {
+    //- Dissimilar types: non-reusable
+    static bool reusable
+    (
+        const tmp<FieldField<Field, Type1>>& tf1,
+        const tmp<FieldField<Field, Type2>>& tf2
+    )
+    {
+        return false;
+    }
+
     //- Dissimilar types: just use size
     static tmp<FieldField<Field, TypeR>> New
     (
@@ -128,6 +158,16 @@ struct reuseTmpTmpFieldField
 template<template<class> class Field, class TypeR, class Type1, class Type12>
 struct reuseTmpTmpFieldField<Field, TypeR, Type1, Type12, TypeR>
 {
+    //- Second input has return type: possibly reusable
+    static bool movable
+    (
+        const tmp<FieldField<Field, Type1>>& tf1,
+        const tmp<FieldField<Field, TypeR>>& tf2
+    )
+    {
+        return tf2.movable();
+    }
+
     //- Second input has return type
     static tmp<FieldField<Field, TypeR>> New
     (
@@ -148,6 +188,16 @@ struct reuseTmpTmpFieldField<Field, TypeR, Type1, Type12, TypeR>
 template<template<class> class Field, class TypeR, class Type2>
 struct reuseTmpTmpFieldField<Field, TypeR, TypeR, TypeR, Type2>
 {
+    //- First input has return type: possibly reusable
+    static bool reusable
+    (
+        const tmp<FieldField<Field, TypeR>>& tf1,
+        const tmp<FieldField<Field, Type2>>& tf2
+    )
+    {
+        return tf1.movable();
+    }
+
     //- First input has return type
     static tmp<FieldField<Field, TypeR>> New
     (
@@ -168,6 +218,16 @@ struct reuseTmpTmpFieldField<Field, TypeR, TypeR, TypeR, Type2>
 template<template<class> class Field, class TypeR>
 struct reuseTmpTmpFieldField<Field, TypeR, TypeR, TypeR, TypeR>
 {
+    //- Both inputs have return type: possibly reusable
+    static bool reusable
+    (
+        const tmp<FieldField<Field, TypeR>>& tf1,
+        const tmp<FieldField<Field, TypeR>>& tf2
+    )
+    {
+        return (tf1.movable() || tf2.movable());
+    }
+
     //- Both inputs have return type
     static tmp<FieldField<Field, TypeR>> New
     (
diff --git a/src/OpenFOAM/fields/Fields/Field/Field.C b/src/OpenFOAM/fields/Fields/Field/Field.C
index e116fa2cb45be36dd67aba81145da67b243a63e8..77d1fd238b52c7a6489c3c154745c06309460d3a 100644
--- a/src/OpenFOAM/fields/Fields/Field/Field.C
+++ b/src/OpenFOAM/fields/Fields/Field/Field.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2023 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -599,7 +599,8 @@ void Foam::Field<Type>::rmap
 template<class Type>
 void Foam::Field<Type>::negate()
 {
-    TFOR_ALL_F_OP_OP_F(Type, *this, =, -, Type, *this)
+    // std::for_each
+    TSEQ_FORALL_F_OP_OP_F_inplace(*this, =, -, *this)
 }
 
 
@@ -629,8 +630,16 @@ void Foam::Field<Type>::replace
     const UList<cmptType>& sf
 )
 {
-    TFOR_ALL_F_OP_FUNC_S_F(Type, *this, ., replace, const direction, d,
-        cmptType, sf)
+    if (this->cdata_bytes() == sf.cdata_bytes())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_S_F_inplace(*this, ., replace, d, sf)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_S_F(*this, ., replace, d, sf)
+    }
 }
 
 
@@ -653,8 +662,7 @@ void Foam::Field<Type>::replace
     const cmptType& c
 )
 {
-    TFOR_ALL_F_OP_FUNC_S_S(Type, *this, ., replace, const direction, d,
-        cmptType, c)
+    TSEQ_FORALL_F_OP_FUNC_S_S(*this, ., replace, d, c)
 }
 
 
@@ -778,7 +786,7 @@ template<class Type>
 template<class Form, class Cmpt, Foam::direction nCmpt>
 void Foam::Field<Type>::operator=(const VectorSpace<Form,Cmpt,nCmpt>& vs)
 {
-    TFOR_ALL_F_OP_S(Type, *this, =, VSType, vs)
+    TSEQ_FORALL_F_OP_S(*this, =, vs)
 }
 
 
@@ -787,7 +795,16 @@ void Foam::Field<Type>::operator=(const VectorSpace<Form,Cmpt,nCmpt>& vs)
 template<class Type>                                                           \
 void Foam::Field<Type>::operator op(const UList<TYPE>& f)                      \
 {                                                                              \
-    TFOR_ALL_F_OP_F(Type, *this, op, TYPE, f)                                  \
+    if (this->cdata_bytes() == f.cdata_bytes())                                \
+    {                                                                          \
+        /* std::for_each */                                                    \
+        TSEQ_FORALL_F_OP_F_inplace(*this, op, f)                               \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        /* std::transform */                                                   \
+        TSEQ_FORALL_F_OP_F(*this, op, f)                                       \
+    }                                                                          \
 }                                                                              \
                                                                                \
 template<class Type>                                                           \
@@ -800,7 +817,7 @@ void Foam::Field<Type>::operator op(const tmp<Field<TYPE>>& tf)                \
 template<class Type>                                                           \
 void Foam::Field<Type>::operator op(const TYPE& t)                             \
 {                                                                              \
-    TFOR_ALL_F_OP_S(Type, *this, op, TYPE, t)                                  \
+    TSEQ_FORALL_F_OP_S(*this, op, t)                                           \
 }
 
 COMPUTED_ASSIGNMENT(Type, +=)
diff --git a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
index 583c4439710246ebc4abfc18b44a55a9cd8d08fc..21fb8b9cb67fd5db8d67931af0976f18d7b6ac32 100644
--- a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
+++ b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,19 +47,32 @@ void component
     const direction d
 )
 {
-    typedef typename Field<Type>::cmptType resultType;
-
-    TFOR_ALL_F_OP_F_FUNC_S
-    (
-        resultType, result, =, Type, f1, .component, const direction, d
-    )
+    if (result.cdata_bytes() == f1.cdata_bytes())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_F_FUNC_S_inplace(result, =, f1, .component, d)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_F_FUNC_S(result, =, f1, .component, d)
+    }
 }
 
 
 template<class Type>
 void T(Field<Type>& result, const UList<Type>& f1)
 {
-    TFOR_ALL_F_OP_F_FUNC(Type, result, =, Type, f1, T)
+    if (result.cdata() == f1.cdata())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_F_FUNC_inplace(result, =, f1, T)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_F_FUNC(result, =, f1, T)
+    }
 }
 
 
@@ -72,11 +85,19 @@ void pow
 {
     typedef typename powProduct<Type, r>::type resultType;
 
-    TFOR_ALL_F_OP_FUNC_F_S
-    (
-        resultType, result, =, pow, Type, f1, resultType,
-        pTraits<resultType>::zero
-    )
+    if (result.cdata_bytes() == f1.cdata_bytes())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_F_S_inplace
+        (result, =, pow, f1, pTraits<resultType>::zero)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_F_S
+        (result, =, pow, f1, pTraits<resultType>::zero)
+    }
+
 }
 
 template<class Type, direction r>
@@ -116,9 +137,16 @@ void sqr
     const UList<Type>& f1
 )
 {
-    typedef typename outerProduct<Type, Type>::type resultType;
-
-    TFOR_ALL_F_OP_FUNC_F(resultType, result, =, sqr, Type, f1)
+    if (result.cdata_bytes() == f1.cdata_bytes())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_F_inplace(result, =, sqr, f1)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_F(result, =, sqr, f1)
+    }
 }
 
 template<class Type>
@@ -150,9 +178,16 @@ void magSqr
     const UList<Type>& f1
 )
 {
-    typedef typename typeOfMag<Type>::type resultType;
-
-    TFOR_ALL_F_OP_FUNC_F(resultType, result, =, magSqr, Type, f1)
+    if (result.cdata_bytes() == f1.cdata_bytes())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_F_inplace(result, =, magSqr, f1)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_F(result, =, magSqr, f1)
+    }
 }
 
 template<class Type>
@@ -186,9 +221,16 @@ void mag
     const UList<Type>& f1
 )
 {
-    typedef typename typeOfMag<Type>::type resultType;
-
-    TFOR_ALL_F_OP_FUNC_F(resultType, result, =, mag, Type, f1)
+    if (result.cdata_bytes() == f1.cdata_bytes())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_F_inplace(result, =, mag, f1)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_F(result, =, mag, f1)
+    }
 }
 
 template<class Type>
@@ -222,9 +264,16 @@ void cmptMax
     const UList<Type>& f1
 )
 {
-    typedef typename Field<Type>::cmptType resultType;
-
-    TFOR_ALL_F_OP_FUNC_F(resultType, result, =, cmptMax, Type, f1)
+    if (result.cdata_bytes() == f1.cdata_bytes())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_F_inplace(result, =, cmptMax, f1)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_F(result, =, cmptMax, f1)
+    }
 }
 
 template<class Type>
@@ -254,9 +303,16 @@ void cmptMin
     const UList<Type>& f1
 )
 {
-    typedef typename Field<Type>::cmptType resultType;
-
-    TFOR_ALL_F_OP_FUNC_F(resultType, result, =, cmptMin, Type, f1)
+    if (result.cdata_bytes() == f1.cdata_bytes())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_F_inplace(result, =, cmptMin, f1)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_F(result, =, cmptMin, f1)
+    }
 }
 
 template<class Type>
@@ -286,9 +342,16 @@ void cmptAv
     const UList<Type>& f1
 )
 {
-    typedef typename Field<Type>::cmptType resultType;
-
-    TFOR_ALL_F_OP_FUNC_F(resultType, result, =, cmptAv, Type, f1)
+    if (result.cdata_bytes() == f1.cdata_bytes())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_F_inplace(result, =, cmptAv, f1)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_F(result, =, cmptAv, f1)
+    }
 }
 
 template<class Type>
@@ -314,7 +377,16 @@ tmp<Field<typename Field<Type>::cmptType>> cmptAv(const tmp<Field<Type>>& tf1)
 template<class Type>
 void cmptMag(Field<Type>& result, const UList<Type>& f1)
 {
-    TFOR_ALL_F_OP_FUNC_F(Type, result, =, cmptMag, Type, f1)
+    if (result.cdata_bytes() == f1.cdata_bytes())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_F_inplace(result, =, cmptMag, f1)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_F(result, =, cmptMag, f1)
+    }
 }
 
 template<class Type>
@@ -338,7 +410,16 @@ tmp<Field<Type>> cmptMag(const tmp<Field<Type>>& tf1)
 template<class Type>
 void cmptMagSqr(Field<Type>& result, const UList<Type>& f1)
 {
-    TFOR_ALL_F_OP_FUNC_F(Type, result, =, cmptMagSqr, Type, f1)
+    if (result.cdata_bytes() == f1.cdata_bytes())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_F_inplace(result, =, cmptMagSqr, f1)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_F(result, =, cmptMagSqr, f1)
+    }
 }
 
 template<class Type>
@@ -375,7 +456,7 @@ Type max(const UList<Type>& f1)
     if (f1.size())
     {
         Type result(f1[0]);
-        TFOR_ALL_S_OP_FUNC_F_S(Type, result, =, max, Type, f1, Type, result)
+        TSEQ_FORALL_S_OP_FUNC_F_S(result, =, max, f1, result)
         return result;
     }
 
@@ -390,7 +471,7 @@ Type min(const UList<Type>& f1)
     if (f1.size())
     {
         Type result(f1[0]);
-        TFOR_ALL_S_OP_FUNC_F_S(Type, result, =, min, Type, f1, Type, result)
+        TSEQ_FORALL_S_OP_FUNC_F_S(result, =, min, f1, result)
         return result;
     }
 
@@ -409,7 +490,7 @@ Type sum(const UList<Type>& f1)
     if (f1.size())
     {
         // Use resultType() as functional cast
-        TFOR_ALL_S_OP_FUNC_F(resultType, result, +=, resultType, Type, f1)
+        TSEQ_FORALL_S_OP_FUNC_F(result, +=, resultType, f1)
     }
 
     return Type(result);
@@ -432,15 +513,12 @@ Type maxMagSqr(const UList<Type>& f1)
     if (f1.size())
     {
         Type result(f1[0]);
-        TFOR_ALL_S_OP_FUNC_F_S
+        TSEQ_FORALL_S_OP_FUNC_F_S
         (
-            Type,
             result,
             =,
             maxMagSqrOp<Type>(),
-            Type,
             f1,
-            Type,
             result
         )
         return result;
@@ -457,15 +535,12 @@ Type minMagSqr(const UList<Type>& f1)
     if (f1.size())
     {
         Type result(f1[0]);
-        TFOR_ALL_S_OP_FUNC_F_S
+        TSEQ_FORALL_S_OP_FUNC_F_S
         (
-            Type,
             result,
             =,
             minMagSqrOp<Type>(),
-            Type,
             f1,
-            Type,
             result
         )
         return result;
@@ -485,7 +560,8 @@ sumProd(const UList<Type>& f1, const UList<Type>& f2)
     resultType result = Zero;
     if (f1.size() && (f1.size() == f2.size()))
     {
-        TFOR_ALL_S_OP_F_OP_F(resultType, result, +=, Type, f1, &&, Type, f2)
+        // std::transform
+        TSEQ_FORALL_S_OP_F_OP_F(result, +=, f1, &&, f2)
     }
     return result;
 }
@@ -497,15 +573,12 @@ Type sumCmptProd(const UList<Type>& f1, const UList<Type>& f2)
     Type result = Zero;
     if (f1.size() && (f1.size() == f2.size()))
     {
-        TFOR_ALL_S_OP_FUNC_F_F
+        TSEQ_FORALL_S_OP_FUNC_F_F
         (
-            Type,
             result,
             +=,
             cmptMultiply,
-            Type,
             f1,
-            Type,
             f2
         )
     }
@@ -522,7 +595,7 @@ sumSqr(const UList<Type>& f1)
     resultType result = Zero;
     if (f1.size())
     {
-        TFOR_ALL_S_OP_FUNC_F(resultType, result, +=, sqr, Type, f1)
+        TSEQ_FORALL_S_OP_FUNC_F(result, +=, sqr, f1)
     }
     return result;
 }
@@ -547,7 +620,7 @@ sumMag(const UList<Type>& f1)
     resultType result = Zero;
     if (f1.size())
     {
-        TFOR_ALL_S_OP_FUNC_F(resultType, result, +=, mag, Type, f1)
+        TSEQ_FORALL_S_OP_FUNC_F(result, +=, mag, f1)
     }
     return result;
 }
@@ -561,7 +634,7 @@ Type sumCmptMag(const UList<Type>& f1)
     Type result = Zero;
     if (f1.size())
     {
-        TFOR_ALL_S_OP_FUNC_F(Type, result, +=, cmptMag, Type, f1)
+        TSEQ_FORALL_S_OP_FUNC_F(result, +=, cmptMag, f1)
     }
     return result;
 }
@@ -773,8 +846,20 @@ void OpFunc                                                                    \
     const UList<Type2>& f2                                                     \
 )                                                                              \
 {                                                                              \
-    typedef typename product<Type1, Type2>::type resultType;                   \
-    TFOR_ALL_F_OP_F_OP_F(resultType, result, =, Type1, f1, Op, Type2, f2)      \
+    if                                                                         \
+    (                                                                          \
+        result.cdata_bytes() == f1.cdata_bytes()                               \
+     || result.cdata_bytes() == f2.cdata_bytes()                               \
+    )                                                                          \
+    {                                                                          \
+        /* std::for_each */                                                    \
+        TSEQ_FORALL_F_OP_F_OP_F_inplace(result, =, f1, Op, f2)                 \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        /* std::transform */                                                   \
+        TSEQ_FORALL_F_OP_F_OP_F(result, =, f1, Op, f2)                         \
+    }                                                                          \
 }                                                                              \
                                                                                \
 template<class Type1, class Type2>                                             \
@@ -829,9 +914,21 @@ void OpFunc                                                                    \
     const VectorSpace<Form,Cmpt,nCmpt>& vs                                     \
 )                                                                              \
 {                                                                              \
-    typedef typename product<Type, Form>::type resultType;                     \
-    TFOR_ALL_F_OP_F_OP_S                                                       \
-        (resultType, result, =,Type, f1, Op, Form, static_cast<const Form&>(vs))\
+    if                                                                         \
+    (                                                                          \
+        result.cdata_bytes() == f1.cdata_bytes()                               \
+    )                                                                          \
+    {                                                                          \
+        /* std::for_each */                                                    \
+        TSEQ_FORALL_F_OP_F_OP_S_inplace                                        \
+        (result, =, f1, Op, static_cast<const Form&>(vs))                      \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        /* std::transform */                                                   \
+        TSEQ_FORALL_F_OP_F_OP_S                                                \
+        (result, =, f1, Op, static_cast<const Form&>(vs))                      \
+    }                                                                          \
 }                                                                              \
                                                                                \
 template<class Type, class Form, class Cmpt, direction nCmpt>                  \
@@ -867,9 +964,21 @@ void OpFunc                                                                    \
     const UList<Type>& f1                                                      \
 )                                                                              \
 {                                                                              \
-    typedef typename product<Form, Type>::type resultType;                     \
-    TFOR_ALL_F_OP_S_OP_F                                                       \
-        (resultType, result, =,Form,static_cast<const Form&>(vs), Op, Type, f1)\
+    if                                                                         \
+    (                                                                          \
+        result.cdata_bytes() == f1.cdata_bytes()                               \
+    )                                                                          \
+    {                                                                          \
+        /* std::for_each */                                                    \
+        TSEQ_FORALL_F_OP_S_OP_F_inplace                                        \
+        (result, =, static_cast<const Form&>(vs), Op, f1)                      \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        /* std::transform */                                                   \
+        TSEQ_FORALL_F_OP_S_OP_F                                                \
+        (result, =, static_cast<const Form&>(vs), Op, f1)                      \
+    }                                                                          \
 }                                                                              \
                                                                                \
 template<class Form, class Cmpt, direction nCmpt, class Type>                  \
diff --git a/src/OpenFOAM/fields/Fields/Field/FieldFunctionsM.C b/src/OpenFOAM/fields/Fields/Field/FieldFunctionsM.C
index 39829843ab153be3f4ae541d6c92df8ae1cce532..47838432238f8224a060716758c449e34f4254ce 100644
--- a/src/OpenFOAM/fields/Fields/Field/FieldFunctionsM.C
+++ b/src/OpenFOAM/fields/Fields/Field/FieldFunctionsM.C
@@ -40,7 +40,14 @@ void Func                                                                      \
     const UList<Type1>& f1                                                     \
 )                                                                              \
 {                                                                              \
-    TFOR_ALL_F_OP_FUNC_F(ReturnType, result, =, ::Foam::Func, Type1, f1)       \
+    if (result.cdata_bytes() == f1.cdata_bytes())                              \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_FUNC_F_inplace(result, =, ::Foam::Func, f1)           \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_FUNC_F(result, =, ::Foam::Func, f1)                   \
+    }                                                                          \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -78,7 +85,14 @@ void OpFunc                                                                    \
     const UList<Type1>& f1                                                     \
 )                                                                              \
 {                                                                              \
-    TFOR_ALL_F_OP_OP_F(ReturnType, result, =, Op, Type1, f1)                   \
+    if (result.cdata_bytes() == f1.cdata_bytes())                              \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_OP_F_inplace(result, =, Op, f1)                       \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_OP_F(result, =, Op, f1)                               \
+    }                                                                          \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -117,10 +131,18 @@ void Func                                                                      \
     const UList<Type2>& f2                                                     \
 )                                                                              \
 {                                                                              \
-    TFOR_ALL_F_OP_FUNC_F_F                                                     \
+    if                                                                         \
     (                                                                          \
-        ReturnType, result, =, ::Foam::Func, Type1, f1, Type2, f2              \
+        result.cdata_bytes() == f1.cdata_bytes()                               \
+     || result.cdata_bytes() == f2.cdata_bytes()                               \
     )                                                                          \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_FUNC_F_F_inplace(result, =, ::Foam::Func, f1, f2)     \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_FUNC_F_F(result, =, ::Foam::Func, f1, f2)             \
+    }                                                                          \
 }
 
 #define BINARY_FUNCTION_INTERFACE(ReturnType, Type1, Type2, Func)              \
@@ -194,10 +216,14 @@ void Func                                                                      \
     const UList<Type2>& f2                                                     \
 )                                                                              \
 {                                                                              \
-    TFOR_ALL_F_OP_FUNC_S_F                                                     \
-    (                                                                          \
-        ReturnType, result, =, ::Foam::Func, Type1, s1, Type2, f2              \
-    )                                                                          \
+    if (result.cdata_bytes() == f2.cdata_bytes())                              \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_FUNC_S_F_inplace(result, =, ::Foam::Func, s1, f2)     \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_FUNC_S_F(result, =, ::Foam::Func, s1, f2)             \
+    }                                                                          \
 }
 
 #define BINARY_FUNCTION_INTERFACE_SF(ReturnType, Type1, Type2, Func)           \
@@ -242,10 +268,14 @@ void Func                                                                      \
     const Type2& s2                                                            \
 )                                                                              \
 {                                                                              \
-    TFOR_ALL_F_OP_FUNC_F_S                                                     \
-    (                                                                          \
-        ReturnType, result, =, ::Foam::Func, Type1, f1, Type2, s2              \
-    )                                                                          \
+    if (result.cdata_bytes() == f1.cdata_bytes())                              \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_FUNC_F_S_inplace(result, =, ::Foam::Func, f1, s2)     \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_FUNC_F_S(result, =, ::Foam::Func, f1, s2)             \
+    }                                                                          \
 }
 
 #define BINARY_FUNCTION_INTERFACE_FS(ReturnType, Type1, Type2, Func)           \
@@ -297,7 +327,18 @@ void OpFunc                                                                    \
     const UList<Type2>& f2                                                     \
 )                                                                              \
 {                                                                              \
-    TFOR_ALL_F_OP_F_OP_F(ReturnType, result, =, Type1, f1, Op, Type2, f2)      \
+    if                                                                         \
+    (                                                                          \
+        result.cdata_bytes() == f1.cdata_bytes()                               \
+     || result.cdata_bytes() == f2.cdata_bytes()                               \
+    )                                                                          \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_F_OP_F_inplace(result, =, f1, Op, f2)                 \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_F_OP_F(result, =, f1, Op, f2)                         \
+    }                                                                          \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -365,7 +406,14 @@ void OpFunc                                                                    \
     const UList<Type2>& f2                                                     \
 )                                                                              \
 {                                                                              \
-    TFOR_ALL_F_OP_S_OP_F(ReturnType, result, =, Type1, s1, Op, Type2, f2)      \
+    if (result.cdata_bytes() == f2.cdata_bytes())                              \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_S_OP_F_inplace(result, =, s1, Op, f2)                 \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_S_OP_F(result, =, s1, Op, f2)                         \
+    }                                                                          \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -404,7 +452,14 @@ void OpFunc                                                                    \
     const Type2& s2                                                            \
 )                                                                              \
 {                                                                              \
-    TFOR_ALL_F_OP_F_OP_S(ReturnType, result, =, Type1, f1, Op, Type2, s2)      \
+    if (result.cdata_bytes() == f1.cdata_bytes())                              \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_F_OP_S_inplace(result, =, f1, Op, s2)                 \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_F_OP_S(result, =, f1, Op, s2)                         \
+    }                                                                          \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -451,10 +506,21 @@ void Func                                                                      \
     const UList<Type3>& f3                                                     \
 )                                                                              \
 {                                                                              \
-    TFOR_ALL_F_OP_FUNC_F_F_F                                                   \
+    if                                                                         \
     (                                                                          \
-        ReturnType, result, =, ::Foam::Func, Type1, f1, Type2, f2, Type3, f3   \
+        result.cdata_bytes() == f1.cdata_bytes()                               \
+     || result.cdata_bytes() == f2.cdata_bytes()                               \
+     || result.cdata_bytes() == f3.cdata_bytes()                               \
     )                                                                          \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_FUNC_F_F_F_inplace                                    \
+        (result, =, ::Foam::Func, f1, f2, f3)                                  \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_FUNC_F_F_F                                            \
+        (result, =, ::Foam::Func, f1, f2, f3)                                  \
+    }                                                                          \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -586,10 +652,20 @@ void Func                                                                      \
     const Type3& s3                                                            \
 )                                                                              \
 {                                                                              \
-    TFOR_ALL_F_OP_FUNC_F_F_S                                                   \
+    if                                                                         \
     (                                                                          \
-        ReturnType, result, =, ::Foam::Func, Type1, f1, Type2, f2, Type3, s3   \
+        result.cdata_bytes() == f1.cdata_bytes()                               \
+     || result.cdata_bytes() == f2.cdata_bytes()                               \
     )                                                                          \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_FUNC_F_F_S_inplace                                    \
+        (result, =, ::Foam::Func, f1, f2, s3)                                  \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        TSEQ_FORALL_F_OP_FUNC_F_F_S                                            \
+        (result, =, ::Foam::Func, f1, f2, s3)                                  \
+    }                                                                          \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
diff --git a/src/OpenFOAM/fields/Fields/Field/FieldM.H b/src/OpenFOAM/fields/Fields/Field/FieldM.H
index 06e24b1574a3bd94c76fc8742efd869909ec6768..1cc2efe01b57c7cb804bd8e05a0369a8961f7fa4 100644
--- a/src/OpenFOAM/fields/Fields/Field/FieldM.H
+++ b/src/OpenFOAM/fields/Fields/Field/FieldM.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2022-2023 OpenCFD Ltd.
+    Copyright (C) 2022-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,8 +32,8 @@ Description
 #ifndef Foam_FieldM_H
 #define Foam_FieldM_H
 
-#include "error.H"
-#include "ListLoopM.H"  // For list access macros
+#include "errorCheckFields.H"   // Field size checks (fulldebug mode)
+#include "ListLoopM.H"          // List access macros
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -42,126 +42,20 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef FULLDEBUG
-
-template<class Type1, class Type2>
-void checkFields
-(
-    const UList<Type1>& f1,
-    const UList<Type2>& f2,
-    const char* op
-)
-{
-    if (f1.size() != f2.size())
-    {
-        FatalErrorInFunction
-            << " Field<"<<pTraits<Type1>::typeName<<"> f1("<<f1.size()<<')'
-            << " and Field<"<<pTraits<Type2>::typeName<<"> f2("<<f2.size()<<')'
-            << endl
-            << " for operation " << op
-            << abort(FatalError);
-    }
-}
-
-template<class Type1, class Type2, class Type3>
-void checkFields
-(
-    const UList<Type1>& f1,
-    const UList<Type2>& f2,
-    const UList<Type3>& f3,
-    const char* op
-)
-{
-    if (f1.size() != f2.size() || f1.size() != f3.size())
-    {
-        FatalErrorInFunction
-            << " Field<"<<pTraits<Type1>::typeName<<"> f1("<<f1.size()<<')'
-            << ", Field<"<<pTraits<Type2>::typeName<<"> f2("<<f2.size()<<')'
-            << " and Field<"<<pTraits<Type3>::typeName<<"> f3("<<f3.size()<<')'
-            << endl
-            << "    for operation " << op
-            << abort(FatalError);
-    }
-}
-
-template<class Type1, class Type2, class Type3, class Type4>
-void checkFields
-(
-    const UList<Type1>& f1,
-    const UList<Type2>& f2,
-    const UList<Type3>& f3,
-    const UList<Type4>& f4,
-    const char* op
-)
-{
-    if
-    (
-        f1.size() != f2.size()
-     || f1.size() != f3.size()
-     || f1.size() != f4.size()
-    )
-    {
-        FatalErrorInFunction
-            << " Field<"<<pTraits<Type1>::typeName<<"> f1("<<f1.size()<<')'
-            << ", Field<"<<pTraits<Type2>::typeName<<"> f2("<<f2.size()<<')'
-            << ", Field<"<<pTraits<Type3>::typeName<<"> f3("<<f3.size()<<')'
-            << " and Field<"<<pTraits<Type4>::typeName<<"> f4("<<f4.size()<<')'
-            << endl
-            << "    for operation " << op
-            << abort(FatalError);
-    }
-}
-
-#else
-
-template<class Type1, class Type2>
-void checkFields
-(
-    const UList<Type1>&,
-    const UList<Type2>&,
-    const char*
-)
-{}
-
-template<class Type1, class Type2, class Type3>
-void checkFields
-(
-    const UList<Type1>&,
-    const UList<Type2>&,
-    const UList<Type3>&,
-    const char*
-)
-{}
-
-template<class Type1, class Type2, class Type3, class Type4>
-void checkFields
-(
-    const UList<Type1>&,
-    const UList<Type2>&,
-    const UList<Type3>&,
-    const UList<Type4>&,
-    const char*
-)
-{}
-
-#endif
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
+// ------------
 // Unary Free Function : f1 OP Func(f2)
 
-#define TFOR_ALL_F_OP_FUNC_F(typeF1, f1, OP, FUNC, typeF2, f2)                 \
+#define TSEQ_FORALL_F_OP_FUNC_F_impl(f1, OP, FUNC, f2, Restrict)               \
 {                                                                              \
     /* Check fields have same size */                                          \
     checkFields(f1, f2, "f1 " #OP " " #FUNC "(f2)");                           \
                                                                                \
     /* Field access */                                                         \
-    List_ACCESS(typeF1, f1, f1P);                                              \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
+    TFOR_List_ACCESS(f1, f1P, Restrict);                                       \
+    TFOR_List_CONST_ACCESS(f2, f2P, Restrict);                                 \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: f1 OP FUNC(f2) */                                                 \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -170,20 +64,32 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_FUNC_F(f1, OP, FUNC, f2) \
+        TSEQ_FORALL_F_OP_FUNC_F_impl(f1, OP, FUNC, f2, __restrict__)
+
+#define TSEQ_FORALL_F_OP_FUNC_F_inplace(f1, OP, FUNC, f2) \
+        TSEQ_FORALL_F_OP_FUNC_F_impl(f1, OP, FUNC, f2,)
+
+// old-style
+#define TFOR_ALL_F_OP_FUNC_F(typeF1, f1, OP, FUNC, typeF2, f2) \
+     TSEQ_FORALL_F_OP_FUNC_F_impl(f1, OP, FUNC, f2, __restrict__)
 
+
+// ------------
 // Nullary Member Function : f1 OP f2.FUNC()
+// NB: the '.' is not passed in the func parameter
 
-#define TFOR_ALL_F_OP_F_FUNC(typeF1, f1, OP, typeF2, f2, FUNC)                 \
+#define TSEQ_FORALL_F_OP_F_FUNC_impl(f1, OP, f2, FUNC, Restrict)               \
 {                                                                              \
     /* Check fields have same size */                                          \
-    checkFields(f1, f2, "f1 " #OP " f2" #FUNC);                                \
+    checkFields(f1, f2, "f1 " #OP " f2 ." #FUNC "()");                         \
                                                                                \
     /* Field access */                                                         \
-    List_ACCESS(typeF1, f1, f1P);                                              \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
+    TFOR_List_ACCESS(f1, f1P, Restrict);                                       \
+    TFOR_List_CONST_ACCESS(f2, f2P, Restrict);                                 \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: f1 OP f2.FUNC() */                                                \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -192,21 +98,32 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_F_FUNC(f1, OP, f2, FUNC) \
+        TSEQ_FORALL_F_OP_F_FUNC_impl(f1, OP, f2, FUNC, __restrict__)
+
+#define TSEQ_FORALL_F_OP_F_FUNC_inplace(f1, OP, f2, FUNC) \
+        TSEQ_FORALL_F_OP_F_FUNC_impl(f1, OP, f2, FUNC,)
 
+// old-style
+#define TFOR_ALL_F_OP_F_FUNC(typeF1, f1, OP, typeF2, f2, FUNC) \
+     TSEQ_FORALL_F_OP_F_FUNC_impl(f1, OP, f2, FUNC, __restrict__)
+
+
+// ------------
 // Binary Free Function : f1 OP FUNC(f2, f3)
 
-#define TFOR_ALL_F_OP_FUNC_F_F(typeF1, f1, OP, FUNC, typeF2, f2, typeF3, f3)   \
+#define TSEQ_FORALL_F_OP_FUNC_F_F_impl(f1, OP, FUNC, f2, f3, Restrict)         \
 {                                                                              \
     /* Check fields have same size */                                          \
     checkFields(f1, f2, f3, "f1 " #OP " " #FUNC "(f2, f3)");                   \
                                                                                \
     /* Field access */                                                         \
-    List_ACCESS(typeF1, f1, f1P);                                              \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
-    List_CONST_ACCESS(typeF3, f3, f3P);                                        \
+    TFOR_List_ACCESS(f1, f1P, Restrict);                                       \
+    TFOR_List_CONST_ACCESS(f2, f2P, Restrict);                                 \
+    TFOR_List_CONST_ACCESS(f3, f3P, Restrict);                                 \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: f1 OP FUNC(f2, f3) */                                             \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -215,20 +132,31 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_FUNC_F_F(f1, OP, FUNC, f2, f3) \
+        TSEQ_FORALL_F_OP_FUNC_F_F_impl(f1, OP, FUNC, f2, f3, __restrict__)
+
+#define TSEQ_FORALL_F_OP_FUNC_F_F_inplace(f1, OP, FUNC, f2, f3) \
+        TSEQ_FORALL_F_OP_FUNC_F_F_impl(f1, OP, FUNC, f2, f3,)
+
+// old-style
+#define TFOR_ALL_F_OP_FUNC_F_F(typeF1, f1, OP, FUNC, typeF2, f2, typeF3, f3) \
+     TSEQ_FORALL_F_OP_FUNC_F_F_impl(f1, OP, FUNC, f2, f3, __restrict__)
+
 
+// ------------
 // [reduction] Binary Free Function : s OP FUNC(f1, f2)
 
-#define TFOR_ALL_S_OP_FUNC_F_F(typeS, s, OP, FUNC, typeF1, f1, typeF2, f2)     \
+#define TSEQ_FORALL_S_OP_FUNC_F_F_impl(s, OP, FUNC, f1, f2, Unused)            \
 {                                                                              \
     /* Check fields have same size */                                          \
     checkFields(f1, f2, "s " #OP " " #FUNC "(f1, f2)");                        \
                                                                                \
     /* Field access */                                                         \
-    List_CONST_ACCESS(typeF1, f1, f1P);                                        \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
+    TFOR_List_CONST_ACCESS(f1, f1P, __restrict__);                             \
+    TFOR_List_CONST_ACCESS(f2, f2P, __restrict__);                             \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: s OP FUNC(f1, f2) */                                              \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas, reduction... */                                                \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -237,20 +165,31 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_S_OP_FUNC_F_F(s, OP, FUNC, f1, f2) \
+        TSEQ_FORALL_S_OP_FUNC_F_F_impl(s, OP, FUNC, f1, f2, __restrict__)
 
+#define TSEQ_FORALL_S_OP_FUNC_F_F_inplace(s, OP, FUNC, f1, f2) \
+        TSEQ_FORALL_S_OP_FUNC_F_F_impl(s, OP, FUNC, f1, f2,)
+
+// old-style
+#define TFOR_ALL_S_OP_FUNC_F_F(typeS, s, OP, FUNC, typeF1, f1, typeF2, f2) \
+     TSEQ_FORALL_S_OP_FUNC_F_F_impl(s, OP, FUNC, f1, f2, __restrict__)
+
+
+// ------------
 // Binary Free Function : f1 OP FUNC(f2, s)
 
-#define TFOR_ALL_F_OP_FUNC_F_S(typeF1, f1, OP, FUNC, typeF2, f2, typeS, s)     \
+#define TSEQ_FORALL_F_OP_FUNC_F_S_impl(f1, OP, FUNC, f2, s, Restrict)          \
 {                                                                              \
     /* Check fields have same size */                                          \
     checkFields(f1, f2, "f1 " #OP " " #FUNC "(f2, s)");                        \
                                                                                \
     /* Field access */                                                         \
-    List_ACCESS(typeF1, f1, f1P);                                              \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
+    TFOR_List_ACCESS(f1, f1P, Restrict);                                       \
+    TFOR_List_CONST_ACCESS(f2, f2P, Restrict);                                 \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: f1 OP FUNC(f2, s) */                                              \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -259,16 +198,27 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_FUNC_F_S(f1, OP, FUNC, f2, s) \
+        TSEQ_FORALL_F_OP_FUNC_F_S_impl(f1, OP, FUNC, f2, s, __restrict__)
+
+#define TSEQ_FORALL_F_OP_FUNC_F_S_inplace(f1, OP, FUNC, f2, s) \
+        TSEQ_FORALL_F_OP_FUNC_F_S_impl(f1, OP, FUNC, f2, s,)
+
+// old-style
+#define TFOR_ALL_F_OP_FUNC_F_S(typeF1, f1, OP, FUNC, typeF2, f2, typeS, s) \
+     TSEQ_FORALL_F_OP_FUNC_F_S_impl(f1, OP, FUNC, f2, s, __restrict__)
 
+
+// ------------
 // [reduction] Binary Free Function : s1 OP FUNC(f, s2)
 
-#define TFOR_ALL_S_OP_FUNC_F_S(typeS1, s1, OP, FUNC, typeF, f, typeS2, s2)     \
+#define TSEQ_FORALL_S_OP_FUNC_F_S_impl(s1, OP, FUNC, f, s2, Unused)            \
 {                                                                              \
     /* Field access */                                                         \
-    List_CONST_ACCESS(typeF, f, fP);                                           \
+    TFOR_List_CONST_ACCESS(f, fP, __restrict__);                               \
+    TFOR_List_LENGTH(f, loop_len);                                             \
                                                                                \
     /* Loop: s1 OP FUNC(f, s2) */                                              \
-    const label loop_len = (f).size();                                         \
                                                                                \
     /* pragmas, reduction... */                                                \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -277,20 +227,31 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_S_OP_FUNC_F_S(s1, OP, FUNC, f, s2) \
+        TSEQ_FORALL_S_OP_FUNC_F_S_impl(s1, OP, FUNC, f, s2, __restrict__)
+
+#define TSEQ_FORALL_S_OP_FUNC_F_S_inplace(s1, OP, FUNC, f, s2) \
+        TSEQ_FORALL_S_OP_FUNC_F_S_impl(s1, OP, FUNC, f, s2,)
+
+// old-style
+#define TFOR_ALL_S_OP_FUNC_F_S(typeS1, s1, OP, FUNC, typeF, f, typeS2, s2) \
+     TSEQ_FORALL_S_OP_FUNC_F_S_impl(s1, OP, FUNC, f, s2, __restrict__)
 
+
+// ------------
 // Binary Free Function : f1 OP FUNC(s, f2)
 
-#define TFOR_ALL_F_OP_FUNC_S_F(typeF1, f1, OP, FUNC, typeS, s, typeF2, f2)     \
+#define TSEQ_FORALL_F_OP_FUNC_S_F_impl(f1, OP, FUNC, s, f2, Restrict)          \
 {                                                                              \
     /* Check fields have same size */                                          \
     checkFields(f1, f2, "f1 " #OP " " #FUNC "(s, f2)");                        \
                                                                                \
     /* Field access */                                                         \
-    List_ACCESS(typeF1, f1, f1P);                                              \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
+    TFOR_List_ACCESS(f1, f1P, Restrict);                                       \
+    TFOR_List_CONST_ACCESS(f2, f2P, Restrict);                                 \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: f1 OP1 f2 OP2 f3 */                                               \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -299,16 +260,27 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_FUNC_S_F(f1, OP, FUNC, s, f2) \
+        TSEQ_FORALL_F_OP_FUNC_S_F_impl(f1, OP, FUNC, s, f2, __restrict__)
+
+#define TSEQ_FORALL_F_OP_FUNC_S_F_inplace(f1, OP, FUNC, s, f2) \
+        TSEQ_FORALL_F_OP_FUNC_S_F_impl(f1, OP, FUNC, s, f2,)
 
+// old-style
+#define TFOR_ALL_F_OP_FUNC_S_F(typeF1, f1, OP, FUNC, typeS, s, typeF2, f2) \
+     TSEQ_FORALL_F_OP_FUNC_S_F_impl(f1, OP, FUNC, s, f2, __restrict__)
+
+
+// ------------
 // Binary Free Function : f1 OP FUNC(s1, s2)
 
-#define TFOR_ALL_F_OP_FUNC_S_S(typeF1, f1, OP, FUNC, typeS1, s1, typeS2, s2)   \
+#define TSEQ_FORALL_F_OP_FUNC_S_S_impl(f1, OP, FUNC, s1, s2, Unused)           \
 {                                                                              \
     /* Field access */                                                         \
-    List_ACCESS(typeF1, f1, f1P);                                              \
+    TFOR_List_ACCESS(f1, f1P, __restrict__);                                   \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: f1 OP FUNC(s1, s2) */                                             \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -317,20 +289,32 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_FUNC_S_S(f1, OP, FUNC, s1, s2) \
+        TSEQ_FORALL_F_OP_FUNC_S_S_impl(f1, OP, FUNC, s1, s2, __restrict__)
+
+#define TSEQ_FORALL_F_OP_FUNC_S_S_inplace(f1, OP, FUNC, s1, s2) \
+        TSEQ_FORALL_F_OP_FUNC_S_S_impl(f1, OP, FUNC, s1, s2,)
+
+// old-style
+#define TFOR_ALL_F_OP_FUNC_S_S(typeF1, f1, OP, FUNC, typeS1, s1, typeS2, s2) \
+     TSEQ_FORALL_F_OP_FUNC_S_S_impl(f1, OP, FUNC, s1, s2, __restrict__)
 
-// Unary Member Function : f1 OP f2 FUNC(s)
 
-#define TFOR_ALL_F_OP_F_FUNC_S(typeF1, f1, OP, typeF2, f2, FUNC, typeS, s)     \
+// ------------
+// Unary Member Function : f1 OP f2.FUNC(s)
+// NB: the '.' is passed in the func parameter
+
+#define TSEQ_FORALL_F_OP_F_FUNC_S_impl(f1, OP, f2, FUNC, s, Restrict)          \
 {                                                                              \
     /* Check fields have same size */                                          \
     checkFields(f1, f2, "f1 " #OP " f2 " #FUNC "(s)");                         \
                                                                                \
     /* Field access */                                                         \
-    List_ACCESS(typeF1, f1, f1P);                                              \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
+    TFOR_List_ACCESS(f1, f1P, Restrict);                                       \
+    TFOR_List_CONST_ACCESS(f2, f2P, Restrict);                                 \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: f1 OP f2 FUNC(s) */                                               \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -339,25 +323,35 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_F_FUNC_S(f1, OP, f2, FUNC, s) \
+        TSEQ_FORALL_F_OP_F_FUNC_S_impl(f1, OP, f2, FUNC, s, __restrict__)
+
+#define TSEQ_FORALL_F_OP_F_FUNC_S_inplace(f1, OP, f2, FUNC, s) \
+        TSEQ_FORALL_F_OP_F_FUNC_S_impl(f1, OP, f2, FUNC, s,)
+
+// old-style
+#define TFOR_ALL_F_OP_F_FUNC_S(typeF1, f1, OP, typeF2, f2, FUNC, typeS, s) \
+     TSEQ_FORALL_F_OP_F_FUNC_S_impl(f1, OP, f2, FUNC, s, __restrict__)
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// ------------
 // Ternary Free Function : f1 OP FUNC(f2, f3, f4)
 
-#define TFOR_ALL_F_OP_FUNC_F_F_F\
-(typeF1, f1, OP, FUNC, typeF2, f2, typeF3, f3, typeF4, f4)                     \
+#define TSEQ_FORALL_F_OP_FUNC_F_F_F_impl(f1, OP, FUNC, f2, f3, f4, Restrict)   \
 {                                                                              \
     /* Check fields have same size */                                          \
     checkFields(f1, f2, f3, f4, "f1 " #OP " " #FUNC "(f2, f3, f4)");           \
                                                                                \
     /* Field access */                                                         \
-    List_ACCESS(typeF1, f1, f1P);                                              \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
-    List_CONST_ACCESS(typeF3, f3, f3P);                                        \
-    List_CONST_ACCESS(typeF4, f4, f4P);                                        \
+    TFOR_List_ACCESS(f1, f1P, Restrict);                                       \
+    TFOR_List_CONST_ACCESS(f2, f2P, Restrict);                                 \
+    TFOR_List_CONST_ACCESS(f3, f3P, Restrict);                                 \
+    TFOR_List_CONST_ACCESS(f4, f4P, Restrict);                                 \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: f1 OP FUNC(f2, f3, f4) */                                         \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -366,21 +360,32 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_FUNC_F_F_F(f1, OP, FUNC, f2, f3, f4) \
+        TSEQ_FORALL_F_OP_FUNC_F_F_F_impl(f1, OP, FUNC, f2, f3, f4, __restrict__)
+
+#define TSEQ_FORALL_F_OP_FUNC_F_F_F_inplace(f1, OP, FUNC, f2, f3, f4) \
+        TSEQ_FORALL_F_OP_FUNC_F_F_F_impl(f1, OP, FUNC, f2, f3, f4,)
+
+// old-style
+#define TFOR_ALL_F_OP_FUNC_F_F_F(type1, f1, OP, FUNC, type2, f2, type3, f3, type4, f4) \
+     TSEQ_FORALL_F_OP_FUNC_F_F_F_impl(f1, OP, FUNC, f2, f3, f4, __restrict__)
+
+
+// ------------
 // Ternary Free Function : f1 OP FUNC(f2, f3, s4)
 
-#define TFOR_ALL_F_OP_FUNC_F_F_S\
-(typeF1, f1, OP, FUNC, typeF2, f2, typeF3, f3, typeF4, s4)                     \
+#define TSEQ_FORALL_F_OP_FUNC_F_F_S_impl(f1, OP, FUNC, f2, f3, s4, Restrict)   \
 {                                                                              \
     /* Check fields have same size */                                          \
     checkFields(f1, f2, f3, "f1 " #OP " " #FUNC "(f2, f3, s)");                \
                                                                                \
     /* Field access */                                                         \
-    List_ACCESS(typeF1, f1, f1P);                                              \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
-    List_CONST_ACCESS(typeF3, f3, f3P);                                        \
+    TFOR_List_ACCESS(f1, f1P, Restrict);                                       \
+    TFOR_List_CONST_ACCESS(f2, f2P, Restrict);                                 \
+    TFOR_List_CONST_ACCESS(f3, f3P, Restrict);                                 \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: f1 OP FUNC(f2, f3, s4) */                                         \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -389,23 +394,34 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_FUNC_F_F_S(f1, OP, FUNC, f2, f3, s4) \
+        TSEQ_FORALL_F_OP_FUNC_F_F_S_impl(f1, OP, FUNC, f2, f3, s4, __restrict__)
+
+#define TSEQ_FORALL_F_OP_FUNC_F_F_S_inplace(f1, OP, FUNC, f2, f3, s4) \
+        TSEQ_FORALL_F_OP_FUNC_F_F_S_impl(f1, OP, FUNC, f2, f3, s4,)
+
+// old-style
+#define TFOR_ALL_F_OP_FUNC_F_F_S(type1, f1, OP, FUNC, type2, f2, type3, f3, type4, s4) \
+     TSEQ_FORALL_F_OP_FUNC_F_F_S_impl(f1, OP, FUNC, f2, f3, s4, __restrict__)
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// ------------
 // Member operator : this field f1 OP1 f2 OP2 f3
 
-#define TFOR_ALL_F_OP_F_OP_F(typeF1, f1, OP1, typeF2, f2, OP2, typeF3, f3)     \
+#define TSEQ_FORALL_F_OP_F_OP_F_impl(f1, OP1, f2, OP2, f3, Restrict)           \
 {                                                                              \
     /* Check fields have same size */                                          \
     checkFields(f1, f2, f3, "f1 " #OP1 " f2 " #OP2 " f3");                     \
                                                                                \
     /* Field access */                                                         \
-    List_ACCESS(typeF1, f1, f1P);                                              \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
-    List_CONST_ACCESS(typeF3, f3, f3P);                                        \
+    TFOR_List_ACCESS(f1, f1P, Restrict);                                       \
+    TFOR_List_CONST_ACCESS(f2, f2P, Restrict);                                 \
+    TFOR_List_CONST_ACCESS(f3, f3P, Restrict);                                 \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: f1 OP1 f2 OP2 f3 */                                               \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -414,20 +430,31 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_F_OP_F(f1, OP1, f2, OP2, f3) \
+        TSEQ_FORALL_F_OP_F_OP_F_impl(f1, OP1, f2, OP2, f3, __restrict__)
+
+#define TSEQ_FORALL_F_OP_F_OP_F_inplace(f1, OP1, f2, OP2, f3) \
+        TSEQ_FORALL_F_OP_F_OP_F_impl(f1, OP1, f2, OP2, f3,)
+
+// old-style
+#define TFOR_ALL_F_OP_F_OP_F(type1, f1, OP1, type2, f2, OP2, type3, f3) \
+     TSEQ_FORALL_F_OP_F_OP_F_impl(f1, OP1, f2, OP2, f3, __restrict__)
 
+
+// ------------
 // Member operator : this field f1 OP1 s OP2 f2
 
-#define TFOR_ALL_F_OP_S_OP_F(typeF1, f1, OP1, typeS, s, OP2, typeF2, f2)       \
+#define TSEQ_FORALL_F_OP_S_OP_F_impl(f1, OP1, s, OP2, f2, Restrict)            \
 {                                                                              \
     /* Check fields have same size */                                          \
     checkFields(f1, f2, "f1 " #OP1 " s " #OP2 " f2");                          \
                                                                                \
     /* Field access */                                                         \
-    List_ACCESS(typeF1, f1, f1P);                                              \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
+    TFOR_List_ACCESS(f1, f1P, Restrict);                                       \
+    TFOR_List_CONST_ACCESS(f2, f2P, Restrict);                                 \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: f1 OP1 s OP2 f2 */                                                \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -436,20 +463,31 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_S_OP_F(f1, OP1, s, OP2, f2) \
+        TSEQ_FORALL_F_OP_S_OP_F_impl(f1, OP1, s, OP2, f2, __restrict__)
+
+#define TSEQ_FORALL_F_OP_S_OP_F_inplace(f1, OP1, s, OP2, f2) \
+        TSEQ_FORALL_F_OP_S_OP_F_impl(f1, OP1, s, OP2, f2,)
 
+// old-style
+#define TFOR_ALL_F_OP_S_OP_F(type1, f1, OP1, typeS, s, OP2, type2, f2) \
+    TSEQ_FORALL_F_OP_S_OP_F_impl(f1, OP1, s, OP2, f2, __restrict__)
+
+
+// ------------
 // Member operator : this field f1 OP1 f2 OP2 s
 
-#define TFOR_ALL_F_OP_F_OP_S(typeF1, f1, OP1, typeF2, f2, OP2, typeS, s)       \
+#define TSEQ_FORALL_F_OP_F_OP_S_impl(f1, OP1, f2, OP2, s, Restrict)            \
 {                                                                              \
     /* Check fields have same size */                                          \
     checkFields(f1, f2, "f1 " #OP1 " f2 " #OP2 " s");                          \
                                                                                \
     /* Field access */                                                         \
-    List_ACCESS(typeF1, f1, f1P);                                              \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
+    TFOR_List_ACCESS(f1, f1P, Restrict);                                       \
+    TFOR_List_CONST_ACCESS(f2, f2P, Restrict);                                 \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop f1 OP1 s OP2 f2 */                                                 \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -458,20 +496,31 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_F_OP_S(f1, OP1, f2, OP2, s) \
+        TSEQ_FORALL_F_OP_F_OP_S_impl(f1, OP1, f2, OP2, s, __restrict__)
+
+#define TSEQ_FORALL_F_OP_F_OP_S_inplace(f1, OP1, f2, OP2, s) \
+        TSEQ_FORALL_F_OP_F_OP_S_impl(f1, OP1, f2, OP2, s,)
+
+// old-style
+#define TFOR_ALL_F_OP_F_OP_S(type1, f1, OP1, type2, f2, OP2, typeS, s) \
+     TSEQ_FORALL_F_OP_F_OP_S_impl(f1, OP1, f2, OP2, s, __restrict__)
+
 
+// ------------
 // Member operator : this field f1 OP f2
 
-#define TFOR_ALL_F_OP_F(typeF1, f1, OP, typeF2, f2)                            \
+#define TSEQ_FORALL_F_OP_F_impl(f1, OP, f2, Restrict)                          \
 {                                                                              \
     /* Check fields have same size */                                          \
     checkFields(f1, f2, "f1 " #OP " f2");                                      \
                                                                                \
     /* Field access */                                                         \
-    List_ACCESS(typeF1, f1, f1P);                                              \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
+    TFOR_List_ACCESS(f1, f1P, Restrict);                                       \
+    TFOR_List_CONST_ACCESS(f2, f2P, Restrict);                                 \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: f1 OP f2 */                                                       \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -480,20 +529,31 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_F(f1, OP, f2) \
+        TSEQ_FORALL_F_OP_F_impl(f1, OP, f2, __restrict__)
 
+#define TSEQ_FORALL_F_OP_F_inplace(f1, OP, f2) \
+        TSEQ_FORALL_F_OP_F_impl(f1, OP, f2,)
+
+// old-style
+#define TFOR_ALL_F_OP_F(type1, f1, OP, type2, f2) \
+     TSEQ_FORALL_F_OP_F_impl(f1, OP, f2, __restrict__)
+
+
+// ------------
 // Member operator : this field f1 OP1 OP2 f2
 
-#define TFOR_ALL_F_OP_OP_F(typeF1, f1, OP1, OP2, typeF2, f2)                   \
+#define TSEQ_FORALL_F_OP_OP_F_impl(f1, OP1, OP2, f2, Restrict)                 \
 {                                                                              \
     /* Check fields have same size */                                          \
     checkFields(f1, f2, #OP1 " " #OP2 " f2");                                  \
                                                                                \
     /* Field access */                                                         \
-    List_ACCESS(typeF1, f1, f1P);                                              \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
+    TFOR_List_ACCESS(f1, f1P, Restrict);                                       \
+    TFOR_List_CONST_ACCESS(f2, f2P, Restrict);                                 \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: f1 OP1 OP2 f2 */                                                  \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -502,16 +562,27 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_OP_F(f1, OP1, OP2, f2) \
+        TSEQ_FORALL_F_OP_OP_F_impl(f1, OP1, OP2, f2, __restrict__)
+
+#define TSEQ_FORALL_F_OP_OP_F_inplace(f1, OP1, OP2, f2) \
+        TSEQ_FORALL_F_OP_OP_F_impl(f1, OP1, OP2, f2,)
+
+// old-style
+#define TFOR_ALL_F_OP_OP_F(typeF1, f1, OP1, OP2, typeF2, f2) \
+     TSEQ_FORALL_F_OP_OP_F_impl(f1, OP1, OP2, f2, __restrict__)
 
+
+// ------------
 // Member operator : this field f OP s
 
-#define TFOR_ALL_F_OP_S(typeF, f, OP, typeS, s)                                \
+#define TSEQ_FORALL_F_OP_S_impl(f, OP, s, Unused)                              \
 {                                                                              \
     /* Field access */                                                         \
-    List_ACCESS(typeF, f, fP);                                                 \
+    TFOR_List_ACCESS(f, fP, __restrict__);                                     \
+    TFOR_List_LENGTH(f, loop_len);                                             \
                                                                                \
     /* Loop: f OP s */                                                         \
-    const label loop_len = (f).size();                                         \
                                                                                \
     /* pragmas... */                                                           \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -520,18 +591,29 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_F_OP_S(f, OP, s) \
+        TSEQ_FORALL_F_OP_S_impl(f, OP, s, __restrict__)
+
+#define TSEQ_FORALL_F_OP_S_inplace(f, OP, s) \
+        TSEQ_FORALL_F_OP_S_impl(f, OP, s,)
+
+// old-style
+#define TFOR_ALL_F_OP_S(typeF, f, OP, typeS, s) \
+     TSEQ_FORALL_F_OP_S_impl(f, OP, s, __restrict__)
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Friend operator function : s OP f, allocates storage for s
+// ------------
+// Friend operator function : s OP f
 
-#define TFOR_ALL_S_OP_F(typeS, s, OP, typeF, f)                                \
+#define TSEQ_FORALL_S_OP_F_impl(s, OP, f, Unused)                              \
 {                                                                              \
     /* Field access */                                                         \
-    List_CONST_ACCESS(typeF, f, fP);                                           \
+    TFOR_List_CONST_ACCESS(f, fP, __restrict__);                               \
+    TFOR_List_LENGTH(f, loop_len);                                             \
                                                                                \
     /* Loop: s OP f */                                                         \
-    const label loop_len = (f).size();                                         \
                                                                                \
     /* pragmas, reduction... */                                                \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -540,17 +622,28 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_S_OP_F(s, OP, f) \
+        TSEQ_FORALL_S_OP_F_impl(s, OP, f, __restrict__)
+
+#define TSEQ_FORALL_S_OP_F_inplace(s, OP, f) \
+        TSEQ_FORALL_S_OP_F_impl(s, OP, f,)
 
-// Friend operator function : s OP1 f1 OP2 f2, allocates storage for s
+// old-style
+#define TFOR_ALL_S_OP_F(typeS, s, OP, typeF, f) \
+     TSEQ_FORALL_S_OP_F_impl(s, OP, f, __restrict__)
 
-#define TFOR_ALL_S_OP_F_OP_F(typeS, s, OP1, typeF1, f1, OP2, typeF2, f2)       \
+
+// ------------
+// Friend operator function : s OP1 f1 OP2 f2
+
+#define TSEQ_FORALL_S_OP_F_OP_F_impl(s, OP1, f1, OP2, f2, Unused)              \
 {                                                                              \
     /* Field access */                                                         \
-    List_CONST_ACCESS(typeF1, f1, f1P);                                        \
-    List_CONST_ACCESS(typeF2, f2, f2P);                                        \
+    TFOR_List_ACCESS(f1, f1P, __restrict__);                                   \
+    TFOR_List_CONST_ACCESS(f2, f2P, __restrict__);                             \
+    TFOR_List_LENGTH(f1, loop_len);                                            \
                                                                                \
     /* Loop: s OP1 f1 OP2 f2 */                                                \
-    const label loop_len = (f1).size();                                        \
                                                                                \
     /* pragmas, reduction... */                                                \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -559,15 +652,27 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_S_OP_F_OP_F(s, OP1, f1, OP2, f2) \
+        TSEQ_FORALL_S_OP_F_OP_F_impl(s, OP1, f1, OP2, f2, __restrict__)
+
+#define TSEQ_FORALL_S_OP_F_OP_F_inplace(s, OP1, f1, OP2, f2) \
+        TSEQ_FORALL_S_OP_F_OP_F_impl(s, OP1, f1, OP2, f2,)
+
+// old-style
+#define TFOR_ALL_S_OP_F_OP_F(typeS, s, OP1, type1, f1, OP2, type2, f2) \
+     TSEQ_FORALL_S_OP_F_OP_F_impl(s, OP1, f1, OP2, f2, __restrict__)
+
 
-// Friend operator function : s OP FUNC(f), allocates storage for s
-#define TFOR_ALL_S_OP_FUNC_F(typeS, s, OP, FUNC, typeF, f)                     \
+// ------------
+// Friend operator function : s OP FUNC(f)
+
+#define TSEQ_FORALL_S_OP_FUNC_F_impl(s, OP, FUNC, f, Unused)                   \
 {                                                                              \
     /* Field access */                                                         \
-    List_CONST_ACCESS(typeF, f, fP);                                           \
+    TFOR_List_CONST_ACCESS(f, fP, __restrict__);                               \
+    TFOR_List_LENGTH(f, loop_len);                                             \
                                                                                \
     /* Loop: s OP FUNC(f) */                                                   \
-    const label loop_len = (f).size();                                         \
                                                                                \
     /* pragmas, reduction... */                                                \
     for (label i = 0; i < loop_len; ++i)                                       \
@@ -576,6 +681,16 @@ void checkFields
     }                                                                          \
 }
 
+#define TSEQ_FORALL_S_OP_FUNC_F(s, OP, FUNC, f) \
+        TSEQ_FORALL_S_OP_FUNC_F_impl(s, OP, FUNC, f, __restrict__)
+
+#define TSEQ_FORALL_S_OP_FUNC_F_inplace(s, OP, FUNC, f) \
+        TSEQ_FORALL_S_OP_FUNC_F_impl(s, OP, FUNC, f,)
+
+// old-style
+#define TFOR_ALL_S_OP_FUNC_F(typeS, s, OP, FUNC, typeF, f) \
+     TSEQ_FORALL_S_OP_FUNC_F_impl(s, OP, FUNC, f, __restrict__)
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/fields/Fields/Field/FieldReuseFunctions.H b/src/OpenFOAM/fields/Fields/Field/FieldReuseFunctions.H
index 0d23585fe97e036719900484295feb43d51405ea..c5e03634cd348604d9ed780f9403d06d9e5d724c 100644
--- a/src/OpenFOAM/fields/Fields/Field/FieldReuseFunctions.H
+++ b/src/OpenFOAM/fields/Fields/Field/FieldReuseFunctions.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,6 +40,15 @@ namespace Foam
 template<class TypeR, class Type1>
 struct reuseTmp
 {
+    //- Dissimilar types: non-reusable
+    static bool reusable
+    (
+        const tmp<Field<Type1>>& tf1
+    )
+    {
+        return false;
+    }
+
     //- Pass-through to tmp New
     static tmp<Field<TypeR>> New(const Field<Type1>& f1)
     {
@@ -57,6 +66,15 @@ struct reuseTmp
 template<class TypeR>
 struct reuseTmp<TypeR, TypeR>
 {
+    //- Identical types: possibly reusable
+    static bool reusable
+    (
+        const tmp<Field<TypeR>>& tf1
+    )
+    {
+        return tf1.movable();
+    }
+
     //- Identical input and return types:
     //- allow optional copy assignment of the initial content
     static tmp<Field<TypeR>> New
@@ -99,6 +117,16 @@ template<class TypeR> tmp<Field<TypeR>> New
 template<class TypeR, class Type1, class Type12, class Type2>
 struct reuseTmpTmp
 {
+    //- Dissimilar types: non-reusable
+    static bool reusable
+    (
+        const tmp<Field<Type1>>& tf1,
+        const tmp<Field<Type2>>& tf2
+    )
+    {
+        return false;
+    }
+
     //- Dissimilar types: just use size
     static tmp<Field<TypeR>> New
     (
@@ -114,6 +142,16 @@ struct reuseTmpTmp
 template<class TypeR, class Type1, class Type12>
 struct reuseTmpTmp<TypeR, Type1, Type12, TypeR>
 {
+    //- Second input has return type: possibly reusable
+    static bool reusable
+    (
+        const tmp<Field<Type1>>& tf1,
+        const tmp<Field<TypeR>>& tf2
+    )
+    {
+        return tf2.movable();
+    }
+
     //- Second input has return type
     static tmp<Field<TypeR>> New
     (
@@ -134,6 +172,16 @@ struct reuseTmpTmp<TypeR, Type1, Type12, TypeR>
 template<class TypeR, class Type2>
 struct reuseTmpTmp<TypeR, TypeR, TypeR, Type2>
 {
+    //- First input has return type: possibly reusable
+    static bool reusable
+    (
+        const tmp<Field<TypeR>>& tf1,
+        const tmp<Field<Type2>>& tf2
+    )
+    {
+        return tf1.movable();
+    }
+
     //- First input has return type
     static tmp<Field<TypeR>> New
     (
@@ -154,6 +202,16 @@ struct reuseTmpTmp<TypeR, TypeR, TypeR, Type2>
 template<class TypeR>
 struct reuseTmpTmp<TypeR, TypeR, TypeR, TypeR>
 {
+    //- Both inputs have return type: possibly reusable
+    static bool reusable
+    (
+        const tmp<Field<TypeR>>& tf1,
+        const tmp<Field<TypeR>>& tf2
+    )
+    {
+        return (tf1.movable() || tf2.movable());
+    }
+
     //- Both inputs have return type
     static tmp<Field<TypeR>> New
     (
diff --git a/src/OpenFOAM/fields/Fields/Field/ListLoopM.H b/src/OpenFOAM/fields/Fields/Field/ListLoopM.H
index 60479fdda5913d5543fc0c0b630e79ccf00fcbfe..689144c7a73b953fb6b453529de1a505f771cd1b 100644
--- a/src/OpenFOAM/fields/Fields/Field/ListLoopM.H
+++ b/src/OpenFOAM/fields/Fields/Field/ListLoopM.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,6 +32,22 @@ Description
 #ifndef Foam_ListLoopM_H
 #define Foam_ListLoopM_H
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Element access looping
+
+// The List size
+#define TFOR_List_LENGTH(f, n)                  \
+    const Foam::label n = (f).size()
+
+// List begin, non-const access with __restrict__ (or not)
+#define TFOR_List_ACCESS(f, fp, Restrict)       \
+    auto* const Restrict fp = (f).begin()
+
+// List begin, const access with __restrict__ (or not)
+#define TFOR_List_CONST_ACCESS(f, fp, Restrict) \
+    const auto* const Restrict fp = (f).begin()
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Element access looping
 
@@ -45,8 +61,8 @@ Description
 
 // Loop over all elements
 #define List_FOR_ALL(f, i)                      \
-        const label _n##i = (f).size();         \
-        for (label i=0; i<_n##i; ++i)
+        const Foam::label _n##i = (f).size();   \
+        for (Foam::label i = 0; i < _n##i; ++i)
 
 // Current element (non-const access)
 #define List_ELEM(fp, i)  (fp[i])
diff --git a/src/OpenFOAM/fields/Fields/Field/errorCheckFields.H b/src/OpenFOAM/fields/Fields/Field/errorCheckFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..fbc6cc21e5959903b4de06ac6d28c05d679d1690
--- /dev/null
+++ b/src/OpenFOAM/fields/Fields/Field/errorCheckFields.H
@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2024 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/>.
+
+Description
+    Field size checks for field operations
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_errorCheckFields_H
+#define Foam_errorCheckFields_H
+
+#include "error.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type1, class Type2>
+void checkFields
+(
+    const UList<Type1>& f1,
+    const UList<Type2>& f2,
+    const char* op
+)
+{
+    #ifdef FULLDEBUG
+    if (f1.size() != f2.size())
+    {
+        FatalErrorInFunction
+            << " Field<"<<pTraits<Type1>::typeName<<"> f1("<<f1.size()<<')'
+            << " and Field<"<<pTraits<Type2>::typeName<<"> f2("<<f2.size()<<')'
+            << nl
+            << " for operation " << op << endl
+            << abort(FatalError);
+    }
+    #endif
+}
+
+
+template<class Type1, class Type2, class Type3>
+void checkFields
+(
+    const UList<Type1>& f1,
+    const UList<Type2>& f2,
+    const UList<Type3>& f3,
+    const char* op
+)
+{
+    #ifdef FULLDEBUG
+    if
+    (
+        f1.size() != f2.size()
+     || f1.size() != f3.size()
+    )
+    {
+        FatalErrorInFunction
+            << " Field<"<<pTraits<Type1>::typeName<<"> f1("<<f1.size()<<')'
+            << ", Field<"<<pTraits<Type2>::typeName<<"> f2("<<f2.size()<<')'
+            << " and Field<"<<pTraits<Type3>::typeName<<"> f3("<<f3.size()<<')'
+            << nl
+            << " for operation " << op << endl
+            << abort(FatalError);
+    }
+    #endif
+}
+
+
+template<class Type1, class Type2, class Type3, class Type4>
+void checkFields
+(
+    const UList<Type1>& f1,
+    const UList<Type2>& f2,
+    const UList<Type3>& f3,
+    const UList<Type4>& f4,
+    const char* op
+)
+{
+    #ifdef FULLDEBUG
+    if
+    (
+        f1.size() != f2.size()
+     || f1.size() != f3.size()
+     || f1.size() != f4.size()
+    )
+    {
+        FatalErrorInFunction
+            << " Field<"<<pTraits<Type1>::typeName<<"> f1("<<f1.size()<<')'
+            << ", Field<"<<pTraits<Type2>::typeName<<"> f2("<<f2.size()<<')'
+            << ", Field<"<<pTraits<Type3>::typeName<<"> f3("<<f3.size()<<')'
+            << " and Field<"<<pTraits<Type4>::typeName<<"> f4("<<f4.size()<<')'
+            << nl
+            << " for operation " << op << endl
+            << abort(FatalError);
+    }
+    #endif
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/Fields/boolField/boolField.C b/src/OpenFOAM/fields/Fields/boolField/boolField.C
index 6d4caf0408860ac96e936cff7accfdef2b16c0f2..b161fbe29d4b2bb2c4251a5c8fbd4663de4e7d44 100644
--- a/src/OpenFOAM/fields/Fields/boolField/boolField.C
+++ b/src/OpenFOAM/fields/Fields/boolField/boolField.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,7 +45,7 @@ template<>
 void boolField::negate()
 {
     // std::for_each, std::logical_not
-    TFOR_ALL_F_OP_OP_F(bool, *this, =, !, bool, *this)
+    TSEQ_FORALL_F_OP_OP_F_inplace(*this, =, !, *this)
 }
 
 
diff --git a/src/OpenFOAM/fields/Fields/complex/complexField.C b/src/OpenFOAM/fields/Fields/complex/complexField.C
index 9614215088c565ea8818091edb79ab9b9b16ddd6..abb3fe8d0a0cd5fff826bac3c9943bf40392cf6b 100644
--- a/src/OpenFOAM/fields/Fields/complex/complexField.C
+++ b/src/OpenFOAM/fields/Fields/complex/complexField.C
@@ -255,7 +255,7 @@ complex sumProd(const UList<complex>& f1, const UList<complex>& f2)
     if (f1.size() && (f1.size() == f2.size()))
     {
         // std::inner_product
-        TFOR_ALL_S_OP_F_OP_F(complex, result, +=, complex, f1, *, complex, f2)
+        TSEQ_FORALL_S_OP_F_OP_F(result, +=, f1, *, f2)
     }
     return result;
 }
diff --git a/src/OpenFOAM/fields/Fields/scalarField/scalarField.C b/src/OpenFOAM/fields/Fields/scalarField/scalarField.C
index d1f00a8c8cccd0aa326878620e1ccee32a77b795..82c2cb7d62d42170e2df0a1ea6c48a5df4a79249 100644
--- a/src/OpenFOAM/fields/Fields/scalarField/scalarField.C
+++ b/src/OpenFOAM/fields/Fields/scalarField/scalarField.C
@@ -68,11 +68,16 @@ void scalarField::replace(const direction, const scalar& s)
 
 void stabilise(scalarField& res, const UList<scalar>& sf, const scalar s)
 {
-    // std::transform
-    TFOR_ALL_F_OP_FUNC_S_F
-    (
-        scalar, res, =, ::Foam::stabilise, scalar, s, scalar, sf
-    )
+    if (res.cdata() == sf.cdata())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_S_F_inplace(res, =, ::Foam::stabilise, s, sf)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_S_F(res, =, ::Foam::stabilise, s, sf)
+    }
 }
 
 tmp<scalarField> stabilise(const UList<scalar>& sf, const scalar s)
@@ -96,11 +101,11 @@ tmp<scalarField> stabilise(const tmp<scalarField>& tsf, const scalar s)
 template<>
 float sumProd(const UList<float>& f1, const UList<float>& f2)
 {
-    float result = 0.0;
+    float result(0);
     if (f1.size() && (f1.size() == f2.size()))
     {
         // std::inner_product
-        TFOR_ALL_S_OP_F_OP_F(float, result, +=, float, f1, *, float, f2)
+        TSEQ_FORALL_S_OP_F_OP_F(result, +=, f1, *, f2)
     }
     return result;
 }
@@ -109,11 +114,11 @@ float sumProd(const UList<float>& f1, const UList<float>& f2)
 template<>
 double sumProd(const UList<double>& f1, const UList<double>& f2)
 {
-    double result = 0.0;
+    double result(0);
     if (f1.size() && (f1.size() == f2.size()))
     {
         // std::inner_product
-        TFOR_ALL_S_OP_F_OP_F(double, result, +=, double, f1, *, double, f2)
+        TSEQ_FORALL_S_OP_F_OP_F(result, +=, f1, *, f2)
     }
     return result;
 }
@@ -188,7 +193,16 @@ UNARY_FUNCTION(scalar, scalar, paToBar)
 #define BesselFunc(func)                                                       \
 void func(scalarField& res, const int n, const UList<scalar>& sf)              \
 {                                                                              \
-    TFOR_ALL_F_OP_FUNC_S_F(scalar, res, =, ::Foam::func, int, n, scalar, sf)   \
+    if (res.cdata() == sf.cdata())                                             \
+    {                                                                          \
+        /* std::for_each */                                                    \
+        TSEQ_FORALL_F_OP_FUNC_S_F_inplace(res, =, ::Foam::func, n, sf)         \
+    }                                                                          \
+    else                                                                       \
+    {                                                                          \
+        /* std::transform */                                                   \
+        TSEQ_FORALL_F_OP_FUNC_S_F(res, =, ::Foam::func, n, sf)                 \
+    }                                                                          \
 }                                                                              \
                                                                                \
 tmp<scalarField> func(const int n, const UList<scalar>& sf)                    \
diff --git a/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C b/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
index 1894e8961fb95064f85f368966c6a36428d82ffe..902739a657af7ac71e3e14dc78cab86dc9687aa7 100644
--- a/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
+++ b/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,8 +54,16 @@ UNARY_FUNCTION(symmTensor, symmTensor, cof)
 void inv(Field<symmTensor>& result, const UList<symmTensor>& f1)
 {
     // With 'failsafe' invert
-    // std::transform
-    TFOR_ALL_F_OP_F_FUNC(symmTensor, result, =, symmTensor, f1, safeInv)
+    if (result.cdata() == f1.cdata())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_F_FUNC_inplace(result, =, f1, safeInv)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_F_FUNC(result, =, f1, safeInv)
+    }
 }
 
 tmp<symmTensorField> inv(const UList<symmTensor>& tf)
diff --git a/src/OpenFOAM/fields/Fields/symmTransformField/symmTransformField.C b/src/OpenFOAM/fields/Fields/symmTransformField/symmTransformField.C
index eee4cb2869c8f9bcf95b85be65fe60685cced61e..cca7336e79eca4436652c76b900f55793fcbff5d 100644
--- a/src/OpenFOAM/fields/Fields/symmTransformField/symmTransformField.C
+++ b/src/OpenFOAM/fields/Fields/symmTransformField/symmTransformField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2023 OpenCFD Ltd.
+    Copyright (C) 2023-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,11 +39,16 @@ void Foam::transform
     const Field<Type>& fld
 )
 {
-    // std::transform
-    TFOR_ALL_F_OP_FUNC_S_F
-    (
-        Type, result, =, transform, symmTensor, rot, Type, fld
-    );
+    if (result.cdata() == fld.cdata())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_S_F_inplace(result, =, transform, rot, fld)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_S_F(result, =, transform, rot, fld)
+    }
 }
 
 
@@ -55,16 +60,23 @@ void Foam::transform
     const Field<Type>& fld
 )
 {
+    // Does not handle corner case where result and rot are identical...
+
     if (rot.size() == 1)
     {
         return transform(result, rot.front(), fld);
     }
 
-    // std::transform
-    TFOR_ALL_F_OP_FUNC_F_F
-    (
-        Type, result, =, transform, symmTensor, rot, Type, fld
-    );
+    if (result.cdata() == fld.cdata())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_F_F_inplace(result, =, transform, rot, fld)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_F_F(result, =, transform, rot, fld)
+    }
 }
 
 
diff --git a/src/OpenFOAM/fields/Fields/tensorField/tensorField.C b/src/OpenFOAM/fields/Fields/tensorField/tensorField.C
index 4591719aaac58c9e7216861d5a198c4bd4b60512..f5bd683b3eb29443c7b3639a3d529e31c9ac4d88 100644
--- a/src/OpenFOAM/fields/Fields/tensorField/tensorField.C
+++ b/src/OpenFOAM/fields/Fields/tensorField/tensorField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,8 +55,16 @@ UNARY_FUNCTION(tensor, tensor, cof)
 void inv(Field<tensor>& result, const UList<tensor>& f1)
 {
     // With 'failsafe' invert
-    // std::transform
-    TFOR_ALL_F_OP_F_FUNC(tensor, result, =, tensor, f1, safeInv)
+    if (result.cdata() == f1.cdata())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_F_FUNC_inplace(result, =, f1, safeInv)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_F_FUNC(result, =, f1, safeInv)
+    }
 }
 
 tmp<tensorField> inv(const UList<tensor>& tf)
@@ -88,7 +96,7 @@ tmp<Field<tensor>> transformFieldMask<tensor>
 {
     auto tres = tmp<tensorField>::New(stf.size());
     auto& res = tres.ref();
-    TFOR_ALL_F_OP_F(tensor, res, =, symmTensor, stf)
+    TSEQ_FORALL_F_OP_F(res, =, stf)
     return tres;
 }
 
diff --git a/src/OpenFOAM/fields/Fields/transformField/transformField.C b/src/OpenFOAM/fields/Fields/transformField/transformField.C
index e069e324909f00ae26369e4c5f180fece27a95ef..aaaf1dc6037cbacb8886e17352d3ad85bf8bb736 100644
--- a/src/OpenFOAM/fields/Fields/transformField/transformField.C
+++ b/src/OpenFOAM/fields/Fields/transformField/transformField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,14 +34,22 @@ License
 
 void Foam::transform
 (
-    vectorField& rtf,
+    vectorField& result,
     const quaternion& q,
-    const vectorField& tf
+    const vectorField& fld
 )
 {
     const tensor rot = q.R();
-    // std::transform
-    TFOR_ALL_F_OP_FUNC_S_F(vector, rtf, =, transform, tensor, rot, vector, tf)
+    if (result.cdata() == fld.cdata())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_S_F_inplace(result, =, transform, rot, fld)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_S_F(result, =, transform, rot, fld)
+    }
 }
 
 
@@ -82,8 +90,16 @@ void Foam::transformPoints
     // Check if any translation
     if (mag(trans) > VSMALL)
     {
-        // std::transform
-        TFOR_ALL_F_OP_F_OP_S(vector, result, =, vector, fld, -, vector, trans);
+        if (result.cdata() == fld.cdata())
+        {
+            // std::for_each
+            TSEQ_FORALL_F_OP_F_OP_S_inplace(result, =, fld, -, trans)
+        }
+        else
+        {
+            // std::transform
+            TSEQ_FORALL_F_OP_F_OP_S(result, =, fld, -, trans)
+        }
     }
     else
     {
diff --git a/src/OpenFOAM/fields/Fields/transformField/transformFieldTemplates.C b/src/OpenFOAM/fields/Fields/transformField/transformFieldTemplates.C
index af032f213d9bf899045ce958bebef4715edca0d1..63e0e0460f2bbc49136f5ca1933636bda091dd1b 100644
--- a/src/OpenFOAM/fields/Fields/transformField/transformFieldTemplates.C
+++ b/src/OpenFOAM/fields/Fields/transformField/transformFieldTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,11 +39,16 @@ void Foam::transform
     const Field<Type>& fld
 )
 {
-    // std::transform
-    TFOR_ALL_F_OP_FUNC_S_F
-    (
-        Type, result, =, transform, tensor, rot, Type, fld
-    );
+    if (result.cdata() == fld.cdata())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_S_F_inplace(result, =, transform, rot, fld)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_S_F(result, =, transform, rot, fld)
+    }
 }
 
 
@@ -55,16 +60,23 @@ void Foam::transform
     const Field<Type>& fld
 )
 {
+    // Does not handle corner case where result and rot are identical...
+
     if (rot.size() == 1)
     {
         return transform(result, rot.front(), fld);
     }
 
-    // std::transform
-    TFOR_ALL_F_OP_FUNC_F_F
-    (
-        Type, result, =, transform, tensor, rot, Type, fld
-    );
+    if (result.cdata() == fld.cdata())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_F_F_inplace(result, =, transform, rot, fld)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_F_F(result, =, transform, rot, fld)
+    }
 }
 
 
@@ -165,11 +177,16 @@ void Foam::invTransform
     const Field<Type>& fld
 )
 {
-    // std::transform
-    TFOR_ALL_F_OP_FUNC_S_F
-    (
-        Type, result, =, invTransform, tensor, rot, Type, fld
-    );
+    if (result.cdata() == fld.cdata())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_S_F_inplace(result, =, invTransform, rot, fld)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_S_F(result, =, invTransform, rot, fld)
+    }
 }
 
 
@@ -186,11 +203,16 @@ void Foam::invTransform
         return invTransform(result, rot.front(), fld);
     }
 
-    // std::transform
-    TFOR_ALL_F_OP_FUNC_F_F
-    (
-        Type, result, =, invTransform, tensor, rot, Type, fld
-    );
+    if (result.cdata() == fld.cdata())
+    {
+        // std::for_each
+        TSEQ_FORALL_F_OP_FUNC_F_F_inplace(result, =, invTransform, rot, fld)
+    }
+    else
+    {
+        // std::transform
+        TSEQ_FORALL_F_OP_FUNC_F_F(result, =, invTransform, rot, fld)
+    }
 }
 
 
diff --git a/src/OpenFOAM/fields/Fields/transformList/transformList.C b/src/OpenFOAM/fields/Fields/transformList/transformList.C
index 4ed8482718434824f66c554db74bdfb932c24a90..739d2ae4536c3aa43a5d915ba314dcd63cd9bd32 100644
--- a/src/OpenFOAM/fields/Fields/transformList/transformList.C
+++ b/src/OpenFOAM/fields/Fields/transformList/transformList.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,12 +37,12 @@ Foam::List<T> Foam::transform
     const UList<T>& field
 )
 {
-    const label loopLen = field.size();
+    const label loop_len = field.size();
 
-    List<T> result(loopLen);
+    List<T> result(loop_len);
 
     /* pragmas... */
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         result[i] = transform(rotTensor, field[i]);
     }
@@ -54,10 +54,10 @@ Foam::List<T> Foam::transform
 template<class T>
 void Foam::transformList(const tensor& rotTensor, UList<T>& field)
 {
-    const label loopLen = field.size();
+    const label loop_len = field.size();
 
     /* pragmas... */
-    for (label i = 0; i < loopLen; ++i)
+    for (label i = 0; i < loop_len; ++i)
     {
         field[i] = transform(rotTensor, field[i]);
     }
@@ -73,10 +73,10 @@ void Foam::transformList(const tensorField& rotTensor, UList<T>& field)
     }
     else if (rotTensor.size() == field.size())
     {
-        const label loopLen = field.size();
+        const label loop_len = field.size();
 
         /* pragmas... */
-        for (label i = 0; i < loopLen; ++i)
+        for (label i = 0; i < loop_len; ++i)
         {
             field[i] = transform(rotTensor[i], field[i]);
         }
@@ -96,7 +96,7 @@ void Foam::transformList(const tensor& rotTensor, Map<T>& field)
 {
     forAllIters(field, iter)
     {
-        T& value = iter.val();
+        auto& value = iter.val();
         value = transform(rotTensor, value);
     }
 }
@@ -124,7 +124,7 @@ void Foam::transformList(const tensor& rotTensor, EdgeMap<T>& field)
 {
     forAllIters(field, iter)
     {
-        T& value = iter.val();
+        auto& value = iter.val();
         value = transform(rotTensor, value);
     }
 }
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldReuseFunctions.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldReuseFunctions.H
index edc15dea5a12db0092d00749d20d845fa41b49e6..40c3a4f8308d9f44d99007da70f475da2fb7ccc9 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldReuseFunctions.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldReuseFunctions.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -87,6 +87,15 @@ template
 >
 struct reuseTmpGeometricField
 {
+    //- Dissimilar types: non-reusable
+    static bool reusable
+    (
+        const tmp<GeometricField<Type1, PatchField, GeoMesh>>& tf1
+    )
+    {
+        return false;
+    }
+
     //- Pass-through to New GeometricField
     static tmp<GeometricField<TypeR, PatchField, GeoMesh>> New
     (
@@ -126,6 +135,15 @@ struct reuseTmpGeometricField
 template<class TypeR, template<class> class PatchField, class GeoMesh>
 struct reuseTmpGeometricField<TypeR, TypeR, PatchField, GeoMesh>
 {
+    //- Identical types: possibly reusable
+    static bool reusable
+    (
+        const tmp<GeometricField<TypeR, PatchField, GeoMesh>>& tf1
+    )
+    {
+        return Detail::reusable(tf1);
+    }
+
     //- Allow optional copy assignment of the initial content
     //- for identical input and output types
     static tmp<GeometricField<TypeR, PatchField, GeoMesh>> New
@@ -156,7 +174,7 @@ struct reuseTmpGeometricField<TypeR, TypeR, PatchField, GeoMesh>
 
         if (initCopy)
         {
-            tresult.ref() == tf1();
+            tresult.ref() == f1;
         }
 
         return tresult;
@@ -187,6 +205,9 @@ tmp
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Two-parameter versions
+
 template
 <
     class TypeR,
@@ -198,6 +219,16 @@ template
 >
 struct reuseTmpTmpGeometricField
 {
+    //- Dissimilar types: non-reusable
+    static bool reusable
+    (
+        const tmp<GeometricField<Type1, PatchField, GeoMesh>>& tf1,
+        const tmp<GeometricField<Type2, PatchField, GeoMesh>>& tf2
+    )
+    {
+        return false;
+    }
+
     static tmp<GeometricField<TypeR, PatchField, GeoMesh>> New
     (
         const tmp<GeometricField<Type1, PatchField, GeoMesh>>& tf1,
@@ -231,6 +262,17 @@ struct reuseTmpTmpGeometricField
     TypeR, Type1, Type12, TypeR, PatchField, GeoMesh
 >
 {
+    //- Second input has return type: possibly reusable
+    static bool reusable
+    (
+        const tmp<GeometricField<Type1, PatchField, GeoMesh>>& tf1,
+        const tmp<GeometricField<TypeR, PatchField, GeoMesh>>& tf2
+    )
+    {
+        return Detail::reusable(tf2);
+    }
+
+    //- Second input has return type
     static tmp<GeometricField<TypeR, PatchField, GeoMesh>> New
     (
         const tmp<GeometricField<Type1, PatchField, GeoMesh>>& tf1,
@@ -272,6 +314,17 @@ struct reuseTmpTmpGeometricField
     TypeR, TypeR, TypeR, Type2, PatchField, GeoMesh
 >
 {
+    //- First input has return type: possibly reusable
+    static bool reusable
+    (
+        const tmp<GeometricField<TypeR, PatchField, GeoMesh>>& tf1,
+        const tmp<GeometricField<Type2, PatchField, GeoMesh>>& tf2
+    )
+    {
+        return Detail::reusable(tf2);
+    }
+
+    //- First input has return type
     static tmp<GeometricField<TypeR, PatchField, GeoMesh>> New
     (
         const tmp<GeometricField<TypeR, PatchField, GeoMesh>>& tf1,
@@ -307,6 +360,16 @@ struct reuseTmpTmpGeometricField
     TypeR, TypeR, TypeR, TypeR, PatchField, GeoMesh
 >
 {
+    //- Both inputs have return type: possibly reusable
+    static bool reusable
+    (
+        const tmp<GeometricField<TypeR, PatchField, GeoMesh>>& tf1,
+        const tmp<GeometricField<TypeR, PatchField, GeoMesh>>& tf2
+    )
+    {
+        return (Detail::reusable(tf1) || Detail::reusable(tf2));
+    }
+
     static tmp<GeometricField<TypeR, PatchField, GeoMesh>> New
     (
         const tmp<GeometricField<TypeR, PatchField, GeoMesh>>& tf1,