diff --git a/applications/test/expressionTemplates/Make/files b/applications/test/expressionTemplates/Make/files
new file mode 100644
index 0000000000..9d866b294b
--- /dev/null
+++ b/applications/test/expressionTemplates/Make/files
@@ -0,0 +1,5 @@
+Test-expressionTemplates-volFields.C
+/* Test-expressionTemplates-pointFields.C */
+/* Test-expressionTemplates-Fields.C */
+
+EXE = $(FOAM_USER_APPBIN)/Test-expressionTemplates
diff --git a/applications/test/expressionTemplates/Make/options b/applications/test/expressionTemplates/Make/options
new file mode 100644
index 0000000000..ede08caa60
--- /dev/null
+++ b/applications/test/expressionTemplates/Make/options
@@ -0,0 +1,11 @@
+EXE_INC = \
+    -O0 \
+    -I$(FOAM_SRC)/meshTools/lnInclude \
+    -I$(FOAM_SRC)/finiteVolume/lnInclude \
+
+EXE_LIBS = \
+/*    -rpath /Volumes/case_sensitive/OpenFOAM/work/feature-et/platforms */ \
+/*    -rpath /Volumes/case_sensitive/OpenFOAM/work/feature-et/platforms/darwin64GccDPInt32Opt/lib */ \
+/*    -rpath /Volumes/case_sensitive/OpenFOAM/work/feature-et/platforms/darwin64GccDPInt32Opt/lib/$(FOAM_MPI)  */\
+    -L/Volumes/case_sensitive/OpenFOAM/work/feature-et/platforms/darwin64GccDPInt32Opt/lib \
+    -lfiniteVolume
diff --git a/applications/test/expressionTemplates/Test-expressionTemplates-pointFields.C b/applications/test/expressionTemplates/Test-expressionTemplates-pointFields.C
new file mode 100644
index 0000000000..47b0017c11
--- /dev/null
+++ b/applications/test/expressionTemplates/Test-expressionTemplates-pointFields.C
@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 M. Janssens
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-expressionTemplates
+
+\*---------------------------------------------------------------------------*/
+
+#include "Time.H"
+#include "argList.H"
+#include "polyMesh.H"
+#include "pointMesh.H"
+#include "pointFields.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createPolyMesh.H"
+
+    {
+        scalarField vals0({1.0, 2.0, 3.0});
+        const Expression::ListConstRefWrap<scalar> wvals0(vals0);
+
+        scalarField vals1;
+        Expression::ListRefWrap<scalar> wvals1(vals1.size(), vals1);
+        wvals1 = wvals0;
+        return 0;
+    }
+
+
+
+
+    const pointMesh& pMesh = pointMesh::New(mesh);
+
+    // Field, dimensionedScalar as List expression
+    {
+        scalarField vals({1.0, 2.0, 3.0});
+        dimensionedScalar d(dimless, 4.0);
+
+        scalarField result(d.expr(vals.size())*vals.expr());
+
+        Pout<< "result:" << result << endl;
+    }
+
+    Info<< "Reading field p\n" << endl;
+    pointScalarField p
+    (
+        IOobject
+        (
+            "pointDisplacement",
+            runTime.timeName(),
+            pMesh.thisDb(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pMesh
+    );
+
+    pointScalarField result
+    (
+        IOobject
+        (
+            "result",
+            runTime.timeName(),
+            pMesh.thisDb(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            IOobject::NO_REGISTER
+        ),
+        pMesh,
+        dimensionedScalar("zero", Foam::sqr(p().dimensions()), 0)
+    );
+
+    // Construct value + dimensions with correct sizing
+    const auto oneField(dimensionedScalar(dimless, 1.0).expr(p));
+
+    // Make expression
+    auto expression = oneField * p.expr() + p.expr();
+
+    // Combine expressions
+    auto newExpression = sqr(expression);
+
+    // Assign values
+    result = newExpression;
+
+    DebugVar(result);
+
+
+    // Construct from expression
+    pointScalarField result2("result2", pMesh, sqrt(p.expr()));
+
+    DebugVar(result2);
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/expressionTemplates/Test-expressionTemplates-volFields.C b/applications/test/expressionTemplates/Test-expressionTemplates-volFields.C
new file mode 100644
index 0000000000..bd5b2a3219
--- /dev/null
+++ b/applications/test/expressionTemplates/Test-expressionTemplates-volFields.C
@@ -0,0 +1,451 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    volPointInterpolationTest
+
+\*---------------------------------------------------------------------------*/
+
+#include "Time.H"
+#include "argList.H"
+#include "fvMesh.H"
+#include "ListExpression.H"
+#include "GeometricFieldExpression.H"
+#include "fvCFD.H"
+#include "fvMatrixExpression.H"
+#include <ratio>
+#include <chrono>
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+tmp<volScalarField> someFunction(const volScalarField& fld)
+{
+    return fld*1.0001;
+}
+
+
+template<class Type>
+void fusedGaussFvmLaplacian
+(
+    fvMatrix<Type>& fvm,
+    const surfaceInterpolationScheme<scalar>& interpGammaScheme,
+    const fv::snGradScheme<Type>& snGradScheme,
+    const GeometricField<scalar, fvPatchField, volMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    // Replacement for gaussLaplacianScheme::fvmLaplacian with scalar gamma
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> surfaceType;
+
+    const auto& mesh = vf.mesh();
+
+    // Expression for weights
+    const auto weights = interpGammaScheme.weights(gamma).expr();
+
+    // Expression for gamma_face * magSf
+    const auto gammaMagSf =
+        Expression::interpolate(gamma.expr(), weights, mesh)
+      * mesh.magSf().expr();
+
+    // Expression for deltaCoeffs
+    const auto deltaCoeffs = snGradScheme.deltaCoeffs(vf).expr();
+
+    // Construct matrix
+    Expression::fvmLaplacianUncorrected(fvm, gammaMagSf, deltaCoeffs);
+
+    if (snGradScheme.corrected())
+    {
+        // Wrap correction
+        const auto corr(snGradScheme.correction(vf).expr());
+        const auto V = mesh.V().expr();
+
+        if (mesh.fluxRequired(vf.name()))
+        {
+            fvm.faceFluxCorrectionPtr() = std::make_unique<surfaceType>
+            (
+                IOobject
+                (
+                    "faceFluxCorr",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    IOobject::NO_REGISTER
+                ),
+                mesh,
+                gamma.dimensions()
+               *mesh.magSf().dimensions()
+               *corr.data().dimensions()
+            );
+            auto& faceFluxCorr = *fvm.faceFluxCorrectionPtr();
+            faceFluxCorr = gammaMagSf*corr;
+
+            fvm.source() =
+                fvm.source().expr()
+              - (
+                    V * fvc::div
+                    (
+                        faceFluxCorr
+                    )().primitiveField().expr()
+                );
+        }
+        else
+        {
+            // Temporary field
+            surfaceType faceFluxCorr
+            (
+                IOobject
+                (
+                    "faceFluxCorr",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    IOobject::NO_REGISTER
+                ),
+                mesh,
+                gamma.dimensions()
+               *mesh.magSf().dimensions()
+               *corr.data().dimensions()
+            );
+            faceFluxCorr = gammaMagSf*corr;
+
+            fvm.source() =
+                fvm.source().expr()
+              - (
+                    V * fvc::div
+                    (
+                        faceFluxCorr
+                    )().primitiveField().expr()
+                );
+        }
+    }
+}
+
+
+using namespace std::chrono;
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh.thisDb(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+    
+    {
+        volScalarField p2("p2", p);
+        DebugVar(p2.boundaryFieldRef());
+    
+        for (auto& pf : p.boundaryFieldRef())
+        {
+            scalarField newVals(pf.size());
+            forAll(pf, i)
+            {
+                newVals[i] = scalar(i);
+            }
+            pf == newVals;
+        }
+
+        //std::vector<const scalarField*> ptrs;
+        // List<const scalarField*> ptrs;
+        // for (const auto& pp : p.boundaryField())
+        // {
+        //     ptrs.push_back(&pp);
+        // }
+        DebugVar(sizeof(long));
+        DebugVar(p.boundaryField());
+        Expression::ListsConstRefWrap<scalarField> expr(p.boundaryField());
+
+        const auto twoA = expr + expr;
+        Expression::ListsRefWrap<scalarField> expr2(p2.boundaryFieldRef());
+        Pout<< "**before assignment twoA:" << twoA.size()
+            << " expr2:" << expr2.size() << endl;
+        expr2 = twoA;
+        Pout<< "**after assignment twoA:" << twoA.size()
+            << " expr2:" << expr2.size() << endl;
+        DebugVar(p2.boundaryField());
+        // forAll(expr, i)
+        // {
+        //     Pout<< "i:" << i
+        //         //<< " expr:" << expr[i]
+        //         //<< " twoA:" << twoA[i]
+        //         << " expr2:" << expr2[i]
+        //         << endl;
+        // }
+        return 0;
+    }
+
+    {
+        //DebugVar(linearInterpolate(p));
+        //auto tweights = linear<scalar>(mesh).weights(p);
+        //DebugVar(tweights);
+
+        surfaceScalarField result
+        (
+            IOobject
+            (
+                "result",
+                runTime.timeName(),
+                mesh.thisDb(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            dimensionedScalar(p.dimensions(), 0)
+        );
+        //result = Expression::interpolate
+        //(
+        //    p.expr(),
+        //    tweights().expr(),
+        //    mesh
+        //);
+
+        result = Expression::linearInterpolate(p.expr(), mesh);
+
+
+        DebugVar(result);
+
+        return 0;
+    }
+
+
+
+    volScalarField p2
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh.thisDb(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE,
+            IOobject::NO_REGISTER
+        ),
+        mesh
+    );
+
+    // Expresions of volFields
+    {
+        volScalarField result
+        (
+            IOobject
+            (
+                "result",
+                runTime.timeName(),
+                mesh.thisDb(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            dimensionedScalar(p.dimensions(), 0)
+        );
+
+        {
+            Pout<< "No expression templates:" << endl;
+            const high_resolution_clock::time_point t1 =
+                high_resolution_clock::now();
+            result = p + p2;
+            const high_resolution_clock::time_point t2 =
+                high_resolution_clock::now();
+            const duration<double> time_span = t2 - t1;
+            Pout<< "Operation time:" << time_span.count() << endl;
+        }
+        {
+            Pout<< "With expression templates:" << endl;
+            const high_resolution_clock::time_point t1 =
+                high_resolution_clock::now();
+            result = p.expr() + p2.expr();
+            const high_resolution_clock::time_point t2 =
+                high_resolution_clock::now();
+            const duration<double> time_span = t2 - t1;
+            Pout<< "Operation time:" << time_span.count() << endl;
+        }
+
+//        const auto oldDimensions = p.dimensions();
+//        p.dimensions().reset(dimless);
+//        p2.dimensions().reset(dimless);
+//        result.dimensions().reset(dimless);
+//        {
+//
+//            Pout<< "Complex expression : No expression templates:" << endl;
+//            const high_resolution_clock::time_point t1 =
+//                high_resolution_clock::now();
+//            result = cos(p + 0.5*sqrt(p2-sin(p)));
+//            const high_resolution_clock::time_point t2 =
+//                high_resolution_clock::now();
+//            const duration<double> time_span = t2 - t1;
+//            Pout<< "Operation time:" << time_span.count() << endl;
+//        }
+//        {
+//            Pout<< "Complex expression : With expression templates:" << endl;
+//            const high_resolution_clock::time_point t1 =
+//                high_resolution_clock::now();
+//            const auto zeroDotFive
+//            (
+//                dimensionedScalar(dimless, 0.5).expr(p)
+//            );
+//            result = cos(p.expr() + zeroDotFive*sqrt(p2.expr()-sin(p.expr())));
+//            const high_resolution_clock::time_point t2 =
+//                high_resolution_clock::now();
+//            const duration<double> time_span = t2 - t1;
+//            Pout<< "Operation time:" << time_span.count() << endl;
+//        }
+//        p.dimensions().reset(oldDimensions);
+//        p2.dimensions().reset(oldDimensions);
+//        result.dimensions().reset(oldDimensions);
+
+        return 0;
+//        auto expression = someFunction(p).expr() + someFunction(p).expr();
+//        result = expression;
+//        DebugVar(result);
+    }
+    {
+        volScalarField result
+        (
+            IOobject
+            (
+                "result",
+                runTime.timeName(),
+                mesh.thisDb(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            dimensionedScalar(p.dimensions(), 0)
+        );
+        auto expression = someFunction(p).expr() + someFunction(p).expr();
+        result = expression;
+        DebugVar(result);
+    }
+
+    // Expresions of volFields
+    {
+        volScalarField result
+        (
+            "result",
+            mesh,
+            sqr(p.expr() + p.expr())
+        );
+        DebugVar(result);
+    }
+    {
+        // Fill p with some values
+        forAll(p, celli)
+        {
+            p[celli] = celli;
+        }
+        p.correctBoundaryConditions();
+
+        // Interpolate to surface field
+        surfaceScalarField result
+        (
+            "result",
+            mesh,
+            Expression::interpolate    //<volExpr, surfaceExpr>
+            (
+                p.expr(),
+                mesh.surfaceInterpolation::weights().expr(),
+                mesh
+            )
+        );
+        DebugVar(result);
+    }
+    {
+        // For testing as a replacement of laplacian weights
+        const volScalarField gamma
+        (
+            IOobject
+            (
+                "gamma",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            dimensionedScalar(dimless, 1.0)
+        );
+
+        fvMatrix<scalar> fvm
+        (
+            p,
+            gamma.dimensions()*mesh.magSf().dimensions()*p.dimensions()
+        );
+
+        const linear<scalar> interpGammaScheme(mesh);
+        const fv::correctedSnGrad<scalar> snGradScheme(mesh);
+
+        fusedGaussFvmLaplacian
+        (
+            fvm,
+            interpGammaScheme,
+            snGradScheme,
+            gamma,
+            p
+        );
+
+        DebugVar(fvm.source());
+    }
+
+
+
+    // Expressions of fvMatrix
+    {
+        tmp<fvMatrix<scalar>> tm0(fvm::laplacian(p));
+        const fvMatrix<scalar>& m0 = tm0();
+        DebugVar(m0.dimensions());
+
+        tmp<fvMatrix<scalar>> tm1(fvm::laplacian(p));
+        const fvMatrix<scalar>& m1 = tm1();
+        DebugVar(m1.dimensions());
+
+        fvMatrix<scalar> m2(p, m0.expr() + m1.expr());
+        DebugVar(m2.dimensions());
+    }
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/etc/config.sh/settings b/etc/config.sh/settings
index a615345631..ae5e714e64 100644
--- a/etc/config.sh/settings
+++ b/etc/config.sh/settings
@@ -81,8 +81,8 @@ Darwin)
         case "$WM_COMPILE_CONTROL" in
         (*version=*) ;;
         (*)
-            WM_COMPILER=Clang
-            echo "openfoam (darwin): using clang instead of gcc" 1>&2
+            #WM_COMPILER=Clang
+            #echo "openfoam (darwin): using clang instead of gcc" 1>&2
             ;;
         esac
     fi
diff --git a/src/Allwmake b/src/Allwmake
index 04f82ffe97..8e87d20311 100755
--- a/src/Allwmake
+++ b/src/Allwmake
@@ -55,8 +55,6 @@ wmake $targetType fileFormats
 wmake $targetType surfMesh
 wmake $targetType meshTools
 
-#------------------------------------------------------------------------------
-
 wmake $targetType finiteArea
 wmake $targetType finiteVolume
 wmake $targetType fused/finiteVolume
diff --git a/src/OpenFOAM/containers/Lists/List/List.H b/src/OpenFOAM/containers/Lists/List/List.H
index eaaaaea6a9..fd2b176638 100644
--- a/src/OpenFOAM/containers/Lists/List/List.H
+++ b/src/OpenFOAM/containers/Lists/List/List.H
@@ -45,6 +45,7 @@ SourceFiles
 
 #include "autoPtr.H"
 #include "UList.H"
+#include "SLListFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,7 +66,6 @@ typedef List<bool> boolList;        //!< A List of bools
 typedef List<char> charList;        //!< A List of chars
 typedef List<label> labelList;      //!< A List of labels
 
-
 /*---------------------------------------------------------------------------*\
                            Class List Declaration
 \*---------------------------------------------------------------------------*/
@@ -395,6 +395,37 @@ public:
 
         //- Alias for resize()
         void setSize(label n, const T& val) { this->resize(n, val); }
+
+
+    // Expression templates
+
+        //- Construct from value expression
+        template<typename E>
+        explicit List
+        (
+            const Expression::ListExpression
+            <
+                E,
+                typename E::value_type
+            >& expr
+        )
+        {
+            expr.evaluate(*this);
+        }
+
+        //- Assign values from expression
+        template<typename E>
+        void operator=
+        (
+            const Expression::ListExpression
+            <
+                E,
+                typename E::value_type
+            >& expr
+        )
+        {
+            expr.evaluate(*this);
+        }
 };
 
 
diff --git a/src/OpenFOAM/containers/Lists/List/UList.H b/src/OpenFOAM/containers/Lists/List/UList.H
index 7614575dc8..8982132e9c 100644
--- a/src/OpenFOAM/containers/Lists/List/UList.H
+++ b/src/OpenFOAM/containers/Lists/List/UList.H
@@ -80,6 +80,11 @@ typedef UList<bool> boolUList;      //!< A UList of bools
 typedef UList<char> charUList;      //!< A UList of chars
 typedef UList<label> labelUList;    //!< A UList of labels
 
+namespace Expression
+{
+template<class E, typename ValType> class ListExpression;
+template<class T> class ListConstRefWrap;
+}
 
 /*---------------------------------------------------------------------------*\
                            Class UList Declaration
@@ -697,6 +702,39 @@ public:
         {
             return this->contains(val, pos);
         }
+
+
+    // Expression templates
+
+        //- Have unique tag
+        using is_List = void;
+
+        //- Wrap value as expression
+        auto expr() const
+        {
+            return Expression::ListConstRefWrap<T>
+            (
+                reinterpret_cast<const List<T>&>(*this)
+            );
+        }
+
+        //- Assign values from expression
+        template<typename E>
+        void operator=
+        (
+            const Expression::ListExpression
+            <
+                E,
+                typename E::value_type
+            >& expr
+        )
+        {
+            // bypass expr.evaluate to avoid resize ...
+            for (label i = 0; i < expr.size(); ++i)
+            {
+                operator[](i) = expr[i];
+            }
+        }
 };
 
 
diff --git a/src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H b/src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H
index 9f840e2698..f166c6eae6 100644
--- a/src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H
+++ b/src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H
@@ -60,6 +60,12 @@ template<class Type> class dimensioned;
 template<class Type>
 Istream& operator>>(Istream& is, dimensioned<Type>& dt);
 
+namespace Expression
+{
+template<class T> class UniformListWrap;
+template<class GeoField> class UniformGeometricFieldWrap;
+}
+
 /*---------------------------------------------------------------------------*\
                            Class dimensioned Declaration
 \*---------------------------------------------------------------------------*/
@@ -449,6 +455,28 @@ public:
         {
             return getOrAddToDict(name, dict, deflt);
         }
+
+
+    // Expression templates
+
+        //- Wrap value as constant-value List expression
+        auto expr(const label size) const
+        {
+            return Expression::UniformListWrap<Type>(size, value_);
+        }
+
+        //- Wrap value as constant-value GeometricField expression. Supplied
+        //  fields is used for sizing only and needs to be valid at time of
+        //  evaluation.
+        template<class GeoField>
+        auto expr(const GeoField& fld) const
+        {
+            return Expression::UniformGeometricFieldWrap<GeoField>
+            (
+                fld,
+                *this
+            );
+        }
 };
 
 
diff --git a/src/OpenFOAM/expressionTemplates/GenericExpression.H b/src/OpenFOAM/expressionTemplates/GenericExpression.H
new file mode 100644
index 0000000000..24d83bc631
--- /dev/null
+++ b/src/OpenFOAM/expressionTemplates/GenericExpression.H
@@ -0,0 +1,252 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 M. Janssens
+-------------------------------------------------------------------------------
+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/>.
+
+Namespace
+    Foam::GenericExpression
+
+Description
+    Expression templates.
+
+SourceFiles
+    GenericExpression.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_GenericExpression_H
+#define Foam_GenericExpression_H
+
+#include <cassert>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace Expression
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class GenericExpression Declarations
+\*---------------------------------------------------------------------------*/
+
+template<typename E>
+class GenericExpression
+{
+public:
+    static constexpr bool is_leaf = false;
+
+    auto evaluate() const
+    {
+        // Delegation to the actual expression type. This avoids dynamic
+        // polymorphism (a.k.a. virtual functions in C++)
+        return static_cast<E const&>(*this).evaluate();
+    }
+};
+
+
+// Macros
+// ~~~~~~
+
+#define EXPRESSION_FUNCTION1(BaseClass, Func, BaseFunc, OpFunc)                \
+template<typename E1>                                                          \
+class OpFunc                                                                   \
+:                                                                              \
+    public BaseClass<OpFunc<E1>>                                               \
+{                                                                              \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    OpFunc(const E1& u)                                                        \
+    :                                                                          \
+        u_(u)                                                                  \
+    {}                                                                         \
+                                                                               \
+    auto evaluate() const                                                      \
+    {                                                                          \
+        return BaseFunc(u_.evaluate());                                        \
+    }                                                                          \
+};                                                                             \
+template<typename E1>                                                          \
+OpFunc<E1> Func                                                                \
+(                                                                              \
+    const BaseClass<E1>& u                                                     \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1>(static_cast<const E1&>(u));                              \
+}
+
+
+#define EXPRESSION_FUNCTION2(BaseClass, Func, BaseFunc, OpFunc)                \
+template<typename E1, typename E2>                                             \
+class OpFunc                                                                   \
+:                                                                              \
+    public BaseClass<OpFunc<E1, E2>>                                           \
+{                                                                              \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+    typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    OpFunc(const E1& u, const E2& v)                                           \
+    :                                                                          \
+     u_(u), v_(v)                                                              \
+    {}                                                                         \
+    auto evaluate() const                                                      \
+    {                                                                          \
+        return BaseFunc(u_.evaluate(), v_.evaluate());                         \
+    }                                                                          \
+};                                                                             \
+template<typename E1, typename E2>                                             \
+OpFunc<E1, E2> Func                                                            \
+(                                                                              \
+    const BaseClass<E1>& u,                                                    \
+    const BaseClass<E2>& v                                                     \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1, E2>                                                      \
+    (                                                                          \
+        static_cast<const E1&>(u),                                             \
+        static_cast<const E2&>(v)                                              \
+    );                                                                         \
+}
+
+
+#define EXPRESSION_OPERATOR(BaseClass, Op, BaseOp, OpFunc)                     \
+template<typename E1, typename E2>                                             \
+class OpFunc                                                                   \
+:                                                                              \
+    public BaseClass                                                           \
+    <                                                                          \
+        OpFunc<E1, E2>                                                         \
+    >                                                                          \
+{                                                                              \
+    /* cref if leaf, copy otherwise */                                         \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+    typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    OpFunc(const E1& u, const E2& v)                                           \
+    :                                                                          \
+     u_(u), v_(v)                                                              \
+    {}                                                                         \
+    auto evaluate() const                                                      \
+    {                                                                          \
+        return u_.evaluate() BaseOp v_.evaluate();                             \
+    }                                                                          \
+};                                                                             \
+template<typename E1, typename E2>                                             \
+OpFunc<E1, E2>                                                                 \
+operator Op                                                                    \
+(                                                                              \
+    BaseClass<E1> const& u,                                                    \
+    BaseClass<E2> const& v                                                     \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1, E2>                                                      \
+    (                                                                          \
+        static_cast<const E1&>(u),                                             \
+        static_cast<const E2&>(v)                                              \
+    );                                                                         \
+}
+
+
+// Do '-' separately until we work out macro expansion...
+template<typename E1>
+class g_negate
+:
+    public GenericExpression
+    <
+        g_negate<E1>
+    >
+{
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;
+
+public:
+    static constexpr bool is_leaf = false;
+
+    g_negate(const E1& u)
+    :
+        u_(u)
+    {}
+
+    auto evaluate() const
+    {
+        return -u_.evaluate();
+    }
+};
+template<typename E1>
+g_negate<E1> operator-
+(
+    const GenericExpression<E1>& u
+)
+{
+    return g_negate<E1>(static_cast<const E1&>(u));
+}
+
+
+EXPRESSION_FUNCTION1(GenericExpression, sqr, Foam::sqr, g_sqr)
+EXPRESSION_FUNCTION1(GenericExpression, sqrt, Foam::sqrt, g_sqrt)
+EXPRESSION_FUNCTION1(GenericExpression, magSqr, Foam::magSqr, g_magSqr)
+EXPRESSION_FUNCTION1(GenericExpression, symm, Foam::symm, g_symm)
+EXPRESSION_FUNCTION1(GenericExpression, pow2, Foam::pow2, g_pow2)
+EXPRESSION_FUNCTION1(GenericExpression, pow3, Foam::pow3, g_pow3)
+EXPRESSION_FUNCTION1(GenericExpression, pow4, Foam::pow4, g_pow4)
+EXPRESSION_FUNCTION1(GenericExpression, operator~, Foam::operator~, g_compl)
+
+//TBD. Parse problem
+//EXPRESSION_FUNCTION1(GenericExpression, operator!, Foam::operator!, g_not)
+
+#undef EXPRESSION_FUNCTION1
+
+EXPRESSION_FUNCTION2(GenericExpression, operator+, Foam::operator+, g_add);
+EXPRESSION_FUNCTION2(GenericExpression, operator-, Foam::operator-, g_subtract);
+EXPRESSION_FUNCTION2(GenericExpression, operator*, Foam::operator*, g_multiply);
+EXPRESSION_FUNCTION2(GenericExpression, operator/, Foam::operator/, g_divide);
+
+#undef EXPRESSION_FUNCTION2
+
+EXPRESSION_OPERATOR(GenericExpression, ||, ||, g_or)
+EXPRESSION_OPERATOR(GenericExpression, &&, &&, g_and)
+EXPRESSION_OPERATOR(GenericExpression, &, &, g_bitand)
+EXPRESSION_OPERATOR(GenericExpression, |, |, g_bitor)
+EXPRESSION_OPERATOR(GenericExpression, ^, ^, g_xor)
+
+#undef EXPRESSION_OPERATOR
+
+} // End namespace Expression
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/expressionTemplates/GeometricFieldExpression.H b/src/OpenFOAM/expressionTemplates/GeometricFieldExpression.H
new file mode 100644
index 0000000000..b748c33f24
--- /dev/null
+++ b/src/OpenFOAM/expressionTemplates/GeometricFieldExpression.H
@@ -0,0 +1,1879 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2024,2025 M. Janssens
+    Copyright (C) 2025 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/>.
+
+Namespace
+    Foam::Expression
+
+Description
+    Expression templates for GeometricFields
+
+SourceFiles
+    GeometricFieldExpression.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_GeometricFieldExpression_H
+#define Foam_GeometricFieldExpression_H
+
+#include "ListExpression.H"
+#include "pointMesh.H"
+//#include "fixedValuePointPatchField.H"  // no assignable yet. tbd.
+//#include <boost/core/demangle.hpp>
+//#include <typeinfo>
+//#include <iostream>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace Expression
+{
+
+/*---------------------------------------------------------------------------*\
+                   Namespace Expression Declarations
+\*---------------------------------------------------------------------------*/
+
+
+// Expression of GeometricField
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+template
+<
+    typename E,
+    typename IntExpr,
+    typename UncoupledPatchExpr,
+    typename CoupledPatchExpr,
+    typename Type
+>
+class GeometricFieldExpression
+{
+protected:
+
+    const dimensionSet dimensions_;
+    const orientedType oriented_;
+
+public:
+    static constexpr bool is_leaf = false;
+
+    GeometricFieldExpression
+    (
+        const dimensionSet& dimensions,
+        const orientedType oriented
+    )
+    :
+        dimensions_(dimensions),
+        oriented_(oriented)
+    {
+        //Pout<< "GeometricFieldExpression :" << nl
+        //    << "    E          : "
+        //    << boost::core::demangle(typeid(E).name())
+        //    << nl
+        //    << "    IntExpr   : "
+        //    << boost::core::demangle(typeid(IntExpr).name()) << nl
+        //    << "    dimensions : " << dimensions_ << nl
+        //    << "    oriented   : " << oriented_ << nl
+        //    << endl;
+    }
+
+    // TBD. Do we still need this? Is faster than always going through
+    // internalField() though.
+    Type operator[](const label i) const
+    {
+        // Delegation to the actual expression type. This avoids dynamic
+        // polymorphism (a.k.a. virtual functions in C++)
+        return static_cast<E const&>(*this)[i];
+    }
+    auto size() const noexcept { return static_cast<const E&>(*this).size(); }
+
+    const dimensionSet& dimensions() const noexcept
+    {
+        return dimensions_;
+    }
+
+    const orientedType& oriented() const noexcept
+    {
+        return oriented_;
+    }
+
+    IntExpr internalField() const
+    {
+        return static_cast<const E&>(*this).internalField();
+    }
+
+    IntExpr internalField()
+    {
+        return static_cast<E&>(*this).internalField();
+    }
+
+    UncoupledPatchExpr patchField(const label i)
+    {
+        return static_cast<E&>(*this).patchField(i);
+    }
+
+    UncoupledPatchExpr patchField(const label i) const
+    {
+        return static_cast<const E&>(*this).patchField(i);
+    }
+
+    CoupledPatchExpr coupledPatchField(const label i)
+    {
+        return static_cast<E&>(*this).coupledPatchField(i);
+    }
+
+    CoupledPatchExpr coupledPatchField(const label i) const
+    {
+        return static_cast<const E&>(*this).coupledPatchField(i);
+    }
+
+    template<class Op>
+    auto access(const Op& cop, const label i) const
+    {
+        return static_cast<const E&>(*this).access(cop, i);
+    }
+
+    //- Helper to evaluate a GeometricField
+    template<class GeoField>
+    GeoField& evaluate(GeoField& fld, const bool force = false) const
+    {
+        fld.dimensions() = this->dimensions();
+        fld.oriented() = this->oriented();
+
+        assert(fld.size() == this->size());
+        using ActualType = typename IntExpr::value_type;
+        ListRefWrap<ActualType> wfld(fld.internalFieldRef(), internalField());
+
+        // Do boundary field
+        auto& bfld = fld.boundaryFieldRef();
+        const label n = bfld.size();
+        for (label i = 0; i < n; ++i)
+        {
+            auto& pfld = bfld[i];
+            if (force || pfld.assignable())
+            {
+                if (pfld.coupled())
+                {
+                    // Get patch expression
+                    const auto patchExpr = this->coupledPatchField(i);
+
+                    // Note: check can be removed once assignable also on
+                    //       pointPatchFields. TBD.
+                    if constexpr
+                    (
+                        std::is_same_v<pointMesh, typename GeoField::Mesh>
+                    )
+                    {
+                        const auto* vp = isA<Field<Type>>(pfld);
+                        if (vp)
+                        {
+                            patchExpr.evaluate(const_cast<Field<Type>&>(*vp));
+                        }
+                    }
+                    else
+                    {
+                        patchExpr.evaluate(pfld);
+                    }
+                }
+                else
+                {
+                    // Get patch expression
+                    const auto patchExpr = this->patchField(i);
+
+                    // Note: check can be removed once assignable also on
+                    //       pointPatchFields. TBD.
+                    if constexpr
+                    (
+                        std::is_same_v<pointMesh, typename GeoField::Mesh>
+                    )
+                    {
+                        const auto* vp = isA<Field<Type>>(pfld);
+                        if (vp)
+                        {
+                            patchExpr.evaluate(const_cast<Field<Type>&>(*vp));
+                        }
+                    }
+                    else
+                    {
+                        patchExpr.evaluate(pfld);
+                    }
+                }
+            }
+        }
+        fld.correctLocalBoundaryConditions();
+        return fld;
+    }
+};
+
+
+template<class GeoField>
+class GeometricFieldRefWrap
+:
+    public GeometricFieldExpression
+    <
+        GeometricFieldRefWrap<GeoField>,
+        ListRefWrap<typename GeoField::value_type>,
+        ListRefWrap<typename GeoField::value_type>,
+        ListTmpWrap<Field<typename GeoField::value_type>>,
+        typename GeoField::value_type
+    >
+{
+public:
+
+    static constexpr bool is_leaf = false;  //true;
+
+    //- The GeometricField type
+    typedef GeoField this_type;
+
+    //- Type to return for internal field
+    typedef typename GeoField::value_type value_type;
+
+    //- Type to return for patchField
+    typedef ListRefWrap<value_type> IntExpr;
+    typedef ListRefWrap<value_type> UncoupledPatchExpr;
+    typedef ListTmpWrap<Field<value_type>> CoupledPatchExpr;
+
+
+private:
+
+    this_type& elems_;
+
+
+public:
+
+    //- Copy construct
+    GeometricFieldRefWrap(this_type& elems)
+    :
+        GeometricFieldExpression
+        <
+            GeometricFieldRefWrap<GeoField>,
+            IntExpr,
+            UncoupledPatchExpr,
+            CoupledPatchExpr,
+            value_type
+        >
+        (
+            elems.dimensions(),
+            elems.oriented()
+        ),
+        elems_(elems)
+    {}
+
+    //- Move construct
+    GeometricFieldRefWrap(this_type&& elems)
+    :
+        GeometricFieldExpression
+        <
+            GeometricFieldRefWrap<GeoField>,
+            IntExpr,
+            UncoupledPatchExpr,
+            CoupledPatchExpr,
+            value_type
+        >
+        (
+            elems.dimensions(),
+            elems.oriented()
+        ),
+        elems_(elems)
+    {}
+
+    // Construct from ListExpression, forcing its evaluation.
+    template<typename E>
+    GeometricFieldRefWrap
+    (
+        this_type& elems,
+        const GeometricFieldExpression
+        <
+            E,
+            typename E::IntExpr,
+            typename E::UncoupledPatchExpr,
+            typename E::CoupledPatchExpr,
+            typename E::value_type
+        >& expr
+    )
+    :
+        GeometricFieldExpression
+        <
+            GeometricFieldRefWrap<GeoField>,
+            IntExpr,
+            UncoupledPatchExpr,
+            CoupledPatchExpr,
+            value_type
+        >
+        (
+            elems.dimensions(),
+            elems.oriented()
+        ),
+        elems_(elems)
+    {
+        expr.evaluate(elems_);
+    }
+
+    //- Assignment
+    template<typename E>
+    void operator=
+    (
+        const GeometricFieldExpression
+        <
+            E,
+            typename E::IntExpr,
+            typename E::UncoupledPatchExpr,
+            typename E::CoupledPatchExpr,
+            typename E::value_type
+        >& expr
+    )
+    {
+        expr.evaluate(elems_);
+    }
+
+    //- Evaluate and return as GeoField. Rename to evaluate to make it clear
+    //- it takes time? Or leave as indexing for convenience?
+    template<typename E>
+    this_type& evaluate
+    (
+        const GeometricFieldExpression
+        <
+            E,
+            typename E::IntExpr,
+            typename E::UncoupledPatchExpr,
+            typename E::CoupledPatchExpr,
+            typename E::value_type
+        >& expr
+    )
+    {
+        return expr.evaluate(elems_);
+    }
+
+    value_type operator[](const label i) const
+    {
+        return elems_[i];
+    }
+
+    value_type& operator[](const label  i)
+    {
+        return elems_[i];
+    }
+
+    auto size() const noexcept
+    {
+        return elems_.size();
+    }
+
+    IntExpr internalField()
+    {
+        auto& fld = elems_.internalFieldRef();
+        return IntExpr(fld.size(), fld);
+    }
+
+    IntExpr internalField() const
+    {
+        auto& fld = elems_.internalField();
+        return IntExpr(fld.size(), fld);
+    }
+
+    UncoupledPatchExpr patchField(const label i)
+    {
+        auto& fld = elems_.boundaryFieldRef()[i];
+        if constexpr
+        (
+            std::is_same_v<pointMesh, typename GeoField::Mesh>
+        )
+        {
+            const auto* vp = isA<Field<value_type>>(fld);
+            if (vp)
+            {
+                auto& v = *vp;
+                return UncoupledPatchExpr
+                (
+                    v.size(),
+                    const_cast<Field<value_type>&>(v)
+                );
+            }
+            else
+            {
+                return UncoupledPatchExpr(0, List<value_type>::null());
+            }
+        }
+        else
+        {
+            return UncoupledPatchExpr(fld.size(), fld);
+        }
+    }
+
+    UncoupledPatchExpr patchField(const label i) const
+    {
+        const auto& fld = elems_.boundaryField()[i];
+        if constexpr
+        (
+            std::is_same_v<pointMesh, typename GeoField::Mesh>
+        )
+        {
+            const auto* vp = isA<Field<value_type>>(fld);
+            if (vp)
+            {
+                return UncoupledPatchExpr(vp->size(), *vp);
+            }
+            else
+            {
+                return UncoupledPatchExpr(0, List<value_type>::null());
+            }
+        }
+        else
+        {
+            return UncoupledPatchExpr(fld.size(), fld);
+        }
+    }
+
+    CoupledPatchExpr coupledPatchField(const label i)
+    {
+        auto& fld = elems_.boundaryFieldRef()[i];
+        if constexpr
+        (
+            std::is_same_v<pointMesh, typename GeoField::Mesh>
+        )
+        {
+            const auto* vp = isA<Field<value_type>>(fld);
+            if (vp)
+            {
+                auto& v = *vp;
+                return CoupledPatchExpr
+                (
+                    //v.size(),
+                    const_cast<Field<value_type>&>(v)
+                );
+            }
+            else
+            {
+                return CoupledPatchExpr(List<value_type>::null());
+            }
+        }
+        else
+        {
+            return CoupledPatchExpr(fld);
+        }
+    }
+
+    CoupledPatchExpr coupledPatchField(const label i) const
+    {
+        const auto& fld = elems_.boundaryField()[i];
+        if constexpr
+        (
+            std::is_same_v<pointMesh, typename GeoField::Mesh>
+        )
+        {
+            const auto* vp = isA<Field<value_type>>(fld);
+            if (vp)
+            {
+                //return CoupledPatchExpr(vp->size(), *vp);
+                return CoupledPatchExpr(*vp);
+            }
+            else
+            {
+                //return CoupledPatchExpr(0, List<value_type>::null());
+                return CoupledPatchExpr(List<value_type>::null());
+            }
+        }
+        else
+        {
+            //return UncoupledPatchExpr(fld.size(), fld);
+            return CoupledPatchExpr(fld);
+        }
+    }
+
+    template<class AccessOp>
+    auto access(const AccessOp& cop, const label i) const
+    {
+        return cop(elems_, i);
+    }
+};
+
+
+template<class GeoField>
+class GeometricFieldConstRefWrap
+:
+    public GeometricFieldExpression
+    <
+        GeometricFieldConstRefWrap<GeoField>,
+        ListConstRefWrap<typename GeoField::value_type>,
+        ListConstRefWrap<typename GeoField::value_type>,
+        ListConstTmpWrap<Field<typename GeoField::value_type>>,
+        typename GeoField::value_type
+    >
+{
+public:
+
+    static constexpr bool is_leaf = false;  //true;
+
+    //- The GeometricField type
+    typedef GeoField this_type;
+
+    //- Type to return for internal field
+    typedef typename GeoField::value_type value_type;
+
+    //- Type to return for patchField
+    typedef ListConstRefWrap<value_type> IntExpr;
+    typedef ListConstRefWrap<value_type> UncoupledPatchExpr;
+    typedef ListConstTmpWrap<Field<value_type>> CoupledPatchExpr;
+
+
+private:
+
+    const this_type& elems_;
+
+
+public:
+
+    // Construct from components
+    GeometricFieldConstRefWrap(const this_type& elems)
+    :
+        GeometricFieldExpression
+        <
+            GeometricFieldConstRefWrap<GeoField>,
+            IntExpr,
+            UncoupledPatchExpr,
+            CoupledPatchExpr,
+            value_type
+        >
+        (
+            elems.dimensions(),
+            elems.oriented()
+        ),
+        elems_(elems)
+    {}
+
+    // return the underlying data
+    const auto& data() const
+    {
+        return elems_;
+    }
+
+    value_type operator[](const label i) const
+    {
+        return elems_[i];
+    }
+
+    auto size() const noexcept
+    {
+        return elems_.size();
+    }
+
+    IntExpr internalField() const
+    {
+        return elems_.internalField();
+    }
+
+    UncoupledPatchExpr patchField(const label i) const
+    {
+        //return elems_.boundaryField()[i];
+        const auto& fld = elems_.boundaryField()[i];
+        if constexpr
+        (
+            std::is_same_v<pointMesh, typename GeoField::Mesh>
+        )
+        {
+            const auto* vp = isA<Field<value_type>>(fld);
+            if (vp)
+            {
+                return UncoupledPatchExpr(*vp);
+            }
+            else
+            {
+                return UncoupledPatchExpr(List<value_type>::null());
+            }
+        }
+        else
+        {
+            return UncoupledPatchExpr(fld);
+        }
+    }
+
+    CoupledPatchExpr coupledPatchField(const label i) const
+    {
+        const auto& fld = elems_.boundaryField()[i];
+        if constexpr
+        (
+            std::is_same_v<pointMesh, typename GeoField::Mesh>
+        )
+        {
+            const auto* vp = isA<Field<value_type>>(fld);
+            if (vp)
+            {
+                return CoupledPatchExpr(*vp);
+            }
+            else
+            {
+                return CoupledPatchExpr(tmp<Field<value_type>>(nullptr));
+            }
+        }
+        else
+        {
+            return CoupledPatchExpr(fld);
+        }
+    }
+
+    template<class AccessOp>
+    auto access(const AccessOp& cop, const label i) const
+    {
+        return cop(elems_, i);
+    }
+};
+
+
+template<class GeoField>
+class GeometricFieldConstTmpWrap
+:
+    public GeometricFieldExpression
+    <
+        GeometricFieldConstTmpWrap<GeoField>,
+        ListConstRefWrap<typename GeoField::value_type>,
+        ListConstRefWrap<typename GeoField::value_type>,
+        ListConstTmpWrap<Field<typename GeoField::value_type>>,
+        typename GeoField::value_type
+    >
+{
+public:
+
+    //- Have expressions use copy to maintain tmp refCount
+    static constexpr bool is_leaf = false;
+
+    //- The GeometricField type
+    typedef GeoField this_type;
+
+    //- Type to return for internal field
+    typedef typename GeoField::value_type value_type;
+
+    //- Type to return for patchField
+    typedef ListConstRefWrap<value_type> IntExpr;
+    typedef ListConstRefWrap<value_type> UncoupledPatchExpr;
+    typedef ListConstTmpWrap<Field<value_type>> CoupledPatchExpr;
+
+
+private:
+
+    const tmp<this_type> elems_;
+
+
+public:
+
+    // Construct from components
+    GeometricFieldConstTmpWrap(const tmp<this_type>& elems)
+    :
+        GeometricFieldExpression
+        <
+            GeometricFieldConstTmpWrap<GeoField>,
+            IntExpr,
+            UncoupledPatchExpr,
+            CoupledPatchExpr,
+            value_type
+        >
+        (
+            elems().dimensions(),
+            elems().oriented()
+        ),
+        elems_(elems, true)
+    {
+        //Pout<< "GeometricFieldConstTmpWrap :" << nl
+        //    << "    this_type  : "
+        //    << boost::core::demangle(typeid(this_type).name())
+        //    << nl
+        //    << "    PatchExpr  : "
+        //    << boost::core::demangle(typeid(IntExpr).name()) << nl
+        //    << endl;
+    }
+
+    GeometricFieldConstTmpWrap
+    (
+        const GeometricFieldConstTmpWrap<GeoField>& w
+    )
+    :
+        GeometricFieldExpression
+        <
+            GeometricFieldConstTmpWrap<GeoField>,
+            IntExpr,
+            UncoupledPatchExpr,
+            CoupledPatchExpr,
+            value_type
+        >
+        (
+            w.dimensions(),
+            w.oriented()
+        ),
+        //elems_(std::move(w.elems_))
+        elems_(w.elems_)
+    {
+        //Pout<< "GeometricFieldConstTmpWrap copy :" << nl
+        //    << "    this_type  : "
+        //    << boost::core::demangle(typeid(this_type).name())
+        //    << nl
+        //    << "    PatchExpr  : "
+        //    << boost::core::demangle(typeid(IntExpr).name()) << nl
+        //    << endl;
+    }
+
+    ~GeometricFieldConstTmpWrap()
+    {
+        //Pout<< "~GeometricFieldConstTmpWrap destructor :" << nl
+        //    << "    this_type  : "
+        //    << boost::core::demangle(typeid(this_type).name())
+        //    << nl
+        //    << endl;
+    }
+
+    // return the underlying data
+    const auto& data() const
+    {
+        return elems_();
+    }
+
+    value_type operator[](const label i) const
+    {
+        // Redirection for every access. Slow?
+        return elems_()[i];
+    }
+
+    auto size() const noexcept
+    {
+        return elems_().size();
+    }
+
+    IntExpr internalField() const
+    {
+        return elems_().internalField();
+    }
+
+    UncoupledPatchExpr patchField(const label i) const
+    {
+        const auto& fld = elems_().boundaryField()[i];
+        if constexpr
+        (
+            std::is_same_v<pointMesh, typename GeoField::Mesh>
+        )
+        {
+            const auto* vp = isA<Field<value_type>>(fld);
+            if (vp)
+            {
+                return UncoupledPatchExpr(*vp);
+            }
+            else
+            {
+                return UncoupledPatchExpr(List<value_type>::null());
+            }
+        }
+        else
+        {
+            return UncoupledPatchExpr(fld);
+        }
+
+    }
+
+    CoupledPatchExpr coupledPatchField(const label i) const
+    {
+        const auto& fld = elems_().boundaryField()[i];
+        if constexpr
+        (
+            std::is_same_v<pointMesh, typename GeoField::Mesh>
+        )
+        {
+            const auto* vp = isA<Field<value_type>>(fld);
+            if (vp)
+            {
+                return CoupledPatchExpr(*vp);
+            }
+            else
+            {
+                return CoupledPatchExpr(List<value_type>::null());
+            }
+        }
+        else
+        {
+            return CoupledPatchExpr(fld);
+        }
+    }
+
+    template<class AccessOp>
+    auto access(const AccessOp& cop, const label i) const
+    {
+        return cop(elems_(), i);
+    }
+};
+
+
+// Expressions on GeometricFields
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+// Do '+' ourselves - just to show. Others done with macros below.
+template<typename E1, typename E2>
+class GF_add
+:
+    public GeometricFieldExpression
+    <
+        GF_add<E1, E2>,
+        List_add
+        <
+            typename E1::IntExpr,
+            typename E2::IntExpr
+        >,
+        List_add
+        <
+            typename E1::UncoupledPatchExpr,
+            typename E2::UncoupledPatchExpr
+        >,
+        List_add
+        <
+            typename E1::CoupledPatchExpr,
+            typename E2::CoupledPatchExpr
+        >,
+        typename E1::value_type
+    >
+{
+    // cref if leaf, copy otherwise
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;
+    typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_;
+
+public:
+    static constexpr bool is_leaf = false;
+
+    //- Type to return for internal field
+    typedef typename E1::value_type value_type;
+
+    //- Type to return for patchField
+    typedef List_add
+    <
+        typename E1::IntExpr,
+        typename E2::IntExpr
+    > IntExpr;
+    typedef List_add
+    <
+        typename E1::UncoupledPatchExpr,
+        typename E2::UncoupledPatchExpr
+    > UncoupledPatchExpr;
+    typedef List_add
+    <
+        typename E1::CoupledPatchExpr,
+        typename E2::CoupledPatchExpr
+    > CoupledPatchExpr;
+
+    GF_add(const E1& u, const E2& v)
+    :
+        GeometricFieldExpression
+        <
+            GF_add<E1, E2>,
+            IntExpr,
+            UncoupledPatchExpr,
+            CoupledPatchExpr,
+            value_type
+        >
+        (
+            u.dimensions() + v.dimensions(),
+            u.oriented()    // + v.oriented()
+        ),
+        u_(u),
+        v_(v)
+    {
+        assert(u.size() == v.size());
+        static_assert
+        (
+            std::is_same_v
+            <
+                typename E1::value_type,
+                typename E2::value_type
+            > == true
+        );
+    }
+    auto operator[](const label i) const
+    {
+        return u_[i] + v_[i];
+    }
+    auto size() const noexcept { return u_.size(); }
+
+    IntExpr internalField() const
+    {
+        return IntExpr(u_.internalField(), v_.internalField());
+    }
+
+    UncoupledPatchExpr patchField(const label i) const
+    {
+        return UncoupledPatchExpr(u_.patchField(i), v_.patchField(i));
+    }
+
+    CoupledPatchExpr coupledPatchField(const label i) const
+    {
+        return CoupledPatchExpr
+        (
+            u_.coupledPatchField(i),
+            v_.coupledPatchField(i)
+        );
+    }
+
+    template<class AccessOp>
+    auto access(const AccessOp& cop, const label i) const
+    {
+        return UncoupledPatchExpr(u_.access(cop, i), v_.access(cop, i));
+    }
+};
+template<typename E1, typename E2>
+GF_add<E1, E2>
+operator+
+(
+    GeometricFieldExpression
+    <
+        E1,
+        typename E1::IntExpr,
+        typename E1::UncoupledPatchExpr,
+        typename E1::CoupledPatchExpr,
+        typename E1::value_type
+    > const& u,
+    GeometricFieldExpression
+    <
+        E2,
+        typename E2::IntExpr,
+        typename E2::UncoupledPatchExpr,
+        typename E2::CoupledPatchExpr,
+        typename E2::value_type
+    > const& v
+)
+{
+    return GF_add<E1, E2>
+    (
+        static_cast<const E1&>(u),
+        static_cast<const E2&>(v)
+    );
+}
+
+
+#define EXPRESSION_GF_FUNCTION1_DIMLESS(Func, BaseFunc, WrapType, OpFunc)      \
+template<typename E1>                                                          \
+class OpFunc                                                                   \
+:                                                                              \
+    public GeometricFieldExpression                                            \
+    <                                                                          \
+        OpFunc<E1>,                                                            \
+        WrapType<typename E1::IntExpr>,                                        \
+        WrapType<typename E1::UncoupledPatchExpr>,                             \
+        WrapType<typename E1::CoupledPatchExpr>,                               \
+        typename E1::value_type                                                \
+    >                                                                          \
+{                                                                              \
+    /* cref if leaf, copy otherwise */                                         \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    /* Type to return for internal field */                                    \
+    typedef typename E1::value_type value_type;                                \
+                                                                               \
+    /* Type to return for patchField */                                        \
+    typedef WrapType<typename E1::IntExpr> IntExpr;                            \
+    typedef WrapType<typename E1::UncoupledPatchExpr> UncoupledPatchExpr;      \
+    typedef WrapType<typename E1::CoupledPatchExpr> CoupledPatchExpr;          \
+                                                                               \
+    OpFunc(const E1& u)                                                        \
+    :                                                                          \
+        GeometricFieldExpression                                               \
+        <                                                                      \
+            OpFunc<E1>,                                                        \
+            IntExpr,                                                           \
+            UncoupledPatchExpr,                                                \
+            CoupledPatchExpr,                                                  \
+            value_type                                                         \
+        >                                                                      \
+        (                                                                      \
+            u.dimensions(),                                                    \
+            u.oriented()                                                       \
+        ),                                                                     \
+        u_(u)                                                                  \
+    {}                                                                         \
+    auto operator[](const label i) const                                       \
+    {                                                                          \
+        return BaseFunc(u_[i]);                                                \
+    }                                                                          \
+    auto size() const noexcept { return u_.size(); }                           \
+                                                                               \
+    IntExpr internalField() const                                              \
+    {                                                                          \
+        return IntExpr(u_.internalField());                                    \
+    }                                                                          \
+                                                                               \
+    UncoupledPatchExpr patchField(const label i) const                         \
+    {                                                                          \
+        return UncoupledPatchExpr(u_.patchField(i));                           \
+    }                                                                          \
+                                                                               \
+    CoupledPatchExpr coupledPatchField(const label i) const                    \
+    {                                                                          \
+        return CoupledPatchExpr(u_.coupledPatchField(i));                      \
+    }                                                                          \
+                                                                               \
+    template<class AccessOpOp>                                                 \
+    auto access(const AccessOpOp& cop, const label i) const                    \
+    {                                                                          \
+        return UncoupledPatchExpr(u_.access(cop, i));                          \
+    }                                                                          \
+};                                                                             \
+template<typename E1>                                                          \
+OpFunc<E1> Func                                                                \
+(                                                                              \
+    GeometricFieldExpression                                                   \
+    <                                                                          \
+        E1,                                                                    \
+        typename E1::IntExpr,                                                  \
+        typename E1::UncoupledPatchExpr,                                       \
+        typename E1::CoupledPatchExpr,                                         \
+        typename E1::value_type                                                \
+    > const& u                                                                 \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1>(static_cast<const E1&>(u));                              \
+}
+
+
+#define EXPRESSION_GF_FUNCTION1(Func, BaseFunc, WrapType, OpFunc)              \
+template<typename E1>                                                          \
+class OpFunc                                                                   \
+:                                                                              \
+    public GeometricFieldExpression                                            \
+    <                                                                          \
+        OpFunc<E1>,                                                            \
+        WrapType<typename E1::IntExpr>,                                        \
+        WrapType<typename E1::UncoupledPatchExpr>,                             \
+        WrapType<typename E1::CoupledPatchExpr>,                               \
+        typename E1::value_type                                                \
+    >                                                                          \
+{                                                                              \
+    /* cref if leaf, copy otherwise */                                         \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    /* Type to return for internal field */                                    \
+    typedef typename E1::value_type value_type;                                \
+                                                                               \
+    /* Type to return for patchField */                                        \
+    typedef WrapType<typename E1::IntExpr> IntExpr;                            \
+    typedef WrapType<typename E1::UncoupledPatchExpr> UncoupledPatchExpr;      \
+    typedef WrapType<typename E1::CoupledPatchExpr> CoupledPatchExpr;          \
+                                                                               \
+    OpFunc(const E1& u)                                                        \
+    :                                                                          \
+        GeometricFieldExpression                                               \
+        <                                                                      \
+            OpFunc<E1>,                                                        \
+            IntExpr,                                                           \
+            UncoupledPatchExpr,                                                \
+            CoupledPatchExpr,                                                  \
+            value_type                                                         \
+        >                                                                      \
+        (                                                                      \
+            BaseFunc(u.dimensions()),                                          \
+            BaseFunc(u.oriented())                                             \
+        ),                                                                     \
+        u_(u)                                                                  \
+    {}                                                                         \
+    auto operator[](const label i) const                                       \
+    {                                                                          \
+        return BaseFunc(u_[i]);                                                \
+    }                                                                          \
+    auto size() const noexcept { return u_.size(); }                           \
+                                                                               \
+    IntExpr internalField() const                                              \
+    {                                                                          \
+        return IntExpr(u_.internalField());                                    \
+    }                                                                          \
+                                                                               \
+    UncoupledPatchExpr patchField(const label i) const                         \
+    {                                                                          \
+        return UncoupledPatchExpr(u_.patchField(i));                           \
+    }                                                                          \
+                                                                               \
+    CoupledPatchExpr coupledPatchField(const label i) const                    \
+    {                                                                          \
+        return CoupledPatchExpr(u_.coupledPatchField(i));                      \
+    }                                                                          \
+                                                                               \
+    template<class AccessOpOp>                                                 \
+    auto access(const AccessOpOp& cop, const label i) const                    \
+    {                                                                          \
+        return UncoupledPatchExpr(u_.access(cop, i));                          \
+    }                                                                          \
+};                                                                             \
+template<typename E1>                                                          \
+OpFunc<E1> Func                                                                \
+(                                                                              \
+    GeometricFieldExpression                                                   \
+    <                                                                          \
+        E1,                                                                    \
+        typename E1::IntExpr,                                                  \
+        typename E1::UncoupledPatchExpr,                                       \
+        typename E1::CoupledPatchExpr,                                         \
+        typename E1::value_type                                                \
+    > const& u                                                                 \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1>(static_cast<const E1&>(u));                              \
+}
+
+
+//// Do negation separately - problem with macro expansion. Tbd.
+//template<typename E1>
+//class GeometricFieldNegate
+//:
+//    public GeometricFieldExpression
+//    <
+//        GeometricFieldNegate<E1>,
+//        List_negate<typename E1::IntExpr>,
+//        List_negate<typename E1::UncoupledPatchExpr>,
+//        List_negate<typename E1::CoupledPatchExpr>,
+//        typename E1::value_type
+//    >
+//{
+//    /* cref if leaf, copy otherwise */ 
+//    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;
+//
+//public:
+//    static constexpr bool is_leaf = false;
+//
+//    /* Type to return for internal field */
+//    typedef typename E1::value_type value_type;
+//
+//    /* Type to return for patchField */
+//    typedef List_negate<typename E1::IntExpr> IntExpr;
+//    typedef List_negate<typename E1::UncoupledPatchExpr> UncoupledPatchExpr;
+//    typedef List_negate<typename E1::CoupledPatchExpr> CoupledPatchExpr;
+//
+//    GeometricFieldNegate(const E1& u)
+//    :
+//        GeometricFieldExpression
+//        <
+//            GeometricFieldNegate<E1>,
+//            IntExpr,
+//            UncoupledPatchExpr,
+//            CoupledPatchExpr,
+//            value_type
+//        >
+//        (
+//            -u.dimensions(),
+//            -u.oriented()
+//        ),
+//        u_(u)
+//    {}
+//    auto operator[](const label i) const
+//    {
+//        return -u_[i];
+//    }
+//    label size() const { return u_.size(); }
+//
+//    IntExpr internalField() const
+//    {
+//        return PatchExpr(u_.internalField());
+//    }
+//
+//    UncoupledPatchExpr patchField(const label i) const
+//    {
+//        return UncoupledPatchExpr(u_.patchField(i));
+//    }
+//};
+//template<typename E1>
+//GeometricFieldNegate<E1> operator-
+//(
+//    GeometricFieldExpression
+//    <
+//        E1,
+//        typename E1::IntExpr,
+//        typename E1::UncoupledPatchExpr,
+//        typename E1::CoupledPatchExpr,
+//        typename E1::value_type
+//    > const& u
+//)
+//{
+//    return GeometricFieldNegate<E1>(static_cast<const E1&>(u));
+//}
+
+#define EXPRESSION_GF_OPERATOR(Op, WrapType, OpFunc)                           \
+template<typename E1, typename E2>                                             \
+class OpFunc                                                                   \
+:                                                                              \
+    public GeometricFieldExpression                                            \
+    <                                                                          \
+        OpFunc<E1, E2>,                                                        \
+        WrapType                                                               \
+        <                                                                      \
+            typename E1::IntExpr,                                              \
+            typename E2::IntExpr                                               \
+        >,                                                                     \
+        WrapType                                                               \
+        <                                                                      \
+            typename E1::UncoupledPatchExpr,                                   \
+            typename E2::UncoupledPatchExpr                                    \
+        >,                                                                     \
+        WrapType                                                               \
+        <                                                                      \
+            typename E1::CoupledPatchExpr,                                     \
+            typename E2::CoupledPatchExpr                                      \
+        >,                                                                     \
+        typename E1::value_type                                                \
+    >                                                                          \
+{                                                                              \
+    /* cref if leaf, copy otherwise */                                         \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+    typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    /* Type to return for internal field */                                    \
+    typedef typename E1::value_type value_type;                                \
+                                                                               \
+    /* Type to return for patchField */                                        \
+    typedef WrapType                                                           \
+    <                                                                          \
+        typename E1::IntExpr,                                                  \
+        typename E2::IntExpr                                                   \
+    > IntExpr;                                                                 \
+    typedef WrapType                                                           \
+    <                                                                          \
+        typename E1::UncoupledPatchExpr,                                       \
+        typename E2::UncoupledPatchExpr                                        \
+    > UncoupledPatchExpr;                                                      \
+    typedef WrapType                                                           \
+    <                                                                          \
+        typename E1::CoupledPatchExpr,                                         \
+        typename E2::CoupledPatchExpr                                          \
+    > CoupledPatchExpr;                                                        \
+                                                                               \
+    OpFunc(const E1& u, const E2& v)                                           \
+    :                                                                          \
+        GeometricFieldExpression                                               \
+        <                                                                      \
+            OpFunc<E1, E2>,                                                    \
+            IntExpr,                                                           \
+            UncoupledPatchExpr,                                                \
+            CoupledPatchExpr,                                                  \
+            typename E1::value_type                                            \
+        >                                                                      \
+        (                                                                      \
+            u.dimensions() Op v.dimensions(),                                  \
+            u.oriented() Op v.oriented()                                       \
+        ),                                                                     \
+        u_(u),                                                                 \
+        v_(v)                                                                  \
+    {                                                                          \
+      /*static_assert                                                          \
+        (                                                                      \
+            std::is_same_v                                                     \
+            <                                                                  \
+                typename E1::value_type,                                       \
+                typename E2::value_type                                        \
+            > == true                                                          \
+        );*/                                                                   \
+        assert(u.size() == v.size());                                          \
+    }                                                                          \
+    auto operator[](const label i) const                                       \
+    {                                                                          \
+        return u_[i] Op v_[i];                                                 \
+    }                                                                          \
+    auto size() const noexcept { return u_.size(); }                           \
+                                                                               \
+    IntExpr internalField() const                                              \
+    {                                                                          \
+        return IntExpr(u_.internalField(), v_.internalField());                \
+    }                                                                          \
+                                                                               \
+    UncoupledPatchExpr patchField(const label i) const                         \
+    {                                                                          \
+        return UncoupledPatchExpr(u_.patchField(i), v_.patchField(i));         \
+    }                                                                          \
+                                                                               \
+    CoupledPatchExpr coupledPatchField(const label i) const                    \
+    {                                                                          \
+        return CoupledPatchExpr                                                \
+        (                                                                      \
+            u_.coupledPatchField(i),                                           \
+            v_.coupledPatchField(i)                                            \
+        );                                                                     \
+    }                                                                          \
+                                                                               \
+    template<class AccessOp>                                                   \
+    auto access(const AccessOp& cop, const label i) const                      \
+    {                                                                          \
+        return UncoupledPatchExpr                                              \
+        (                                                                      \
+            u_.access(cop, i),                                                 \
+            v_.access(cop, i)                                                  \
+        );                                                                     \
+    }                                                                          \
+};                                                                             \
+template<typename E1, typename E2>                                             \
+OpFunc<E1, E2>                                                                 \
+operator Op                                                                    \
+(                                                                              \
+    GeometricFieldExpression                                                   \
+    <                                                                          \
+        E1,                                                                    \
+        typename E1::IntExpr,                                                  \
+        typename E1::UncoupledPatchExpr,                                       \
+        typename E1::CoupledPatchExpr,                                         \
+        typename E1::value_type                                                \
+    > const& u,                                                                \
+    GeometricFieldExpression                                                   \
+    <                                                                          \
+        E2,                                                                    \
+        typename E2::IntExpr,                                                  \
+        typename E2::UncoupledPatchExpr,                                       \
+        typename E2::CoupledPatchExpr,                                         \
+        typename E2::value_type                                                \
+    > const& v                                                                 \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1, E2>                                                      \
+    (                                                                          \
+        static_cast<const E1&>(u),                                             \
+        static_cast<const E2&>(v)                                              \
+    );                                                                         \
+}
+
+#define EXPRESSION_GF_FUNCTION2(Func, BaseFunc, WrapType, OpFunc)              \
+template<typename E1, typename E2>                                             \
+class OpFunc                                                                   \
+:                                                                              \
+    public GeometricFieldExpression                                            \
+    <                                                                          \
+        OpFunc<E1, E2>,                                                        \
+        WrapType                                                               \
+        <                                                                      \
+            typename E1::IntExpr,                                              \
+            typename E2::IntExpr                                               \
+        >,                                                                     \
+        WrapType                                                               \
+        <                                                                      \
+            typename E1::UncoupledPatchExpr,                                   \
+            typename E2::UncoupledPatchExpr                                    \
+        >,                                                                     \
+        WrapType                                                               \
+        <                                                                      \
+            typename E1::CoupledPatchExpr,                                     \
+            typename E2::CoupledPatchExpr                                      \
+        >,                                                                     \
+        typename E1::value_type                                                \
+    >                                                                          \
+{                                                                              \
+    /* cref if leaf, copy otherwise */                                         \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+    typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    /* Type to return for internal field */                                    \
+    typedef typename E1::value_type value_type;                                \
+                                                                               \
+    /* Type to return for patchField */                                        \
+    typedef WrapType                                                           \
+    <                                                                          \
+        typename E1::IntExpr,                                                  \
+        typename E2::IntExpr                                                   \
+    > IntExpr;                                                                 \
+    typedef WrapType                                                           \
+    <                                                                          \
+        typename E1::UncoupledPatchExpr,                                       \
+        typename E2::UncoupledPatchExpr                                        \
+    > UncoupledPatchExpr;                                                      \
+    typedef WrapType                                                           \
+    <                                                                          \
+        typename E1::CoupledPatchExpr,                                         \
+        typename E2::CoupledPatchExpr                                          \
+    > CoupledPatchExpr;                                                        \
+                                                                               \
+    OpFunc(const E1& u, const E2& v)                                           \
+    :                                                                          \
+        GeometricFieldExpression                                               \
+        <                                                                      \
+            OpFunc<E1, E2>,                                                    \
+            IntExpr,                                                           \
+            UncoupledPatchExpr,                                                \
+            CoupledPatchExpr,                                                  \
+            value_type                                                         \
+        >                                                                      \
+        (                                                                      \
+            BaseFunc(u.dimensions(), v.dimensions()),                          \
+            BaseFunc(u.oriented(), v.oriented())                               \
+        ),                                                                     \
+        u_(u),                                                                 \
+        v_(v)                                                                  \
+    {                                                                          \
+        /*static_assert                                                        \
+        (                                                                      \
+            std::is_same_v                                                     \
+            <                                                                  \
+                typename E1::value_type,                                       \
+                typename E2::value_type                                        \
+            > == true                                                          \
+        );*/                                                                   \
+    }                                                                          \
+    auto operator[](const label i) const                                       \
+    {                                                                          \
+        return BaseFunc(u_[i], v_[i]);                                         \
+    }                                                                          \
+    auto size() const noexcept {return u_.size(); }                            \
+                                                                               \
+    IntExpr internalField() const                                              \
+    {                                                                          \
+        return IntExpr(u_.internalField(), v_.internalField());                \
+    }                                                                          \
+                                                                               \
+    UncoupledPatchExpr patchField(const label i) const                         \
+    {                                                                          \
+        return UncoupledPatchExpr(u_.patchField(i), v_.patchField(i));         \
+    }                                                                          \
+                                                                               \
+    CoupledPatchExpr coupledPatchField(const label i) const                    \
+    {                                                                          \
+        return CoupledPatchExpr                                                \
+        (                                                                      \
+            u_.coupledPatchField(i),                                           \
+            v_.coupledPatchField(i)                                            \
+        );                                                                     \
+    }                                                                          \
+                                                                               \
+    template<class AccessOp>                                                   \
+    auto access(const AccessOp& cop, const label i) const                      \
+    {                                                                          \
+        return UncoupledPatchExpr                                              \
+        (                                                                      \
+            u_.access(cop, i),                                                 \
+            v_.access(cop, i)                                                  \
+        );                                                                     \
+    }                                                                          \
+};                                                                             \
+template<typename E1, typename E2>                                             \
+OpFunc<E1, E2> Func                                                            \
+(                                                                              \
+    GeometricFieldExpression                                                   \
+    <                                                                          \
+        E1,                                                                    \
+        typename E1::IntExpr,                                                  \
+        typename E1::UncoupledPatchExpr,                                       \
+        typename E1::CoupledPatchExpr,                                         \
+        typename E1::value_type                                                \
+    > const& u,                                                                \
+    GeometricFieldExpression                                                   \
+    <                                                                          \
+        E2,                                                                    \
+        typename E2::IntExpr,                                                  \
+        typename E2::UncoupledPatchExpr,                                       \
+        typename E2::CoupledPatchExpr,                                         \
+        typename E2::value_type                                                \
+    > const& v                                                                 \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1, E2>                                                      \
+    (                                                                          \
+        static_cast<const E1&>(u),                                             \
+        static_cast<const E2&>(v)                                              \
+    );                                                                         \
+}
+
+
+//- GeometricField expressions.
+
+// macro arguments:
+//  - 'sin' : name of new function (in Expression namespace)
+//  - '::sin' : per-element function to call
+//  - 'List_sin' : helper class to handle patchFields
+//  - 'GF_sin' : generated expression helper class
+EXPRESSION_GF_FUNCTION1_DIMLESS(sin, ::sin, List_sin, GF_sin)
+EXPRESSION_GF_FUNCTION1_DIMLESS(cos, ::cos, List_cos, GF_cos)
+EXPRESSION_GF_FUNCTION1_DIMLESS(tan, ::tan, List_tan, GF_tan)
+EXPRESSION_GF_FUNCTION1_DIMLESS(sinh, ::sinh, List_sinh, GF_sinh)
+EXPRESSION_GF_FUNCTION1_DIMLESS(cosh, ::cosh, List_cosh, GF_cosh)
+EXPRESSION_GF_FUNCTION1_DIMLESS(tanh, ::tanh, List_tanh, GF_tanh)
+
+#undef EXPRESSION_GF_FUNCTION1_DIMLESS
+
+EXPRESSION_GF_FUNCTION1(sqr, Foam::sqr, List_sqr, GF_sqr)
+EXPRESSION_GF_FUNCTION1(sqrt, Foam::sqrt, List_sqrt, GF_sqrt)
+EXPRESSION_GF_FUNCTION1(magSqr, Foam::magSqr, List_magSqr, GF_magSqr)
+
+EXPRESSION_GF_FUNCTION1(mag, Foam::mag, List_mag, GF_mag)
+EXPRESSION_GF_FUNCTION1(symm, Foam::symm, List_symm, GF_symm)
+EXPRESSION_GF_FUNCTION1(pow2, Foam::pow2, List_pow2, GF_pow2)
+EXPRESSION_GF_FUNCTION1(pow3, Foam::pow3, List_pow3, GF_pow3)
+EXPRESSION_GF_FUNCTION1(pow4, Foam::pow4, List_pow4, GF_pow4)
+
+#undef EXPRESSION_GF_FUNCTION1
+
+EXPRESSION_GF_FUNCTION2(min, Foam::min, List_min, GF_min)
+EXPRESSION_GF_FUNCTION2(max, Foam::max, List_max, GF_max)
+
+#undef EXPRESSION_GF_FUNCTION2
+
+EXPRESSION_GF_OPERATOR(-, List_subtract, GF_subtract)
+EXPRESSION_GF_OPERATOR(*, List_multiply, GF_multiply)
+EXPRESSION_GF_OPERATOR(/, List_divide, GF_divide)
+EXPRESSION_GF_OPERATOR(&, List_dot, GF_dot)
+
+#undef EXPRESSION_GF_OPERATOR
+
+
+
+// Expressions on constants
+// ~~~~~~~~~~~~~~~~~~~~~~~~
+
+template<class GeoField>
+class UniformGeometricFieldWrap
+:
+    public GeometricFieldExpression
+    <
+        UniformGeometricFieldWrap<GeoField>,
+        UniformListWrap<typename GeoField::value_type>,
+        UniformListWrap<typename GeoField::value_type>,
+        UniformListWrap<typename GeoField::value_type>,
+        typename GeoField::value_type
+    >
+{
+public:
+
+    static constexpr bool is_leaf = false;  //true;
+
+    //- Type to return for internal field
+    typedef typename GeoField::value_type value_type;
+
+    //- Type to return for patchField
+    typedef UniformListWrap<value_type> IntExpr;
+    typedef UniformListWrap<value_type> UncoupledPatchExpr;
+    typedef UniformListWrap<value_type> CoupledPatchExpr;
+
+
+private:
+
+    const GeoField& elems_;
+
+    const value_type val_;
+
+public:
+
+    // Construct from field (sizes only) and value
+    UniformGeometricFieldWrap(const GeoField& elems, const value_type val)
+    :
+        GeometricFieldExpression
+        <
+            UniformGeometricFieldWrap<GeoField>,
+            IntExpr,
+            UncoupledPatchExpr,
+            CoupledPatchExpr,
+            value_type
+        >
+        (
+            dimless,
+            orientedType()
+        ),
+        elems_(elems),
+        val_(val)
+    {}
+
+    // Construct from field (sizes only) and value, dimension
+    UniformGeometricFieldWrap
+    (
+        const GeoField& elems,
+        const dimensioned<value_type> val
+    )
+    :
+        GeometricFieldExpression
+        <
+            UniformGeometricFieldWrap<GeoField>,
+            IntExpr,
+            UncoupledPatchExpr,
+            CoupledPatchExpr,
+            value_type
+        >
+        (
+            val.dimensions(),
+            orientedType()
+        ),
+        elems_(elems),
+        val_(val.value())
+    {}
+
+    value_type operator[](const label i) const
+    {
+        return val_;
+    }
+
+    auto size() const noexcept
+    {
+        return elems_.size();
+    }
+
+    IntExpr internalField() const
+    {
+        return IntExpr(elems_.internalField().size(), val_);
+    }
+
+    UncoupledPatchExpr patchField(const label i) const
+    {
+        return UncoupledPatchExpr(elems_.boundaryField()[i].size(), val_);
+    }
+
+    CoupledPatchExpr coupledPatchField(const label i) const
+    {
+        return CoupledPatchExpr(elems_.boundaryField()[i].size(), val_);
+    }
+
+    template<class AccessOp>
+    auto access(const AccessOp& cop, const label i) const
+    {
+        // Size+value only
+        return patchField(i);
+    }
+};
+
+
+//- Fully self-contained constant field wrapper. Not needed?
+//template<class GeoField>
+//class UniformGeometricFieldWrap2
+//:
+//    public GeometricFieldExpression
+//    <
+//        UniformGeometricFieldWrap2<GeoField>,
+//        UniformListWrap<typename GeoField::value_type>,
+//        typename GeoField::value_type
+//    >
+//{
+//public:
+//
+//    static constexpr bool is_leaf = true;
+//
+//    //- Type to return for internal field
+//    typedef typename GeoField::value_type value_type;
+//
+//    //- Type to return for patchField
+//    typedef UniformListWrap<value_type> PatchExpr;
+//
+//
+//private:
+//
+//    const value_type val_;
+//
+//    const label nInternal_;
+//
+//    labelList patchSizes_;
+//
+//    //const dimensionSet& dimensions_;
+//
+//public:
+//
+//    // Construct from field (sizes only) and value
+//    UniformGeometricFieldWrap2(const GeoField& elems, const value_type val)
+//    :
+//        val_(val),
+//        nInternal_(elems.size()),
+//        patchSizes_(elems.boundaryField().size())
+//        //dimensions_(elems.dimensions()),
+//    {
+//        for (const auto& pfld : elems.boundaryField())
+//        {
+//            patchSizes_[pfld.patch().index()] = pfld.size();
+//        }
+//    }
+//
+//    // Construct from field (sizes only) and value, dimension
+//    UniformGeometricFieldWrap2
+//    (
+//        const GeoField& elems,
+//        const dimensioned<value_type> val
+//    )
+//    :
+//        val_(val.value()),
+//        nInternal_(elems.size()),
+//        patchSizes_(elems.boundaryField().size())
+//        //dimensions_(val.dimensions()),
+//    {
+//        for (const auto& pfld : elems.boundaryField())
+//        {
+//            patchSizes_[pfld.patch().index()] = pfld.size();
+//        }
+//    }
+//
+//    value_type operator[](const label i) const
+//    {
+//        return val_;
+//    }
+//
+//    auto size() const noexcept
+//    {
+//        return nInternal_;
+//    }
+//
+//    PatchExpr internalField() const
+//    {
+//        return PatchExpr(nInternal_, val_);
+//    }
+//
+//    PatchExpr patchField(const label i) const
+//    {
+//        return PatchExpr(patchSizes_[i], val_);
+//    }
+//};
+
+
+/*
+// Some fvm functions using expressions
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+// Supply a multiplier, either 1 or -1
+template<class Expr, class Expr2>
+void Su
+(
+    fvMatrix<typename Expr::value_type>& m,
+    const Expr2& mult,
+    const Expr& expression
+)
+{
+    //- Wrap of List as an expression
+    typedef ListConstRefWrap<typename Expr::value_type> expr;
+    //- Evaluator of an expression
+    typedef ListRefWrap<typename Expr::value_type> evaluator;
+
+    // Wrap mesh volume
+    const expr V(m.psi().mesh().V());
+
+    // Add expression to source
+    auto& s = m.source();
+    evaluator(s, expr(s) + mult*V*expression);
+}
+
+
+// Always add
+template<class Expr>
+void Su
+(
+    fvMatrix<typename Expr::value_type>& m,
+    const Expr& expression
+)
+{
+    //- Wrap of List as an expression
+    typedef ListConstRefWrap<typename Expr::value_type> expr;
+    //- Evaluator of an expression
+    typedef ListRefWrap<typename Expr::value_type> evaluator;
+
+    // Wrap mesh volume
+    const expr V(m.psi().mesh().V());
+
+    // Add expression to source
+    auto& s = m.source();
+    evaluator(s, expr(s) + V*expression);
+}
+
+
+template<class Expr>
+void rhs
+(
+    fvMatrix<typename Expr::value_type>& m,
+    const Expr& expression
+)
+{
+    Su(m, expression);
+}
+
+
+template<class Expr, class Expr2>
+void Sp
+(
+    fvMatrix<typename Expr::value_type>& m,
+    const Expr2& mult,
+    const Expr& expression
+)
+{
+    //- Wrap of List as an expression
+    typedef ListConstRefWrap<typename Expr::value_type> expr;
+    //- Evaluator of an expression
+    typedef ListRefWrap<typename Expr::value_type> evaluator;
+
+    // Wrap mesh volume
+    const expr V(m.psi().mesh().V());
+
+    // Add expression onto diag
+    auto& d = m.diag();
+    evaluator(d, expr(d) - mult*V*expression);
+}
+
+
+template<class Expr, class Expr2>
+void SuSp
+(
+    fvMatrix<typename Expr::value_type>& m,
+    const Expr2& mult,
+    const Expr& expression
+)
+{
+    //- Wrap of constant as a list expression
+    typedef UniformListWrap<scalar> constant;
+    //- Wrap of List as an expression
+    typedef ListConstRefWrap<typename Expr::value_type> expr;
+    //- Evaluator of an expression
+    typedef ListRefWrap<typename Expr::value_type> evaluator;
+
+    const constant constantZero(expression.size(), 0.0);
+
+    // Wrap mesh volume
+    const expr V(m.psi().mesh().V());
+
+    // Linearise expression
+    auto& diag = m.diag();
+    auto& source = m.source();
+
+    // diag() += domain*max(susp.field(), scalar(0));
+    evaluator
+    (
+        diag,
+        expr(diag) - mult*V*max(expression, constantZero)
+    );
+    // source() -= domain*min(susp.field(), scalar(0))*fld.primitiveField();
+    evaluator
+    (
+        source,
+        expr(source)
+      + (
+            mult
+           *V
+           *min(expression, constantZero)
+           *expr(m.psi().internalField())
+        )
+    );
+}
+*/
+
+
+} // End namespace Expression
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/expressionTemplates/ListExpression.H b/src/OpenFOAM/expressionTemplates/ListExpression.H
new file mode 100644
index 0000000000..07f372c7e3
--- /dev/null
+++ b/src/OpenFOAM/expressionTemplates/ListExpression.H
@@ -0,0 +1,2411 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2024,2025 M. Janssens
+    Copyright (C) 2025 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/>.
+
+Namespace
+    Foam::ListExpression
+
+Description
+    Expression templates for List
+
+SourceFiles
+    ListExpression.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_ListExpression_H
+#define Foam_ListExpression_H
+
+#include <cassert>
+#include <execution>
+#include "List.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace Expression
+{
+
+/*---------------------------------------------------------------------------*\
+                       Namespace Expression Declarations
+\*---------------------------------------------------------------------------*/
+
+template<typename E, typename ValType>
+class ListExpression
+{
+public:
+    static constexpr bool is_leaf = false;
+
+    ValType operator[](const label i) const
+    {
+        // Delegation to the actual expression type. This avoids dynamic
+        // polymorphism (a.k.a. virtual functions in C++)
+        return static_cast<E const&>(*this)[i];
+    }
+    auto size() const noexcept { return static_cast<E const&>(*this).size(); }
+
+    //- Helper: assign to passed in list
+    template<class Container>
+    Container& evaluate(Container& lst) const
+    {
+        lst.resize_nocopy(this->size());
+
+#ifdef _OPENMP
+        int loopLen = this->size();
+        _Pragma("omp target teams distribute parallel for if(loopLen > 1000)")
+        for(int i = 0; i < loopLen; ++i)
+        {
+            auto src = static_cast<E const&>(*this).cbegin() + i;
+            lst[i] = *src;
+        }
+#else
+    std::copy
+    (
+        #if defined(__APPLE__) && defined(__clang__)
+        // Currently not working (2025-10-06) with par_unseq. Works with
+        // sequential.
+        #else
+        std::execution::par_unseq,
+        #endif
+        static_cast<E const&>(*this).cbegin(),
+        static_cast<E const&>(*this).cend(),
+        lst.begin()
+    );
+#endif
+        return lst;
+    }
+};
+
+
+// Wrap of List
+// ~~~~~~~~~~~~
+
+template<class T>
+class ListWrap
+:
+    public ListExpression<ListWrap<T>, T>
+{
+    List<T> elems_;
+
+public:
+    static constexpr bool is_leaf = true;
+
+    typedef T value_type;
+
+    //- Copy construct
+    ListWrap(const ListWrap<T>& elems)
+    :
+        elems_(elems.data())
+    {}
+
+    //- Move construct
+    ListWrap(ListWrap<T>&& elems)
+    :
+        elems_(std::move(elems.data()))
+    {}
+
+    // Construct from initializer list 
+    ListWrap(std::initializer_list<T> elems)
+    :
+        elems_(elems)
+    {}
+
+    ListWrap(const UList<T>& elems)
+    :
+        elems_(elems)
+    {}
+
+    // Construct from ListExpression, forcing its evaluation.
+    template<typename E>
+    ListWrap(const ListExpression<E, typename E::value_type>& expr)
+    :
+        elems_(expr.size())
+    {
+        expr.evaluate(elems_);
+    }
+
+    // Access underlying data
+    const auto& data() const
+    { 
+        return elems_; 
+    }
+
+    // Access underlying data
+    auto& data()
+    { 
+        return elems_; 
+    }
+
+    value_type operator[](const label i) const
+    {
+        return elems_[i];
+    }
+
+    value_type& operator[](const label i)
+    {
+        return elems_[i];
+    }
+
+    auto size() const noexcept
+    {
+        return elems_.size();
+    }
+
+    typedef typename List<T>::const_iterator const_iterator;
+    auto cbegin() const
+    {
+        return elems_.cbegin();
+    }
+    auto cend() const
+    {
+        return elems_.cend();
+    }
+
+    typedef typename List<T>::iterator iterator;
+    auto begin()
+    {
+        return elems_.begin();
+    }
+    auto end()
+    {
+        return elems_.end();
+    }
+};
+
+
+// Wrap of non-const reference to List
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+template<class T>
+class ListRefWrap
+:
+    public ListExpression<ListRefWrap<T>, T>
+{
+    List<T>& elems_;
+
+public:
+    static constexpr bool is_leaf = false;  //true;
+
+    typedef T value_type;
+
+    //- Copy construct
+    ListRefWrap(ListRefWrap<T>& w)
+    :
+        elems_(w.elems_)
+    {}
+
+    //- Move construct
+    ListRefWrap(ListRefWrap<T>&& elems)
+    :
+        elems_(elems)
+    {}
+
+    //- Construct from List and expected size. This is to enforce that
+    //- our size() member function returns the correct size. Use the
+    //- construct-from-expression otherwise.
+    ListRefWrap(const label size, List<T>& elems)
+    :
+        elems_(elems)
+    {
+        elems_.resize_nocopy(size);
+    }
+
+    // Construct from ListExpression, forcing its evaluation.
+    template<typename E>
+    ListRefWrap
+    (
+        List<T>& elems,
+        const ListExpression<E, typename E::value_type>& expr
+    )
+    :
+        elems_(elems)
+    {
+        expr.evaluate(elems_);
+    }
+
+    //- Assignment
+    template<typename E>
+    void operator=(const ListExpression<E, typename E::value_type>& expr)
+    {
+        expr.evaluate(elems_);
+    }
+
+    // Access underlying data
+    const auto& data() const
+    { 
+        return elems_; 
+    }
+
+    // Access underlying data
+    auto& data()
+    { 
+        return elems_; 
+    }
+
+    value_type operator[](const label i) const
+    {
+        return elems_[i];
+    }
+
+    value_type& operator[](const label i)
+    {
+        return elems_[i];
+    }
+
+    auto size() const noexcept
+    {
+        return elems_.size();
+    }
+
+    typedef typename List<T>::const_iterator const_iterator;
+    auto cbegin() const
+    {
+        return elems_.cbegin();
+    }
+    auto cend() const
+    {
+        return elems_.cend();
+    }
+
+    typedef typename List<T>::iterator iterator;
+    auto begin()
+    {
+        return elems_.begin();
+    }
+    auto end()
+    {
+        return elems_.end();
+    }
+};
+
+
+// Wrap of const reference to List
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+template<class T>
+class ListConstRefWrap
+:
+    public ListExpression<ListConstRefWrap<T>, T>
+{
+    const List<T>& elems_;
+
+public:
+
+    // ! Store as copy (since holds reference)
+    static constexpr bool is_leaf = false;
+
+    typedef T value_type;
+
+    // construct from components
+    ListConstRefWrap(const List<T>& elems)
+    :
+        elems_(elems)
+    {}
+
+    // returns the underlying data
+    const auto& data() const
+    { 
+        return elems_; 
+    }
+
+    value_type operator[](const label i) const
+    {
+        return elems_[i];
+    }
+
+    auto size() const noexcept
+    {
+        return elems_.size();
+    }
+
+    //typedef typename List<T>::const_iterator const_iterator;
+    struct const_iterator
+    {
+        // Minimalistic const_iter implementation
+        typedef typename List<T>::const_iterator ConstIter;
+        typedef typename std::iterator_traits<ConstIter> ConstIterTraits;
+        using difference_type = typename ConstIterTraits::difference_type;
+        using value_type = typename ConstIterTraits::value_type;
+        using pointer = typename ConstIterTraits::pointer;
+        using reference = typename ConstIterTraits::reference;
+//using iterator_category = typename ConstIterTraits::iterator_category;
+        using iterator_category = std::random_access_iterator_tag;
+
+        ConstIter base_;
+
+        const_iterator(ConstIter it) : base_(it) {}
+        const_iterator() : base_() {}
+        ConstIter& operator=(ConstIter it)
+        {
+            base_ = it;
+            return base_;
+        }
+
+    // Indexing operators
+
+        value_type operator*() const
+        {
+            return *base_;
+        }
+        //value_type operator->() const
+        //{
+        //    return *base_;
+        //}
+        value_type operator[](const difference_type i) const
+        {
+           return *(base_+i);
+        }
+
+    // Increment/Decrement
+
+        const_iterator& operator++()
+        {
+            ++base_;
+            return *this;
+        }
+        const_iterator operator++(int)
+        {
+            const_iterator old(*this);
+            ++base_;
+            return old;
+        }
+        const_iterator& operator--()
+        {
+            --base_;
+            return *this;
+        }
+        const_iterator operator--(int)
+        {
+            const_iterator old(*this);
+            --base_;
+            return old;
+        }
+
+    // Arithmetic operators
+
+        const_iterator operator+(const difference_type i) const
+        {
+            return const_iterator(base_+i);
+        }
+        const_iterator operator-(const difference_type i) const
+        {
+            return const_iterator(base_-i);
+        }
+        difference_type operator-(const const_iterator& it) const
+        {
+            return (base_ - it.base_);
+        }
+
+    // Comparison operators
+
+        bool operator==(const const_iterator& it) const
+        {
+            return (base_ == it.base_);
+        }
+        bool operator!=(const const_iterator& it) const
+        {
+            return (base_ != it.base_);
+        }
+        bool operator<(const const_iterator& it) const
+        {
+            return (base_ < it.base_);
+        }
+        bool operator<=(const const_iterator& it) const
+        {
+            return (base_ <= it.base_);
+        }
+        bool operator>(const const_iterator& it) const
+        {
+            return (base_ > it.base_);
+        }
+        bool operator>=(const const_iterator& it) const
+        {
+            return (base_ >= it.base_);
+        }
+
+    // Compound assignment 
+
+        const_iterator& operator+=(const difference_type n)
+        {
+            base_ += n;
+            return *this;
+        }
+        const_iterator& operator-=(const difference_type n)
+        {
+            base_ -= n;
+            return *this;
+        }
+    };
+    auto cbegin() const
+    {
+        return const_iterator(elems_.cbegin());
+    }
+    auto cend() const
+    {
+        return const_iterator(elems_.cend());
+    }
+};
+
+
+// Wrap of tmp to List
+// ~~~~~~~~~~~~~~~~~~~
+
+template<class T>
+class ListTmpWrap
+:
+    public ListExpression<ListTmpWrap<T>, typename T::value_type>
+{
+    tmp<T> elems_;
+
+public:
+
+    //- Have expressions use copy to maintain tmp refCount
+    static constexpr bool is_leaf = false;
+
+    typedef typename T::value_type value_type;
+
+    // Construct from components
+    ListTmpWrap(const tmp<T>& elems)
+    :
+        elems_(elems)
+    {}
+
+    // Construct from components
+    ListTmpWrap(tmp<T>& elems)
+    :
+        elems_(elems)
+    {}
+
+    // Construct from ListExpression, forcing its evaluation.
+    template<typename E>
+    ListTmpWrap
+    (
+        tmp<T>& elems,
+        const ListExpression<E, typename E::value_type>& expr
+    )
+    :
+        elems_(elems)
+    {
+        expr.evaluate(elems_.ref());
+    }
+
+    //- Assignment
+    template<typename E>
+    void operator=(const ListExpression<E, typename E::value_type>& expr)
+    {
+        expr.evaluate(elems_.ref());
+    }
+
+    // returns the underlying data
+    const T& data() const
+    {
+        return elems_();
+    }
+
+    value_type operator[](const label i) const
+    {
+        return elems_()[i];
+    }
+
+    auto size() const noexcept
+    {
+        return elems_().size();
+    }
+
+    typedef typename T::const_iterator const_iterator;
+    auto cbegin() const
+    {
+        return elems_().cbegin();
+    }
+    auto cend() const
+    {
+        return elems_().cend();
+    }
+    typedef typename List<T>::iterator iterator;
+    auto begin()
+    {
+        return elems_().begin();
+    }
+    auto end()
+    {
+        return elems_().end();
+    }
+};
+
+
+template<class T>
+class ListConstTmpWrap
+:
+    public ListExpression<ListConstTmpWrap<T>, typename T::value_type>
+{
+    const tmp<T> elems_;
+
+public:
+
+    //- Have expressions use copy to maintain tmp refCount
+    static constexpr bool is_leaf = false;
+
+    typedef typename T::value_type value_type;
+
+    // Construct from components
+    ListConstTmpWrap(const tmp<T>& elems)
+    :
+        elems_(elems)
+    {}
+
+    // returns the underlying data
+    const T& data() const
+    {
+        return elems_();
+    }
+
+    value_type operator[](const label i) const
+    {
+        return elems_()[i];
+    }
+
+    auto size() const noexcept
+    {
+        return elems_().size();
+    }
+
+    //- Define our own iterator so we can refer to it. 'normal', e.g. double*
+    //  iterator will not work we'd have to use the iterator traits to get
+    //  at its e.g. difference_type from other classes. TBD.
+    struct const_iterator
+    {
+        // Minimalistic const_iter implementation
+        typedef typename T::const_iterator ConstIter;
+        typedef typename std::iterator_traits<ConstIter> ConstIterTraits;
+        using difference_type = typename ConstIterTraits::difference_type;
+        using value_type = typename ConstIterTraits::value_type;
+        using pointer = typename ConstIterTraits::pointer;
+        using reference = typename ConstIterTraits::reference;
+        using iterator_category = typename ConstIterTraits::iterator_category;
+
+        ConstIter base_;
+
+        const_iterator(ConstIter it) : base_(it) {}
+        value_type operator*() const
+        {
+            return *base_;
+        }
+        value_type operator[](const difference_type i) const
+        {
+            return base_[i];
+        }
+        const_iterator& operator++()
+        {
+            ++base_;
+            return *this;
+        }
+        const_iterator operator++(int)
+        {
+            const_iterator old(*this);
+            ++base_;
+            return old;
+        }
+        const_iterator& operator--()
+        {
+            --base_;
+            return *this;
+        }
+        const_iterator operator--(int)
+        {
+            const_iterator old(*this);
+            --base_;
+            return old;
+        }
+        const_iterator operator+(const difference_type i) const
+        {
+            return const_iterator(base_+i);
+        }
+        const_iterator operator-(const difference_type i) const
+        {
+            return const_iterator(base_-i);
+        }
+        difference_type operator-(const const_iterator& it) const
+        {
+            return (base_ - it.base_);
+        }
+        bool operator==(const const_iterator& it) const
+        {
+            return (base_ == it.base_);
+        }
+        bool operator!=(const const_iterator& it) const
+        {
+            return (base_ != it.base_);
+        }
+        bool operator<(const const_iterator& it) const
+        {
+            return (base_ < it.base_);
+        }
+        bool operator<=(const const_iterator& it) const
+        {
+            return (base_ <= it.base_);
+        }
+        bool operator>(const const_iterator& it) const
+        {
+            return (base_ > it.base_);
+        }
+        bool operator>=(const const_iterator& it) const
+        {
+            return (base_ >= it.base_);
+        }
+        const_iterator& operator+=(const difference_type n)
+        {
+            base_ += n;
+            return base_;
+        }
+        const_iterator& operator-=(const difference_type n)
+        {
+            base_ -= n;
+            return base_;
+        }
+    };
+    //typedef typename T::const_iterator const_iterator;
+    auto cbegin() const
+    {
+        return const_iterator(elems_().cbegin());
+    }
+    auto cend() const
+    {
+        return const_iterator(elems_().cend());
+    }
+};
+
+
+// Wrap of indirection
+// ~~~~~~~~~~~~~~~~~~~
+
+// E1 : expression to get value of element i
+template<typename E1>
+class IndirectConstWrap
+:
+    public ListExpression
+    <
+        IndirectConstWrap<E1>,
+        typename E1::value_type
+    >
+{
+    //- Expression for cell values
+    typename std::conditional
+    <
+        E1::is_leaf,
+        const E1&,
+        const E1
+    >::type values_;
+
+    const labelUList& addr_;
+
+public:
+
+    // ! Store as copy (since holds reference)
+    static constexpr bool is_leaf = false;
+
+    typedef typename E1::value_type value_type;
+
+    // construct from components
+    IndirectConstWrap
+    (
+        const E1& values,
+        const labelUList& addr
+    )
+    :
+        values_(values),
+        addr_(addr)
+    {}
+
+    value_type operator[](const label i) const
+    {
+        return values_[addr_[i]];
+    }
+    label size() const
+    {
+        return addr_.size();
+    }
+
+    struct const_iterator
+    {
+        // See IndirectListBase::const_iterator
+
+        //- Step type for addressing
+        using difference_type = label;
+
+        //- Iterator definition for E1
+        typedef typename E1::const_iterator ConstIter;
+        //using difference_type = typename ConstIter::difference_type;
+        using value_type = typename ConstIter::value_type;
+        using pointer = typename ConstIter::pointer;
+        using reference = typename ConstIter::reference;
+        using iterator_category = std::random_access_iterator_tag;
+
+        ConstIter begin_;
+        typename labelUList::const_iterator iter_;
+
+        const_iterator
+        (
+            ConstIter begin,
+            typename labelUList::const_iterator addrIter
+        )
+        :
+            begin_(begin),
+            iter_(addrIter)
+        {}
+        auto operator*() const
+        {
+            return *(begin_ + *iter_);
+        }
+        value_type operator[](const difference_type i) const
+        {
+            return *(begin_+iter_[i]);
+        }
+        const_iterator& operator++()
+        {
+            ++iter_;
+            return *this;
+        }
+        const_iterator operator++(int)
+        {
+            const_iterator old(*this);
+            ++iter_;
+            return old;
+        }
+        const_iterator& operator--()
+        {
+            --iter_;
+            return *this;
+        }
+        const_iterator operator--(int)
+        {
+            const_iterator old(*this);
+            --iter_;
+            return old;
+        }
+        const_iterator operator+(const difference_type i) const
+        {
+            return const_iterator(begin_, iter_+i);
+        }
+        const_iterator operator-(const difference_type i) const
+        {
+            return const_iterator(begin_, iter_-i);
+        }
+        difference_type operator-(const const_iterator& it) const
+        {
+            return (iter_ - it.iter_);
+        }
+        bool operator==(const const_iterator& it) const
+        {
+            return iter_ == it.iter_;
+        }
+        bool operator!=(const const_iterator& it) const
+        {
+            return iter_ != it.iter_;
+        }
+        bool operator<(const const_iterator& it) const
+        {
+            return (iter_ < it.iter_);
+        }
+        bool operator<=(const const_iterator& it) const
+        {
+            return (iter_ <= it.iter_);
+        }
+        bool operator>(const const_iterator& it) const
+        {
+            return (iter_ > it.iter_);
+        }
+        bool operator>=(const const_iterator& it) const
+        {
+            return (iter_ >= it.iter_);
+        }
+        const_iterator& operator+=(const difference_type n)
+        {
+            iter_ += n;
+            return *this;
+        }
+        const_iterator& operator-=(const difference_type n)
+        {
+            iter_ -= n;
+            return *this;
+        }
+    };
+    auto cbegin() const
+    {
+        return const_iterator(values_.cbegin(), addr_.cbegin());
+    }
+    auto cend() const
+    {
+        return const_iterator(values_.cbegin(), addr_.cend());
+    }
+};
+
+
+// Expresssions on List
+// ~~~~~~~~~~~~~~~~~~~~
+
+#define EXPRESSION_FUNCTION1(Func, BaseFunc, OpFunc)                           \
+template<typename E1>                                                          \
+class OpFunc                                                                   \
+:                                                                              \
+    public ListExpression<OpFunc<E1>, typename E1::value_type>                 \
+{                                                                              \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    typedef typename E1::value_type value_type;                                \
+                                                                               \
+    OpFunc(const E1& u)                                                        \
+    :                                                                          \
+        u_(u)                                                                  \
+    {}                                                                         \
+                                                                               \
+    auto operator[](const label i) const                                       \
+    {                                                                          \
+        return BaseFunc(u_[i]);                                                \
+    }                                                                          \
+                                                                               \
+    auto size() const noexcept { return u_.size(); }                           \
+                                                                               \
+    struct const_iterator                                                      \
+    {                                                                          \
+        typedef typename E1::const_iterator ConstIter;                         \
+        using difference_type = typename ConstIter::difference_type;           \
+        using value_type = typename ConstIter::value_type;                     \
+        using pointer = typename ConstIter::pointer;                           \
+        using reference = typename ConstIter::reference;                       \
+        using iterator_category = typename ConstIter::iterator_category;       \
+                                                                               \
+        ConstIter base_;                                                       \
+                                                                               \
+        const_iterator(ConstIter it) : base_(it) {}                            \
+        auto operator*() const                                           \
+        {                                                                      \
+            return BaseFunc(*base_);                                           \
+        }                                                                      \
+        auto operator[](const difference_type i) const                   \
+        {                                                                      \
+            return BaseFunc(base_[i]);                                         \
+        }                                                                      \
+        const_iterator& operator++()                                           \
+        {                                                                      \
+            ++base_;                                                           \
+            return *this;                                                      \
+        }                                                                      \
+        const_iterator operator++(int)                                         \
+        {                                                                      \
+            const_iterator old(*this);                                         \
+            ++base_;                                                           \
+            return old;                                                        \
+        }                                                                      \
+        const_iterator& operator--()                                           \
+        {                                                                      \
+            --base_;                                                           \
+            return *this;                                                      \
+        }                                                                      \
+        const_iterator operator--(int)                                         \
+        {                                                                      \
+            const_iterator old(*this);                                         \
+            --base_;                                                           \
+            return old;                                                        \
+        }                                                                      \
+        const_iterator operator+(const difference_type i) const                \
+        {                                                                      \
+            return const_iterator(base_+i);                                    \
+        }                                                                      \
+        const_iterator operator-(const difference_type i) const                \
+        {                                                                      \
+            return const_iterator(base_-i);                                    \
+        }                                                                      \
+        difference_type operator-(const const_iterator& it) const              \
+        {                                                                      \
+            return (base_ - it.base_);                                         \
+        }                                                                      \
+        bool operator==(const const_iterator& it) const                        \
+        {                                                                      \
+            return base_ == it.base_;                                          \
+        }                                                                      \
+        bool operator!=(const const_iterator& it) const                        \
+        {                                                                      \
+            return base_ != it.base_;                                          \
+        }                                                                      \
+        bool operator<(const const_iterator& it) const                         \
+        {                                                                      \
+            return (base_ < it.base_);                                         \
+        }                                                                      \
+        bool operator<=(const const_iterator& it) const                        \
+        {                                                                      \
+            return (base_ <= it.base_);                                        \
+        }                                                                      \
+        bool operator>(const const_iterator& it) const                         \
+        {                                                                      \
+            return (base_ > it.base_);                                         \
+        }                                                                      \
+        bool operator>=(const const_iterator& it) const                        \
+        {                                                                      \
+            return (base_ >= it.base_);                                        \
+        }                                                                      \
+        const_iterator& operator+=(const difference_type n)                    \
+        {                                                                      \
+            base_ += n;                                                        \
+            return *this;                                                      \
+        }                                                                      \
+        const_iterator& operator-=(const difference_type n)                    \
+        {                                                                      \
+            base_ -= n;                                                        \
+            return *this;                                                      \
+        }                                                                      \
+    };                                                                         \
+    const_iterator cbegin() const                                              \
+    {                                                                          \
+        return const_iterator(u_.cbegin());                                    \
+    }                                                                          \
+    const_iterator cend() const                                                \
+    {                                                                          \
+        return const_iterator(u_.cend());                                      \
+    }                                                                          \
+};                                                                             \
+template<typename E1>                                                          \
+OpFunc<E1> Func                                                                \
+(                                                                              \
+    const ListExpression<E1, typename E1::value_type>& u                       \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1>(static_cast<const E1&>(u));                              \
+}
+
+
+// Do '-' separately until we work out macro expansion...
+template<typename E1>
+class List_negate
+:
+    public ListExpression<List_negate<E1>, typename E1::value_type>
+{
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;
+
+public:
+    static constexpr bool is_leaf = false;
+
+    typedef typename E1::value_type value_type;
+
+    List_negate(const E1& u)
+    :
+        u_(u)
+    {}
+
+    auto operator[](const label i) const
+    {
+        return -u_[i];
+    }
+
+    auto size() const noexcept { return u_.size(); }
+
+    struct const_iterator
+    {
+        // Minimalistic const_iter implementation
+
+        typedef typename E1::const_iterator ConstIter;
+        using difference_type = typename ConstIter::difference_type;
+        using value_type = typename ConstIter::value_type;
+        using pointer = typename ConstIter::pointer;
+        using reference = typename ConstIter::reference;
+        using iterator_category = typename ConstIter::iterator_category;
+
+        ConstIter base_;
+
+        const_iterator(ConstIter it) : base_(it) {}
+        value_type operator*() const
+        {
+            return -(*base_);
+        }
+        value_type operator[](const difference_type i) const
+        {
+           return -*(base_+i);
+        }
+        const_iterator& operator++()
+        {
+            ++base_;
+            return *this;
+        }
+        const_iterator operator++(int)
+        {
+            const_iterator old(*this);
+            ++base_;
+            return old;
+        }
+        const_iterator& operator--()
+        {
+            --base_;
+            return *this;
+        }
+        const_iterator operator--(int)
+        {
+            const_iterator old(*this);
+            --base_;
+            return old;
+        }
+        const_iterator operator+(const difference_type i) const
+        {
+            return const_iterator(base_+i);
+        }
+        const_iterator operator-(const difference_type i) const
+        {
+            return const_iterator(base_-i);
+        }
+        difference_type operator-(const const_iterator& it) const
+        {
+            return (base_ - it.base_);
+        }
+        bool operator==(const const_iterator& it) const
+        {
+            return base_ == it.base_;
+        }
+        bool operator!=(const const_iterator& it) const
+        {
+            return base_ != it.base_;
+        }
+        bool operator<(const const_iterator& it) const
+        {
+            return (base_ < it.base_);
+        }
+        bool operator<=(const const_iterator& it) const
+        {
+            return (base_ <= it.base_);
+        }
+        bool operator>(const const_iterator& it) const
+        {
+            return (base_ > it.base_);
+        }
+        bool operator>=(const const_iterator& it) const
+        {
+            return (base_ >= it.base_);
+        }
+        const_iterator& operator+=(const difference_type n)
+        {
+            base_ += n;
+            return base_;
+        }
+        const_iterator& operator-=(const difference_type n)
+        {
+            base_ -= n;
+            return base_;
+        }
+    };
+    const_iterator cbegin() const
+    {
+        return const_iterator(u_.cbegin());
+    }
+    const_iterator cend() const
+    {
+        return const_iterator(u_.cend());
+    }
+};
+template<typename E1>
+List_negate<E1> operator-
+(
+    const ListExpression<E1, typename E1::value_type>& u
+)
+{
+    return List_negate<E1>(static_cast<const E1&>(u));
+}
+
+
+#define EXPRESSION_FUNCTION2(Func, BaseFunc, OpFunc)                           \
+template<typename E1, typename E2>                                             \
+class OpFunc                                                                   \
+:                                                                              \
+    public ListExpression                                                      \
+    <                                                                          \
+        OpFunc<E1, E2>,                                                        \
+        typename E1::value_type                                                \
+    >                                                                          \
+{                                                                              \
+    /* cref if leaf, copy otherwise */                                         \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+    typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    /* Type to return for internal field */                                    \
+    typedef typename E1::value_type value_type;                                \
+                                                                               \
+    OpFunc(const E1& u, const E2& v)                                           \
+    :                                                                          \
+     u_(u), v_(v)                                                              \
+    {                                                                          \
+        assert(u.size() == v.size());                                          \
+    }                                                                          \
+    auto operator[](const label i) const                                       \
+    {                                                                          \
+        return BaseFunc(u_[i], v_[i]);                                         \
+    }                                                                          \
+    auto size() const noexcept {return u_.size(); }                            \
+                                                                               \
+    struct const_iterator                                                      \
+    {                                                                          \
+        typedef typename E1::const_iterator E1Iter;                            \
+        typedef typename E2::const_iterator E2Iter;                            \
+        using difference_type = typename E1Iter::difference_type;              \
+        using value_type = typename E1Iter::value_type;                        \
+        using pointer = typename E1Iter::pointer;                              \
+        using reference = typename E1Iter::reference;                          \
+        using iterator_category = typename E1Iter::iterator_category;          \
+                                                                               \
+        E1Iter uIter_;                                                         \
+        E2Iter vIter_;                                                         \
+                                                                               \
+        const_iterator(E1Iter uIter, E2Iter vIter)                             \
+        :                                                                      \
+            uIter_(uIter),                                                     \
+            vIter_(vIter)                                                      \
+        {}                                                                     \
+        auto operator*() const                                                   \
+        {                                                                      \
+            return BaseFunc(*uIter_, *vIter_);                                 \
+        }                                                                      \
+        const_iterator& operator++()                                           \
+        {                                                                      \
+            ++uIter_;                                                          \
+            ++vIter_;                                                          \
+            return *this;                                                      \
+        }                                                                      \
+        difference_type operator-(const const_iterator& it) const              \
+        {                                                                      \
+            return (uIter_ - it.uIter_);                                       \
+        }                                                                      \
+        const_iterator operator+(const label i) const                          \
+        {                                                                      \
+            return const_iterator(uIter_+i, vIter_+i);                         \
+        }                                                                      \
+        bool operator==(const const_iterator& it) const                        \
+        {                                                                      \
+            return uIter_ == it.uIter_ && vIter_ == it.vIter_;                 \
+        }                                                                      \
+        bool operator!=(const const_iterator& it) const                        \
+        {                                                                      \
+            return !operator==(it);                                            \
+        }                                                                      \
+        bool operator<(const const_iterator& it) const                         \
+        {                                                                      \
+            return (uIter_ < it.uIter_);                                       \
+        }                                                                      \
+    };                                                                         \
+    const_iterator cbegin() const                                              \
+    {                                                                          \
+        return const_iterator(u_.cbegin(), v_.cbegin());                       \
+    }                                                                          \
+    const_iterator cend() const                                                \
+    {                                                                          \
+        return const_iterator(u_.cend(), v_.cend());                           \
+    }                                                                          \
+};                                                                             \
+template<typename E1, typename E2>                                             \
+OpFunc<E1, E2> Func                                                            \
+(                                                                              \
+    ListExpression                                                             \
+    <                                                                          \
+        E1,                                                                    \
+        typename E1::value_type                                                \
+    > const& u,                                                                \
+    ListExpression                                                             \
+    <                                                                          \
+        E2,                                                                    \
+        typename E2::value_type                                                \
+    > const& v                                                                 \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1, E2>                                                      \
+    (                                                                          \
+        static_cast<const E1&>(u),                                             \
+        static_cast<const E2&>(v)                                              \
+    );                                                                         \
+}
+#define EXPRESSION_OPERATOR(Op, OpFunc)                                        \
+template<typename E1, typename E2>                                             \
+class OpFunc                                                                   \
+:                                                                              \
+    public ListExpression                                                      \
+    <                                                                          \
+        OpFunc<E1, E2>,                                                        \
+        typename E1::value_type                                                \
+    >                                                                          \
+{                                                                              \
+    /* cref if leaf, copy otherwise */                                         \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+    typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    /* Type to return for internal field */                                    \
+    typedef typename E1::value_type value_type;                                \
+                                                                               \
+    OpFunc(const E1& u, const E2& v)                                           \
+    :                                                                          \
+     u_(u), v_(v)                                                              \
+    {                                                                          \
+        assert(u.size() == v.size());                                          \
+    }                                                                          \
+    auto operator[](const label i) const                                       \
+    {                                                                          \
+        return u_[i] Op v_[i];                                                 \
+    }                                                                          \
+    auto size() const noexcept { return u_.size(); }                           \
+                                                                               \
+    struct const_iterator                                                      \
+    {                                                                          \
+        typedef typename E1::const_iterator E1Iter;                            \
+        typedef typename E2::const_iterator E2Iter;                            \
+        using difference_type = typename E1Iter::difference_type;              \
+        using value_type = typename E1Iter::value_type;                        \
+        using pointer = typename E1Iter::pointer;                              \
+        using reference = typename E1Iter::reference;                          \
+        using iterator_category = typename E1Iter::iterator_category;          \
+                                                                               \
+        E1Iter uIter_;                                                         \
+        E2Iter vIter_;                                                         \
+                                                                               \
+        const_iterator(E1Iter uIter, E2Iter vIter)                             \
+        :                                                                      \
+            uIter_(uIter),                                                     \
+            vIter_(vIter)                                                      \
+        {}                                                                     \
+        value_type operator*() const                                           \
+        {                                                                      \
+            return *uIter_ Op *vIter_;                                         \
+        }                                                                      \
+        value_type operator[](const difference_type i) const                   \
+        {                                                                      \
+           return uIter_[i] Op vIter_[i];                                      \
+        }                                                                      \
+        const_iterator& operator++()                                           \
+        {                                                                      \
+            ++uIter_;                                                          \
+            ++vIter_;                                                          \
+            return *this;                                                      \
+        }                                                                      \
+        const_iterator operator++(int)                                         \
+        {                                                                      \
+            const_iterator old(*this);                                         \
+            ++uIter_;                                                          \
+            ++vIter_;                                                          \
+            return old;                                                        \
+        }                                                                      \
+        const_iterator& operator--()                                           \
+        {                                                                      \
+            --uIter_;                                                          \
+            --vIter_;                                                          \
+            return *this;                                                      \
+        }                                                                      \
+        const_iterator operator--(int)                                         \
+        {                                                                      \
+            const_iterator old(*this);                                         \
+            --uIter_;                                                          \
+            --vIter_;                                                          \
+            return old;                                                        \
+        }                                                                      \
+        const_iterator operator+(const difference_type i) const                \
+        {                                                                      \
+            return const_iterator(uIter_+i, vIter_+i);                         \
+        }                                                                      \
+        const_iterator operator-(const difference_type i) const                \
+        {                                                                      \
+            return const_iterator(uIter_-i, vIter_-i);                         \
+        }                                                                      \
+        difference_type operator-(const const_iterator& it) const              \
+        {                                                                      \
+            return (uIter_ - it.uIter_);                                       \
+        }                                                                      \
+        bool operator==(const const_iterator& it) const                        \
+        {                                                                      \
+            return uIter_ == it.uIter_ && vIter_ == it.vIter_;                 \
+        }                                                                      \
+        bool operator!=(const const_iterator& it) const                        \
+        {                                                                      \
+            return !operator==(it);                                            \
+        }                                                                      \
+        bool operator<(const const_iterator& it) const                         \
+        {                                                                      \
+            return (uIter_ < it.uIter_);                                       \
+        }                                                                      \
+        bool operator<=(const const_iterator& it) const                        \
+        {                                                                      \
+            return (uIter_ <= it.uIter_);                                      \
+        }                                                                      \
+        bool operator>(const const_iterator& it) const                         \
+        {                                                                      \
+            return (uIter_ > it.uIter_);                                       \
+        }                                                                      \
+        bool operator>=(const const_iterator& it) const                        \
+        {                                                                      \
+            return (uIter_ >= it.uIter_);                                      \
+        }                                                                      \
+        const_iterator& operator+=(const difference_type n)                    \
+        {                                                                      \
+            uIter_ += n;                                                       \
+            vIter_ += n;                                                       \
+            return *this;                                                      \
+        }                                                                      \
+        const_iterator& operator-=(const difference_type n)                    \
+        {                                                                      \
+            uIter_ -= n;                                                       \
+            vIter_ -= n;                                                       \
+            return *this;                                                      \
+        }                                                                      \
+    };                                                                         \
+    const_iterator cbegin() const                                              \
+    {                                                                          \
+        return const_iterator(u_.cbegin(), v_.cbegin());                       \
+    }                                                                          \
+    const_iterator cend() const                                                \
+    {                                                                          \
+        return const_iterator(u_.cend(), v_.cend());                           \
+    }                                                                          \
+};                                                                             \
+template<typename E1, typename E2>                                             \
+OpFunc<E1, E2>                                                                 \
+operator Op                                                                    \
+(                                                                              \
+    ListExpression                                                             \
+    <                                                                          \
+        E1,                                                                    \
+        typename E1::value_type                                                \
+    > const& u,                                                                \
+    ListExpression                                                             \
+    <                                                                          \
+        E2,                                                                    \
+        typename E2::value_type                                                \
+    > const& v                                                                 \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1, E2>                                                      \
+    (                                                                          \
+        static_cast<const E1&>(u),                                             \
+        static_cast<const E2&>(v)                                              \
+    );                                                                         \
+}
+
+
+//- List functions
+// macro arguments:
+//  - 'sin' : name of new function (in ListExpression namespace)
+//  - '::sin' : per-element function to call
+//  - 'List_sin' : generated expression helper class
+EXPRESSION_FUNCTION1(sin, ::sin, List_sin)
+EXPRESSION_FUNCTION1(cos, ::cos, List_cos)
+EXPRESSION_FUNCTION1(tan, ::tan, List_tan)
+EXPRESSION_FUNCTION1(sinh, ::sinh, List_sinh)
+EXPRESSION_FUNCTION1(cosh, ::cosh, List_cosh)
+EXPRESSION_FUNCTION1(tanh, ::tanh, List_tanh)
+EXPRESSION_FUNCTION1(sqr, Foam::sqr, List_sqr)
+EXPRESSION_FUNCTION1(sqrt, Foam::sqrt, List_sqrt)
+EXPRESSION_FUNCTION1(magSqr, Foam::magSqr, List_magSqr)
+EXPRESSION_FUNCTION1(mag, Foam::mag, List_mag)
+EXPRESSION_FUNCTION1(symm, Foam::symm, List_symm)
+//EXPRESSION_FUNCTION1(-, -, List_negate)
+EXPRESSION_FUNCTION1(pow2, Foam::pow2, List_pow2)
+EXPRESSION_FUNCTION1(pow3, Foam::pow3, List_pow3)
+EXPRESSION_FUNCTION1(pow4, Foam::pow4, List_pow4)
+
+#undef EXPRESSION_FUNCTION1
+
+EXPRESSION_FUNCTION2(min, Foam::min, List_min)
+EXPRESSION_FUNCTION2(max, Foam::max, List_max)
+
+#undef EXPRESSION_FUNCTION2
+
+//- List operators
+EXPRESSION_OPERATOR(+, List_add)
+EXPRESSION_OPERATOR(-, List_subtract)
+EXPRESSION_OPERATOR(*, List_multiply)
+EXPRESSION_OPERATOR(/, List_divide)
+EXPRESSION_OPERATOR(&, List_dot)
+
+#undef EXPRESSION_OPERATOR
+
+//- Not possible: a *= b has to return a new type so it can never be in-place
+//EXPRESSION_LIST_OPERATOR(+=, List_inplace_add)
+//EXPRESSION_LIST_OPERATOR(-=, List_inplace_subtract)
+//EXPRESSION_LIST_OPERATOR(*=, List_inplace_multiply)
+//EXPRESSION_LIST_OPERATOR(/=, List_inplace_divide)
+
+
+// Interpolation
+// ~~~~~~~~~~~~~
+
+// E1 : expression to get value of element i
+// E2 : expression to get value of element i
+// E3 : expression to get interpolation weight
+template<typename E1, typename E2, typename E3>
+class List_interpolate
+:
+    public ListExpression
+    <
+        List_interpolate<E1, E2, E3>,
+        typename E1::value_type
+    >
+{
+    //- Expression for cell values1
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;
+
+    //- Expression for cell values2
+    typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_;
+
+    //- Expression for interpolation weight
+    typename std::conditional<E3::is_leaf, const E3&, const E3>::type weight_;
+
+
+public:
+
+    // ! Store as copy (since holds reference)
+    static constexpr bool is_leaf = false;
+
+    typedef typename E1::value_type value_type;
+
+    // construct from components
+    List_interpolate(const E1& u, const E2& v, const E3& weight)
+    :
+        u_(u),
+        v_(v),
+        weight_(weight)
+    {}
+
+    value_type operator[](const label i) const
+    {
+        const auto ownVal = u_[i];
+        const auto neiVal = v_[i];
+        return weight_[i]*(ownVal-neiVal) + neiVal;
+    }
+    label size() const
+    {
+        return weight_.size();
+    }
+
+    struct const_iterator
+    {
+        typedef typename E1::const_iterator E1Iter;
+        typedef typename E2::const_iterator E2Iter;
+        typedef typename E3::const_iterator E3Iter;
+        using difference_type = typename E3Iter::difference_type;
+        using value_type = typename E1Iter::value_type;
+        using pointer = typename E1Iter::pointer;
+        using reference = typename E1Iter::reference;
+        using iterator_category = typename E3Iter::iterator_category;
+        E1Iter uIter_;
+        E2Iter vIter_;
+        E3Iter weightIter_;
+
+        const_iterator(E1Iter uIter, E2Iter vIter, E3Iter weightIter)
+        :
+            uIter_(uIter),
+            vIter_(vIter),
+            weightIter_(weightIter)
+        {}
+        value_type operator*() const
+        {
+            const auto ownVal = *uIter_;
+            const auto neiVal = *vIter_;
+            return *weightIter_*(ownVal-neiVal) + neiVal;
+        }
+        value_type operator[](const difference_type i) const
+        {
+            const auto ownVal = uIter_[i];
+            const auto neiVal = vIter_[i];
+            return weightIter_[i]*(ownVal-neiVal) + neiVal;
+        }
+        const_iterator& operator++()
+        {
+            ++uIter_;
+            ++vIter_;
+            ++weightIter_;
+            return *this;
+        }
+        const_iterator operator++(int)
+        {
+            const_iterator old(*this);
+            ++uIter_;
+            ++vIter_;
+            ++weightIter_;
+            return old;
+        }
+        const_iterator& operator--()
+        {
+            --uIter_;
+            --vIter_;
+            --weightIter_;
+            return *this;
+        }
+        const_iterator operator--(int)
+        {
+            const_iterator old(*this);
+            --uIter_;
+            --vIter_;
+            --weightIter_;
+            return old;
+        }
+        const_iterator operator+(const label i) const
+        {
+            return const_iterator(uIter_+i, vIter_+i, weightIter_+i);
+        }
+        const_iterator operator-(const label i) const
+        {
+            return const_iterator(uIter_-i, vIter_-i, weightIter_-i);
+        }
+        difference_type operator-(const const_iterator& it) const
+        {
+            return (uIter_ - it.uIter_);
+        }
+        bool operator==(const const_iterator& it) const
+        {
+            return
+                uIter_ == it.uIter_
+             && vIter_ == it.vIter_
+             && weightIter_ == it.weightIter_;
+        }
+        bool operator!=(const const_iterator& it) const
+        {
+            return !operator==(it);
+        }
+        bool operator<(const const_iterator& it) const
+        {
+            return (uIter_ < it.uIter_);
+        }
+        bool operator<=(const const_iterator& it) const
+        {
+            return (uIter_ <= it.uIter_);
+        }
+        bool operator>(const const_iterator& it) const
+        {
+            return (uIter_ > it.uIter_);
+        }
+        bool operator>=(const const_iterator& it) const
+        {
+            return (uIter_ >= it.uIter_);
+        }
+        const_iterator& operator+=(const difference_type n)
+        {
+            uIter_ += n;
+            vIter_ += n;
+            weightIter_ += n;
+            return *this;
+        }
+        const_iterator& operator-=(const difference_type n)
+        {
+            uIter_ -= n;
+            vIter_ -= n;
+            weightIter_ -= n;
+            return *this;
+        }
+    };
+    const_iterator cbegin() const
+    {
+        return const_iterator(u_.cbegin(), v_.cbegin(), weight_.cbegin());
+    }
+    const_iterator cend() const
+    {
+        return const_iterator(u_.cend(), v_.cend(), weight_.cend());
+    }
+};
+
+
+// Expressions on constants
+// ~~~~~~~~~~~~~~~~~~~~~~~~
+
+template<class T>
+class UniformListWrap
+:
+    public ListExpression<UniformListWrap<T>, T>
+{
+private:
+
+    const label size_;
+
+    const T val_;
+
+
+public:
+
+    static constexpr bool is_leaf = false;
+
+    typedef T value_type;
+
+    // construct from components
+    UniformListWrap(const label size, const T val)
+    :
+        size_(size),
+        val_(val)
+    {}
+
+    // construct from components
+    UniformListWrap(const label size, const dimensioned<T>& val)
+    :
+        size_(size),
+        val_(val.value())
+    {}
+
+    T operator[](const label i) const
+    {
+        return val_;
+    }
+    auto size() const noexcept
+    {
+        return size_;
+    }
+
+    struct const_iterator
+    {
+        // Minimalistic const_iter implementation
+        using difference_type = label;
+        using value_type = const T;
+        using pointer = const T*;
+        using reference = const T&;
+        using iterator_category = std::random_access_iterator_tag;
+
+        difference_type count_;
+
+        const T val_;
+
+        const_iterator(const difference_type count, const T val)
+        :
+            count_(count),
+            val_(val)
+        {}
+        value_type operator*() const
+        {
+            return val_;
+        }
+        value_type operator[](const difference_type i) const
+        {
+           return val_;
+        }
+        const_iterator& operator++()
+        {
+            count_++;
+            return *this;
+        }
+        const_iterator operator++(int)
+        {
+            const_iterator old(*this);
+            ++count_;
+            return old;
+        }
+        const_iterator& operator--()
+        {
+            count_--;
+            return *this;
+        }
+        const_iterator operator--(int)
+        {
+            const_iterator old(*this);
+            --count_;
+            return old;
+        }
+        const_iterator operator+(const difference_type i) const
+        {
+            return const_iterator(count_+i, val_);
+        }
+        const_iterator operator-(const difference_type i) const
+        {
+            return const_iterator(count_-i, val_);
+        }
+        difference_type operator-(const const_iterator& it) const
+        {
+            return (count_ - it.count_);
+        }
+        bool operator==(const const_iterator& it) const
+        {
+            return val_ == it.val_ && count_ == it.count_;
+        }
+        bool operator!=(const const_iterator& it) const
+        {
+            return !operator==(it);
+        }
+        bool operator<(const const_iterator& it) const
+        {
+            return (count_ < it.count_);
+        }
+        bool operator<=(const const_iterator& it) const
+        {
+            return (count_ <= it.count_);
+        }
+        bool operator>(const const_iterator& it) const
+        {
+            return (count_ > it.count_);
+        }
+        bool operator>=(const const_iterator& it) const
+        {
+            return (count_ >= it.count_);
+        }
+        const_iterator& operator+=(const difference_type n)
+        {
+            count_ += n;
+            return *this;
+        }
+        const_iterator& operator-=(const difference_type n)
+        {
+            count_ -= n;
+            return *this;
+        }
+    };
+    auto cbegin() const
+    {
+        return const_iterator(0, val_);
+    }
+    auto cend() const
+    {
+        return const_iterator(size_, val_);
+    }
+};
+
+
+//- Wrapper of multiple lists (e.g. FieldField)
+template <class V>
+class ListsConstRefWrap
+:
+    public ListExpression<ListsConstRefWrap<V>, typename V::value_type>
+{
+    List<const V*> data;
+    List<size_t> offsets;
+    size_t total_size = 0;
+
+    void init()
+    {
+        offsets.resize_nocopy(data.size());
+        total_size = 0;
+        for (label i = 0; i < data.size(); ++i)
+        {
+            offsets[i] = total_size;
+            total_size += data[i]->size();
+        }
+    }
+
+public:
+    static constexpr bool is_leaf = false;  //true;
+
+    using value_type = typename V::value_type;  //T;
+    typedef ListsConstRefWrap<V> this_type;
+
+    ListsConstRefWrap() = delete;
+
+    ListsConstRefWrap(const UList<const V*>& d)
+    :
+        data(d)
+    {
+        init();
+    }
+    ListsConstRefWrap(List<const V*>&& d)
+    :
+        data(d)
+    {
+        init();
+    }
+    template<class Container>
+    ListsConstRefWrap(const Container& ds)
+    {
+        data.resize_nocopy(ds.size());
+        forAll(ds, i)
+        {
+            const V& vals = static_cast<const V&>(ds[i]);
+            data[i] = &vals;
+        }
+        init();
+    }
+    template<class Container>
+    ListsConstRefWrap(const Container& ds, const labelUList& elems)
+    {
+        data.resize_nocopy(elems.size());
+        forAll(ds, i)
+        {
+            const V& vals = static_cast<const V&>(ds[elems[i]]);
+            data[i] = &vals;
+        }
+        init();
+    }
+
+    void push_back(const V* vec)
+    {
+        data.push_back(vec);
+        offsets.push_back(total_size);
+        total_size += vec->size();
+    }
+
+    value_type operator[](label i) const
+    {    
+        const label vec_index =
+            std::upper_bound
+            (
+                offsets.begin(),
+                offsets.end(),
+                i
+            )
+         - offsets.begin() - 1;
+        const label elem_index = i - offsets[vec_index];
+        return (*data[vec_index])[elem_index];
+    }
+    label size() const { return total_size; }
+
+
+    typedef value_type* pointer;
+    typedef value_type& reference;
+    typedef const value_type* const_pointer;
+    typedef const value_type& const_reference;
+    typedef label difference_type;
+
+    struct const_iterator
+    {
+        using iterator_category = std::random_access_iterator_tag;
+        using difference_type   = this_type::difference_type;
+        using value_type  = this_type::value_type;
+        using pointer     = this_type::const_pointer;
+        using reference   = this_type::const_reference;
+
+        const this_type& base_; 
+        label index_;
+
+        const_iterator(const this_type& base, const label index)
+        :
+            base_(base),
+            index_(index)
+        {}
+        const_iterator(const const_iterator& it)
+        :
+            base_(it.base_),
+            index_(it.index_)
+        {}
+        const_iterator() = delete;  //: base_(0) {}
+        const_iterator& operator=(const_iterator& it) = default;
+
+    // Indexing operators
+
+        value_type operator*() const
+        {
+            return base_[index_];
+        }
+        //value_type operator->() const
+        //{
+        //    return *base_;
+        //}
+        value_type operator[](const difference_type i) const
+        {
+           return base_[index_];
+        }
+
+    // Increment/Decrement
+
+        const_iterator& operator++()
+        {
+            ++index_;
+            return *this;
+        }
+        const_iterator operator++(int)
+        {
+            const_iterator old(*this);
+            ++index_;
+            return old;
+        }
+        const_iterator& operator--()
+        {
+            --index_;
+            return *this;
+        }
+        const_iterator operator--(int)
+        {
+            const_iterator old(*this);
+            --index_;
+            return old;
+        }
+
+    // Arithmetic operators
+
+        const_iterator operator+(const difference_type i) const
+        {
+            return const_iterator(base_, index_+i);
+        }
+        const_iterator operator-(const difference_type i) const
+        {
+            return const_iterator(base_, index_-i);
+        }
+        difference_type operator-(const const_iterator& it) const
+        {
+            return (index_ - it.index_);
+        }
+
+    // Comparison operators
+
+        bool operator==(const const_iterator& it) const
+        {
+            return (index_ == it.index_);
+        }
+        bool operator!=(const const_iterator& it) const
+        {
+            return (index_ != it.index_);
+        }
+        bool operator<(const const_iterator& it) const
+        {
+            return (index_ < it.index_);
+        }
+        bool operator<=(const const_iterator& it) const
+        {
+            return (index_ <= it.index_);
+        }
+        bool operator>(const const_iterator& it) const
+        {
+            return (index_ > it.index_);
+        }
+        bool operator>=(const const_iterator& it) const
+        {
+            return (index_ >= it.index_);
+        }
+
+    // Compound assignment 
+
+        const_iterator& operator+=(const difference_type n)
+        {
+            index_ += n;
+            return *this;
+        }
+        const_iterator& operator-=(const difference_type n)
+        {
+            index_ -= n;
+            return *this;
+        }
+    };
+    auto cbegin() const
+    {
+        return const_iterator(*this, 0);
+    }
+    auto cend() const
+    {
+        return const_iterator(*this, this->size());
+    }
+};
+
+
+template <class V>
+class ListsRefWrap
+:
+    public ListExpression<ListsRefWrap<V>, typename V::value_type>//,T>, T>
+{
+    List<V*> data;
+    List<size_t> offsets;
+    size_t total_size = 0;
+
+    void init()
+    {
+        offsets.resize_nocopy(data.size());
+        total_size = 0;
+        for (label i = 0; i < data.size(); ++i)
+        {
+            offsets[i] = total_size;
+            total_size += data[i]->size();
+        }
+    }
+    labelPair whichVecElem(const label i) const
+    {
+        const label vec_index =
+        std::upper_bound
+        (
+            offsets.begin(),
+            offsets.end(),
+            i
+        )
+        - offsets.begin() - 1;
+        const label elem_index = i - offsets[vec_index];
+        return labelPair(vec_index, elem_index);
+    }
+
+public:
+    static constexpr bool is_leaf = false;  //true;
+    
+    using value_type = typename V::value_type;  //T;
+    typedef ListsRefWrap<V> this_type;
+
+    ListsRefWrap() = delete;
+    
+    ListsRefWrap(UList<V*>& d)
+    :
+        data(d)
+    {
+        init();
+    }
+    ListsRefWrap(List<V*>&& d)
+    :
+        data(d)
+    {
+        init();
+    }
+    template<class Container>
+    ListsRefWrap(Container& ds)
+    {
+        data.resize_nocopy(ds.size());
+        forAll(ds, i)
+        {
+            V& vals = static_cast<V&>(ds[i]);
+            data[i] = &vals;
+        }
+        init();
+    }
+    template<class Container>
+    ListsRefWrap(Container& ds, const labelUList& elems)
+    {
+        data.resize_nocopy(elems.size());
+        forAll(ds, i)
+        {
+            V& vals = static_cast<V&>(ds[i]);
+            data[i] = &vals;
+        }
+        init();
+    }
+
+    void push_back(V* vec)
+    {
+        data.push_back(vec);
+        offsets.push_back(total_size);
+        total_size += vec->size();
+    }
+
+    //- Dummy resizing so we can pass *this to
+    //- ListExpression::evaluate. Maybe do resizing externally?
+    void resize_nocopy(const label)
+    {}
+    //- Assignment
+    template<typename E>
+    void operator=(const ListExpression<E, typename E::value_type>& expr)
+    {
+        expr.evaluate(*this);
+    }
+
+    value_type operator[](const label i) const
+    {
+        const auto ve(whichVecElem(i));
+        return (*data[ve.first()])[ve.second()];
+    }
+    value_type& operator[](const label i)
+    {    
+        const auto ve(whichVecElem(i));
+        auto& elems = *data[ve.first()];
+        return elems[ve.second()];
+    }
+    label size() const { return total_size; }
+
+
+    typedef value_type* pointer;
+    typedef value_type& reference;
+    typedef const value_type* const_pointer;
+    typedef const value_type& const_reference;
+    typedef label difference_type;
+
+    struct const_iterator
+    {
+        using iterator_category = std::random_access_iterator_tag;
+        using difference_type   = this_type::difference_type;
+        using value_type  = this_type::value_type;
+        using pointer     = this_type::const_pointer;
+        using reference   = this_type::const_reference;
+
+        const this_type& base_; 
+        label index_;
+
+        const_iterator(const this_type& base, const label index)
+        :
+            base_(base),
+            index_(index)
+        {}
+        const_iterator(const const_iterator& it)
+        :
+            base_(it.base_),
+            index_(it.index_)
+        {}
+        const_iterator(const_iterator& it)
+        :
+            base_(it.base_),
+            index_(it.index_)
+        {}
+        const_iterator() = delete;  //: base_(0) {}
+        const_iterator& operator=(const_iterator& it) = default;
+
+    // Indexing operators
+
+        value_type operator*() const
+        {
+             return base_[index_];
+        }
+        //value_type operator->() const
+        //{
+        //    return *base_;
+        //}
+        value_type operator[](const difference_type i) const
+        {
+           return base_[index_+i];
+        }
+
+    // Increment/Decrement
+
+        const_iterator& operator++()
+        {
+            ++index_;
+            return *this;
+        }
+        const_iterator operator++(int)
+        {
+            const_iterator old(*this);
+            ++index_;
+            return old;
+        }
+        const_iterator& operator--()
+        {
+            --index_;
+            return *this;
+        }
+        const_iterator operator--(int)
+        {
+            const_iterator old(*this);
+            --index_;
+            return old;
+        }
+
+    // Arithmetic operators
+
+        const_iterator operator+(const difference_type i) const
+        {
+            return const_iterator(base_, index_+i);
+        }
+        const_iterator operator-(const difference_type i) const
+        {
+            return const_iterator(base_, index_-i);
+        }
+        difference_type operator-(const const_iterator& it) const
+        {
+            return (index_ - it.index_);
+        }
+
+    // Comparison operators
+
+        bool operator==(const const_iterator& it) const
+        {
+            return (index_ == it.index_);
+        }
+        bool operator!=(const const_iterator& it) const
+        {
+            return (index_ != it.index_);
+        }
+        bool operator<(const const_iterator& it) const
+        {
+            return (index_ < it.index_);
+        }
+        bool operator<=(const const_iterator& it) const
+        {
+            return (index_ <= it.index_);
+        }
+        bool operator>(const const_iterator& it) const
+        {
+            return (index_ > it.index_);
+        }
+        bool operator>=(const const_iterator& it) const
+        {
+            return (index_ >= it.index_);
+        }
+
+    // Compound assignment 
+
+        const_iterator& operator+=(const difference_type n)
+        {
+            index_ += n;
+            return *this;
+        }
+        const_iterator& operator-=(const difference_type n)
+        {
+            index_ -= n;
+            return *this;
+        }
+    };
+    auto cbegin() const
+    {
+        return const_iterator(*this, 0);
+    }
+    auto cend() const
+    {
+        return const_iterator(*this, this->size());
+    }
+
+    struct iterator
+    {
+        using iterator_category = std::random_access_iterator_tag;
+        using difference_type   = this_type::difference_type;
+        using value_type  = this_type::value_type;
+        using pointer     = this_type::pointer;
+        using reference   = this_type::reference;
+
+        this_type& base_; 
+        label index_;
+
+        iterator(this_type& base, const label index)
+        :
+            base_(base),
+            index_(index)
+        {}
+        iterator(const iterator& it) = default;
+        iterator() = delete;
+        iterator& operator=(const iterator& it)
+        {
+            base_ = it.base_;
+            index_ = it.index_;
+            return *this;
+        };
+
+    // Indexing operators
+
+        value_type& operator*()
+        {
+            return base_[index_];
+        }
+        value_type operator*() const
+        {
+            return base_[index_];
+        }
+        value_type operator->() const
+        {
+            return base_[index_];
+        }
+        value_type& operator->()
+        {
+            return base_[index_];
+        }
+        value_type operator[](const difference_type i) const
+        {
+           return base_[index_+i];
+        }
+        value_type& operator[](const difference_type i)
+        {
+           return base_[index_+i];
+        }
+
+    // Increment/Decrement
+
+        iterator& operator++()
+        {
+            ++index_;
+            return *this;
+        }
+        iterator operator++(int)
+        {
+            iterator old(*this);
+            ++index_;
+            return old;
+        }
+        iterator& operator--()
+        {
+            --index_;
+            return *this;
+        }
+        iterator operator--(int)
+        {
+            iterator old(*this);
+            --index_;
+            return old;
+        }
+
+    // Arithmetic operators
+
+        iterator operator+(const difference_type i) const
+        {
+            return iterator(base_, index_+i);
+        }
+        iterator operator-(const difference_type i) const
+        {
+            return iterator(base_, index_-i);
+        }
+        difference_type operator-(const iterator& it) const
+        {
+            return (index_ - it.index_);
+        }
+
+    // Comparison operators
+
+        bool operator==(const iterator& it) const
+        {
+            return (index_ == it.index_);
+        }
+        bool operator!=(const iterator& it) const
+        {
+            return (index_ != it.index_);
+        }
+        bool operator<(const iterator& it) const
+        {
+            return (index_ < it.index_);
+        }
+        bool operator<=(const iterator& it) const
+        {
+            return (index_ <= it.index_);
+        }
+        bool operator>(const iterator& it) const
+        {
+            return (index_ > it.index_);
+        }
+        bool operator>=(const iterator& it) const
+        {
+            return (index_ >= it.index_);
+        }
+
+    // Compound assignment 
+
+        iterator& operator+=(const difference_type n)
+        {
+            index_ += n;
+            return *this;
+        }
+        iterator& operator-=(const difference_type n)
+        {
+            index_ -= n;
+            return *this;
+        }
+    };
+    auto begin()
+    {
+        auto& this_vals = *this;
+        return iterator(this_vals, 0);
+    }
+    auto end()
+    {
+        auto& this_vals = *this;
+        return iterator(this_vals, this->size());
+    }
+};
+
+
+} // End namespace Expression
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/expressionTemplates/boolExpression.H b/src/OpenFOAM/expressionTemplates/boolExpression.H
new file mode 100644
index 0000000000..cceb35d444
--- /dev/null
+++ b/src/OpenFOAM/expressionTemplates/boolExpression.H
@@ -0,0 +1,145 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 M. Janssens
+-------------------------------------------------------------------------------
+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/>.
+
+Namespace
+    Foam::boolExpression
+
+Description
+    Expression templates for bool
+
+SourceFiles
+    boolExpression.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_boolExpression_H
+#define Foam_boolExpression_H
+
+#include "GenericExpression.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace Expression
+{
+
+/*---------------------------------------------------------------------------*\
+                       Namespace Expression Declarations
+\*---------------------------------------------------------------------------*/
+
+// Wrap of non-const bool
+// ~~~~~~~~~~~~~~~~~~~~~~
+
+class boolWrap
+:
+    public GenericExpression<boolWrap>
+{
+    bool elems_;
+
+public:
+    static constexpr bool is_leaf = false;  //true;
+
+    //- Construct from components
+    boolWrap(const bool elems)
+    :
+        elems_(elems)
+    {}
+
+    // Construct from GenericExpression, forcing its evaluation.
+    template<typename E>
+    boolWrap
+    (
+        const GenericExpression<E>& expr
+    )
+    :
+        elems_(expr.evaluate())
+    {}
+
+    auto evaluate() const
+    {
+        return elems_;
+    }
+
+    template<typename E>
+    auto evaluate
+    (
+        const GenericExpression<E>& expr
+    )
+    {
+        elems_ = expr.evaluate();
+        return elems_;
+    }
+
+    //- Assignment
+    template<typename E>
+    void operator=
+    (
+        const GenericExpression<E>& expr
+    )
+    {
+        elems_ = expr.evaluate();
+    }
+};
+
+
+// Wrap of const bool
+// ~~~~~~~~~~~~~~~~~~
+
+class boolConstWrap
+:
+    public GenericExpression<boolConstWrap>
+{
+    const bool elems_;
+
+public:
+
+    // ! Store as copy (since holds reference)
+    static constexpr bool is_leaf = false;
+
+    // construct from components
+    boolConstWrap(const bool elems)
+    :
+        elems_(elems)
+    {}
+
+    auto evaluate() const
+    {
+        return elems_;
+    }
+};
+
+
+} // End namespace Expression
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/expressionTemplates/dimensionSetExpression.H b/src/OpenFOAM/expressionTemplates/dimensionSetExpression.H
new file mode 100644
index 0000000000..09f395860a
--- /dev/null
+++ b/src/OpenFOAM/expressionTemplates/dimensionSetExpression.H
@@ -0,0 +1,148 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 M. Janssens
+-------------------------------------------------------------------------------
+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/>.
+
+Namespace
+    Foam::dimensionSetExpression
+
+Description
+    Expression templates for dimensionSet
+
+SourceFiles
+    dimensionSetExpression.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_dimensionSetExpression_H
+#define Foam_dimensionSetExpression_H
+
+#include "GenericExpression.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace Expression
+{
+
+/*---------------------------------------------------------------------------*\
+                       Namespace Expression Declarations
+\*---------------------------------------------------------------------------*/
+
+// Wrap of non-const reference to dimensionSet
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+class dimensionSetRefWrap
+:
+    public GenericExpression<dimensionSetRefWrap>
+{
+    dimensionSet& elems_;
+
+public:
+    static constexpr bool is_leaf = false;  //true;
+
+    //- Construct from components
+    dimensionSetRefWrap(dimensionSet& elems)
+    :
+        elems_(elems)
+    {}
+
+    // Construct from GenericExpression, forcing its evaluation.
+    template<typename E>
+    dimensionSetRefWrap
+    (
+        dimensionSet& elems,
+        const GenericExpression<E>& expr
+    )
+    :
+        elems_(elems)
+    {
+        elems_ = expr.evaluate();
+    }
+
+    auto evaluate() const
+    {
+        return elems_;
+    }
+
+    template<typename E>
+    auto& evaluate
+    (
+        const GenericExpression<E>& expr
+    )
+    {
+        elems_ = expr.evaluate();
+        return elems_;
+    }
+
+    //- Assignment
+    template<typename E>
+    void operator=
+    (
+        const GenericExpression<E>& expr
+    )
+    {
+        elems_ = expr.evaluate();
+    }
+};
+
+
+// Wrap of const reference to dimensionSet
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+class dimensionSetConstRefWrap
+:
+    public GenericExpression<dimensionSetConstRefWrap>
+{
+    const dimensionSet& elems_;
+
+public:
+
+    // ! Store as copy (since holds reference)
+    static constexpr bool is_leaf = false;
+
+    // construct from components
+    dimensionSetConstRefWrap(const dimensionSet& elems)
+    :
+        elems_(elems)
+    {}
+
+    auto evaluate() const
+    {
+        return elems_;
+    }
+};
+
+
+} // End namespace Expression
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/Fields/Field/Field.H b/src/OpenFOAM/fields/Fields/Field/Field.H
index 198ee76dd2..988b36faee 100644
--- a/src/OpenFOAM/fields/Fields/Field/Field.H
+++ b/src/OpenFOAM/fields/Fields/Field/Field.H
@@ -582,6 +582,37 @@ public:
 
         friend Ostream& operator<< <Type>
         (Ostream&, const tmp<Field<Type>>&);
+
+
+    // Expression templates
+
+        //- Construct from value expression
+        template<typename E>
+        explicit Field
+        (
+            const Expression::ListExpression
+            <
+                E,
+                typename E::value_type
+            >& expr
+        )
+        {
+            expr.evaluate(*this);
+        }
+
+        //- Assign values from expression
+        template<typename E>
+        void operator=
+        (
+            const Expression::ListExpression
+            <
+                E,
+                typename E::value_type
+            >& expr
+        )
+        {
+            expr.evaluate(*this);
+        }
 };
 
 
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
index 69d9dc30f3..a8bc58a252 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
@@ -1675,6 +1675,124 @@ Foam::Ostream& Foam::operator<<
 }
 
 
+// * * * * * * * * * * * * * Expression Templates  * * * * * * * * * * * * * //
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+template<typename E>
+Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
+(
+    const IOobject& io,
+    const Mesh& mesh,
+    const Expression::GeometricFieldExpression
+    <
+        E,
+        typename E::IntExpr,
+        typename E::UncoupledPatchExpr,
+        typename E::CoupledPatchExpr,
+        typename E::value_type
+    >& expr
+)
+:
+    Internal(io, mesh, expr.dimensions(), false),
+    timeIndex_(this->time().timeIndex()),
+    boundaryField_(mesh.boundary(), *this, PatchField<Type>::calculatedType())
+{
+    DebugInFunction
+        << "Creating from expression " << nl << this->info() << endl;
+
+    bool hasRead = readIfPresent();
+    if (!hasRead)
+    {
+        expr.evaluate(*this);
+    }
+}
+
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+template<typename E>
+Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
+(
+    const word& name,
+    const Mesh& mesh,
+    const Expression::GeometricFieldExpression
+    <
+        E,
+        typename E::IntExpr,
+        typename E::UncoupledPatchExpr,
+        typename E::CoupledPatchExpr,
+        typename E::value_type
+    >& expr
+)
+:
+    Internal
+    (
+        IOobject
+        (
+            name,
+            mesh.time().timeName(),
+            mesh.thisDb()
+        ),
+        mesh,
+        expr.dimensions(),
+        false
+    ),
+    timeIndex_(this->time().timeIndex()),
+    boundaryField_(mesh.boundary(), *this, PatchField<Type>::calculatedType())
+{
+    DebugInFunction
+        << "Creating from expression " << nl << this->info() << endl;
+
+    expr.evaluate(*this);
+}
+
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+Foam::Expression::GeometricFieldConstRefWrap
+<
+    Foam::GeometricField<Type, PatchField, GeoMesh>
+>
+Foam::GeometricField<Type, PatchField, GeoMesh>::expr() const
+{
+    return Expression::GeometricFieldConstRefWrap<this_type>(*this);
+}
+
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+template<typename E>
+void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=
+(
+    const Expression::GeometricFieldExpression
+    <
+        E,
+        typename E::IntExpr,
+        typename E::UncoupledPatchExpr,
+        typename E::CoupledPatchExpr,
+        typename E::value_type
+    >& expr
+)
+{
+    Expression::GeometricFieldRefWrap<this_type>(*this, expr);
+}
+
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+template<typename E>
+void Foam::GeometricField<Type, PatchField, GeoMesh>::operator==
+(
+    const Expression::GeometricFieldExpression
+    <
+        E,
+        typename E::IntExpr,
+        typename E::UncoupledPatchExpr,
+        typename E::CoupledPatchExpr,
+        typename E::value_type
+    >& expr
+)
+{
+    expr.evaluate(*this, true);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #undef checkField
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
index 9f42d3f79f..05e882ee0c 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
@@ -42,6 +42,7 @@ SourceFiles
 #define Foam_GeometricField_H
 
 #include "GeometricBoundaryField.H"
+#include "GeometricFieldExpression.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -1142,6 +1143,79 @@ public:
         {
             return this->clamp_range(lo.value(), hi.value());
         }
+
+
+    // Expression templates
+
+        //- Have unique tag
+        using is_GeometricField = void;
+
+        //- Construct from IOobject and value expression. Supports
+        //  read_if_present
+        template<typename E>
+        GeometricField
+        (
+            const IOobject& io,
+            const Mesh& mesh,
+            const Expression::GeometricFieldExpression
+            <
+                E,
+                typename E::IntExpr,
+                typename E::UncoupledPatchExpr,
+                typename E::CoupledPatchExpr,
+                typename E::value_type
+            >& expr
+        );
+
+        //- Construct from name, mesh and value expression.
+        template<typename E>
+        GeometricField
+        (
+            const word& name,
+            const Mesh& mesh,
+            const Expression::GeometricFieldExpression
+            <
+                E,
+                typename E::IntExpr,
+                typename E::UncoupledPatchExpr,
+                typename E::CoupledPatchExpr,
+                typename E::value_type
+            >& expr
+        );
+
+        //- Wrap value as expression
+        typename Expression::GeometricFieldConstRefWrap
+        <
+            this_type
+        > expr() const;
+
+        //- Assign values from expression
+        template<typename E>
+        void operator=
+        (
+            const Expression::GeometricFieldExpression
+            <
+                E,
+                typename E::IntExpr,
+                typename E::UncoupledPatchExpr,
+                typename E::CoupledPatchExpr,
+                typename E::value_type
+            >& expr
+        );
+
+        //- Assign values from expression. Override asssignable.
+        template<typename E>
+        void operator==
+        (
+            const Expression::GeometricFieldExpression
+            <
+                E,
+                typename E::IntExpr,
+                typename E::UncoupledPatchExpr,
+                typename E::CoupledPatchExpr,
+                typename E::value_type
+            >& expr
+        );
 };
 
 
diff --git a/src/OpenFOAM/memory/tmp/tmp.H b/src/OpenFOAM/memory/tmp/tmp.H
index b35587d2e8..bc0b8c9337 100644
--- a/src/OpenFOAM/memory/tmp/tmp.H
+++ b/src/OpenFOAM/memory/tmp/tmp.H
@@ -59,8 +59,17 @@ See also
 namespace Foam
 {
 
+namespace Expression
+{
+template<class GeoField> class GeometricFieldConstTmpWrap;
+template<class T> class ListConstTmpWrap;
+}
+
 // Forward Declarations
 template<class T> class refPtr;
+template<class T> class tmp;
+template<class Type, template<class> class PatchField, class GeoMesh>
+class GeometricField;
 
 /*---------------------------------------------------------------------------*\
                              Class tmp Declaration
@@ -352,6 +361,33 @@ public:
         //  \deprecated(2020-07) - use bool operator
         FOAM_DEPRECATED_FOR(2020-07, "bool operator")
         bool empty() const noexcept { return !ptr_; }
+
+
+    // Expression templates
+
+        //- Wrap value as expression
+        auto expr() const
+        {
+            // Add 'using is_GeometricField = void' to GeometricField
+            // Add 'using is_List = void' to List
+
+            if constexpr
+            (
+                std::is_void_v<typename T::is_GeometricField>
+            )
+            {
+                return Expression::GeometricFieldConstTmpWrap<T>(*this);
+            }
+            else //if constexpr (std::is_void_v<T::is_List>)
+            {
+                return Expression::ListConstTmpWrap<T>(*this);
+            }
+            //else
+            //{
+            //    //static_assert("XXX",YYY);
+            //    return;
+            //}
+        }
 };
 
 
diff --git a/src/OpenFOAM/memory/tmp/tmpI.H b/src/OpenFOAM/memory/tmp/tmpI.H
index 19744903e5..37a1d05d23 100644
--- a/src/OpenFOAM/memory/tmp/tmpI.H
+++ b/src/OpenFOAM/memory/tmp/tmpI.H
@@ -131,7 +131,7 @@ inline Foam::tmp<T>::tmp(const tmp<T>& rhs)
         if (ptr_)
         {
             ptr_->refCount::operator++();
-            this->checkUseCount();
+            //this->checkUseCount();
         }
         else
         {
@@ -162,7 +162,7 @@ inline Foam::tmp<T>::tmp(const tmp<T>& rhs, bool reuse)
             else
             {
                 ptr_->refCount::operator++();
-                this->checkUseCount();
+                //this->checkUseCount();
             }
         }
         else
diff --git a/src/finiteVolume/expressionTemplates/fvMatrixExpression.H b/src/finiteVolume/expressionTemplates/fvMatrixExpression.H
new file mode 100644
index 0000000000..35d41c61f6
--- /dev/null
+++ b/src/finiteVolume/expressionTemplates/fvMatrixExpression.H
@@ -0,0 +1,1747 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 M. Janssens
+-------------------------------------------------------------------------------
+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/>.
+
+Namespace
+    Foam::Expression
+
+Description
+    Expression templates for fvMatrix.
+
+    TBD:
+    - useImplicit_ etc.
+    - lduMatrix::lowerCSR ?
+
+SourceFiles
+    fvMatrixExpression.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_fvMatrixExpression_H
+#define Foam_fvMatrixExpression_H
+
+#include "ListExpression.H"
+#include "GeometricFieldExpression.H"
+//#include <boost/core/demangle.hpp>
+//#include <typeinfo>
+//#include <iostream>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace Expression
+{
+
+/*---------------------------------------------------------------------------*\
+                   Namespace Expression Declarations
+\*---------------------------------------------------------------------------*/
+
+// Expression of fvMatrix
+// ~~~~~~~~~~~~~~~~~~~~~~
+// Bit different:
+// - immediately stores state (hasDiag, hasUpper etc.)
+// - delays execution of building coefficients
+// - big problem: individual components can be built differently so
+//   different template argument ... Needs to be templated on all outputs
+//   separately!
+template
+<
+    typename E,
+    typename DiagExpr,
+    typename UpperExpr,
+    typename LowerExpr,
+    typename FaceFluxExpr,
+    typename SourceExpr
+>
+class fvMatrixExpression
+{
+protected:
+
+    const bool hasDiag_;
+    const bool hasUpper_;
+    const bool hasLower_;
+    const bool hasFaceFluxCorrection_;
+    dimensionSet dimensions_;
+
+public:
+    static constexpr bool is_leaf = false;
+
+    // Constructors
+
+        fvMatrixExpression
+        (
+            const bool hasDiag,
+            const bool hasUpper,
+            const bool hasLower,
+            const bool hasFaceFluxCorrection,
+            const dimensionSet& dimensions
+        )
+        :
+            hasDiag_(hasDiag),
+            hasUpper_(hasUpper),
+            hasLower_(hasLower),
+            hasFaceFluxCorrection_(hasFaceFluxCorrection),
+            dimensions_(dimensions)
+        {
+            //Pout<< "fvMatrixExpression :" << nl
+            //    << "    E        : "
+            //    << boost::core::demangle(typeid(E).name())
+            //    << nl
+            //    << "    DiagExpr : "
+            //    << boost::core::demangle(typeid(DiagExpr).name())
+            //    << endl;
+        }
+
+
+    // lduMatrix part
+
+        DiagExpr diag() const
+        {
+            return static_cast<const E&>(*this).diag();
+        }
+
+        DiagExpr diag()
+        {
+            return static_cast<E&>(*this).diag();
+        }
+
+        UpperExpr upper() const
+        {
+            return static_cast<const E&>(*this).upper();
+        }
+
+        UpperExpr upper()
+        {
+            return static_cast<E&>(*this).upper();
+        }
+
+        LowerExpr lower() const
+        {
+            return static_cast<const E&>(*this).lower();
+        }
+
+        LowerExpr lower()
+        {
+            return static_cast<E&>(*this).lower();
+        }
+
+        bool hasDiag() const
+        {
+             return hasDiag_;
+        }
+        bool hasUpper() const
+        {
+             return hasUpper_;
+        }
+        bool hasLower() const
+        {
+             return hasLower_;
+        }
+        bool diagonal() const
+        {
+             return (hasDiag() && !hasLower() && !hasUpper());
+        }
+        bool symmetric() const
+        {
+             return (hasDiag() && !hasLower() && hasUpper());
+        }
+        bool asymmetric() const
+        {
+             return (hasDiag() && hasLower() && hasUpper());
+        }
+        bool hasFaceFluxCorrection() const
+        {
+             return hasFaceFluxCorrection_;
+        }
+
+
+    // fvMatrix
+
+        SourceExpr source()
+        {
+            return static_cast<E&>(*this).source();
+        }
+
+        SourceExpr source() const
+        {
+            return static_cast<const E&>(*this).source();
+        }
+
+        UpperExpr internalCoeffs(const label i)
+        {
+            return static_cast<E&>(*this).internalCoeffs(i);
+        }
+
+        UpperExpr internalCoeffs(const label i) const
+        {
+            return static_cast<const E&>(*this).internalCoeffs(i);
+        }
+
+        UpperExpr boundaryCoeffs(const label i)
+        {
+            return static_cast<E&>(*this).boundaryCoeffs(i);
+        }
+
+        UpperExpr boundaryCoeffs(const label i) const
+        {
+            return static_cast<const E&>(*this).boundaryCoeffs(i);
+        }
+
+        FaceFluxExpr faceFluxCorrection() const
+        {
+            return static_cast<const E&>(*this).faceFluxCorrection();
+        }
+
+        FaceFluxExpr faceFluxCorrection()
+        {
+            return static_cast<E&>(*this).faceFluxCorrection();
+        }
+
+        const dimensionSet& dimensions() const noexcept
+        {
+            return dimensions_;
+        }
+
+
+    // Other
+
+        //- Helper to evaluate (=construct) an fvMatrix
+        template<class Matrix>
+        Matrix& evaluate(Matrix& m) const
+        {
+            if (hasDiag())
+            {
+                // Assign to diag. Also sets hasDiag.
+                diag().evaluate(m.diag());
+            }
+            if (hasUpper())
+            {
+                upper().evaluate(m.upper());
+            }
+            if (hasLower())
+            {
+                lower().evaluate(m.lower());
+            }
+            if (hasFaceFluxCorrection())
+            {
+                faceFluxCorrection().evaluate(*m.faceFluxCorrectionPtr());
+            }
+
+            // Do boundary
+            auto& intCoeffs = m.internalCoeffs();
+            auto& bouCoeffs = m.boundaryCoeffs();
+            const label n = intCoeffs.size();
+            for (label i = 0; i < n; ++i)
+            {
+                if (intCoeffs.set(i) && intCoeffs[i].size())
+                {
+                    const auto intExpr = internalCoeffs(i);
+                    intExpr.evaluate(intCoeffs[i]);
+                }
+                if (bouCoeffs.set(i) && bouCoeffs[i].size())
+                {
+                    const auto bouExpr = boundaryCoeffs(i);
+                    bouExpr.evaluate(bouCoeffs[i]);
+                }
+            }
+
+            const_cast<dimensionSet&>(m.dimensions()) = dimensions();
+
+            return m;
+        }
+};
+
+
+template<class Matrix>
+class fvMatrixRefWrap
+:
+    public fvMatrixExpression
+    <
+        fvMatrixRefWrap<Matrix>,
+        ListRefWrap<typename Matrix::psiFieldType::value_type>, // diag
+        ListRefWrap<typename Matrix::psiFieldType::value_type>, // upper
+        ListRefWrap<typename Matrix::psiFieldType::value_type>, // lower
+        //ListRefWrap<typename Matrix::psiFieldType::value_type>, // faceflux
+        GeometricFieldRefWrap<typename Matrix::faceFluxFieldType>,
+        ListRefWrap<typename Matrix::psiFieldType::value_type>  // source
+    >
+{
+public:
+
+    static constexpr bool is_leaf = false;  //true;
+
+    //- The fvMatrix type
+    typedef Matrix this_type;
+
+    //- Type to return for internal field
+    typedef typename Matrix::psiFieldType::value_type value_type;
+
+    //- Type to return for containers
+    typedef ListRefWrap<value_type> DiagExpr;
+    typedef ListRefWrap<value_type> UpperExpr;
+    typedef ListRefWrap<value_type> LowerExpr;
+    //typedef ListRefWrap<value_type> FaceFluxExpr;
+    typedef GeometricFieldRefWrap
+    <
+        typename Matrix::faceFluxFieldType
+    > FaceFluxExpr;
+    typedef ListRefWrap<value_type> SourceExpr;
+
+
+private:
+
+    this_type& elems_;
+
+
+public:
+
+    //- Copy construct
+    fvMatrixRefWrap(this_type& elems)
+    :
+        fvMatrixExpression
+        <
+            fvMatrixRefWrap<Matrix>,
+            DiagExpr,
+            UpperExpr,
+            LowerExpr,
+            FaceFluxExpr,
+            SourceExpr
+        >
+        (
+            elems.hasDiag(),
+            elems.hasUpper(),
+            elems.hasLower(),
+            elems.hasFaceFluxCorrection(),
+            elems.dimensions()
+        ),
+        elems_(elems)
+    {}
+
+    //- Move construct
+    fvMatrixRefWrap(this_type&& elems)
+    :
+        fvMatrixExpression
+        <
+            fvMatrixRefWrap<Matrix>,
+            DiagExpr,
+            UpperExpr,
+            LowerExpr,
+            FaceFluxExpr,
+            SourceExpr
+        >
+        (
+            elems.hasDiag(),
+            elems.hasUpper(),
+            elems.hasLower(),
+            elems.hasFaceFluxCorrection(),
+            elems.dimensions()
+        ),
+        elems_(elems)
+    {}
+
+    // Construct from Matrix, forcing its evaluation.
+    template<typename E>
+    fvMatrixRefWrap
+    (
+        this_type& elems,
+        const fvMatrixExpression
+        <
+            E,
+            typename E::DiagExpr,
+            typename E::UpperExpr,
+            typename E::LowerExpr,
+            typename E::FaceFluxExpr,
+            typename E::SourceExpr
+        >& expr
+    )
+    :
+        fvMatrixExpression
+        <
+            fvMatrixRefWrap<Matrix>,
+            DiagExpr,
+            UpperExpr,
+            LowerExpr,
+            FaceFluxExpr,
+            SourceExpr
+        >
+        (
+            expr.hasDiag(),
+            expr.hasUpper(),
+            expr.hasLower(),
+            elems.hasFaceFluxCorrection(),
+            elems.dimensions()
+        ),
+        elems_(elems)
+    {
+        expr.evaluate(elems_);
+    }
+
+    //- Assignment
+    template<typename E>
+    void operator=
+    (
+        const fvMatrixExpression
+        <
+            E,
+            typename E::DiagExpr,
+            typename E::UpperExpr,
+            typename E::LowerExpr,
+            typename E::FaceFluxExpr,
+            typename E::SourceExpr
+        >& expr
+    )
+    {
+        expr.evaluate(elems_);
+    }
+
+    //- Evaluate and return as GeoField. Rename to evaluate to make it clear
+    //- it takes time? Or leave as indexing for convenience?
+    template<typename E>
+    this_type& evaluate
+    (
+        const fvMatrixExpression
+        <
+            E,
+            typename E::DiagExpr,
+            typename E::UpperExpr,
+            typename E::LowerExpr,
+            typename E::FaceFluxExpr,
+            typename E::SourceExpr
+        >& expr
+    )
+    {
+        return expr.evaluate(elems_);
+    }
+
+
+    // lduMatrix part
+
+        DiagExpr diag() const
+        {
+            return elems_.diag();
+        }
+
+        DiagExpr diag()
+        {
+            return elems_.diag();
+        }
+
+        UpperExpr upper() const
+        {
+            return elems_.upper();
+        }
+
+        UpperExpr upper()
+        {
+            return elems_.upper();
+        }
+
+        LowerExpr lower() const
+        {
+            return elems_.lower();
+        }
+
+        LowerExpr lower()
+        {
+            return elems_.lower();
+        }
+
+
+    // fvMatrix
+
+        SourceExpr source() const
+        {
+            return elems_.source();
+        }
+
+        SourceExpr source()
+        {
+            return elems_.source();
+        }
+
+        UpperExpr internalCoeffs(const label i)
+        {
+            return
+            (
+                elems_.internalCoeffs().set(i)
+              ? elems_.internalCoeffs()[i]
+              : List<value_type>::null()
+            );
+        }
+
+        UpperExpr internalCoeffs(const label i) const
+        {
+            return
+            (
+                elems_.internalCoeffs().set(i)
+              ? elems_.internalCoeffs()[i]
+              : List<value_type>::null()
+            );
+        }
+
+        UpperExpr boundaryCoeffs(const label i)
+        {
+            return
+            (
+                elems_.boundaryCoeffs().set(i)
+              ? elems_.boundaryCoeffs()[i]
+              : List<value_type>::null()
+            );
+        }
+
+        UpperExpr boundaryCoeffs(const label i) const
+        {
+            return
+            (
+                elems_.boundaryCoeffs().set(i)
+              ? elems_.boundaryCoeffs()[i]
+              : List<value_type>::null()
+            );
+        }
+
+        FaceFluxExpr faceFluxCorrection() const
+        {
+            return *const_cast<this_type&>
+            (
+                elems_
+            ).faceFluxCorrectionPtr().get();
+        }
+
+        FaceFluxExpr faceFluxCorrection()
+        {
+            return *elems_.faceFluxCorrectionPtr().get();
+        }
+};
+
+
+template<class Matrix>
+class fvMatrixConstRefWrap
+:
+    public fvMatrixExpression
+    <
+        fvMatrixConstRefWrap<Matrix>,
+        ListConstRefWrap<typename Matrix::psiFieldType::value_type>, // diag
+        ListConstRefWrap<typename Matrix::psiFieldType::value_type>, // upper
+        ListConstRefWrap<typename Matrix::psiFieldType::value_type>, // lower
+        //ListConstRefWrap<typename Matrix::psiFieldType::value_type>, // faceflux
+        GeometricFieldConstRefWrap<typename Matrix::faceFluxFieldType>,
+        ListConstRefWrap<typename Matrix::psiFieldType::value_type>  // source
+    >
+{
+public:
+
+    static constexpr bool is_leaf = false;  //true;
+
+    //- The fvMatrix type
+    typedef Matrix this_type;
+
+    //- Type to return for internal field
+    typedef typename Matrix::psiFieldType::value_type value_type;
+
+    //- Type to return for containers
+    typedef ListConstRefWrap<value_type> DiagExpr;
+    typedef ListConstRefWrap<value_type> UpperExpr;
+    typedef ListConstRefWrap<value_type> LowerExpr;
+    //typedef ListConstRefWrap<value_type> FaceFluxExpr;
+    typedef GeometricFieldConstRefWrap
+    <
+        typename Matrix::faceFluxFieldType
+    > FaceFluxExpr;
+    typedef ListConstRefWrap<value_type> SourceExpr;
+
+
+private:
+
+    const this_type& elems_;
+
+
+public:
+
+    // Construct from components
+    fvMatrixConstRefWrap(const this_type& elems)
+    :
+        fvMatrixExpression
+        <
+            fvMatrixConstRefWrap<Matrix>,
+            DiagExpr,
+            UpperExpr,
+            LowerExpr,
+            FaceFluxExpr,
+            SourceExpr
+        >
+        (
+            elems.hasDiag(),
+            elems.hasUpper(),
+            elems.hasLower(),
+            elems.hasFaceFluxCorrection(),
+            elems.dimensions()
+        ),
+        elems_(elems)
+    {}
+
+
+    // lduMatrix part
+
+        DiagExpr diag() const
+        {
+            return elems_.diag();
+        }
+
+        UpperExpr upper() const
+        {
+            return elems_.upper();
+        }
+
+        LowerExpr lower() const
+        {
+            return elems_.lower();
+        }
+
+
+    // fvMatrix
+
+        SourceExpr source() const
+        {
+            return elems_.source();
+        }
+
+        UpperExpr internalCoeffs(const label i) const
+        {
+            return
+            (
+                elems_.internalCoeffs().set(i)
+              ? elems_.internalCoeffs()[i]
+              : List<value_type>::null()
+            );
+        }
+
+        UpperExpr boundaryCoeffs(const label i) const
+        {
+            return
+            (
+                elems_.boundaryCoeffs().set(i)
+              ? elems_.boundaryCoeffs()[i]
+              : List<value_type>::null()
+            );
+        }
+
+        FaceFluxExpr faceFluxCorrection() const
+        {
+            return *const_cast<this_type&>
+            (
+                elems_
+            ).faceFluxCorrectionPtr().get();
+        }
+};
+
+
+// Expressions on fvMatrices
+// ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+// operator+
+template<typename E1, typename E2>
+class fvm_add
+:
+    public fvMatrixExpression
+    <
+        fvm_add<E1, E2>,
+        List_add                    // DiagExpr
+        <
+            typename E1::DiagExpr,
+            typename E2::DiagExpr
+        >,
+        List_add                    // UpperExpr
+        <
+            typename E1::UpperExpr,
+            typename E2::UpperExpr
+        >,
+        List_add                    // LowerExpr
+        <
+            typename E1::LowerExpr,
+            typename E2::LowerExpr
+        >,
+        GF_add                    // FaceFluxExpr
+        <
+            typename E1::FaceFluxExpr,
+            typename E2::FaceFluxExpr
+        >,
+        List_add                    // SourceExpr
+        <
+            typename E1::SourceExpr,
+            typename E2::SourceExpr
+        >
+    >
+{
+    // cref if leaf, copy otherwise
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;
+    typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_;
+
+public:
+    static constexpr bool is_leaf = false;
+
+    //- Type to return for internal field
+    typedef typename E1::value_type value_type;
+
+    //- Type to return for containers
+    typedef List_add
+    <
+        typename E1::DiagExpr,
+        typename E2::DiagExpr
+    > DiagExpr;
+    typedef List_add
+    <
+        typename E1::UpperExpr,
+        typename E2::UpperExpr
+    > UpperExpr;
+    typedef List_add
+    <
+        typename E1::LowerExpr,
+        typename E2::LowerExpr
+    > LowerExpr;
+    typedef GF_add
+    <
+        typename E1::FaceFluxExpr,
+        typename E2::FaceFluxExpr
+    > FaceFluxExpr;
+    typedef List_add
+    <
+        typename E1::SourceExpr,
+        typename E2::SourceExpr
+    > SourceExpr;
+
+    fvm_add(const E1& u, const E2& v)
+    :
+        fvMatrixExpression
+        <
+            fvm_add<E1, E2>,
+            DiagExpr,
+            UpperExpr,
+            LowerExpr,
+            FaceFluxExpr,
+            SourceExpr
+        >
+        (
+            u.hasDiag() || v.hasDiag(),
+            u.hasUpper() || v.hasUpper(),
+            u.hasLower() || v.hasLower(),
+            u.hasFaceFluxCorrection() || v.hasFaceFluxCorrection(),
+            u.dimensions()
+        ),
+        u_(u),
+        v_(v)
+    {
+        //Pout<< "fvm_add :" << nl
+        //    << "    u_ : " << boost::core::demangle(typeid(u_).name())
+        //    << nl
+        //    << "    v_ : " << boost::core::demangle(typeid(v_).name())
+        //    << endl;
+    }
+
+
+    // lduMatrix part
+
+        DiagExpr diag() const
+        {
+            return u_.diag() + v_.diag();
+        }
+
+        UpperExpr upper() const
+        {
+            return u_.upper() + v_.upper();
+        }
+
+        LowerExpr lower() const
+        {
+            return u_.lower() + v_.lower();
+        }
+
+    // fvMatrix
+
+        SourceExpr source() const
+        {
+            return u_.source() + v_.source();
+        }
+
+        UpperExpr boundaryCoeffs(const label i) const
+        {
+            return u_.boundaryCoeffs(i) + v_.boundaryCoeffs(i);
+        }
+
+        UpperExpr internalCoeffs(const label i) const
+        {
+            return u_.internalCoeffs(i) + v_.internalCoeffs(i);
+        }
+
+        FaceFluxExpr faceFluxCorrection() const
+        {
+            return u_.faceFluxCorrection() + v_.faceFluxCorrection();
+        }
+};
+template<typename E1, typename E2>
+fvm_add<E1, E2>
+operator+
+(
+    fvMatrixExpression
+    <
+        E1,
+        typename E1::DiagExpr,
+        typename E1::UpperExpr,
+        typename E1::LowerExpr,
+        typename E1::FaceFluxExpr,
+        typename E1::SourceExpr
+    > const& u,
+    fvMatrixExpression
+    <
+        E2,
+        typename E2::DiagExpr,
+        typename E2::UpperExpr,
+        typename E2::LowerExpr,
+        typename E2::FaceFluxExpr,
+        typename E2::SourceExpr
+    > const& v
+)
+{
+    return fvm_add<E1, E2>
+    (
+        static_cast<const E1&>(u),
+        static_cast<const E2&>(v)
+    );
+}
+
+
+// operator-, operator==
+template<typename E1, typename E2>
+class fvm_subtract
+:
+    public fvMatrixExpression
+    <
+        fvm_subtract<E1, E2>,
+        List_subtract                   // diag
+        <
+            typename E1::DiagExpr,
+            typename E2::DiagExpr
+        >,
+        List_subtract                   // upper
+        <
+            typename E1::UpperExpr,
+            typename E2::UpperExpr
+        >,
+        List_subtract                   // lower
+        <
+            typename E1::LowerExpr,
+            typename E2::LowerExpr
+        >,
+        GF_subtract                   // faceflux
+        <
+            typename E1::FaceFluxExpr,
+            typename E2::FaceFluxExpr
+        >,
+        List_subtract                   // source
+        <
+            typename E1::SourceExpr,
+            typename E2::SourceExpr
+        >
+    >
+{
+    // cref if leaf, copy otherwise
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;
+    typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_;
+
+public:
+    static constexpr bool is_leaf = false;
+
+    //- Type to return for internal field
+    typedef typename E1::value_type value_type;
+
+    //- Type to return for containers
+    typedef List_subtract
+    <
+        typename E1::DiagExpr,
+        typename E2::DiagExpr
+    > DiagExpr;
+    typedef List_subtract
+    <
+        typename E1::UpperExpr,
+        typename E2::UpperExpr
+    > UpperExpr;
+    typedef List_subtract
+    <
+        typename E1::LowerExpr,
+        typename E2::LowerExpr
+    > LowerExpr;
+    typedef GF_subtract
+    <
+        typename E1::FaceFluxExpr,
+        typename E2::FaceFluxExpr
+    > FaceFluxExpr;
+    typedef List_subtract
+    <
+        typename E1::SourceExpr,
+        typename E2::SourceExpr
+    > SourceExpr;
+
+    fvm_subtract(const E1& u, const E2& v)
+    :
+        fvMatrixExpression
+        <
+            fvm_subtract<E1, E2>,
+            DiagExpr,
+            UpperExpr,
+            LowerExpr,
+            FaceFluxExpr,
+            SourceExpr
+        >
+        (
+            u.hasDiag() || v.hasDiag(),
+            u.hasUpper() || v.hasUpper(),
+            u.hasLower() || v.hasLower(),
+            u.hasFaceFluxCorrection() || v.hasFaceFluxCorrection(),
+            u.dimensions()
+        ),
+        u_(u),
+        v_(v)
+    {}
+
+
+    // lduMatrix part
+
+        DiagExpr diag() const
+        {
+            return u_.diag() - v_.diag();
+        }
+
+        UpperExpr upper() const
+        {
+            return u_.upper() - v_.upper();
+        }
+
+        LowerExpr lower() const
+        {
+            return u_.lower() - v_.lower();
+        }
+
+
+    // fvMatrix
+
+        SourceExpr source() const
+        {
+            return u_.source() - v_.source();
+        }
+
+        UpperExpr internalCoeffs(const label i) const
+        {
+            return u_.internalCoeffs(i) - v_.internalCoeffs(i);
+        }
+
+        UpperExpr boundaryCoeffs(const label i) const
+        {
+            return u_.boundaryCoeffs(i) - v_.boundaryCoeffs(i);
+        }
+
+        FaceFluxExpr faceFluxCorrection() const
+        {
+            return u_.faceFluxCorrection() - v_.faceFluxCorrection();
+        }
+};
+template<typename E1, typename E2>
+fvm_subtract<E1, E2>
+operator-
+(
+    fvMatrixExpression
+    <
+        E1,
+        typename E1::DiagExpr,
+        typename E1::UpperExpr,
+        typename E1::LowerExpr,
+        typename E1::FaceFluxExpr,
+        typename E1::SourceExpr
+    > const& u,
+    fvMatrixExpression
+    <
+        E2,
+        typename E2::DiagExpr,
+        typename E2::UpperExpr,
+        typename E2::LowerExpr,
+        typename E2::FaceFluxExpr,
+        typename E2::SourceExpr
+    > const& v
+)
+{
+    return fvm_subtract<E1, E2>
+    (
+        static_cast<const E1&>(u),
+        static_cast<const E2&>(v)
+    );
+}
+template<typename E1, typename E2>
+fvm_subtract<E1, E2>
+operator==
+(
+    fvMatrixExpression
+    <
+        E1,
+        typename E1::DiagExpr,
+        typename E1::UpperExpr,
+        typename E1::LowerExpr,
+        typename E1::FaceFluxExpr,
+        typename E1::SourceExpr
+    > const& u,
+    fvMatrixExpression
+    <
+        E2,
+        typename E2::DiagExpr,
+        typename E2::UpperExpr,
+        typename E2::LowerExpr,
+        typename E2::FaceFluxExpr,
+        typename E2::SourceExpr
+    > const& v
+)
+{
+    return fvm_subtract<E1, E2>
+    (
+        static_cast<const E1&>(u),
+        static_cast<const E2&>(v)
+    );
+}
+
+
+// Negation : operator-
+template<typename E1>
+class fvm_negate
+:
+    public fvMatrixExpression
+    <
+        fvm_negate<E1>,
+        List_negate<typename E1::DiagExpr>,
+        List_negate<typename E1::UpperExpr>,
+        List_negate<typename E1::LowerExpr>,
+        List_negate<typename E1::FaceFluxExpr>,
+        List_negate<typename E1::SourceExpr>
+    >
+{
+    // cref if leaf, copy otherwise
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;
+
+public:
+    static constexpr bool is_leaf = false;
+
+    //- Type to return for internal field
+    typedef typename E1::value_type value_type;
+
+    //- Type to return for containers
+    typedef List_negate<typename E1::DiagExpr> DiagExpr;
+    typedef List_negate<typename E1::UpperExpr> UpperExpr;
+    typedef List_negate<typename E1::LowerExpr> LowerExpr;
+    typedef List_negate<typename E1::FaceFluxExpr> FaceFluxExpr;
+    typedef List_negate<typename E1::SourceExpr> SourceExpr;
+
+    fvm_negate(const E1& u)
+    :
+        fvMatrixExpression
+        <
+            fvm_negate<E1>,
+            DiagExpr,
+            UpperExpr,
+            LowerExpr,
+            FaceFluxExpr,
+            SourceExpr
+        >
+        (
+            u.hasDiag(),
+            u.hasUpper(),
+            u.hasLower(),
+            u.hasFaceFluxCorrection(),
+            u.dimensions()
+        ),
+        u_(u)
+    {}
+
+
+    // lduMatrix part
+
+        DiagExpr diag() const
+        {
+            return -u_.diag();
+        }
+
+        UpperExpr upper() const
+        {
+            return -u_.upper();
+        }
+
+        LowerExpr lower() const
+        {
+            return -u_.lower();
+        }
+
+
+    // fvMatrix
+
+        SourceExpr source() const
+        {
+            return -u_.source();
+        }
+
+        UpperExpr internalCoeffs(const label i) const
+        {
+            return -u_.internalCoeffs(i);
+        }
+
+        UpperExpr boundaryCoeffs(const label i) const
+        {
+            return -u_.boundaryCoeffs(i);
+        }
+
+        FaceFluxExpr faceFluxCorrection() const
+        {
+            return -u_.faceFluxCorrection();
+        }
+};
+template<typename E1>
+fvm_negate<E1>
+operator-
+(
+    fvMatrixExpression
+    <
+        E1,
+        typename E1::DiagExpr,
+        typename E1::UpperExpr,
+        typename E1::LowerExpr,
+        typename E1::FaceFluxExpr,
+        typename E1::SourceExpr
+    > const& u
+)
+{
+    return fvm_negate<E1>(static_cast<const E1&>(u));
+}
+
+
+//- fvMatrix, internal parts and dimensions only. No upper/lower/boundaryCoeffs
+//  etc. Used for source manipulation.
+template<class E1, class E2>    //, class FaceFluxExpr>
+class fvMatrixInternal
+:
+    public fvMatrixExpression
+    <
+        fvMatrixInternal<E1, E2>,
+        E1, // Diag
+        E1, // Upper, Not used
+        E1, // Lower, Not used
+        GeometricField<typename E1::value_type, fvPatchField, volMesh>, // FaceFlux, Not used
+        E2  // Source
+    >
+{
+public:
+
+    static constexpr bool is_leaf = false;  //true;
+
+    //- Type to return for internal field
+    typedef typename E1::value_type value_type;
+
+    //- Type to return for containers
+    typedef E1 DiagExpr;
+    typedef E1 UpperExpr;   // not yet used
+    typedef E1 LowerExpr;   // not yet used
+    typedef GeometricField<value_type, fvPatchField, volMesh>
+        FaceFluxExpr;   // not yet used
+    typedef E2 SourceExpr;
+
+
+private:
+
+    const E1 diag_;
+    const SourceExpr source_;
+
+
+public:
+
+    // Construct from component expressions
+    fvMatrixInternal
+    (
+        const dimensionSet& dimensions,
+        const DiagExpr& diag,
+        const SourceExpr& source
+    )
+    :
+        fvMatrixExpression
+        <
+            fvMatrixInternal<DiagExpr, SourceExpr>,  //, FaceFluxExpr>,
+            DiagExpr,
+            UpperExpr,       // upper not yet used
+            LowerExpr,       // lower not yet used
+            FaceFluxExpr,    // face flux not used
+            SourceExpr
+        >
+        (
+            true,       //elems.hasDiag(),
+            false,      //elems.hasUpper(),
+            false,      //elems.hasLower(),
+            false,      //elems.hasFaceFluxCorrection(),
+            dimensions  //elems.dimensions()
+        ),
+        diag_(diag),
+        source_(source)
+    {
+        //Pout<< "fvMatrixInternal :" << nl
+        //    << "    diag_  : " << boost::core::demangle(typeid(diag_).name())
+        //    << nl
+        //    << "    source_: "
+        //    << boost::core::demangle(typeid(source_).name())
+        //    << nl
+        //    << "    *this  : " << boost::core::demangle(typeid(*this).name())
+        //    << endl;
+    }
+
+    // Construct from component expressions
+    fvMatrixInternal
+    (
+        const dimensionSet& dimensions,
+        const DiagExpr&& diag,
+        const SourceExpr&& source
+    )
+    :
+        fvMatrixExpression
+        <
+            fvMatrixInternal<DiagExpr, SourceExpr>,  //, FaceFluxExpr>,
+            DiagExpr,
+            UpperExpr,       // upper not yet used
+            LowerExpr,       // lower not yet used
+            FaceFluxExpr,       // faceflux not yet used
+            SourceExpr
+        >
+        (
+            true,       //elems.hasDiag(),
+            false,      //elems.hasUpper(),
+            false,      //elems.hasLower(),
+            false,      //elems.hasFaceFluxCorrection(),
+            dimensions  //elems.dimensions()
+        ),
+        diag_(diag),
+        source_(source)
+    {}
+
+
+    // lduMatrix part
+
+        auto diag() const
+        {
+            return diag_;
+        }
+
+        auto upper() const
+        {
+            return List<value_type>::null();
+        }
+
+        auto lower() const
+        {
+            return List<value_type>::null();
+        }
+
+
+    // fvMatrix
+
+        auto source() const
+        {
+            return source_;
+        }
+
+        auto internalCoeffs(const label i) const
+        {
+            return List<value_type>::null();
+        }
+
+        auto boundaryCoeffs(const label i) const
+        {
+            return List<value_type>::null();
+        }
+
+        auto faceFluxCorrection() const
+        {
+            return GeometricField<value_type, fvPatchField, volMesh>::null();
+        }
+};
+template<typename E1>
+auto SuSp
+(
+    const dimensionSet& dimensions, // dimensions of expression
+    const E1& expression,
+    const GeometricField<typename E1::value_type, fvPatchField, volMesh>& fld
+)
+{
+    //- Wrap of constant as a list expression
+    typedef Expression::UniformListWrap<scalar> constant;
+    //- Wrap of List as an expression
+    typedef ListConstRefWrap<typename E1::value_type> expr;
+
+    // Wrap zero
+    const constant constantZero(fld.size(), 0.0);
+
+    // Wrap mesh volume
+    const expr V(fld.mesh().V());
+
+    return fvMatrixInternal
+    (
+        dimVol*dimensions*fld.dimensions(),
+        V*max(expression, constantZero),
+        V*min(expression, constantZero)*expr(fld.internalField())
+    );
+}
+template<typename E1>
+auto Sp
+(
+    const dimensionSet& dimensions, // dimensions of expression
+    const E1& expression,
+    const GeometricField<typename E1::value_type, fvPatchField, volMesh>& fld
+)
+{
+    //- Wrap of constant as a list expression
+    typedef Expression::UniformListWrap<scalar> constant;
+    //- Wrap of List as an expression
+    typedef ListConstRefWrap<typename E1::value_type> expr;
+
+    // Wrap zero
+    const constant constantZero(fld.size(), 0.0);
+
+    // Wrap mesh volume
+    const expr V(fld.mesh().V());
+
+    return fvMatrixInternal
+    (
+        dimVol*dimensions*fld.dimensions(),
+        V*expr(fld.internalField()),
+        constantZero    // Still provides an expression for source. TBD.
+    );
+}
+template<typename E1>
+auto Su
+(
+    const dimensionSet& dimensions, // dimensions of expression
+    const E1& expression,
+    const GeometricField<typename E1::value_type, fvPatchField, volMesh>& fld
+)
+{
+    //- Wrap of constant as a list expression
+    typedef Expression::UniformListWrap<scalar> constant;
+    //- Wrap of List as an expression
+    typedef ListConstRefWrap<typename E1::value_type> expr;
+
+    // Wrap zero
+    const constant constantZero(fld.size(), 0.0);
+
+    // Wrap mesh volume
+    const expr V(fld.mesh().V());
+
+    return fvMatrixInternal
+    (
+        dimVol*dimensions*fld.dimensions(),
+        constantZero    // Still provides an expression for diag. TBD.
+        -V*expression
+    );
+}
+
+
+//XXXXXXXX
+// Alternative in-place fvMatrix sources. Can be removed once ones above
+// work.
+template<class Expr>
+void Su
+(
+    fvMatrix<typename Expr::value_type>& m,
+    const Expr& expression
+)
+{
+    //- Wrap of List as an expression
+    typedef ListConstRefWrap<typename Expr::value_type> expr;
+    //- Evaluator of an expression
+    typedef ListRefWrap<typename Expr::value_type> evaluator;
+
+    // Wrap mesh volume
+    const expr V(m.psi().mesh().V());
+
+    // Add expression to source
+    auto& s = m.source();
+    evaluator(s, expr(s) + V*expression);
+}
+template<class Expr>
+void rhs
+(
+    fvMatrix<typename Expr::value_type>& m,
+    const Expr& expression
+)
+{
+    Su(m, expression);
+}
+template<class Expr, class Expr2>
+void Sp
+(
+    fvMatrix<typename Expr::value_type>& m,
+    const Expr2& mult,
+    const Expr& expression
+)
+{
+    //- Wrap of List as an expression
+    typedef ListConstRefWrap<typename Expr::value_type> expr;
+    //- Evaluator of an expression
+    typedef ListRefWrap<typename Expr::value_type> evaluator;
+
+    // Wrap mesh volume
+    const expr V(m.psi().mesh().V());
+
+    // Add expression onto diag
+    auto& d = m.diag();
+    evaluator(d, expr(d) - mult*V*expression);
+}
+template<class Expr, class Expr2>
+void SuSp
+(
+    fvMatrix<typename Expr::value_type>& m,
+    const Expr2& mult,
+    const Expr& expression
+)
+{
+    //- Wrap of constant as a list expression
+    typedef UniformListWrap<scalar> constant;
+    //- Wrap of List as an expression
+    typedef ListConstRefWrap<typename Expr::value_type> expr;
+    //- Evaluator of an expression
+    typedef ListRefWrap<typename Expr::value_type> evaluator;
+
+    const constant constantZero(expression.size(), 0.0);
+
+    // Wrap mesh volume
+    const expr V(m.psi().mesh().V());
+
+    // Linearise expression
+    auto& diag = m.diag();
+    auto& source = m.source();
+
+    // diag() += domain*max(susp.field(), scalar(0));
+    evaluator
+    (
+        diag,
+        expr(diag) - mult*V*max(expression, constantZero)
+    );
+    // source() -= domain*min(susp.field(), scalar(0))*fld.primitiveField();
+    evaluator
+    (
+        source,
+        expr(source)
+      + (
+            mult
+           *V
+           *min(expression, constantZero)
+           *expr(m.psi().internalField())
+        )
+    );
+}
+//XXXXXXXX
+
+
+
+// Some discretisation using expressions
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+template<typename E1, typename E2>
+class interpolate
+:
+    public GeometricFieldExpression
+    <
+        interpolate<E1, E2>,
+
+        //- Type to return for internal field
+        List_interpolate
+        <
+            IndirectConstWrap<typename E1::IntExpr>,
+            IndirectConstWrap<typename E1::IntExpr>,
+            typename E2::IntExpr
+        >,
+
+        //- Type to return for uncoupled patch field
+        typename E1::UncoupledPatchExpr,
+
+        //- Type to return for coupled patch field
+        List_interpolate
+        <
+            IndirectConstWrap<typename E1::IntExpr>,
+            typename E1::CoupledPatchExpr,
+            typename E2::CoupledPatchExpr
+        >,
+
+        //- Type of actual values
+        typename E1::value_type
+    >
+{
+    //- Internal values
+    typename std::conditional
+    <
+        E1::is_leaf,
+        const E1&,
+        const E1
+    >::type cellVals_;
+
+    //- Interpolation weights
+    typename std::conditional
+    <
+        E2::is_leaf,
+        const E2&,
+        const E2
+    >::type faceWeights_;
+
+    //- Addressing container
+    const fvMesh& mesh_;
+
+public:
+
+    static constexpr bool is_leaf = false;
+
+    //- Type to return for internal field
+    typedef List_interpolate
+    <
+        IndirectConstWrap<typename E1::IntExpr>,
+        IndirectConstWrap<typename E1::IntExpr>,
+        typename E2::IntExpr
+    > IntExpr;
+
+    //- Type to return for uncoupled patch field
+    typedef typename E1::UncoupledPatchExpr UncoupledPatchExpr;
+
+    //- Type to return for coupled patch field
+    typedef List_interpolate
+    <
+        IndirectConstWrap<typename E1::IntExpr>,
+        typename E1::CoupledPatchExpr,
+        typename E2::CoupledPatchExpr
+    > CoupledPatchExpr;
+
+    //- Type of actual values
+    typedef typename E1::value_type value_type;
+
+
+    // Construct from components
+    interpolate
+    (
+        const E1& cellVals,
+        const E2& faceWeights,
+        const fvMesh& mesh
+    )
+    :
+        GeometricFieldExpression
+        <
+            interpolate<E1, E2>,
+            IntExpr,
+            UncoupledPatchExpr,
+            CoupledPatchExpr,
+            value_type
+        >
+        (
+            cellVals.dimensions()*faceWeights.dimensions(),
+            faceWeights.oriented()
+        ),
+        cellVals_(cellVals),
+        faceWeights_(faceWeights),
+        mesh_(mesh)
+    {}
+
+    value_type operator[](const label i) const
+    {
+        const auto ownVal = cellVals_[mesh_.owner()[i]];
+        const auto neiVal = cellVals_[mesh_.neighbour()[i]];
+        return faceWeights_[i]*(ownVal-neiVal) + neiVal;
+    }
+
+    label size() const
+    {
+        return faceWeights_.size(); // number of internal faces
+    }
+
+    IntExpr internalField() const
+    {
+        return IntExpr
+        (
+            IndirectConstWrap<typename E1::IntExpr>
+            (
+                cellVals_.internalField(),
+                mesh_.owner()
+            ),
+            IndirectConstWrap<typename E1::IntExpr>
+            (
+                cellVals_.internalField(),
+                mesh_.neighbour()
+            ),
+            faceWeights_.internalField()
+        );
+    }
+
+    UncoupledPatchExpr patchField(const label i) const
+    {
+        return UncoupledPatchExpr(cellVals_.patchField(i));
+    }
+
+    CoupledPatchExpr coupledPatchField(const label i) const
+    {
+        return CoupledPatchExpr
+        (
+            IndirectConstWrap<typename E1::IntExpr>
+            (
+                cellVals_.internalField(),
+                mesh_.boundaryMesh()[i].faceCells()
+            ),
+            cellVals_.access
+            (
+                [&](const auto& fld, const label i)
+                {
+                    return ListConstTmpWrap<Field<value_type>>
+                    (
+                        fld.boundaryField()[i].patchNeighbourField()
+                    );
+                },
+                i
+            ),
+            faceWeights_.coupledPatchField(i)
+        );
+    }
+};
+
+
+template<typename E1>
+class linearInterpolate
+:
+    public interpolate
+    <
+        E1,
+        GeometricFieldConstRefWrap<surfaceScalarField>
+    >
+{
+public:
+
+    // Construct from components
+    linearInterpolate
+    (
+        const E1& cellVals,
+        const fvMesh& mesh
+    )
+    :
+        interpolate
+        <
+            E1,
+            GeometricFieldConstRefWrap<surfaceScalarField>
+        >
+        (
+            cellVals,
+            mesh.surfaceInterpolation::weights().expr(),
+            mesh
+        )
+    {}
+};
+
+
+template<class Type, class E1, class E2>
+void fvmLaplacianUncorrected
+(
+    fvMatrix<Type>& fvm,
+    const E1& gammaMagSf,
+    const E2& deltaCoeffs
+)
+{
+    typedef Expression::ListConstRefWrap<Type> expr;
+    typedef Expression::ListRefWrap<Type> evaluator;
+
+    const auto& vf = fvm.psi();
+
+    // Set dimensions. Or should we just check?
+    const_cast<dimensionSet&>(fvm.dimensions()).reset
+    (
+        deltaCoeffs.dimensions()
+      * gammaMagSf.dimensions()
+      * vf.dimensions()
+    );
+
+    Expression::ListRefWrap<scalar>
+    (
+        fvm.upper(),
+        deltaCoeffs.internalField()*gammaMagSf.internalField()
+    );
+    fvm.negSumDiag();
+
+    forAll(vf.boundaryField(), patchi)
+    {
+        const auto& pvf = vf.boundaryField()[patchi];
+        auto& intCoeffs = fvm.internalCoeffs()[patchi];
+        auto& bouCoeffs = fvm.boundaryCoeffs()[patchi];
+
+        if (pvf.coupled())
+        {
+            // Evaluate deltaCoeffs
+            scalarField pDeltaCoeffs;
+            ListRefWrap<scalar>(pDeltaCoeffs, deltaCoeffs.patchField(patchi));
+
+            //intCoeffs = pGamma*pvf.gradientInternalCoeffs(pDeltaCoeffs);
+            evaluator
+            (
+                intCoeffs,
+                expr(pvf.gradientInternalCoeffs(pDeltaCoeffs)())
+               *gammaMagSf.patchField(patchi)
+            );
+            //bouCoeffs = -pGamma*pvf.gradientBoundaryCoeffs(pDeltaCoeffs);
+            evaluator
+            (
+                bouCoeffs,
+               -expr(pvf.gradientBoundaryCoeffs(pDeltaCoeffs)())
+               *gammaMagSf.patchField(patchi)
+            );
+        }
+        else
+        {
+            //intCoeffs = pGamma*pvf.gradientInternalCoeffs();
+            evaluator
+            (
+                intCoeffs,
+                expr(pvf.gradientInternalCoeffs()())
+               *gammaMagSf.patchField(patchi)
+            );
+            //bouCoeffs = -pGamma*pvf.gradientBoundaryCoeffs();
+            evaluator
+            (
+                bouCoeffs,
+               -expr(pvf.gradientBoundaryCoeffs()())
+               *gammaMagSf.patchField(patchi)
+            );
+        }
+    }
+}
+
+
+} // End namespace Expression
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 5c338c4f27..89088f2191 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -2984,6 +2984,77 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const fvMatrix<Type>& fvm)
 }
 
 
+// * * * * * * * * * * * * * Expression Templates  * * * * * * * * * * * * * //
+
+template<class Type>
+template<typename E>
+Foam::fvMatrix<Type>::fvMatrix
+(
+    const GeometricField<Type, fvPatchField, volMesh>& psi,
+    const Expression::fvMatrixExpression
+    <
+        E,
+        typename E::DiagExpr,
+        typename E::UpperExpr,
+        typename E::LowerExpr,
+        typename E::FaceFluxExpr,
+        typename E::SourceExpr
+   >& expr
+)
+:
+    lduMatrix(psi.mesh()),
+    psi_(psi),
+    useImplicit_(false),
+    lduAssemblyName_(),
+    nMatrix_(0),
+    dimensions_(expr.dimensions()),
+    source_(psi.size(), Zero),
+    internalCoeffs_(psi.mesh().boundary().size()),
+    boundaryCoeffs_(psi.mesh().boundary().size())
+{
+    DebugInFunction
+        << "Constructing fvMatrix<Type> from expression for field "
+        << psi_.name() << endl;
+
+    checkImplicit();
+
+    // Fill in diag,upper,lower etc.
+    expr.evaluate(*this);
+
+    auto& psiRef = this->psi(0);
+    const label currentStatePsi = psiRef.eventNo();
+    psiRef.boundaryFieldRef().updateCoeffs();
+    psiRef.eventNo() = currentStatePsi;
+}
+
+
+template<class Type>
+Foam::Expression::fvMatrixConstRefWrap<Foam::fvMatrix<Type>>
+Foam::fvMatrix<Type>::expr() const
+{
+    return Expression::fvMatrixConstRefWrap<Foam::fvMatrix<Type>>(*this);
+}
+
+
+template<class Type>
+template<typename E>
+void Foam::fvMatrix<Type>::operator=
+(
+    const Expression::fvMatrixExpression
+    <
+        E,
+        typename E::DiagExpr,
+        typename E::UpperExpr,
+        typename E::LowerExpr,
+        typename E::FaceFluxExpr,
+        typename E::SourceExpr
+   >& expr
+)
+{
+    expr.evaluate(*this);
+}
+
+
 // * * * * * * * * * * * * * * * * Solvers * * * * * * * * * * * * * * * * * //
 
 #include "fvMatrixSolve.C"
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
index ff78e157b0..b56b7ebced 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
@@ -53,6 +53,8 @@ SourceFiles
 #include "lduPrimitiveMeshAssembly.H"
 #include "lduMesh.H"
 
+#include "fvMatrixExpression.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -741,6 +743,47 @@ public:
             Ostream&,
             const fvMatrix<Type>&
         );
+
+
+    // Expression templates
+
+        //- Construct given a field to solve for and expressions for
+        //- all the components
+        template<typename E>
+        fvMatrix
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& psi,
+            const Expression::fvMatrixExpression
+            <
+                E,
+                typename E::DiagExpr,
+                typename E::UpperExpr,
+                typename E::LowerExpr,
+                typename E::FaceFluxExpr,
+                typename E::SourceExpr
+           >& expr
+        );
+
+        //- Wrap value as expression
+        typename Expression::fvMatrixConstRefWrap
+        <
+            fvMatrix<Type>
+        > expr() const;
+
+        //- Assign values from expression
+        template<typename E>
+        void operator=
+        (
+            const Expression::fvMatrixExpression
+            <
+                E,
+                typename E::DiagExpr,
+                typename E::UpperExpr,
+                typename E::LowerExpr,
+                typename E::FaceFluxExpr,
+                typename E::SourceExpr
+           >& expr
+       );
 };
 
 
diff --git a/src/fused/finiteVolume/fusedGaussLaplacianScheme.C b/src/fused/finiteVolume/fusedGaussLaplacianScheme.C
index 733828511a..8ee1c4c3b3 100644
--- a/src/fused/finiteVolume/fusedGaussLaplacianScheme.C
+++ b/src/fused/finiteVolume/fusedGaussLaplacianScheme.C
@@ -45,6 +45,87 @@ namespace fv
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+template<class Type, class GType>
+template<class E1, class E2>
+void fusedGaussLaplacianScheme<Type, GType>::fvmCorrection
+(
+    fvMatrix<Type>& fvm,
+    const dimensionSet& gammaDim,
+    const E1& gammaMagSf,
+    const E2& corr
+)
+{
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> surfaceType;
+
+    const auto& vf = fvm.psi();
+    const fvMesh& mesh = vf.mesh();
+    const auto V = mesh.V().expr();
+
+    if (mesh.fluxRequired(vf.name()))
+    {
+        fvm.faceFluxCorrectionPtr() = std::make_unique<surfaceType>
+        (
+            IOobject
+            (
+                "faceFluxCorr",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            gammaDim
+           *mesh.magSf().dimensions()
+           *corr.data().dimensions()
+        );
+        auto& faceFluxCorr = *fvm.faceFluxCorrectionPtr();
+
+        faceFluxCorr = corr*gammaMagSf;
+
+        fvm.source() =
+            fvm.source().expr()
+          - (
+                fvc::div
+                (
+                    faceFluxCorr
+                )().primitiveField().expr()
+               *V
+            );
+    }
+    else
+    {
+        surfaceType faceFluxCorr
+        (
+            IOobject
+            (
+                "faceFluxCorr",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            gammaDim
+           *mesh.magSf().dimensions()
+           *corr.data().dimensions()
+        );
+        faceFluxCorr = corr*gammaMagSf;
+
+        fvm.source() =
+            fvm.source().expr()
+          - (
+                fvc::div
+                (
+                    faceFluxCorr
+                )().primitiveField().expr()
+               *V
+            );
+    }
+}
+
+
 template<class Type, class GType>
 tmp<fvMatrix<Type>>
 fusedGaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
diff --git a/src/fused/finiteVolume/fusedGaussLaplacianScheme.H b/src/fused/finiteVolume/fusedGaussLaplacianScheme.H
index 85753bff01..4e636b4470 100644
--- a/src/fused/finiteVolume/fusedGaussLaplacianScheme.H
+++ b/src/fused/finiteVolume/fusedGaussLaplacianScheme.H
@@ -65,6 +65,17 @@ class fusedGaussLaplacianScheme
 {
     // Private Member Functions
 
+        //- Helper for calculate source correction. Partially expression
+        //  templated. Move to fvMatrixExpression.H?
+        template<class E1, class E2>
+        static void fvmCorrection
+        (
+            fvMatrix<Type>& fvm,
+            const dimensionSet& gammaDim,
+            const E1& gammaMagSf,
+            const E2& corr
+        );
+
         //- See gaussLaplacianScheme
         tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> gammaSnGradCorr
         (
@@ -198,6 +209,8 @@ defineFvmLaplacianScalarGamma(sphericalTensor);
 defineFvmLaplacianScalarGamma(symmTensor);
 defineFvmLaplacianScalarGamma(tensor);
 
+#undef defineFvmLaplacianScalarGamma
+
 // Unweighted laplacian
 template<>
 tmp<GeometricField<scalar, fvPatchField, volMesh>>
@@ -227,20 +240,6 @@ fusedGaussLaplacianScheme<vector, scalar>::fvcLaplacian
     const GeometricField<scalar, fvPatchField, volMesh>&,
     const GeometricField<vector, fvPatchField, volMesh>&
 );
-template<>
-tmp<fvMatrix<scalar>>
-fusedGaussLaplacianScheme<scalar, scalar>::fvmLaplacian
-(
-    const GeometricField<scalar, fvPatchField, volMesh>&,
-    const GeometricField<scalar, fvPatchField, volMesh>&
-);
-template<>
-tmp<fvMatrix<vector>>
-fusedGaussLaplacianScheme<vector, scalar>::fvmLaplacian
-(
-    const GeometricField<scalar, fvPatchField, volMesh>&,
-    const GeometricField<vector, fvPatchField, volMesh>&
-);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/fused/finiteVolume/fusedGaussLaplacianSchemes.C b/src/fused/finiteVolume/fusedGaussLaplacianSchemes.C
index 0433c423d3..8b96b25919 100644
--- a/src/fused/finiteVolume/fusedGaussLaplacianSchemes.C
+++ b/src/fused/finiteVolume/fusedGaussLaplacianSchemes.C
@@ -28,6 +28,7 @@ License
 
 #include "fusedGaussLaplacianScheme.H"
 #include "fvMesh.H"
+#include "fvMatrixExpression.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -40,61 +41,81 @@ Foam::tmp<Foam::fvMatrix<Foam::Type>>                                          \
 Foam::fv::fusedGaussLaplacianScheme<Foam::Type, Foam::scalar>::                \
 fvmLaplacian                                                                   \
 (                                                                              \
-    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,           \
+    const GeometricField<scalar, fvPatchField, volMesh>& gamma,                \
     const GeometricField<Type, fvPatchField, volMesh>& vf                      \
 )                                                                              \
 {                                                                              \
     DebugPout<< "fusedGaussLaplacianScheme::fvmLaplacian on " << vf.name()     \
-        << " with scalar gamma " << gamma.name() << endl;                      \
+        << " with scalar gamma " << gamma.name() << " and ET" << endl;         \
                                                                                \
     const fvMesh& mesh = this->mesh();                                         \
                                                                                \
-    GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf              \
-    (                                                                          \
-        gamma*mesh.magSf()                                                     \
-    );                                                                         \
+    const auto weights = this->tinterpGammaScheme_().weights(gamma).expr();    \
+    const auto gammaMagSf =                                                    \
+        Expression::interpolate(gamma.expr(), weights, mesh)                   \
+      * mesh.magSf().expr();                                                   \
+    /* For compatibility with linearInterpolate : avoid orientation. TBD. */   \
+    const_cast<orientedType&>(gammaMagSf.oriented()).setOriented(false);       \
+    const auto deltaCoeffs = this->tsnGradScheme_().deltaCoeffs(vf).expr();    \
                                                                                \
-    tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected                         \
+    tmp<fvMatrix<Type>> tfvm                                                   \
     (                                                                          \
-        gammaMagSf,                                                            \
-        this->tsnGradScheme_().deltaCoeffs(vf),                                \
-        vf                                                                     \
+        new fvMatrix<Type>                                                     \
+        (                                                                      \
+            vf,                                                                \
+            gamma.dimensions()*mesh.magSf().dimensions()*vf.dimensions()       \
+        )                                                                      \
     );                                                                         \
     fvMatrix<Type>& fvm = tfvm.ref();                                          \
                                                                                \
+    Expression::fvmLaplacianUncorrected(fvm, gammaMagSf, deltaCoeffs);         \
+                                                                               \
     if (this->tsnGradScheme_().corrected())                                    \
     {                                                                          \
-        if (mesh.fluxRequired(vf.name()))                                      \
-        {                                                                      \
-            fvm.faceFluxCorrectionPtr() = std::make_unique                     \
-            <                                                                  \
-                GeometricField<Type, fvsPatchField, surfaceMesh>               \
-            >                                                                  \
-            (                                                                  \
-                gammaMagSf*this->tsnGradScheme_().correction(vf)               \
-            );                                                                 \
-                                                                               \
-            fvm.source() -=                                                    \
-                mesh.V()*                                                      \
-                fvc::div                                                       \
-                (                                                              \
-                    *fvm.faceFluxCorrectionPtr()                               \
-                )().primitiveField();                                          \
-        }                                                                      \
-        else                                                                   \
-        {                                                                      \
-            fvm.source() -=                                                    \
-                mesh.V()*                                                      \
-                fvc::div                                                       \
-                (                                                              \
-                    gammaMagSf*this->tsnGradScheme_().correction(vf)           \
-                )().primitiveField();                                          \
-        }                                                                      \
+        const auto corr(this->tsnGradScheme_().correction(vf).expr());         \
+        fvmCorrection(fvm, gamma.dimensions(), gammaMagSf, corr);              \
     }                                                                          \
                                                                                \
     return tfvm;                                                               \
 }                                                                              \
                                                                                \
+template<>                                                                     \
+Foam::tmp<Foam::fvMatrix<Foam::Type>>                                          \
+Foam::fv::fusedGaussLaplacianScheme<Foam::Type, Foam::scalar>::                \
+fvmLaplacian                                                                   \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,           \
+    const GeometricField<Type, fvPatchField, volMesh>& vf                      \
+)                                                                              \
+{                                                                              \
+    DebugPout<< "fusedGaussLaplacianScheme::fvmLaplacian on " << vf.name()     \
+        << " with interpolated gamma " << gamma.name() << " and ET" << endl;   \
+                                                                               \
+    const fvMesh& mesh = this->mesh();                                         \
+                                                                               \
+    const auto gammaMagSf = gamma.expr()* mesh.magSf().expr();                 \
+    const auto deltaCoeffs = this->tsnGradScheme_().deltaCoeffs(vf).expr();    \
+                                                                               \
+    tmp<fvMatrix<Type>> tfvm                                                   \
+    (                                                                          \
+        new fvMatrix<Type>                                                     \
+        (                                                                      \
+            vf,                                                                \
+            gamma.dimensions()*mesh.magSf().dimensions()*vf.dimensions()       \
+        )                                                                      \
+    );                                                                         \
+    fvMatrix<Type>& fvm = tfvm.ref();                                          \
+                                                                               \
+    Expression::fvmLaplacianUncorrected(fvm, gammaMagSf, deltaCoeffs);         \
+                                                                               \
+    if (this->tsnGradScheme_().corrected())                                    \
+    {                                                                          \
+        const auto corr(this->tsnGradScheme_().correction(vf).expr());         \
+        fvmCorrection(fvm, gamma.dimensions(), gammaMagSf, corr);              \
+    }                                                                          \
+                                                                               \
+    return tfvm;                                                               \
+}                                                                              \
                                                                                \
 template<>                                                                     \
 Foam::tmp<Foam::GeometricField<Foam::Type, Foam::fvPatchField, Foam::volMesh>> \
@@ -500,38 +521,6 @@ Foam::fv::fusedGaussLaplacianScheme<Foam::vector, Foam::scalar>::fvcLaplacian
 }
 
 
-template<>
-Foam::tmp<Foam::fvMatrix<Foam::scalar>>
-Foam::fv::fusedGaussLaplacianScheme<Foam::scalar, Foam::scalar>::fvmLaplacian
-(
-    const GeometricField<scalar, fvPatchField, volMesh>& gamma,
-    const GeometricField<scalar, fvPatchField, volMesh>& vf
-)
-{
-    // TBD
-    DebugPout
-        << "fusedGaussLaplacianScheme<scalar, scalar>::fvmLaplacian"
-        << " on " << vf.name() << " with gamma " << gamma.name() << endl;
-    return fvmLaplacian(this->tinterpGammaScheme_().interpolate(gamma)(), vf);
-}
-
-
-template<>
-Foam::tmp<Foam::fvMatrix<Foam::vector>>
-Foam::fv::fusedGaussLaplacianScheme<Foam::vector, Foam::scalar>::fvmLaplacian
-(
-    const GeometricField<scalar, fvPatchField, volMesh>& gamma,
-    const GeometricField<vector, fvPatchField, volMesh>& vf
-)
-{
-    // TBD
-    DebugPout
-        << "fusedGaussLaplacianScheme<vector, scalar>::fvmLaplacian"
-        << " on " << vf.name() << " with gamma " << gamma.name() << endl;
-    return fvmLaplacian(this->tinterpGammaScheme_().interpolate(gamma)(), vf);
-}
-
-
 template<>
 Foam::tmp
 <