STYLE: ListExpression: copyright, indenting

Feature grad update

See merge request Development/openfoam!764

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

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

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;
+}
+
+
+// ************************************************************************* //

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;
+}
+
+
+// ************************************************************************* //

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

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);
+        }
 };
 
 

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];
+            }
+        }
 };
 
 

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
+            );
+        }
 };
 
 

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
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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);
+        }
 };
 
 

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

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
+        );
 };
 
 

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;
+            //}
+        }
 };
 
 

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

src/Pstream/mpi/UPstreamRequest.C

@@ -855,7 +855,7 @@ bool Foam::UPstream::finishedRequest(const label i)
     // This allows MPI to progress behind the scenes if it wishes.
 
     int flag = 0;
-    if (i < 0 || i >= PstreamGlobals::outstandingRequests_.size())
+    if (i >= 0 && i < PstreamGlobals::outstandingRequests_.size())
     {
         auto& request = PstreamGlobals::outstandingRequests_[i];
 

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020,2025 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -189,23 +189,6 @@ public:
                 template<class TrackingData>
                 inline bool equal(const deltaData&, TrackingData& td) const;
 
-                //- Interpolate between two values (lerp). Returns true if
-                //- causes changes. Not sure if needs to be specialised between
-                //- face and cell and what index is needed...
-                template<class TrackingData>
-                inline bool interpolate
-                (
-                    const polyMesh&,
-                    const point& pt,
-                    const label i0,
-                    const deltaData& f0,
-                    const label i1,
-                    const deltaData& f1,
-                    const scalar weight,
-                    const scalar tol,
-                    TrackingData& td
-                );
-
 
         // Member Operators
 
@@ -313,17 +296,6 @@ public:
 template<>
 struct is_contiguous<LESModels::smoothDelta::deltaData> : std::true_type {};
 
-//- Interpolation - used in e.g. cyclicAMI interpolation
-inline LESModels::smoothDelta::deltaData lerp
-(
-    const LESModels::smoothDelta::deltaData& a,
-    const LESModels::smoothDelta::deltaData& b,
-    const scalar t
-)
-{
-    return LESModels::smoothDelta::deltaData(lerp(a.delta(), b.delta(), t));
-}
-
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDeltaDeltaDataI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019,2025 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -191,40 +191,6 @@ inline bool Foam::LESModels::smoothDelta::deltaData::equal
 }
 
 
-template<class TrackingData>
-inline bool Foam::LESModels::smoothDelta::deltaData::interpolate
-(
-    const polyMesh&,
-    const point& pt,
-    const label i0,
-    const deltaData& f0,
-    const label i1,
-    const deltaData& f1,
-    const scalar weight,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    if (f0.valid(td) && f1.valid(td))
-    {
-        const deltaData w2(lerp(f0.delta(), f1.delta(), weight));
-        return update(w2, 1.0, tol, td);
-    }
-    else if (f0.valid(td))
-    {
-        return update(f0, 1.0, tol, td);
-    }
-    else if (f1.valid(td))
-    {
-        return update(f1, 1.0, tol, td);
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::LESModels::smoothDelta::deltaData::operator==

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C

@@ -218,9 +218,11 @@ tmp<vectorField> atmBoundaryLayer::U(const vectorField& pCf) const
     const scalar groundMin = zDir() & ppMin_;
 
     // (YGCJ:Table 1, RH:Eq. 6, HW:Eq. 5)
+    scalarField zEff(max((zDir() & pCf) - groundMin - d + z0, z0));
+
     scalarField Un
     (
-        (Ustar(z0)/kappa_)*log(((zDir() & pCf) - groundMin - d + z0)/z0)
+        (Ustar(z0)/kappa_)*log(zEff/z0)
     );
 
     return flowDir()*Un;
@@ -235,9 +237,9 @@ tmp<scalarField> atmBoundaryLayer::k(const vectorField& pCf) const
     const scalar groundMin = zDir() & ppMin_;
 
     // (YGCJ:Eq. 21; RH:Eq. 7, HW:Eq. 6 when C1=0 and C2=1)
-    return
-        sqr(Ustar(z0))/sqrt(Cmu_)
-       *sqrt(C1_*log(((zDir() & pCf) - groundMin - d + z0)/z0) + C2_);
+    scalarField zEff(max((zDir() & pCf) - groundMin - d + z0, z0));
+
+    return sqr(Ustar(z0))/sqrt(Cmu_)*sqrt(C1_*log(zEff/z0) + C2_);
 }
 
 
@@ -249,9 +251,9 @@ tmp<scalarField> atmBoundaryLayer::epsilon(const vectorField& pCf) const
     const scalar groundMin = zDir() & ppMin_;
 
     // (YGCJ:Eq. 22; RH:Eq. 8, HW:Eq. 7 when C1=0 and C2=1)
-    return
-        pow3(Ustar(z0))/(kappa_*((zDir() & pCf) - groundMin - d + z0))
-       *sqrt(C1_*log(((zDir() & pCf) - groundMin - d + z0)/z0) + C2_);
+    scalarField zEff(max((zDir() & pCf) - groundMin - d + z0, z0));
+
+    return pow3(Ustar(z0))/(kappa_*zEff)*sqrt(C1_*log(zEff/z0) + C2_);
 }
 
 
@@ -263,7 +265,9 @@ tmp<scalarField> atmBoundaryLayer::omega(const vectorField& pCf) const
     const scalar groundMin = zDir() & ppMin_;
 
     // (YGJ:Eq. 13)
-    return Ustar(z0)/(kappa_*sqrt(Cmu_)*((zDir() & pCf) - groundMin - d + z0));
+    scalarField zEff(max((zDir() & pCf) - groundMin - d + z0, z0));
+
+    return Ustar(z0)/(kappa_*sqrt(Cmu_)*zEff);
 }
 
 

src/dynamicMesh/meshCut/directions/directionInfo/directionInfo.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020,2025 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -242,23 +242,6 @@ public:
             template<class TrackingData>
             inline bool equal(const directionInfo&, TrackingData& td) const;
 
-            //- Interpolate between two values (lerp). Returns true if
-            //- causes changes. Not sure if needs to be specialised between
-            //- face and cell and what index is needed...
-            template<class TrackingData>
-            inline bool interpolate
-            (
-                const polyMesh&,
-                const point& pt,
-                const label i0,
-                const directionInfo& f0,
-                const label i1,
-                const directionInfo& f1,
-                const scalar weight,
-                const scalar tol,
-                TrackingData& td
-            );
-
 
     // Member Operators
 

src/dynamicMesh/meshCut/directions/directionInfo/directionInfoI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020,2025 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -287,35 +287,6 @@ inline bool Foam::directionInfo::equal
 }
 
 
-template<class TrackingData>
-inline bool Foam::directionInfo::interpolate
-(
-    const polyMesh& mesh,
-    const point& pt,
-    const label i0,
-    const directionInfo& f0,
-    const label i1,
-    const directionInfo& f1,
-    const scalar weight,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    if (f0.valid(td))
-    {
-        return updateFace(mesh, -1, f0, tol, td);
-    }
-    if (f1.valid(td))
-    {
-        return updateFace(mesh, -1, f1, tol, td);
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::directionInfo::operator==

src/dynamicMesh/meshCut/wallLayerCells/wallNormalInfo/wallNormalInfo.H

@@ -186,23 +186,6 @@ public:
                 TrackingData& td
             );
 
-            //- Interpolate between two values (lerp). Returns true if
-            //- causes changes. Not sure if needs to be specialised between
-            //- face and cell and what index is needed...
-            template<class TrackingData>
-            inline bool interpolate
-            (
-                const polyMesh&,
-                const point& pt,
-                const label i0,
-                const wallNormalInfo& f0,
-                const label i1,
-                const wallNormalInfo& f1,
-                const scalar weight,
-                const scalar tol,
-                TrackingData& td
-            );
-
             //- Test for equality, with TrackingData
             template<class TrackingData>
             inline bool equal(const wallNormalInfo&, TrackingData& td) const;

src/dynamicMesh/meshCut/wallLayerCells/wallNormalInfo/wallNormalInfoI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020,2025 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -191,35 +191,6 @@ inline bool Foam::wallNormalInfo::equal
 }
 
 
-template<class TrackingData>
-inline bool Foam::wallNormalInfo::interpolate
-(
-    const polyMesh&,
-    const point& pt,
-    const label i0,
-    const wallNormalInfo& f0,
-    const label i1,
-    const wallNormalInfo& f1,
-    const scalar weight,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    if (f0.valid(td))
-    {
-        return update(f0, td);
-    }
-    if (f1.valid(td))
-    {
-        return update(f1, td);
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::wallNormalInfo::operator==

src/dynamicMesh/polyTopoChange/polyTopoChange/refinementData.H

@@ -197,23 +197,6 @@ public:
             template<class TrackingData>
             inline bool equal(const refinementData&, TrackingData& td) const;
 
-            //- Interpolate between two values (lerp). Returns true if
-            //- causes changes. Not sure if needs to be specialised between
-            //- face and cell and what index is needed...
-            template<class TrackingData>
-            inline bool interpolate
-            (
-                const polyMesh&,
-                const point& pt,
-                const label i0,
-                const refinementData& f0,
-                const label i1,
-                const refinementData& f1,
-                const scalar weight,
-                const scalar tol,
-                TrackingData& td
-            );
-
 
     // Member Operators
 

src/dynamicMesh/polyTopoChange/polyTopoChange/refinementDataI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020,2025 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -248,35 +248,6 @@ inline bool Foam::refinementData::equal
 }
 
 
-template<class TrackingData>
-inline bool Foam::refinementData::interpolate
-(
-    const polyMesh& mesh,
-    const point& pt,
-    const label i0,
-    const refinementData& f0,
-    const label i1,
-    const refinementData& f1,
-    const scalar weight,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    if (f0.valid(td))
-    {
-        return updateFace(mesh, -1, f0, tol, td);
-    }
-    if (f1.valid(td))
-    {
-        return updateFace(mesh, -1, f1, tol, td);
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::refinementData::operator==

src/dynamicMesh/polyTopoChange/polyTopoChange/refinementDistanceData.H

@@ -229,23 +229,6 @@ public:
             inline bool equal(const refinementDistanceData&, TrackingData&)
             const;
 
-            //- Interpolate between two values (lerp). Returns true if
-            //- causes changes. Not sure if needs to be specialised between
-            //- face and cell and what index is needed...
-            template<class TrackingData>
-            inline bool interpolate
-            (
-                const polyMesh&,
-                const point& pt,
-                const label i0,
-                const refinementDistanceData& f0,
-                const label i1,
-                const refinementDistanceData& f1,
-                const scalar weight,
-                const scalar tol,
-                TrackingData& td
-            );
-
 
     // Member Operators
 

src/dynamicMesh/polyTopoChange/polyTopoChange/refinementDistanceDataI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020,2025 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -278,35 +278,6 @@ inline bool Foam::refinementDistanceData::equal
 }
 
 
-template<class TrackingData>
-inline bool Foam::refinementDistanceData::interpolate
-(
-    const polyMesh&,
-    const point& pt,
-    const label i0,
-    const refinementDistanceData& f0,
-    const label i1,
-    const refinementDistanceData& f1,
-    const scalar weight,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    if (f0.valid(td))
-    {
-        return update(pt, f0, tol, td);
-    }
-    if (f1.valid(td))
-    {
-        return update(pt, f1, tol, td);
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::refinementDistanceData::operator==

src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C

@@ -227,15 +227,9 @@ bool Foam::cyclicACMIFvPatchField<Type>::all_ready() const
             recvRequests_.start(),
             recvRequests_.size()
         )
-     && UPstream::finishedRequests
-        (
-            recvRequests1_.start(),
-            recvRequests1_.size()
-        )
     )
     {
         recvRequests_.clear();
-        recvRequests1_.clear();
         ++done;
     }
 
@@ -246,15 +240,9 @@ bool Foam::cyclicACMIFvPatchField<Type>::all_ready() const
             sendRequests_.start(),
             sendRequests_.size()
         )
-     && UPstream::finishedRequests
-        (
-            sendRequests1_.start(),
-            sendRequests1_.size()
-        )
     )
     {
         sendRequests_.clear();
-        sendRequests1_.clear();
         ++done;
     }
 
@@ -272,15 +260,9 @@ bool Foam::cyclicACMIFvPatchField<Type>::ready() const
             recvRequests_.start(),
             recvRequests_.size()
         )
-     && UPstream::finishedRequests
-        (
-            recvRequests1_.start(),
-            recvRequests1_.size()
-        )
     )
     {
         recvRequests_.clear();
-        recvRequests1_.clear();
 
         if
         (
@@ -289,15 +271,9 @@ bool Foam::cyclicACMIFvPatchField<Type>::ready() const
                 sendRequests_.start(),
                 sendRequests_.size()
             )
-         && UPstream::finishedRequests
-            (
-                sendRequests1_.start(),
-                sendRequests1_.size()
-            )
         )
         {
             sendRequests_.clear();
-            sendRequests1_.clear();
         }
 
         return true;
@@ -507,7 +483,7 @@ void Foam::cyclicACMIFvPatchField<Type>::initEvaluate
         const Field<Type> pnf(this->primitiveField(), nbrFaceCells);
 
         // Assert that all receives are known to have finished
-        if (!recvRequests_.empty() || !recvRequests1_.empty())
+        if (!recvRequests_.empty())
         {
             FatalErrorInFunction
                 << "Outstanding recv request(s) on patch "
@@ -518,20 +494,14 @@ void Foam::cyclicACMIFvPatchField<Type>::initEvaluate
 
         // Assume that sends are also OK
         sendRequests_.clear();
-        sendRequests1_.clear();
 
         cyclicACMIPatch_.initInterpolate
         (
             pnf,
             sendRequests_,
-            recvRequests_,
             sendBufs_,
-            recvBufs_,
-
-            sendRequests1_,
-            recvRequests1_,
-            sendBufs1_,
-            recvBufs1_
+            recvRequests_,
+            recvBufs_
         );
     }
 }
@@ -577,15 +547,12 @@ void Foam::cyclicACMIFvPatchField<Type>::evaluate
             (
                 Field<Type>::null(),    // Not used for distributed
                 recvRequests_,
-                recvBufs_,
-                recvRequests1_,
-                recvBufs1_
+                recvBufs_
             ).ptr()
         );
 
         // Receive requests all handled by last function call
         recvRequests_.clear();
-        recvRequests1_.clear();
 
         if (doTransform())
         {
@@ -641,7 +608,7 @@ void Foam::cyclicACMIFvPatchField<Type>::initInterfaceMatrixUpdate
         transformCoupleField(pnf, cmpt);
 
         // Assert that all receives are known to have finished
-        if (!recvRequests_.empty() || !recvRequests1_.empty())
+        if (!recvRequests_.empty())
         {
             FatalErrorInFunction
                 << "Outstanding recv request(s) on patch "
@@ -652,20 +619,14 @@ void Foam::cyclicACMIFvPatchField<Type>::initInterfaceMatrixUpdate
 
         // Assume that sends are also OK
         sendRequests_.clear();
-        sendRequests1_.clear();
 
         cyclicACMIPatch_.initInterpolate
         (
             pnf,
             sendRequests_,
-            recvRequests_,
             scalarSendBufs_,
-            scalarRecvBufs_,
-
-            sendRequests1_,
-            recvRequests1_,
-            scalarSendBufs1_,
-            scalarRecvBufs1_
+            recvRequests_,
+            scalarRecvBufs_
         );
     }
 
@@ -720,14 +681,11 @@ void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
             (
                 solveScalarField::null(),    // Not used for distributed
                 recvRequests_,
-                scalarRecvBufs_,
-                recvRequests1_,
-                scalarRecvBufs1_
+                scalarRecvBufs_
             );
 
         // Receive requests all handled by last function call
         recvRequests_.clear();
-        recvRequests1_.clear();
     }
     else
     {
@@ -780,7 +738,7 @@ void Foam::cyclicACMIFvPatchField<Type>::initInterfaceMatrixUpdate
         transformCoupleField(pnf);
 
         // Assert that all receives are known to have finished
-        if (!recvRequests_.empty() || !recvRequests1_.empty())
+        if (!recvRequests_.empty())
         {
             FatalErrorInFunction
                 << "Outstanding recv request(s) on patch "
@@ -791,20 +749,14 @@ void Foam::cyclicACMIFvPatchField<Type>::initInterfaceMatrixUpdate
 
         // Assume that sends are also OK
         sendRequests_.clear();
-        sendRequests1_.clear();
 
         cyclicACMIPatch_.initInterpolate
         (
             pnf,
             sendRequests_,
-            recvRequests_,
             sendBufs_,
-            recvBufs_,
-
-            sendRequests1_,
-            recvRequests1_,
-            sendBufs1_,
-            recvBufs1_
+            recvRequests_,
+            recvBufs_
         );
     }
 
@@ -846,14 +798,11 @@ void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
             (
                 Field<Type>::null(),    // Not used for distributed
                 recvRequests_,
-                recvBufs_,
-                recvRequests1_,
-                recvBufs1_
+                recvBufs_
             );
 
         // Receive requests all handled by last function call
         recvRequests_.clear();
-        recvRequests1_.clear();
     }
     else
     {

src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H

@@ -102,28 +102,6 @@ class cyclicACMIFvPatchField
             //- Scalar receive buffers
             mutable PtrList<List<solveScalar>> scalarRecvBufs_;
 
-
-            // Only used for AMI caching
-
-                //- Current range of send requests (non-blocking)
-                mutable labelRange sendRequests1_;
-
-                //- Current range of recv requests (non-blocking)
-                mutable labelRange recvRequests1_;
-
-                //- Send buffers
-                mutable PtrList<List<Type>> sendBufs1_;
-
-                //- Receive buffers_
-                mutable PtrList<List<Type>> recvBufs1_;
-
-                //- Scalar send buffers
-                mutable PtrList<List<solveScalar>> scalarSendBufs1_;
-
-                //- Scalar receive buffers
-                mutable PtrList<List<solveScalar>> scalarRecvBufs1_;
-
-
             //- Neighbour coupled internal cell data
             mutable autoPtr<Field<Type>> patchNeighbourFieldPtr_;
 

src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.C

@@ -207,15 +207,9 @@ bool Foam::cyclicAMIFvPatchField<Type>::all_ready() const
             recvRequests_.start(),
             recvRequests_.size()
         )
-    &&  UPstream::finishedRequests
-        (
-            recvRequests1_.start(),
-            recvRequests1_.size()
-        )
     )
     {
         recvRequests_.clear();
-        recvRequests1_.clear();
         ++done;
     }
 
@@ -226,15 +220,9 @@ bool Foam::cyclicAMIFvPatchField<Type>::all_ready() const
             sendRequests_.start(),
             sendRequests_.size()
         )
-    &&  UPstream::finishedRequests
-        (
-            sendRequests1_.start(),
-            sendRequests1_.size()
-        )
     )
     {
         sendRequests_.clear();
-        sendRequests1_.clear();
         ++done;
     }
 
@@ -252,15 +240,9 @@ bool Foam::cyclicAMIFvPatchField<Type>::ready() const
             recvRequests_.start(),
             recvRequests_.size()
         )
-     && UPstream::finishedRequests
-        (
-            recvRequests1_.start(),
-            recvRequests1_.size()
-        )
     )
     {
         recvRequests_.clear();
-        recvRequests1_.clear();
 
         if
         (
@@ -269,15 +251,9 @@ bool Foam::cyclicAMIFvPatchField<Type>::ready() const
                 sendRequests_.start(),
                 sendRequests_.size()
             )
-         && UPstream::finishedRequests
-            (
-                sendRequests1_.start(),
-                sendRequests1_.size()
-            )
         )
         {
             sendRequests_.clear();
-            sendRequests1_.clear();
         }
 
         return true;
@@ -343,18 +319,9 @@ Foam::cyclicAMIFvPatchField<Type>::getNeighbourField
 
 
 template<class Type>
-bool Foam::cyclicAMIFvPatchField<Type>::cacheNeighbourField() const
+bool Foam::cyclicAMIFvPatchField<Type>::cacheNeighbourField()
 {
-    // const auto& AMI = this->ownerAMI();
-
-    // if (AMI.cacheActive())
-    // {
-    //     return false;
-    // }
-    // else
-    {
-        return (FieldBase::localBoundaryConsistency() != 0);
-    }
+    return (FieldBase::localBoundaryConsistency() != 0);
 }
 
 
@@ -383,12 +350,11 @@ Foam::cyclicAMIFvPatchField<Type>::getPatchNeighbourField
         }
 
         const auto& fvp = this->patch();
-        const auto& mesh = fvp.boundaryMesh().mesh();
 
         if
         (
             patchNeighbourFieldPtr_
-        && !mesh.upToDatePoints(this->internalField())
+        && !fvp.boundaryMesh().mesh().upToDatePoints(this->internalField())
         )
         {
             //DebugPout
@@ -452,8 +418,7 @@ template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::cyclicAMIFvPatchField<Type>::patchNeighbourField() const
 {
-    // checkCommunicator = true
-    return this->getPatchNeighbourField(true);
+    return this->getPatchNeighbourField(true);  // checkCommunicator = true
 }
 
 
@@ -526,7 +491,7 @@ void Foam::cyclicAMIFvPatchField<Type>::initEvaluate
         const cyclicAMIPolyPatch& cpp = cyclicAMIPatch_.cyclicAMIPatch();
 
         // Assert that all receives are known to have finished
-        if (!recvRequests_.empty() || !recvRequests1_.empty())
+        if (!recvRequests_.empty())
         {
             FatalErrorInFunction
                 << "Outstanding recv request(s) on patch "
@@ -537,20 +502,14 @@ void Foam::cyclicAMIFvPatchField<Type>::initEvaluate
 
         // Assume that sends are also OK
         sendRequests_.clear();
-        sendRequests1_.clear();
 
         cpp.initInterpolate
         (
             pnf,
             sendRequests_,
-            recvRequests_,
             sendBufs_,
-            recvBufs_,
-
-            sendRequests1_,
-            recvRequests1_,
-            sendBufs1_,
-            recvBufs1_
+            recvRequests_,
+            recvBufs_
         );
     }
 }
@@ -603,15 +562,12 @@ void Foam::cyclicAMIFvPatchField<Type>::evaluate
                 Field<Type>::null(),    // Not used for distributed
                 recvRequests_,
                 recvBufs_,
-                recvRequests1_,
-                recvBufs1_,
                 defaultValues
             ).ptr()
         );
 
         // Receive requests all handled by last function call
         recvRequests_.clear();
-        recvRequests1_.clear();
 
         if (doTransform())
         {
@@ -662,7 +618,7 @@ void Foam::cyclicAMIFvPatchField<Type>::initInterfaceMatrixUpdate
         const cyclicAMIPolyPatch& cpp = cyclicAMIPatch_.cyclicAMIPatch();
 
         // Assert that all receives are known to have finished
-        if (!recvRequests_.empty() || !recvRequests1_.empty())
+        if (!recvRequests_.empty())
         {
             FatalErrorInFunction
                 << "Outstanding recv request(s) on patch "
@@ -673,20 +629,14 @@ void Foam::cyclicAMIFvPatchField<Type>::initInterfaceMatrixUpdate
 
         // Assume that sends are also OK
         sendRequests_.clear();
-        sendRequests1_.clear();
 
         cpp.initInterpolate
         (
             pnf,
             sendRequests_,
-            recvRequests_,
             scalarSendBufs_,
-            scalarRecvBufs_,
-
-            sendRequests1_,
-            recvRequests1_,
-            scalarSendBufs1_,
-            scalarRecvBufs1_
+            recvRequests_,
+            scalarRecvBufs_
         );
     }
 
@@ -741,14 +691,11 @@ void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
                 solveScalarField::null(),   // Not used for distributed
                 recvRequests_,
                 scalarRecvBufs_,
-                recvRequests1_,
-                scalarRecvBufs1_,
                 defaultValues
             );
 
         // Receive requests all handled by last function call
         recvRequests_.clear();
-        recvRequests1_.clear();
     }
     else
     {
@@ -810,7 +757,7 @@ void Foam::cyclicAMIFvPatchField<Type>::initInterfaceMatrixUpdate
         const cyclicAMIPolyPatch& cpp = cyclicAMIPatch_.cyclicAMIPatch();
 
         // Assert that all receives are known to have finished
-        if (!recvRequests_.empty() || !recvRequests1_.empty())
+        if (!recvRequests_.empty())
         {
             FatalErrorInFunction
                 << "Outstanding recv request(s) on patch "
@@ -821,20 +768,14 @@ void Foam::cyclicAMIFvPatchField<Type>::initInterfaceMatrixUpdate
 
         // Assume that sends are also OK
         sendRequests_.clear();
-        sendRequests1_.clear();
 
         cpp.initInterpolate
         (
             pnf,
             sendRequests_,
-            recvRequests_,
             sendBufs_,
-            recvBufs_,
-
-            sendRequests1_,
-            recvRequests1_,
-            sendBufs1_,
-            recvBufs1_
+            recvRequests_,
+            recvBufs_
         );
     }
 
@@ -888,14 +829,11 @@ void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
                 Field<Type>::null(),  // Not used for distributed
                 recvRequests_,
                 recvBufs_,
-                recvRequests1_,
-                recvBufs1_,
                 defaultValues
             );
 
         // Receive requests all handled by last function call
         recvRequests_.clear();
-        recvRequests1_.clear();
     }
     else
     {
@@ -980,7 +918,7 @@ void Foam::cyclicAMIFvPatchField<Type>::manipulateMatrix
         }
 
         // Set internalCoeffs and boundaryCoeffs in the assembly matrix
-        // on cyclicAMI patches to be used in the individual matrix by
+        // on clyclicAMI patches to be used in the individual matrix by
         // matrix.flux()
         if (matrix.psi(mat).mesh().fluxRequired(this->internalField().name()))
         {

src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.H

@@ -113,28 +113,6 @@ class cyclicAMIFvPatchField
             //- Scalar receive buffers
             mutable PtrList<List<solveScalar>> scalarRecvBufs_;
 
-
-            // Only used for AMI caching
-
-                //- Current range of send requests (non-blocking)
-                mutable labelRange sendRequests1_;
-
-                //- Current range of recv requests (non-blocking)
-                mutable labelRange recvRequests1_;
-
-                //- Send buffers
-                mutable PtrList<List<Type>> sendBufs1_;
-
-                //- Receive buffers_
-                mutable PtrList<List<Type>> recvBufs1_;
-
-                //- Scalar send buffers
-                mutable PtrList<List<solveScalar>> scalarSendBufs1_;
-
-                //- Scalar receive buffers
-                mutable PtrList<List<solveScalar>> scalarRecvBufs1_;
-
-
             //- Neighbour coupled internal cell data
             mutable autoPtr<Field<Type>> patchNeighbourFieldPtr_;
 
@@ -156,7 +134,7 @@ class cyclicAMIFvPatchField
         virtual bool all_ready() const;
 
         //- Use neighbour field caching
-        bool cacheNeighbourField() const;
+        static bool cacheNeighbourField();
 
         //- Return neighbour coupled internal cell data
         tmp<Field<Type>> getNeighbourField(const UList<Type>&) const;

src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothData.H

@@ -209,23 +209,6 @@ public:
             template<class TrackingData>
             inline bool equal(const smoothData&, TrackingData& td) const;
 
-            //- Interpolate between two values (lerp). Returns true if
-            //- causes changes. Not sure if needs to be specialised between
-            //- face and cell and what index is needed...
-            template<class TrackingData>
-            inline bool interpolate
-            (
-                const polyMesh&,
-                const point& pt,
-                const label i0,
-                const smoothData& f0,
-                const label i1,
-                const smoothData& f1,
-                const scalar weight,
-                const scalar tol,
-                TrackingData& td
-            );
-
 
     // Member Operators
 

src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothDataI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
-    Copyright (C) 2020,2025 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -191,45 +191,6 @@ inline bool Foam::smoothData::equal
 }
 
 
-template<class TrackingData>
-inline bool Foam::smoothData::interpolate
-(
-    const polyMesh&,
-    const point& pt,
-    const label i0,
-    const smoothData& f0,
-    const label i1,
-    const smoothData& f1,
-    const scalar weight,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    // TBD. What to interpolate? Do we have a position? cell/face centre?
-
-    if (!valid(td))
-    {
-        if (f0.valid(td))
-        {
-            operator=(f0);
-        }
-        else
-        {
-            operator=(f1);
-        }
-    }
-    else if (f0.valid(td))
-    {
-        operator=(f0);
-    }
-    else
-    {
-        operator=(f1);
-    }
-    return true;
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::smoothData::operator==

src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepData.H

@@ -207,23 +207,6 @@ public:
             template<class TrackingData>
             inline bool equal(const sweepData&, TrackingData& td) const;
 
-            //- Interpolate between two values (lerp). Returns true if
-            //- causes changes. Not sure if needs to be specialised between
-            //- face and cell and what index is needed...
-            template<class TrackingData>
-            inline bool interpolate
-            (
-                const polyMesh&,
-                const point& pt,
-                const label i0,
-                const sweepData& f0,
-                const label i1,
-                const sweepData& f1,
-                const scalar weight,
-                const scalar tol,
-                TrackingData& td
-            );
-
 
     // Member Operators
 

src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepDataI.H

@@ -210,45 +210,6 @@ inline bool Foam::sweepData::equal
 }
 
 
-template<class TrackingData>
-inline bool Foam::sweepData::interpolate
-(
-    const polyMesh&,
-    const point& pt,
-    const label i0,
-    const sweepData& f0,
-    const label i1,
-    const sweepData& f1,
-    const scalar weight,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    // TBD. What to interpolate? Do we have a position? cell/face centre?
-
-    if (!valid(td))
-    {
-        if (f0.valid(td))
-        {
-            operator=(f0);
-        }
-        else
-        {
-            operator=(f1);
-        }
-    }
-    else if (f0.valid(td))
-    {
-        operator=(f0);
-    }
-    else
-    {
-        operator=(f1);
-    }
-    return true;
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::sweepData::operator==

src/finiteVolume/finiteVolume/gradSchemes/LeastSquaresGrad/LeastSquaresGrad.C

@@ -1,4 +1,4 @@
-/*---------------------------------------------------------------------------*\
+  /*---------------------------------------------------------------------------*\
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
@@ -36,23 +36,18 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type, class Stencil>
-Foam::tmp
-<
-    Foam::GeometricField
-    <
-        typename Foam::outerProduct<Foam::vector, Type>::type,
-        Foam::fvPatchField,
-    Foam::volMesh
-    >
->
-Foam::fv::LeastSquaresGrad<Type, Stencil>::calcGrad
+void Foam::fv::LeastSquaresGrad<Type, Stencil>::calcGrad
 (
-    const GeometricField<Type, fvPatchField, volMesh>& vtf,
-    const word& name
+    GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& lsGrad,
+    const GeometricField<Type, fvPatchField, volMesh>& vtf
 ) const
 {
     typedef typename outerProduct<vector, Type>::type GradType;
-    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
 
     const fvMesh& mesh = vtf.mesh();
 
@@ -62,24 +57,7 @@ Foam::fv::LeastSquaresGrad<Type, Stencil>::calcGrad
         mesh
     );
 
-    tmp<GradFieldType> tlsGrad
-    (
-        new GradFieldType
-        (
-            IOobject
-            (
-                name,
-                vtf.instance(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensioned<GradType>(vtf.dimensions()/dimLength, Zero),
-            fvPatchFieldBase::extrapolatedCalculatedType()
-        )
-    );
-    GradFieldType& lsGrad = tlsGrad.ref();
+    lsGrad = dimensioned<GradType>(vtf.dimensions()/dimLength, Zero);
     Field<GradType>& lsGradIf = lsGrad;
 
     const extendedCentredCellToCellStencil& stencil = lsv.stencil();
@@ -131,6 +109,50 @@ Foam::fv::LeastSquaresGrad<Type, Stencil>::calcGrad
     // Correct the boundary conditions
     lsGrad.correctBoundaryConditions();
     gaussGrad<Type>::correctBoundaryConditions(vtf, lsGrad);
+}
+
+
+template<class Type, class Stencil>
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+    Foam::volMesh
+    >
+>
+Foam::fv::LeastSquaresGrad<Type, Stencil>::calcGrad
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vtf,
+    const word& name
+) const
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    const fvMesh& mesh = vtf.mesh();
+
+    tmp<GradFieldType> tlsGrad
+    (
+        new GradFieldType
+        (
+            IOobject
+            (
+                name,
+                vtf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            vtf.dimensions()/dimLength,
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+    GradFieldType& lsGrad = tlsGrad.ref();
+
+    calcGrad(lsGrad, vtf);
 
     return tlsGrad;
 }

src/finiteVolume/finiteVolume/gradSchemes/LeastSquaresGrad/LeastSquaresGrad.H

@@ -132,6 +132,14 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& vsf,
             const word& name
         ) const;
+
+        //- Calculate the grad of the given field into supplied field
+        virtual void calcGrad
+        (
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>& res,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
 };
 
 

src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C

@@ -38,19 +38,21 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type>
-Foam::tmp
-<
-    Foam::GeometricField
-    <
-        typename Foam::outerProduct<Foam::vector, Type>::type,
-        Foam::fvPatchField,
-        Foam::volMesh
-    >
->
-Foam::fv::fourthGrad<Type>::calcGrad
+void Foam::fv::fourthGrad<Type>::calcGrad
 (
+    GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& fGrad,
     const GeometricField<Type, fvPatchField, volMesh>& vsf,
-    const word& name
+    const GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& secondfGrad
 ) const
 {
     // The fourth-order gradient is calculated in two passes.  First,
@@ -59,37 +61,11 @@ Foam::fv::fourthGrad<Type>::calcGrad
     // gradient to complete the accuracy.
 
     typedef typename outerProduct<vector, Type>::type GradType;
-    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
 
     const fvMesh& mesh = vsf.mesh();
 
-    // Assemble the second-order least-square gradient
-    // Calculate the second-order least-square gradient
-    tmp<GradFieldType> tsecondfGrad
-      = leastSquaresGrad<Type>(mesh).grad
-        (
-            vsf,
-            "leastSquaresGrad(" + vsf.name() + ")"
-        );
-    const GradFieldType& secondfGrad =
-        tsecondfGrad();
-
-    tmp<GradFieldType> tfGrad
-    (
-        new GradFieldType
-        (
-            IOobject
-            (
-                name,
-                vsf.instance(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            secondfGrad
-        )
-    );
-    GradFieldType& fGrad = tfGrad.ref();
+    // Initialise to gradient
+    fGrad = secondfGrad;
 
     const vectorField& C = mesh.C();
 
@@ -158,9 +134,102 @@ Foam::fv::fourthGrad<Type>::calcGrad
 
     fGrad.correctBoundaryConditions();
     gaussGrad<Type>::correctBoundaryConditions(vsf, fGrad);
+}
+
+
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
+    >
+>
+Foam::fv::fourthGrad<Type>::calcGrad
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vsf,
+    const word& name
+) const
+{
+    // The fourth-order gradient is calculated in two passes.  First,
+    // the standard least-square gradient is assembled.  Then, the
+    // fourth-order correction is added to the second-order accurate
+    // gradient to complete the accuracy.
+
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    const fvMesh& mesh = vsf.mesh();
+
+    // Assemble the second-order least-square gradient
+    // Calculate the second-order least-square gradient
+    tmp<GradFieldType> tsecondfGrad
+      = leastSquaresGrad<Type>(mesh).grad
+        (
+            vsf,
+            "leastSquaresGrad(" + vsf.name() + ")"
+        );
+    const GradFieldType& secondfGrad = tsecondfGrad();
+
+    tmp<GradFieldType> tfGrad
+    (
+        new GradFieldType
+        (
+            IOobject
+            (
+                name,
+                vsf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            vsf.dimensions()/dimLength,
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+
+    calcGrad(tfGrad.ref(), vsf, secondfGrad);
 
     return tfGrad;
 }
 
 
+template<class Type>
+void Foam::fv::fourthGrad<Type>::calcGrad
+(
+    GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& fGrad,
+    const GeometricField<Type, fvPatchField, volMesh>& vsf
+) const
+{
+    // The fourth-order gradient is calculated in two passes.  First,
+    // the standard least-square gradient is assembled.  Then, the
+    // fourth-order correction is added to the second-order accurate
+    // gradient to complete the accuracy.
+
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    const fvMesh& mesh = vsf.mesh();
+
+    // Assemble the second-order least-square gradient
+    // Calculate the second-order least-square gradient
+    tmp<GradFieldType> tsecondfGrad
+      = leastSquaresGrad<Type>(mesh).grad
+        (
+            vsf,
+            "leastSquaresGrad(" + vsf.name() + ")"
+        );
+
+    calcGrad(fGrad, vsf, tsecondfGrad());
+}
+
+
 // ************************************************************************* //

src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.H

@@ -99,6 +99,32 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& vsf,
             const word& name
         ) const;
+
+        //- Calculate the grad of the given field into supplied field
+        virtual void calcGrad
+        (
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>& res,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+        //- Helper : calculate the grad of the given field into supplied field
+        void calcGrad
+        (
+            GeometricField
+            <
+                typename outerProduct<vector, Type>::type,
+                fvPatchField,
+                volMesh
+            >& fGrad,
+            const GeometricField<Type, fvPatchField, volMesh>& vsf,
+            const GeometricField
+            <
+                typename outerProduct<vector, Type>::type,
+                fvPatchField,
+                volMesh
+            >& secondfGrad
+        ) const;
 };
 
 

src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C

@@ -140,6 +140,67 @@ Foam::fv::gaussGrad<Type>::calcGrad
 }
 
 
+template<class Type>
+void Foam::fv::gaussGrad<Type>::calcGrad
+(
+    GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& gGrad,
+    const GeometricField<Type, fvPatchField, volMesh>& vsf
+) const
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+
+    DebugPout<< "gaussGrad<Type>::calcGrad on " << vsf.name()
+        << " into " << gGrad.name() << endl;
+
+    const fvMesh& mesh = vsf.mesh();
+    const labelUList& owner = mesh.owner();
+    const labelUList& neighbour = mesh.neighbour();
+    const vectorField& Sf = mesh.Sf();
+
+    const auto tssf(tinterpScheme_().interpolate(vsf));
+    const auto& ssf = tssf();
+
+    gGrad = dimensioned<GradType>(vsf.dimensions()/dimLength, Zero);
+
+    Field<GradType>& igGrad = gGrad;
+    const Field<Type>& issf = ssf;
+
+    forAll(owner, facei)
+    {
+        const GradType Sfssf = Sf[facei]*issf[facei];
+
+        igGrad[owner[facei]] += Sfssf;
+        igGrad[neighbour[facei]] -= Sfssf;
+    }
+
+    forAll(mesh.boundary(), patchi)
+    {
+        const labelUList& pFaceCells =
+            mesh.boundary()[patchi].faceCells();
+
+        const vectorField& pSf = mesh.Sf().boundaryField()[patchi];
+
+        const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];
+
+        forAll(mesh.boundary()[patchi], facei)
+        {
+            igGrad[pFaceCells[facei]] += pSf[facei]*pssf[facei];
+        }
+    }
+
+    igGrad /= mesh.V();
+
+    gGrad.correctBoundaryConditions();
+
+    correctBoundaryConditions(vsf, gGrad);
+}
+
+
 template<class Type>
 void Foam::fv::gaussGrad<Type>::correctBoundaryConditions
 (

src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.H

@@ -145,6 +145,14 @@ public:
             const word& name
         ) const;
 
+        //- Calculate the grad of the given field into supplied field
+        virtual void calcGrad
+        (
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>& res,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
         //- Correct the boundary values of the gradient using the patchField
         //- snGrad functions
         static void correctBoundaryConditions

src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C

@@ -96,7 +96,7 @@ Foam::fv::gradScheme<Type>::grad
     GradFieldType* pgGrad =
         mesh().objectRegistry::template getObjectPtr<GradFieldType>(name);
 
-    if (!this->mesh().cache(name) || this->mesh().changing())
+    if (!this->mesh().cache(name) || this->mesh().topoChanging())
     {
         // Delete any old occurrences to avoid double registration
         if (pgGrad && pgGrad->ownedByRegistry())
@@ -119,17 +119,14 @@ Foam::fv::gradScheme<Type>::grad
     }
     else
     {
-        if (pgGrad->upToDate(vsf))
+        if (pgGrad->upToDate(vsf) && this->mesh().upToDatePoints(*pgGrad))
         {
             solution::cachePrintMessage("Reusing", name, vsf);
         }
         else
         {
             solution::cachePrintMessage("Updating", name, vsf);
-            delete pgGrad;
-
-            pgGrad = calcGrad(vsf, name).ptr();
-            regIOobject::store(pgGrad);
+            calcGrad(*pgGrad, vsf);
         }
     }
 
@@ -180,4 +177,23 @@ Foam::fv::gradScheme<Type>::grad
 }
 
 
+template<class Type>
+void Foam::fv::gradScheme<Type>::calcGrad
+(
+    GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& result,
+    const GeometricField<Type, fvPatchField, volMesh>& vsf
+) const
+{
+    DebugPout<< "gradScheme<Type>::calcGrad on " << vsf.name()
+        << " into " << result.name() << endl;
+
+    result = calcGrad(vsf, result.name());
+}
+
+
 // ************************************************************************* //

src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H

@@ -142,6 +142,14 @@ public:
             const word& name
         ) const = 0;
 
+        //- Calculate the grad of the given field into supplied field
+        virtual void calcGrad
+        (
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>& res,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
         //- Calculate and return the grad of the given field
         //- which may have been cached
         tmp

src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.C

@@ -62,40 +62,25 @@ void Foam::fv::cellLimitedGrad<Type, Limiter>::limitGradient
 
 
 template<class Type, class Limiter>
-Foam::tmp
-<
-    Foam::GeometricField
-    <
-        typename Foam::outerProduct<Foam::vector, Type>::type,
-        Foam::fvPatchField,
-        Foam::volMesh
-    >
->
-Foam::fv::cellLimitedGrad<Type, Limiter>::calcGrad
+void Foam::fv::cellLimitedGrad<Type, Limiter>::calcGrad
 (
-    const GeometricField<Type, fvPatchField, volMesh>& vsf,
-    const word& name
-) const
-{
-    const fvMesh& mesh = vsf.mesh();
-
-    tmp
-    <
-        GeometricField
-        <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-    > tGrad = basicGradScheme_().calcGrad(vsf, name);
-
-    if (k_ < SMALL)
-    {
-        return tGrad;
-    }
-
     GeometricField
     <
         typename outerProduct<vector, Type>::type,
         fvPatchField,
         volMesh
-    >& g = tGrad.ref();
+    >& g,
+    const GeometricField<Type, fvPatchField, volMesh>& vsf
+) const
+{
+    const fvMesh& mesh = vsf.mesh();
+
+    basicGradScheme_().calcGrad(g, vsf);
+
+    if (k_ < SMALL)
+    {
+        return;
+    }
 
     const labelUList& owner = mesh.owner();
     const labelUList& neighbour = mesh.neighbour();
@@ -227,6 +212,49 @@ Foam::fv::cellLimitedGrad<Type, Limiter>::calcGrad
     limitGradient(limiter, g);
     g.correctBoundaryConditions();
     gaussGrad<Type>::correctBoundaryConditions(vsf, g);
+}
+
+
+template<class Type, class Limiter>
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
+    >
+>
+Foam::fv::cellLimitedGrad<Type, Limiter>::calcGrad
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vsf,
+    const word& name
+) const
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    const fvMesh& mesh = vsf.mesh();
+
+    tmp<GradFieldType> tGrad
+    (
+        new GradFieldType
+        (
+            IOobject
+            (
+                name,
+                vsf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            vsf.dimensions()/dimLength,
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+    
+    calcGrad(tGrad.ref(), vsf);
 
     return tGrad;
 }

src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H

@@ -150,6 +150,14 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& vsf,
             const word& name
         ) const;
+
+        //- Calculate the grad of the given field into supplied field
+        virtual void calcGrad
+        (
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>& res,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
 };
 
 

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"

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
+       );
 };
 
 

src/finiteVolume/fvMesh/fvPatches/constraint/cyclicACMI/cyclicACMIFvPatch.H

@@ -220,14 +220,9 @@ public:
                 (
                     const Field<Type>& fld,
                     labelRange& sendRequests,
-                    labelRange& recvRequests,
                     PtrList<List<Type>>& sendBuffers,
-                    PtrList<List<Type>>& recvBuffers,
-
-                    labelRange& sendRequests1,
-                    labelRange& recvRequests1,
-                    PtrList<List<Type>>& sendBuffers1,
-                    PtrList<List<Type>>& recvBuffers1
+                    labelRange& recvRequests,
+                    PtrList<List<Type>>& recvBuffers
                 ) const
                 {
                     // Make sure areas are up-to-date
@@ -237,14 +232,9 @@ public:
                     (
                         fld,
                         sendRequests,
-                        recvRequests,
                         sendBuffers,
-                        recvBuffers,
-
-                        sendRequests1,
-                        recvRequests1,
-                        sendBuffers1,
-                        recvBuffers1
+                        recvRequests,
+                        recvBuffers
                     );
                 }
 
@@ -253,9 +243,7 @@ public:
                 (
                     const Field<Type>& localFld,
                     const labelRange& requests,  // The receive requests
-                    const PtrList<List<Type>>& recvBuffers,
-                    const labelRange& requests1,  // The receive requests
-                    const PtrList<List<Type>>& recvBuffers1
+                    const PtrList<List<Type>>& recvBuffers
                 ) const
                 {
                     return cyclicACMIPolyPatch_.interpolate
@@ -263,8 +251,6 @@ public:
                         localFld,
                         requests,
                         recvBuffers,
-                        requests1,
-                        recvBuffers1,
                         UList<Type>()
                     );
                 }

src/fused/finiteVolume/fusedGaussGrad.C

@@ -226,6 +226,94 @@ Foam::fv::fusedGaussGrad<Type>::calcGrad
 }
 
 
+template<class Type>
+void Foam::fv::fusedGaussGrad<Type>::calcGrad
+(
+    GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& gGrad,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+
+    const fvMesh& mesh = vf.mesh();
+
+    DebugPout<< "fusedGaussGrad<Type>::calcGrad on " << vf.name()
+        << " into " << gGrad.name() << endl;
+
+    gGrad = dimensioned<GradType>(vf.dimensions()/dimLength, Zero);
+
+    if (this->tinterpScheme_().corrected())
+    {
+        const auto tfaceCorr(this->tinterpScheme_().correction(vf));
+        auto& faceCorr = tfaceCorr();
+
+        DebugPout<< "fusedGaussGrad<Type>::calcGrad corrected interpScheme "
+            << this->tinterpScheme_().type() << endl;
+
+        const auto interpolate = [&]
+        (
+            const vector& area,
+            const scalar lambda,
+
+            const Type& ownVal,
+            const Type& neiVal,
+
+            const Type& correction
+
+        ) -> GradType
+        {
+            return area*((lambda*(ownVal - neiVal) + neiVal) + correction);
+        };
+
+        fvc::surfaceSum
+        (
+            this->tinterpScheme_().weights(vf),
+            vf,
+            faceCorr,
+            interpolate,
+            gGrad,
+            false
+        );
+    }
+    else
+    {
+        DebugPout<< "fusedGaussGrad<Type>::calcGrad uncorrected interpScheme "
+            << this->tinterpScheme_().type() << endl;
+
+        const auto interpolate = [&]
+        (
+            const vector& area,
+            const scalar lambda,
+            const Type& ownVal,
+            const Type& neiVal
+        ) -> GradType
+        {
+            return area*(lambda*(ownVal - neiVal) + neiVal);
+        };
+
+        fvc::surfaceSum
+        (
+            tinterpScheme_().weights(vf),
+            vf,
+            interpolate,
+            gGrad,
+            false
+        );
+    }
+
+    gGrad.primitiveFieldRef() /= mesh.V();
+
+    gGrad.correctBoundaryConditions();
+
+    correctBoundaryConditions(vf, gGrad);
+}
+
+
 template<class Type>
 template<class GradType>
 void Foam::fv::fusedGaussGrad<Type>::correctBoundaryConditions

src/fused/finiteVolume/fusedGaussGrad.H

@@ -145,6 +145,14 @@ public:
             const word& name
         ) const;
 
+        //- Calculate the grad of the given field into supplied field
+        virtual void calcGrad
+        (
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>& res,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
         //- Correct the boundary values of the gradient using the patchField
         //- snGrad functions
         template<class GradType>

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

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>&
-);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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
 <

src/mesh/snappyHexMesh/meshRefinement/wallPoints.H

@@ -251,23 +251,6 @@ public:
             template<class TrackingData>
             inline bool equal(const wallPoints&, TrackingData&) const;
 
-            //- Interpolate between two values (lerp). Returns true if
-            //- causes changes. Not sure if needs to be specialised between
-            //- face and cell and what index is needed...
-            template<class TrackingData>
-            inline bool interpolate
-            (
-                const polyMesh&,
-                const point& pt,
-                const label i0,
-                const wallPoints& f0,
-                const label i1,
-                const wallPoints& f1,
-                const scalar weight,
-                const scalar tol,
-                TrackingData& td
-            );
-
 
     // Member Operators
 

src/mesh/snappyHexMesh/meshRefinement/wallPointsI.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2025 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -443,41 +443,6 @@ inline bool Foam::wallPoints::equal
 }
 
 
-template<class TrackingData>
-inline bool Foam::wallPoints::interpolate
-(
-    const polyMesh&,
-    const point& pt,
-    const label i0,
-    const wallPoints& f0,
-    const label i1,
-    const wallPoints& f1,
-    const scalar weight,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    if (valid(td))
-    {
-        return false;
-    }
-    else if (f0.valid(td))
-    {
-        operator=(f0);
-        return true;
-    }
-    else if (f1.valid(td))
-    {
-        operator=(f1);
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::wallPoints::operator==

src/meshTools/AMIInterpolation/AMIInterpolation/AMICache.C

@@ -1,651 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    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/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "AMICache.H"
-#include "AMIInterpolation.H"
-#include "mathematicalConstants.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-Foam::scalar Foam::AMICache::cacheThetaTolerance_ = 1e-8;
-
-int Foam::AMICache::debug = 0;
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::scalar Foam::AMICache::getRotationAngle(const point& globalPoint) const
-{
-    if (!coordSysPtr_)
-    {
-        FatalErrorInFunction
-            << "No co-ordinate system available for theta evaluation"
-            << abort(FatalError);
-    }
-
-
-    scalar theta = coordSysPtr_->localPosition(globalPoint).y();
-
-    // Ensure 0 < theta < 2pi
-    if (mag(theta) < cacheThetaTolerance_)
-    {
-        theta = 0;
-    }
-    else if (theta < 0)
-    {
-        theta += constant::mathematical::twoPi;
-    }
-
-    return theta;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::AMICache::AMICache(const dictionary& dict, const bool toSource)
-:
-    size_(dict.getOrDefault<label>("cacheSize", 0)),
-    rotationAxis_(dict.getOrDefault<vector>("rotationAxis", Zero)),
-    rotationCentre_(dict.getOrDefault<point>("rotationCentre", Zero)),
-    maxThetaDeg_(dict.getOrDefault<scalar>("maxThetaDeg", 5)),
-    complete_(false),
-    index0_(-1),
-    index1_(-1),
-    interpWeight_(0),
-    coordSysPtr_(nullptr),
-    theta_(),
-    cachedSrcAddress_(),
-    cachedSrcWeights_(),
-    cachedSrcWeightsSum_(),
-    cachedSrcMapPtr_(),
-    cachedTgtAddress_(),
-    cachedTgtWeights_(),
-    cachedTgtWeightsSum_(),
-    cachedTgtMapPtr_(),
-    toSource_(toSource)
-{
-    if (size_ != 0)
-    {
-        theta_.resize(size_, GREAT);
-        cachedSrcAddress_.resize(size_);
-        cachedSrcWeights_.resize(size_);
-        cachedSrcWeightsSum_.resize(size_);
-        cachedSrcMapPtr_.resize(size_);
-        cachedTgtAddress_.resize(size_);
-        cachedTgtWeights_.resize(size_);
-        cachedTgtWeightsSum_.resize(size_);
-        cachedTgtMapPtr_.resize(size_);
-    }
-}
-
-
-Foam::AMICache::AMICache(const bool toSource)
-:
-    size_(0),
-    rotationAxis_(Zero),
-    rotationCentre_(Zero),
-    complete_(false),
-    index0_(-1),
-    index1_(-1),
-    interpWeight_(0),
-    coordSysPtr_(nullptr),
-    theta_(),
-    cachedSrcAddress_(),
-    cachedSrcWeights_(),
-    cachedSrcWeightsSum_(),
-    cachedSrcMapPtr_(),
-    cachedTgtAddress_(),
-    cachedTgtWeights_(),
-    cachedTgtWeightsSum_(),
-    cachedTgtMapPtr_(),
-    toSource_(toSource)
-{}
-
-
-Foam::AMICache::AMICache(const AMICache& cache)
-:
-    size_(cache.size_),
-    rotationAxis_(cache.rotationAxis_),
-    rotationCentre_(cache.rotationCentre_),
-    complete_(cache.complete_),
-    index0_(cache.index0_),
-    index1_(cache.index1_),
-    interpWeight_(cache.interpWeight_),
-    coordSysPtr_(nullptr),  // Need to clone as cylindricalCS
-    theta_(cache.theta_),
-    cachedSrcAddress_(cache.cachedSrcAddress_),
-    cachedSrcWeights_(cache.cachedSrcWeights_),
-    cachedSrcWeightsSum_(cache.cachedSrcWeightsSum_),
-    cachedSrcMapPtr_(cache.cachedSrcMapPtr_.size()),  // Need to clone
-    cachedTgtAddress_(cache.cachedTgtAddress_),
-    cachedTgtWeights_(cache.cachedTgtWeights_),
-    cachedTgtWeightsSum_(cache.cachedTgtWeightsSum_),
-    cachedTgtMapPtr_(cache.cachedTgtMapPtr_.size()),  // Need to clone
-    toSource_(cache.toSource_)
-{
-    if (cache.coordSysPtr_)
-    {
-        coordSysPtr_.reset(new coordSystem::cylindrical(cache.coordSysPtr_()));
-    }
-
-    forAll(cachedSrcMapPtr_, cachei)
-    {
-        cachedSrcMapPtr_[cachei].reset(cache.cachedSrcMapPtr_[cachei].clone());
-    }
-
-    forAll(cachedTgtMapPtr_, cachei)
-    {
-        cachedTgtMapPtr_[cachei].reset(cache.cachedTgtMapPtr_[cachei].clone());
-    }
-}
-
-
-Foam::AMICache::AMICache
-(
-    const AMICache& cache,
-    const AMIInterpolation& fineAMI,
-    const labelList& sourceRestrictAddressing,
-    const labelList& targetRestrictAddressing
-)
-:
-    size_(cache.size_),
-    rotationAxis_(cache.rotationAxis_),
-    rotationCentre_(cache.rotationCentre_),
-    complete_(cache.complete_),
-    index0_(cache.index0_),
-    index1_(cache.index1_),
-    interpWeight_(cache.interpWeight_),
-    coordSysPtr_(nullptr),
-    theta_(cache.theta_),
-    cachedSrcAddress_(cache.size_),
-    cachedSrcWeights_(cache.size_),
-    cachedSrcWeightsSum_(cache.size_),
-    cachedSrcMapPtr_(cache.size_),
-    cachedTgtAddress_(cache.size_),
-    cachedTgtWeights_(cache.size_),
-    cachedTgtWeightsSum_(cache.size_),
-    cachedTgtMapPtr_(cache.size_),
-    toSource_(cache.toSource_)
-{
-    if (size_ > 0 && fineAMI.comm() != -1)
-    {
-        for (label cachei : {index0_, index1_})
-        {
-            if (cachei == -1) continue;
-
-            scalarField dummySrcMagSf;
-            labelListList srcAddress;
-            scalarListList srcWeights;
-            scalarField srcWeightsSum;
-            autoPtr<mapDistribute> tgtMapPtr;
-
-            AMIInterpolation::agglomerate
-            (
-                cache.cachedTgtMapPtr()[cachei],
-                fineAMI.srcMagSf(),
-                cache.cachedSrcAddress()[cachei],
-                cache.cachedSrcWeights()[cachei],
-
-                sourceRestrictAddressing,
-                targetRestrictAddressing,
-
-                dummySrcMagSf,
-                srcAddress,
-                srcWeights,
-                srcWeightsSum,
-                tgtMapPtr,
-                fineAMI.comm()
-            );
-
-            scalarField dummyTgtMagSf;
-            labelListList tgtAddress;
-            scalarListList tgtWeights;
-            scalarField tgtWeightsSum;
-            autoPtr<mapDistribute> srcMapPtr;
-
-            AMIInterpolation::agglomerate
-            (
-                cache.cachedSrcMapPtr()[cachei],
-                fineAMI.tgtMagSf(),
-                cache.cachedTgtAddress()[cachei],
-                cache.cachedTgtWeights()[cachei],
-
-                targetRestrictAddressing,
-                sourceRestrictAddressing,
-
-                dummyTgtMagSf,
-                tgtAddress,
-                tgtWeights,
-                tgtWeightsSum,
-                srcMapPtr,
-                fineAMI.comm()
-            );
-
-            cachedSrcAddress_[cachei] = srcAddress;
-            cachedSrcWeights_[cachei] = srcWeights;
-            cachedSrcWeightsSum_[cachei] = srcWeightsSum;
-            cachedSrcMapPtr_[cachei] = srcMapPtr.clone();
-
-            cachedTgtAddress_[cachei] = tgtAddress;
-            cachedTgtWeights_[cachei] = tgtWeights;
-            cachedTgtWeightsSum_[cachei] = tgtWeightsSum;
-            cachedTgtMapPtr_[cachei] = tgtMapPtr.clone();
-        }
-    }
-}
-
-
-Foam::AMICache::AMICache(Istream& is)
-:
-    size_(readLabel(is)),
-
-    rotationAxis_(is),
-    rotationCentre_(is),
-    maxThetaDeg_(readScalar(is)),
-
-    complete_(readBool(is)),
-
-    index0_(-1),
-    index1_(-1),
-    interpWeight_(0),
-    coordSysPtr_(nullptr),
-    theta_(),
-    cachedSrcAddress_(),
-    cachedSrcWeights_(),
-    cachedSrcWeightsSum_(),
-    cachedSrcMapPtr_(),
-    cachedTgtAddress_(),
-    cachedTgtWeights_(),
-    cachedTgtWeightsSum_(),
-    cachedTgtMapPtr_()
-{
-    const bitSet goodMap(is);
-
-    if (goodMap.size())
-    {
-        is >> index0_
-           >> index1_
-           >> interpWeight_
-           >> theta_;
-
-        const bool goodCoord(readBool(is));
-        if (goodCoord)
-        {
-            coordSysPtr_.reset(new coordSystem::cylindrical(is));
-        }
-
-        is >> cachedSrcAddress_
-           >> cachedSrcWeights_
-           >> cachedSrcWeightsSum_;
-
-        cachedSrcMapPtr_.setSize(goodMap.size());
-        forAll(goodMap, cachei)
-        {
-            if (goodMap[cachei])
-            {
-                cachedSrcMapPtr_[cachei].reset(new mapDistribute(is));
-            }
-        }
-
-        is >> cachedTgtAddress_
-           >> cachedTgtWeights_
-           >> cachedTgtWeightsSum_;
-
-        cachedTgtMapPtr_.setSize(goodMap.size());
-        forAll(goodMap, cachei)
-        {
-            if (goodMap[cachei])
-            {
-                cachedTgtMapPtr_[cachei].reset(new mapDistribute(is));
-            }
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::AMICache::addToCache
-(
-    const AMIInterpolation& ami,
-    const point& globalPoint
-)
-{
-    DebugPout<< "-- addToCache" << endl;
-
-
-    if (!active())
-    {
-        DebugInfo<< "-- addToCache - deactivated" << endl;
-        return;
-    }
-
-    if (!coordSysPtr_)
-    {
-        DebugInfo
-            << "Creating rotation co-ordinate system:"
-            << " rotationCentre:" << rotationCentre_
-            << " rotationAxis:" << rotationAxis_
-            << " p:" << globalPoint
-            << endl;
-
-        coordSysPtr_.reset
-        (
-            new coordSystem::cylindrical(rotationCentre_, rotationAxis_)
-        );
-        DebugPout<< "Coord sys:" << coordSysPtr_() << endl;
-    }
-
-    // Check if cache is complete
-    if (!complete_)
-    {
-        for (const scalar angle : theta_)
-        {
-            // Check against null value; if any are present, cache is incomplete
-            if (angle > constant::mathematical::twoPi)
-            {
-                complete_ = false;
-                break;
-            }
-        }
-    }
-
-    if (!complete_)
-    {
-        const scalar theta = getRotationAngle(globalPoint);
-
-        const label bini = theta/constant::mathematical::twoPi*size_;
-
-        DebugPout<< "  -- bini:" << bini << " for theta:" << theta << endl;
-
-        // Check if already have entry for this bin
-        if (theta_[bini] > constant::mathematical::twoPi)
-        {
-            DebugPout<< "  -- setting cache at index " << bini << endl;
-
-            // New entry
-            theta_[bini] = theta;
-
-            cachedSrcAddress_[bini] = ami.srcAddress();
-            cachedSrcWeights_[bini] = ami.srcWeights();
-            cachedSrcWeightsSum_[bini] = ami.srcWeightsSum();
-
-            if (ami.hasSrcMap())
-            {
-                cachedSrcMapPtr_[bini] = ami.srcMap().clone();
-            }
-
-            cachedTgtAddress_[bini] = ami.tgtAddress();
-            cachedTgtWeights_[bini] = ami.tgtWeights();
-            cachedTgtWeightsSum_[bini] = ami.tgtWeightsSum();
-
-            if (ami.hasTgtMap())
-            {
-                cachedTgtMapPtr_[bini] = ami.tgtMap().clone();
-            }
-        }
-    }
-}
-
-
-bool Foam::AMICache::restoreCache(const point& globalPoint)
-{
-    DebugPout<< "-- restoreCache" << endl;
-
-    index0_ = -1;
-    index1_ = -1;
-    interpWeight_ = -1;
-
-    if (!coordSysPtr_ || size_ == -1)
-    {
-        return false;
-    }
-
-    const scalar theta = getRotationAngle(globalPoint);
-    const scalar twoPi = constant::mathematical::twoPi;
-    const label bini = theta/twoPi*size_;
-
-    DebugPout<< "  -- bini:" << bini << " for theta:" << theta << endl;
-
-    const auto validIndex = [&](const scalar bini)
-    {
-        return theta_[bini] < constant::mathematical::twoPi;
-    };
-
-    // Maximum angle in degrees for which to search for cached bins
-    const scalar maxThetaStencil = maxThetaDeg_*constant::mathematical::pi/180.0;
-
-    if
-    (
-        validIndex(bini)
-     && (
-            mag(theta - theta_[bini]) < cacheThetaTolerance_
-         || mag(theta - twoPi - theta_[bini]) < cacheThetaTolerance_
-        )
-    )
-    {
-        // Hit cached value - no interpolation needed
-        // index1_ = -1 indicates no interpolation
-        index0_ = bini;
-        index1_ = -1;
-        interpWeight_ = 0;
-
-        DebugInfo
-            << "  -- t0:" << theta_[index0_] << " theta:" << theta
-            << " i0:" << index0_ << " i1:" << index1_
-            << " w:" << interpWeight_ << endl;
-        return true;
-    }
-    else
-    {
-        // Find indices and values bracketing theta
-        const label nBin = theta_.size();
-
-        // Participating theta values and bin addresses
-        // - Note we add wrap-around values at start and end
-        DynamicList<scalar> thetap(nBin+2);
-        DynamicList<label> binAddresses(nBin+2);
-
-        // Initialise wrap-around values
-        thetap.push_back(0);
-        binAddresses.push_back(-1);
-        forAll(theta_, thetai)
-        {
-            if (validIndex(thetai))
-            {
-                thetap.push_back(theta_[thetai]);
-                binAddresses.push_back(thetai);
-            }
-        }
-
-        // Check that we have enough data points for interpolation
-        // - We added storage for lower wrap-around value, and we then need
-        //   at least 2 additional values for the interpolation
-        if (thetap.size() < 3)
-        {
-            DebugPout<< "  -- no cache available" << endl;
-            return false;
-        }
-
-        // Set wrap-around values if we have sufficient data
-        thetap[0] = thetap.last() - twoPi;
-        binAddresses[0] = binAddresses.last();
-        thetap.push_back(thetap[1] + twoPi);
-        binAddresses.push_back(binAddresses[1]);
-
-        // Find interpolation indices
-        label loweri = labelMax;
-        label i = 0;
-        while (i < thetap.size())
-        {
-            if (thetap[i] < theta)
-            {
-                loweri = i;
-            }
-            else
-            {
-                break;
-            }
-            ++i;
-        }
-
-        if (loweri == labelMax)
-        {
-            DebugPout<< "  -- no cache available" << endl;
-            return false;
-        }
-
-        label upperi = labelMax;
-        i = thetap.size() - 1;
-        while (i >= 0)
-        {
-            if (thetap[i] > theta)
-            {
-                upperi = i;
-            }
-            else
-            {
-                break;
-            }
-            --i;
-        }
-
-        if (upperi == labelMax)
-        {
-            DebugPout<< "  -- no cache available" << endl;
-            return false;
-        }
-
-        // Ensure distances are valid
-        if (upperi == loweri)
-        {
-            DebugPout
-                << "  -- no cache available: theta:" << theta
-                << " lower:" << loweri << " upper:" << upperi << endl;
-            return false;
-        }
-        if (mag(theta - thetap[loweri]) > maxThetaStencil)
-        {
-            DebugPout
-                << "  -- no cache available: theta:" << theta
-                << " lower:" << thetap[loweri] << endl;
-            return false;
-        }
-        if (mag(theta - thetap[upperi]) > maxThetaStencil)
-        {
-            DebugPout
-                << "  -- no cache available: theta:" << theta
-                << " upper:" << thetap[upperi] << endl;
-            return false;
-        }
-
-        index0_ = binAddresses[loweri];
-        index1_ = binAddresses[upperi];
-        interpWeight_ =
-            (theta - theta_[index0_])/(theta_[index1_] - theta_[index0_]);
-
-        DebugInfo
-            << theta_.size()
-            << "  -- t0:" << theta_[index0_] << " theta:" << theta
-            << " t1:" << theta_[index1_]
-            << " i0:" << index0_ << " i1:" << index1_
-            << " w:" << interpWeight_ << endl;
-
-
-        return true;
-    }
-
-    // If we get here then no valid cache found within stencil
-    DebugPout<< "  -- no cache available" << endl;
-    return false;
-}
-
-
-void Foam::AMICache::write(Ostream& os) const
-{
-    if (size_ > 0)
-    {
-        os.writeEntry("cacheSize", size_);
-        os.writeEntry("rotationAxis", rotationAxis_);
-        os.writeEntry("rotationCentre", rotationCentre_);
-        os.writeEntry("maxThetaDeg", maxThetaDeg_);
-    }
-}
-
-
-bool Foam::AMICache::writeData(Ostream& os) const
-{
-    os  << token::SPACE<< size_
-        << token::SPACE<< rotationAxis_
-        << token::SPACE<< rotationCentre_
-        << token::SPACE<< maxThetaDeg_
-        << token::SPACE<< complete_;
-
-    bitSet goodMap(cachedSrcMapPtr_.size());
-    forAll(goodMap, cachei)
-    {
-        goodMap.set(cachei, cachedSrcMapPtr_[cachei].good());
-    }
-    os  << token::SPACE << goodMap;
-
-    if (goodMap.size())
-    {
-        os  << token::SPACE << index0_
-            << token::SPACE << index1_
-            << token::SPACE << interpWeight_
-            << token::SPACE << theta_;
-
-        os  << token::SPACE << coordSysPtr_.good();
-
-        if (coordSysPtr_.good())
-        {
-            os  << token::SPACE << coordSysPtr_();
-        }
-
-        os  << token::SPACE << cachedSrcAddress_
-            << token::SPACE << cachedSrcWeights_
-            << token::SPACE << cachedSrcWeightsSum_;
-
-        for (const auto& index : goodMap)
-        {
-            os  << token::SPACE << cachedSrcMapPtr_[index]();
-        }
-
-        os  << token::SPACE << cachedTgtAddress_
-            << token::SPACE << cachedTgtWeights_
-            << token::SPACE << cachedTgtWeightsSum_;
-
-        for (const auto& index : goodMap)
-        {
-            os  << token::SPACE << cachedTgtMapPtr_[index]();
-        }
-    }
-
-    return true;
-}
-
-
-// ************************************************************************* //

src/meshTools/AMIInterpolation/AMIInterpolation/AMICache.H

@@ -1,388 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    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/>.
-
-Class
-    Foam::AMICache
-
-Description
-    Provides caching of weights and addressing to AMIInterpolation
-
-SourceFiles
-    AMICache.C
-
-SeeAlso
-    Foam::AMIInterpolation.H
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef Foam_AMICache_H
-#define Foam_AMICache_H
-
-#include "cylindricalCS.H"
-#include "mapDistribute.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-class AMIInterpolation;
-
-/*---------------------------------------------------------------------------*\
-                          Class AMICache Declaration
-\*---------------------------------------------------------------------------*/
-
-class AMICache
-{
-public:
-
-    // Public Data Types
-
-        //- Tolerance used when caching the AMI to identify e.g. if the
-        //- current rotation angle has already been captured
-        static scalar cacheThetaTolerance_;
-
-        static int debug;
-
-private:
-
-    // Private Data
-
-        //- Cache size
-        label size_;
-
-        //- Axis of rotation for rotational cyclics
-        vector rotationAxis_;
-
-        //- Point on axis of rotation for rotational cyclics
-        point rotationCentre_;
-
-        //- Maximum angle in degrees for which to search for cached bins
-        scalar maxThetaDeg_;
-
-        //- Flag to indicate that the cache is complete
-        bool complete_;
-
-        //- Cache index 0 in lists
-        label index0_;
-
-        //- Cache index 1 in lists
-        label index1_;
-
-        //- Interpolation weight between cache indices 0 and 1
-        scalar interpWeight_;
-
-        //- Local co-ordinate system
-        autoPtr<coordSystem::cylindrical> coordSysPtr_;
-
-        //- List of cached angle snapshots
-        List<scalar> theta_;
-
-        //- List of source addresses
-        List<labelListList> cachedSrcAddress_;
-
-        //- List of source weights
-        List<scalarListList> cachedSrcWeights_;
-
-        //- List of source weights sums
-        List<scalarField> cachedSrcWeightsSum_;
-
-        //- List of source parallel maps
-        List<autoPtr<mapDistribute>> cachedSrcMapPtr_;
-
-        //- List of target addresses
-        List<labelListList> cachedTgtAddress_;
-
-        //- List of target weights
-        List<scalarListList> cachedTgtWeights_;
-
-        //- List of target weights sums
-        List<scalarField> cachedTgtWeightsSum_;
-
-        //- List of target parallel maps
-        List<autoPtr<mapDistribute>> cachedTgtMapPtr_;
-
-
-        //- Flag to indicate interpolation direction
-        mutable bool toSource_;
-
-
-    // Private Member Functions
-
-        //- Get rotation angle for point using local co-ordinate system
-        scalar getRotationAngle(const point& globalPoint) const;
-
-
-public:
-
-    // Constructors
-
-        //- Null constructor
-        AMICache(const bool toSource = true);
-
-        //- Construct from dictionary
-        AMICache(const dictionary& dict, const bool toSource = true);
-
-        //- Construct as copy
-        AMICache(const AMICache& cache);
-
-        //- Construct from agglomeration of AMIInterpolation. Agglomeration
-        //- passed in as new coarse size and addressing from fine from coarse
-        AMICache
-        (
-            const AMICache& cache,
-            const AMIInterpolation& fineAMI,
-            const labelList& sourceRestrictAddressing,
-            const labelList& targetRestrictAddressing
-        );
-
-        //- Construct from stream
-        AMICache(Istream& is);
-
-
-    // Member Functions
-
-        //- Return true if cache is active
-        constexpr bool active() const noexcept { return size_ > 0; }
-
-        //- Return cache size
-        constexpr label size() const noexcept { return size_; }
-
-        //- Return true if cache is complete
-        constexpr bool complete() const noexcept { return complete_; }
-
-        //- Return cache lower bound index
-        constexpr label index0() const noexcept { return index0_; }
-
-        //- Return cache upper bound index
-        constexpr label index1() const noexcept { return index1_; }
-
-        //- Return cache interpolation weight
-        constexpr label weight() const noexcept { return interpWeight_; }
-
-        //- Return list of cached rotation angles
-        const List<scalar>& theta() const noexcept { return theta_; }
-
-        //- Return List of source addresses
-        const List<labelListList>& cachedSrcAddress() const noexcept
-        {
-            return cachedSrcAddress_;
-        }
-
-        //- Return List of source weights
-        const List<scalarListList>& cachedSrcWeights() const noexcept
-        {
-            return cachedSrcWeights_;
-        }
-
-        //- Return List of source weights sums
-        const List<scalarField>& cachedSrcWeightsSum() const noexcept
-        {
-            return cachedSrcWeightsSum_;
-        }
-
-        //- Return List of source parallel maps
-        const List<autoPtr<mapDistribute>>& cachedSrcMapPtr() const noexcept
-        {
-            return cachedSrcMapPtr_;
-        }
-
-        //- Return List of target addresses
-        const List<labelListList>& cachedTgtAddress() const noexcept
-        {
-            return cachedTgtAddress_;
-        }
-
-        //- Return List of target weights
-        const List<scalarListList>& cachedTgtWeights() const noexcept
-        {
-            return cachedTgtWeights_;
-        }
-
-        //- Return List of target weights sums
-        const List<scalarField>& cachedTgtWeightsSum() const noexcept
-        {
-            return cachedTgtWeightsSum_;
-        }
-
-        //- Return List of target parallel maps
-        const List<autoPtr<mapDistribute>>& cachedTgtMapPtr() const noexcept
-        {
-            return cachedTgtMapPtr_;
-        }
-
-        //- Apply cached evaluation based on user supplied evaluation function
-        template<class Type, class EvalFunction>
-        bool apply(List<Type>& result, const EvalFunction& eval) const;
-
-        //- Flag that lower bound is applicable
-        constexpr bool applyLower() const noexcept
-        {
-            return index0_ != -1 && index1_ == -1;
-        }
-
-        //- Flag that upper bound is applicable
-        constexpr bool applyUpper() const noexcept
-        {
-            return index0_ == -1 && index1_ != -1;
-        }
-
-        //- Flag that interpolation is applicable
-        constexpr bool applyInterpolate() const noexcept
-        {
-            return index0_ != -1 && index1_ != -1;
-        }
-
-        //- Set the interpolation direction
-        constexpr bool setDirection(bool toSource) const noexcept
-        {
-            toSource_ = toSource;
-            return toSource_;
-        }
-
-        //- Restore AMI weights and addressing from the cache
-        bool restoreCache(const point& globalPoint);
-
-        //- Add AMI weights and addressing to the cache
-        void addToCache
-        (
-            const AMIInterpolation& ami,
-            const point& globalPoint
-        );
-
-        //- Check cache index is valid
-        void checkIndex(const label index) const
-        {
-            if ((index < 0) || (index > size_-1))
-            {
-                FatalErrorInFunction
-                    << "Supplied out of bounds: " << index << "/"
-                    << size_ << abort(FatalError);
-            }
-        }
-
-
-        // Helper functions to retrieve cached values at index0 and index1
-        // Note: uses interpolation direction toSource_
-        #undef  defineMethods01
-        #define defineMethods01(Src, Tgt, idx)                                 \
-            const labelListList& c##Src##Address##idx() const                  \
-            {                                                                  \
-                checkIndex(index##idx##_);                                     \
-                return toSource_ ?                                             \
-                    cached##Src##Address_[index##idx##_]                       \
-                  : cached##Tgt##Address_[index##idx##_];                      \
-            }                                                                  \
-            const scalarListList& c##Src##Weights##idx() const                 \
-            {                                                                  \
-                checkIndex(index##idx##_);                                     \
-                return toSource_ ?                                             \
-                    cached##Src##Weights_[index##idx##_]                       \
-                  : cached##Tgt##Weights_[index##idx##_];                      \
-            }                                                                  \
-            const scalarField& c##Src##WeightsSum##idx() const                 \
-            {                                                                  \
-                checkIndex(index##idx##_);                                     \
-                return toSource_ ?                                             \
-                    cached##Src##WeightsSum_[index##idx##_]                    \
-                  : cached##Tgt##WeightsSum_[index##idx##_];                   \
-            }                                                                  \
-            const autoPtr<mapDistribute>& c##Src##MapPtr##idx() const          \
-            {                                                                  \
-                checkIndex(index##idx##_);                                     \
-                return toSource_ ?                                             \
-                    cached##Src##MapPtr_[index##idx##_]                        \
-                  : cached##Tgt##MapPtr_[index##idx##_];                       \
-            }
-
-        defineMethods01(Src, Tgt, 0)
-        defineMethods01(Src, Tgt, 1)
-        defineMethods01(Tgt, Src, 0)
-        defineMethods01(Tgt, Src, 1)
-
-
-        // Helper functions to retrieve cached values at supplied index
-        // Note: uses interpolation direction toSource_
-        #undef  defineMethodsIndex
-        #define defineMethodsIndex(Src, Tgt)                                   \
-            const labelListList& c##Src##Address(const label index) const      \
-            {                                                                  \
-                checkIndex(index);                                             \
-                return toSource_ ?                                             \
-                    cached##Src##Address_[index]                               \
-                  : cached##Tgt##Address_[index];                              \
-            }                                                                  \
-            const scalarListList& c##Src##Weights(const label index) const     \
-            {                                                                  \
-                checkIndex(index);                                             \
-                return toSource_ ?                                             \
-                    cached##Src##Weights_[index]                               \
-                  : cached##Tgt##Weights_[index];                              \
-            }                                                                  \
-            const scalarField& c##Src##WeightsSum(const label index) const     \
-            {                                                                  \
-                checkIndex(index);                                             \
-                return toSource_ ?                                             \
-                    cached##Src##WeightsSum_[index]                            \
-                  : cached##Tgt##WeightsSum_[index];                           \
-            }                                                                  \
-            const autoPtr<mapDistribute>& c##Src##MapPtr(const label index)    \
-            const                                                              \
-            {                                                                  \
-                checkIndex(index);                                             \
-                return toSource_ ?                                             \
-                    cached##Src##MapPtr_[index]                                \
-                  : cached##Tgt##MapPtr_[index];                               \
-            }
-
-        defineMethodsIndex(Src, Tgt)
-        defineMethodsIndex(Tgt, Src)
-
-
-    // I-O
-
-        //- Write AMI as a dictionary
-        void write(Ostream& os) const;
-
-        //- Write AMI raw
-        bool writeData(Ostream& os) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-    #include "AMICacheTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/meshTools/AMIInterpolation/AMIInterpolation/AMICacheTemplates.C

@@ -1,63 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    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/>.
-
-\*---------------------------------------------------------------------------*/
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type, class EvalFunction>
-bool Foam::AMICache::apply(List<Type>& result, const EvalFunction& eval) const
-{
-    if (applyLower())
-    {
-        eval(result, index0_);
-        return true;
-    }
-    else if (applyUpper())
-    {
-        eval(result, index1_);
-        return true;
-    }
-    else if (applyInterpolate())
-    {
-        List<Type> r0(result);
-        eval(r0, index0_);
-        List<Type> r1(result);
-        eval(r1, index1_);
-
-        //result = (r1 - r0)*interpWeight_ + r0;
-        forAll(result, i)
-        {
-            result[i] = lerp(r0[i], r1[i], interpWeight_);
-        }
-        return true;
-    }
-
-    return false;
-}
-
-
-// ************************************************************************* //

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

@@ -61,26 +61,11 @@ registerOptSwitch
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-void Foam::AMIInterpolation::addToCache(const point& refPt)
-{
-    DebugInfo<< "-- addToCache" << endl;
-
-    cache_.addToCache(*this, refPt);
-}
-
-
-bool Foam::AMIInterpolation::restoreCache(const point& refPt)
-{
-    DebugInfo<< "-- restoreCache" << endl;
-
-    upToDate_ = cache_.restoreCache(refPt);
-
-    return upToDate_;
-}
-
-
 Foam::autoPtr<Foam::indexedOctree<Foam::AMIInterpolation::treeType>>
-Foam::AMIInterpolation::createTree(const primitivePatch &patch) const
+Foam::AMIInterpolation::createTree
+(
+    const primitivePatch& patch
+) const
 {
     treeBoundBox bb(patch.points(), patch.meshPoints());
     bb.inflate(0.01);
@@ -715,8 +700,7 @@ Foam::AMIInterpolation::AMIInterpolation
     tgtWeightsSum_(),
     tgtCentroids_(),
     tgtMapPtr_(nullptr),
-    upToDate_(false),
-    cache_(dict)
+    upToDate_(false)
 {}
 
 
@@ -746,8 +730,7 @@ Foam::AMIInterpolation::AMIInterpolation
     tgtCentroids_(),
     tgtPatchPts_(),
     tgtMapPtr_(nullptr),
-    upToDate_(false),
-    cache_()
+    upToDate_(false)
 {}
 
 
@@ -760,7 +743,7 @@ Foam::AMIInterpolation::AMIInterpolation
 :
     requireMatch_(fineAMI.requireMatch_),
     reverseTarget_(fineAMI.reverseTarget_),
-    lowWeightCorrection_(-1.0), // Deactivated?
+    lowWeightCorrection_(-1.0),
     singlePatchProc_(fineAMI.singlePatchProc_),
     comm_(fineAMI.comm()),  // use fineAMI geomComm if present, comm otherwise
     geomComm_(),
@@ -776,17 +759,16 @@ Foam::AMIInterpolation::AMIInterpolation
     tgtWeightsSum_(),
     tgtPatchPts_(),
     tgtMapPtr_(nullptr),
-    upToDate_(false),
-    cache_(fineAMI.cache(), fineAMI, sourceRestrictAddressing, targetRestrictAddressing)
+    upToDate_(false)
 {
-    const label sourceCoarseSize =
+    label sourceCoarseSize =
     (
         sourceRestrictAddressing.size()
       ? max(sourceRestrictAddressing)+1
       : 0
     );
 
-    const label neighbourCoarseSize =
+    label neighbourCoarseSize =
     (
         targetRestrictAddressing.size()
       ? max(targetRestrictAddressing)+1
@@ -913,15 +895,14 @@ Foam::AMIInterpolation::AMIInterpolation(const AMIInterpolation& ami)
     srcWeights_(ami.srcWeights_),
     srcWeightsSum_(ami.srcWeightsSum_),
     srcCentroids_(ami.srcCentroids_),
-    srcMapPtr_(ami.srcMapPtr_.clone()),
+    srcMapPtr_(nullptr),
     tgtMagSf_(ami.tgtMagSf_),
     tgtAddress_(ami.tgtAddress_),
     tgtWeights_(ami.tgtWeights_),
     tgtWeightsSum_(ami.tgtWeightsSum_),
     tgtCentroids_(ami.tgtCentroids_),
-    tgtMapPtr_(ami.tgtMapPtr_.clone()),
-    upToDate_(ami.upToDate_),
-    cache_(ami.cache_)
+    tgtMapPtr_(nullptr),
+    upToDate_(false)
 {}
 
 
@@ -949,9 +930,7 @@ Foam::AMIInterpolation::AMIInterpolation(Istream& is)
     //tgtPatchPts_(is),
     tgtMapPtr_(nullptr),
 
-    upToDate_(readBool(is)),
-
-    cache_(is)
+    upToDate_(readBool(is))
 {
     // Hopefully no need to stream geomComm_ since only used in processor
     // agglomeration?
@@ -1382,6 +1361,7 @@ const
         }
         else if (ray.distance() < nearest.distance())
         {
+
             nearest = ray;
             nearestFacei = srcFacei;
         }
@@ -1559,8 +1539,6 @@ void Foam::AMIInterpolation::write(Ostream& os) const
     {
         os.writeEntry("lowWeightCorrection", lowWeightCorrection_);
     }
-
-    cache_.write(os);
 }
 
 
@@ -1586,8 +1564,6 @@ bool Foam::AMIInterpolation::writeData(Ostream& os) const
 
         << token::SPACE<< upToDate();
 
-    cache_.writeData(os);
-
     if (distributed() && comm() != -1)
     {
         os  << token::SPACE<< srcMap()

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2025 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,7 +61,7 @@ SourceFiles
 #include "indexedOctree.H"
 #include "treeDataPrimitivePatch.H"
 #include "runTimeSelectionTables.H"
-#include "AMICache.H"
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -78,7 +78,6 @@ public:
 
     // Public Data Types
 
-        //- Flag to store face-to-face intersection triangles; default = false
         static bool cacheIntersections_;
 
         //- Control use of local communicator for AMI communication
@@ -87,10 +86,6 @@ public:
         //  localComm : 2 : like 1 but always include master (for messages)
         static int useLocalComm_;
 
-        //- Tolerance used when caching the AMI to identify e.g. if the
-        //- current rotation angle has already been captured
-        static scalar cacheThetaTolerance_;
-
 
 protected:
 
@@ -150,6 +145,7 @@ protected:
             autoPtr<mapDistribute> srcMapPtr_;
 
 
+
         // Target patch
 
             //- Target face areas
@@ -177,13 +173,9 @@ protected:
             //- Target map pointer - parallel running only
             autoPtr<mapDistribute> tgtMapPtr_;
 
-
         //- Up-to-date flag
         bool upToDate_;
 
-        //- Cache
-        AMICache cache_;
-
 
     // Protected Member Functions
 
@@ -220,10 +212,10 @@ protected:
 
         // Access
 
-            //- Return the original src patch with optionally updated points
+            //- Return the orginal src patch with optionally updated points
             inline const primitivePatch& srcPatch0() const;
 
-            //- Return the original tgt patch with optionally updated points
+            //- Return the orginal tgt patch with optionally updated points
             inline const primitivePatch& tgtPatch0() const;
 
 
@@ -248,6 +240,27 @@ protected:
             );
 
 
+        // Constructor helpers
+
+            static void agglomerate
+            (
+                const autoPtr<mapDistribute>& targetMap,
+                const scalarList& fineSrcMagSf,
+                const labelListList& fineSrcAddress,
+                const scalarListList& fineSrcWeights,
+
+                const labelList& sourceRestrictAddressing,
+                const labelList& targetRestrictAddressing,
+
+                scalarList& srcMagSf,
+                labelListList& srcAddress,
+                scalarListList& srcWeights,
+                scalarField& srcWeightsSum,
+                autoPtr<mapDistribute>& tgtMap,
+                const label comm
+            );
+
+
 public:
 
     //- Runtime type information
@@ -353,26 +366,6 @@ public:
 
     // Member Functions
 
-        // Constructor helpers
-
-            static void agglomerate
-            (
-                const autoPtr<mapDistribute>& targetMap,
-                const scalarList& fineSrcMagSf,
-                const labelListList& fineSrcAddress,
-                const scalarListList& fineSrcWeights,
-
-                const labelList& sourceRestrictAddressing,
-                const labelList& targetRestrictAddressing,
-
-                scalarList& srcMagSf,
-                labelListList& srcAddress,
-                scalarListList& srcWeights,
-                scalarField& srcWeightsSum,
-                autoPtr<mapDistribute>& tgtMap,
-                const label comm
-            );
-
         // Access
 
             //- Is up-to-date?
@@ -386,10 +379,6 @@ public:
                 return old;
             }
 
-            //- Return a reference to the AMI cache
-            const AMICache& cache() const noexcept { return cache_; }
-            AMICache& cache() noexcept { return cache_; }
-
             //- Distributed across processors (singlePatchProc == -1)
             inline bool distributed() const noexcept;
 
@@ -422,8 +411,6 @@ public:
             //  \returns old value
             inline label comm(label communicator) noexcept;
 
-            //- Return true if caching is active
-            inline bool cacheActive() const;
 
             // Source patch
 
@@ -544,8 +531,7 @@ public:
 
             // Low-level
 
-                //- Weighted sum of contributions. Note: cop operates on
-                //- single Type only.
+                //- Weighted sum of contributions
                 template<class Type, class CombineOp>
                 static void weightedSum
                 (
@@ -559,35 +545,20 @@ public:
                     const UList<Type>& defaultValues
                 );
 
-                //- Weighted sum of contributions (includes caching+interp)
-                //- (for primitive types only)
+                //- Weighted sum of contributions
                 template<class Type>
                 void weightedSum
                 (
-                    const bool toSource,
+                    const bool interpolateToSource,
                     const UList<Type>& fld,
                     List<Type>& result,
                     const UList<Type>& defaultValues
                 ) const;
 
-                //- Weighted sum of (potentially distributed) contributions
-                //- and apply caching+interpolation. Note: iop
-                //- operates on single Type only
-                template<class Type, class CombineOp, class InterpolateOp>
-                void interpolate
-                (
-                    const bool toSource,
-                    const UList<Type>& fld,
-                    const CombineOp& cop,
-                    const InterpolateOp& iop,
-                    List<Type>& result,
-                    const UList<Type>& defaultValues
-                ) const;
-
-                //- Interpolate from source to target with supplied op
+                //- Interpolate from target to source with supplied op
                 //- to combine existing value with remote value and weight
                 template<class Type, class CombineOp>
-                void interpolateToTarget
+                void interpolateToSource
                 (
                     const UList<Type>& fld,
                     const CombineOp& cop,
@@ -595,10 +566,10 @@ public:
                     const UList<Type>& defaultValues = UList<Type>::null()
                 ) const;
 
-                //- Interpolate from target to source with supplied op
+                //- Interpolate from source to target with supplied op
                 //- to combine existing value with remote value and weight
                 template<class Type, class CombineOp>
-                void interpolateToSource
+                void interpolateToTarget
                 (
                     const UList<Type>& fld,
                     const CombineOp& cop,
@@ -700,12 +671,6 @@ public:
             )
             const;
 
-            //- Restore AMI weights and addressing from the cache
-            bool restoreCache(const point& refPt);
-
-            //- Add AMI weights and addressing to the cache
-            void addToCache(const point& refPt);
-
 
         // Checks
 

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H

@@ -123,12 +123,6 @@ inline Foam::label Foam::AMIInterpolation::comm(label communicator) noexcept
 }
 
 
-inline bool Foam::AMIInterpolation::cacheActive() const
-{
-    return cache_.active();
-}
-
-
 inline const Foam::List<Foam::scalar>& Foam::AMIInterpolation::srcMagSf() const
 {
     return srcMagSf_;

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationTemplates.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2025 OpenCFD Ltd.
+    Copyright (C) 2015-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -83,19 +83,18 @@ void Foam::AMIInterpolation::weightedSum
 template<class Type>
 void Foam::AMIInterpolation::weightedSum
 (
-    const bool toSource,
+    const bool interpolateToSource,
     const UList<Type>& fld,
     List<Type>& result,
     const UList<Type>& defaultValues
 ) const
 {
-    // Note: using non-caching AMI
     weightedSum
     (
         lowWeightCorrection_,
-        (toSource ? srcAddress_ : tgtAddress_),
-        (toSource ? srcWeights_ : tgtWeights_),
-        (toSource ? srcWeightsSum_ : tgtWeightsSum_),
+        (interpolateToSource ? srcAddress_ : tgtAddress_),
+        (interpolateToSource ? srcWeights_ : tgtWeights_),
+        (interpolateToSource ? srcWeightsSum_ : tgtWeightsSum_),
         fld,
         multiplyWeightedOp<Type, plusEqOp<Type>>(plusEqOp<Type>()),
         result,
@@ -104,223 +103,6 @@ void Foam::AMIInterpolation::weightedSum
 }
 
 
-template<class Type, class CombineOp, class InterpolateOp>
-void Foam::AMIInterpolation::interpolate
-(
-    const bool toSource,
-    const UList<Type>& fld,
-    const CombineOp& cop,
-    const InterpolateOp& iop,
-    List<Type>& result,
-    const UList<Type>& defaultValues
-) const
-{
-    // Note
-    // - behaves as old AMIInterpolation::interpolateToSource if toSource=true
-    // - `result` should be preallocated to correct size and initialized to
-    //   an appropriate value (e.g. Zero)
-
-    // Get data locally and do a weighted sum
-
-    addProfiling(ami, "AMIInterpolation::interpolate");
-
-    cache_.setDirection(toSource);
-
-    auto checkSizes = [&](
-        const UList<Type>& fld,
-        const labelListList& srcAddr,
-        const labelListList& tgtAddr,
-        const UList<Type>& defVals
-    )
-    {
-
-        if (fld.size() != tgtAddr.size())
-        {
-            FatalErrorInFunction
-                << "Supplied field size is not equal to "
-                << (toSource ? "target" : "source") << " patch size" << nl
-                << "    source patch   = " << srcAddr.size() << nl
-                << "    target patch   = " << tgtAddr.size() << nl
-                << "    supplied field = " << fld.size()
-                << abort(FatalError);
-        }
-        else if
-        (
-            (lowWeightCorrection_ > 0) && (defVals.size() != srcAddr.size())
-        )
-        {
-            FatalErrorInFunction
-                << "Employing default values when sum of weights falls below "
-                << lowWeightCorrection_
-                << " but number of default values is not equal to "
-                << (toSource ? "source" : "target") << " patch size" << nl
-                << "    default values = " << defVals.size() << nl
-                << "    source patch   = " << srcAddr.size() << nl
-                << abort(FatalError);
-        }
-    };
-
-
-    // Work space for if distributed
-    List<Type> work;
-
-    List<Type> result0;
-    if (cache_.index0() != -1)
-    {
-        result0 = result;
-
-        const auto& srcAddress = cache_.cSrcAddress0();
-        const auto& srcWeights = cache_.cSrcWeights0();
-        const auto& srcWeightsSum = cache_.cSrcWeightsSum0();
-        const auto& tgtAddress = cache_.cTgtAddress0();
-
-        checkSizes(fld, srcAddress, tgtAddress, defaultValues);
-
-        if (distributed() && cache_.cTgtMapPtr0())
-        {
-            const mapDistribute& map = cache_.cTgtMapPtr0()();
-
-            if (map.comm() == -1)
-            {
-                return;
-            }
-
-            work.resize_nocopy(map.constructSize());
-            SubList<Type>(work, fld.size()) = fld;  // deep copy
-            map.distribute(work);
-        }
-
-        // if constexpr (is_contiguous_scalar<Type>::value)
-        // {
-        //     result0 = Zero;
-        // }
-
-        weightedSum
-        (
-            lowWeightCorrection_,
-            srcAddress,
-            srcWeights,
-            srcWeightsSum,
-            (distributed() ? work : fld),
-            cop,
-            result0,
-            defaultValues
-        );
-    }
-
-    List<Type> result1;
-    if (cache_.index1() != -1)
-    {
-        result1 = result;
-
-        const auto& srcAddress = cache_.cSrcAddress1();
-        const auto& srcWeights = cache_.cSrcWeights1();
-        const auto& srcWeightsSum = cache_.cSrcWeightsSum1();
-        const auto& tgtAddress = cache_.cTgtAddress1();
-
-        checkSizes(fld, srcAddress, tgtAddress, defaultValues);
-
-        if (distributed() && cache_.cTgtMapPtr1())
-        {
-            const mapDistribute& map = cache_.cTgtMapPtr1()();
-
-            if (map.comm() == -1)
-            {
-                return;
-            }
-
-            work.resize_nocopy(map.constructSize());
-            SubList<Type>(work, fld.size()) = fld;  // deep copy
-            map.distribute(work);
-        }
-
-        // if constexpr (is_contiguous_scalar<Type>::value)
-        // {
-        //     result1 = Zero;
-        // }
-
-        weightedSum
-        (
-            lowWeightCorrection_,
-            srcAddress,
-            srcWeights,
-            srcWeightsSum,
-            (distributed() ? work : fld),
-            cop,
-            result1,
-            defaultValues
-        );
-    }
-
-    if (cache_.applyLower())
-    {
-        result = result0;
-    }
-    else if (cache_.applyUpper())
-    {
-        result = result1;
-    }
-    else if (cache_.applyInterpolate())
-    {
-        forAll(result, i)
-        {
-            iop(result[i], i, i, result0[i], i, result1[i], cache_.weight());
-        }
-    }
-    else
-    {
-        // No cache - evaluate the AMI
-
-        const auto& srcAddress = (toSource ? srcAddress_ : tgtAddress_);
-        const auto& srcWeights = (toSource ? srcWeights_ : tgtWeights_);
-        const auto& srcWeightsSum =
-            (toSource ? srcWeightsSum_ : tgtWeightsSum_);
-        const auto& tgtAddress = (toSource ? tgtAddress_ : srcAddress_);
-
-        checkSizes(fld, srcAddress, tgtAddress, defaultValues);
-
-        if (distributed() && tgtMapPtr_)
-        {
-            const mapDistribute& map =
-            (
-                toSource
-              ? tgtMapPtr_()
-              : srcMapPtr_()
-            );
-
-            if (map.comm() == -1)
-            {
-                return;
-            }
-
-            work.resize_nocopy(map.constructSize());
-            SubList<Type>(work, fld.size()) = fld;  // deep copy
-            map.distribute(work);
-        }
-
-        // result.resize_nocopy(srcAddress.size());
-
-        // if constexpr (is_contiguous_scalar<Type>::value)
-        // {
-        //     result = Zero;
-        // }
-
-        weightedSum
-        (
-            lowWeightCorrection_,
-            srcAddress,
-            srcWeights,
-            srcWeightsSum,
-            (distributed() ? work : fld),
-            cop,
-            result,
-            defaultValues
-        );
-    }
-}
-
-
-// Leave API intact below!
 template<class Type, class CombineOp>
 void Foam::AMIInterpolation::interpolateToTarget
 (
@@ -330,31 +112,58 @@ void Foam::AMIInterpolation::interpolateToTarget
     const UList<Type>& defaultValues
 ) const
 {
-    // In-place interpolation
-
     addProfiling(ami, "AMIInterpolation::interpolateToTarget");
 
-    // Wrap lerp operator to operate inplace
-    auto iop = [&]
+    if (fld.size() != srcAddress_.size())
+    {
+        FatalErrorInFunction
+            << "Supplied field size is not equal to source patch size" << nl
+            << "    source patch   = " << srcAddress_.size() << nl
+            << "    target patch   = " << tgtAddress_.size() << nl
+            << "    supplied field = " << fld.size()
+            << abort(FatalError);
+    }
+    else if
     (
-        Type& res,
-        const label i,
-        const label ia,
-        const Type& a,
-        const label ib,
-        const Type& b,
-        const scalar w
+        (lowWeightCorrection_ > 0)
+     && (defaultValues.size() != tgtAddress_.size())
     )
     {
-        res = lerp(a, b, w);
-    };
+        FatalErrorInFunction
+            << "Employing default values when sum of weights falls below "
+            << lowWeightCorrection_
+            << " but supplied default field size is not equal to target "
+            << "patch size" << nl
+            << "    default values = " << defaultValues.size() << nl
+            << "    target patch   = " << tgtAddress_.size() << nl
+            << abort(FatalError);
+    }
 
-    interpolate
+    result.setSize(tgtAddress_.size());
+    List<Type> work;
+
+    if (distributed() && srcMapPtr_)
+    {
+        const mapDistribute& map = srcMapPtr_();
+
+        if (map.comm() == -1)
+        {
+            return;
+        }
+
+        work.resize_nocopy(map.constructSize());
+        SubList<Type>(work, fld.size()) = fld;  // deep copy
+        map.distribute(work);
+    }
+
+    weightedSum
     (
-        false,                  // interpolate to target
-        fld,
+        lowWeightCorrection_,
+        tgtAddress_,
+        tgtWeights_,
+        tgtWeightsSum_,
+        (distributed() ? work : fld),
         cop,
-        iop,
         result,
         defaultValues
     );
@@ -370,32 +179,58 @@ void Foam::AMIInterpolation::interpolateToSource
     const UList<Type>& defaultValues
 ) const
 {
-    // In-place interpolation
-
     addProfiling(ami, "AMIInterpolation::interpolateToSource");
 
-    // Wrap lerp operator to operate inplace
-    auto iop = [&]
+    if (fld.size() != tgtAddress_.size())
+    {
+        FatalErrorInFunction
+            << "Supplied field size is not equal to target patch size" << nl
+            << "    source patch   = " << srcAddress_.size() << nl
+            << "    target patch   = " << tgtAddress_.size() << nl
+            << "    supplied field = " << fld.size()
+            << abort(FatalError);
+    }
+    else if
     (
-        Type& res,
-        const label i,
-        const label ia,
-        const Type& a,
-        const label ib,
-        const Type& b,
-        const scalar w
+        (lowWeightCorrection_ > 0)
+     && (defaultValues.size() != srcAddress_.size())
     )
     {
-        res = lerp(a, b, w);
-    };
+        FatalErrorInFunction
+            << "Employing default values when sum of weights falls below "
+            << lowWeightCorrection_
+            << " but number of default values is not equal to source "
+            << "patch size" << nl
+            << "    default values = " << defaultValues.size() << nl
+            << "    source patch   = " << srcAddress_.size() << nl
+            << abort(FatalError);
+    }
 
+    result.setSize(srcAddress_.size());
+    List<Type> work;
 
-    interpolate
+    if (distributed() && tgtMapPtr_)
+    {
+        const mapDistribute& map = tgtMapPtr_();
+
+        if (map.comm() == -1)
+        {
+            return;
+        }
+
+        work.resize_nocopy(map.constructSize());
+        SubList<Type>(work, fld.size()) = fld;  // deep copy
+        map.distribute(work);
+    }
+
+    weightedSum
     (
-        true,                   // toSource,
-        fld,
+        lowWeightCorrection_,
+        srcAddress_,
+        srcWeights_,
+        srcWeightsSum_,
+        (distributed() ? work : fld),
         cop,
-        iop,
         result,
         defaultValues
     );

src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicACMIGAMGInterfaceField/cyclicACMIGAMGInterfaceField.C

@@ -67,7 +67,9 @@ Foam::cyclicACMIGAMGInterfaceField::cyclicACMIGAMGInterfaceField
     GAMGInterfaceField(GAMGCp, fineInterface),
     cyclicACMIInterface_(refCast<const cyclicACMIGAMGInterface>(GAMGCp)),
     doTransform_(false),
-    rank_(0)
+    rank_(0),
+    sendRequests_(),
+    recvRequests_()
 {
     const cyclicAMILduInterfaceField& p =
         refCast<const cyclicAMILduInterfaceField>(fineInterface);
@@ -87,7 +89,9 @@ Foam::cyclicACMIGAMGInterfaceField::cyclicACMIGAMGInterfaceField
     GAMGInterfaceField(GAMGCp, doTransform, rank),
     cyclicACMIInterface_(refCast<const cyclicACMIGAMGInterface>(GAMGCp)),
     doTransform_(doTransform),
-    rank_(rank)
+    rank_(rank),
+    sendRequests_(),
+    recvRequests_()
 {}
 
 
@@ -100,7 +104,9 @@ Foam::cyclicACMIGAMGInterfaceField::cyclicACMIGAMGInterfaceField
     GAMGInterfaceField(GAMGCp, is),
     cyclicACMIInterface_(refCast<const cyclicACMIGAMGInterface>(GAMGCp)),
     doTransform_(readBool(is)),
-    rank_(readLabel(is))
+    rank_(readLabel(is)),
+    sendRequests_(),
+    recvRequests_()
 {}
 
 
@@ -114,7 +120,9 @@ Foam::cyclicACMIGAMGInterfaceField::cyclicACMIGAMGInterfaceField
     GAMGInterfaceField(GAMGCp, local),
     cyclicACMIInterface_(refCast<const cyclicACMIGAMGInterface>(GAMGCp)),
     doTransform_(false),
-    rank_(0)
+    rank_(0),
+    sendRequests_(),
+    recvRequests_()
 {
     const auto& p = refCast<const cyclicACMILduInterfaceField>(local);
 
@@ -134,15 +142,9 @@ bool Foam::cyclicACMIGAMGInterfaceField::ready() const
             recvRequests_.start(),
             recvRequests_.size()
         )
-     && UPstream::finishedRequests
-        (
-            recvRequests1_.start(),
-            recvRequests1_.size()
-        )
     )
     {
         recvRequests_.clear();
-        recvRequests1_.clear();
 
         if
         (
@@ -151,15 +153,9 @@ bool Foam::cyclicACMIGAMGInterfaceField::ready() const
                 sendRequests_.start(),
                 sendRequests_.size()
             )
-         && UPstream::finishedRequests
-            (
-                sendRequests1_.start(),
-                sendRequests1_.size()
-            )
         )
         {
             sendRequests_.clear();
-            sendRequests1_.clear();
         }
 
         return true;
@@ -214,9 +210,15 @@ void Foam::cyclicACMIGAMGInterfaceField::initInterfaceMatrixUpdate
         // Transform according to the transformation tensors
         transformCoupleField(pnf, cmpt);
 
+        const auto& map =
+        (
+            cyclicACMIInterface_.owner()
+          ? AMI.tgtMap()
+          : AMI.srcMap()
+        );
 
         // Assert that all receives are known to have finished
-        if (!recvRequests_.empty() || !recvRequests1_.empty())
+        if (!recvRequests_.empty())
         {
             FatalErrorInFunction
                 << "Outstanding recv request(s) on patch "
@@ -224,63 +226,22 @@ void Foam::cyclicACMIGAMGInterfaceField::initInterfaceMatrixUpdate
                 << abort(FatalError);
         }
 
-        const auto& cache = AMI.cache();
+        // Assume that sends are also OK
+        sendRequests_.clear();
 
-        if (cache.index0() == -1 && cache.index1() == -1)
-        {
-            const auto& map =
-            (
-                cyclicACMIInterface_.owner()
-              ? AMI.tgtMap()
-              : AMI.srcMap()
-            );
-
-            // Insert send/receive requests (non-blocking). See e.g.
-            // cyclicAMIPolyPatchTemplates.C
-            const label oldWarnComm = UPstream::commWarn(AMI.comm());
-            map.send
-            (
-                pnf,
-                sendRequests_,
-                scalarSendBufs_,
-                recvRequests_,
-                scalarRecvBufs_,
-                19462+cyclicACMIInterface_.index()  // unique offset + patch index
-            );
-            UPstream::commWarn(oldWarnComm);
-        }
-        else
-        {
-            cache.setDirection(cyclicACMIInterface_.owner());
-
-            if (cache.index0() != -1)
-            {
-                const auto& map0 = cache.cTgtMapPtr0()();
-                map0.send
-                (
-                    pnf,
-                    sendRequests_,
-                    scalarSendBufs_,
-                    recvRequests_,
-                    scalarRecvBufs_,
-                    19462+cyclicACMIInterface_.index()   // unique offset + patch index
-                );
-            }
-
-            if (cache.index1() != -1)
-            {
-                const auto& map1 = cache.cTgtMapPtr1()();
-                map1.send
-                (
-                    pnf,
-                    sendRequests1_,
-                    scalarSendBufs1_,
-                    recvRequests1_,
-                    scalarRecvBufs1_,
-                    19463+cyclicACMIInterface_.index()   // unique offset + patch index
-                );
-            }
-        }
+        // Insert send/receive requests (non-blocking). See e.g.
+        // cyclicAMIPolyPatchTemplates.C
+        const label oldWarnComm = UPstream::commWarn(AMI.comm());
+        map.send
+        (
+            pnf,
+            sendRequests_,
+            scalarSendBufs_,
+            recvRequests_,
+            scalarRecvBufs_,
+            19462+cyclicACMIInterface_.index()  // unique offset + patch index
+        );
+        UPstream::commWarn(oldWarnComm);
     }
 
     this->updatedMatrix(false);
@@ -299,8 +260,6 @@ void Foam::cyclicACMIGAMGInterfaceField::updateInterfaceMatrix
     const Pstream::commsTypes
 ) const
 {
-    typedef multiplyWeightedOp<scalar, plusEqOp<scalar>> opType;
-
     const labelUList& faceCells = lduAddr.patchAddr(patchId);
 
     const auto& AMI =
@@ -310,122 +269,48 @@ void Foam::cyclicACMIGAMGInterfaceField::updateInterfaceMatrix
       : cyclicACMIInterface_.neighbPatch().AMI()
     );
 
-    const auto& cache = AMI.cache();
+    DebugPout<< "cyclicACMIGAMGInterfaceField::updateInterfaceMatrix() :"
+        << " interface:" << cyclicACMIInterface_.index()
+        << " size:" << cyclicACMIInterface_.size()
+        << " owner:" << cyclicACMIInterface_.owner()
+        << " AMI distributed:" << AMI.distributed()
+        << endl;
+
 
     if (AMI.distributed() && AMI.comm() != -1)
     {
-        if (cache.index0() == -1 && cache.index1() == -1)
-        {
-            const auto& map =
-            (
-                cyclicACMIInterface_.owner()
-              ? AMI.tgtMap()
-              : AMI.srcMap()
-            );
+        const auto& map =
+        (
+            cyclicACMIInterface_.owner()
+          ? AMI.tgtMap()
+          : AMI.srcMap()
+        );
 
-            // Receive (= copy) data from buffers into work. TBD: receive directly
-            // into slices of work.
-            solveScalarField work;
-            map.receive
-            (
-                recvRequests_,
-                scalarRecvBufs_,
-                work,
-                19462+cyclicACMIInterface_.index()  // unique offset + patch index
-            );
+        // Receive (= copy) data from buffers into work. TBD: receive directly
+        // into slices of work.
+        solveScalarField work;
+        map.receive
+        (
+            recvRequests_,
+            scalarRecvBufs_,
+            work,
+            19462+cyclicACMIInterface_.index()  // unique offset + patch index
+        );
 
-            // Receive requests all handled by last function call
-            recvRequests_.clear();
+        // Receive requests all handled by last function call
+        recvRequests_.clear();
 
-            solveScalarField pnf(faceCells.size(), Zero);
-            AMI.weightedSum
-            (
-                cyclicACMIInterface_.owner(),
-                work,
-                pnf,               // result
-                solveScalarField::null()
-            );
+        solveScalarField pnf(faceCells.size(), Zero);
+        AMI.weightedSum
+        (
+            cyclicACMIInterface_.owner(),
+            work,
+            pnf,               // result
+            solveScalarField::null()
+        );
 
-            // Add result using coefficients
-            this->addToInternalField(result, !add, faceCells, coeffs, pnf);
-        }
-        else
-        {
-            cache.setDirection(cyclicACMIInterface_.owner());
-
-            solveScalarField pnf(faceCells.size());
-            solveScalarField work(faceCells.size());
-
-            if (cache.index0() != -1)
-            {
-                // Receive (= copy) data from buffers into work. TBD: receive directly
-                // into slices of work.
-                work = Zero;
-                cache.cTgtMapPtr0()().receive
-                (
-                    recvRequests_,
-                    scalarRecvBufs_,
-                    work,
-                    19462+cyclicACMIInterface_.index()   // unique offset + patch index
-                );
-
-                // Receive requests all handled by last function call
-                recvRequests_.clear();
-
-                pnf = Zero;
-                AMIInterpolation::weightedSum
-                (
-                    AMI.lowWeightCorrection(),
-                    cache.cSrcAddress0(),
-                    cache.cSrcWeights0(),
-                    cache.cSrcWeightsSum0(),
-                    work,
-                    opType(plusEqOp<scalar>()),
-                    pnf,
-                    solveScalarField::null()
-                );
-
-                pnf *= (1-cache.weight());
-
-                // Add result using coefficients
-                this->addToInternalField(result, !add, faceCells, coeffs, pnf);
-            }
-
-            if (cache.index1() != -1)
-            {
-                // Receive (= copy) data from buffers into work. TBD: receive directly
-                // into slices of work.
-                work = Zero;
-                cache.cTgtMapPtr1()().receive
-                (
-                    recvRequests1_,
-                    scalarRecvBufs1_,
-                    work,
-                    19463+cyclicACMIInterface_.index()   // unique offset + patch index
-                );
-
-                // Receive requests all handled by last function call
-                recvRequests1_.clear();
-
-                pnf = Zero;
-                AMIInterpolation::weightedSum
-                (
-                    AMI.lowWeightCorrection(),
-                    cache.cSrcAddress1(),
-                    cache.cSrcWeights1(),
-                    cache.cSrcWeightsSum1(),
-                    work,
-                    opType(plusEqOp<scalar>()),
-                    pnf,
-                    solveScalarField::null()
-                );
-
-                pnf *= cache.weight();
-
-                // Add result using coefficients
-                this->addToInternalField(result, !add, faceCells, coeffs, pnf);
-            }
-        }
+        // Add result using coefficients
+        this->addToInternalField(result, !add, faceCells, coeffs, pnf);
     }
     else
     {
@@ -433,73 +318,23 @@ void Foam::cyclicACMIGAMGInterfaceField::updateInterfaceMatrix
         const labelUList& nbrFaceCells =
             lduAddr.patchAddr(cyclicACMIInterface_.neighbPatchID());
 
-        solveScalarField work(psiInternal, nbrFaceCells);
+        solveScalarField pnf(psiInternal, nbrFaceCells);
 
         // Transform according to the transformation tensors
-        transformCoupleField(work, cmpt);
+        transformCoupleField(pnf, cmpt);
 
-        if (cache.index0() == -1 && cache.index1() == -1)
+        if (cyclicACMIInterface_.owner())
         {
-            solveScalarField pnf(faceCells.size(), Zero);
-
-            AMI.weightedSum
-            (
-                cyclicACMIInterface_.owner(),
-                work,
-                pnf,                // result
-                solveScalarField::null()
-            );
-
-            const labelUList& faceCells = lduAddr.patchAddr(patchId);
-
-            this->addToInternalField(result, !add, faceCells, coeffs, pnf);
+            pnf = AMI.interpolateToSource(pnf);
         }
         else
         {
-            cache.setDirection(cyclicACMIInterface_.owner());
-
-            solveScalarField pnf(faceCells.size());
-
-            if (cache.index0() != -1)
-            {
-                pnf = Zero;
-                AMIInterpolation::weightedSum
-                (
-                    AMI.lowWeightCorrection(),
-                    cache.cSrcAddress0(),
-                    cache.cSrcWeights0(),
-                    cache.cSrcWeightsSum0(),
-                    work,
-                    opType(plusEqOp<scalar>()),
-                    pnf,
-                    solveScalarField::null()
-                );
-
-                pnf *= (1 - cache.weight());
-
-                this->addToInternalField(result, !add, faceCells, coeffs, pnf);
-            }
-
-            if (cache.index1() != -1)
-            {
-                pnf = Zero;
-                AMIInterpolation::weightedSum
-                (
-                    AMI.lowWeightCorrection(),
-                    cache.cSrcAddress1(),
-                    cache.cSrcWeights1(),
-                    cache.cSrcWeightsSum1(),
-                    work,
-                    opType(plusEqOp<scalar>()),
-                    pnf,
-                    solveScalarField::null()
-                );
-
-                pnf *= cache.weight();
-
-                this->addToInternalField(result, !add, faceCells, coeffs, pnf);
-            }
+            pnf = AMI.interpolateToTarget(pnf);
         }
+
+        const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
+        this->addToInternalField(result, !add, faceCells, coeffs, pnf);
     }
 
     this->updatedMatrix(true);

src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicACMIGAMGInterfaceField/cyclicACMIGAMGInterfaceField.H

@@ -83,20 +83,6 @@ class cyclicACMIGAMGInterfaceField
             //- Scalar receive buffers
             mutable PtrList<List<solveScalar>> scalarRecvBufs_;
 
-            // Only used for AMI caching
-
-                //- Current range of send requests (non-blocking)
-                mutable labelRange sendRequests1_;
-
-                //- Current range of recv requests (non-blocking)
-                mutable labelRange recvRequests1_;
-
-                //- Scalar send buffers
-                mutable PtrList<List<solveScalar>> scalarSendBufs1_;
-
-                //- Scalar receive buffers
-                mutable PtrList<List<solveScalar>> scalarRecvBufs1_;
-
 
     // Private Member Functions
 

src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicAMIGAMGInterfaceField/cyclicAMIGAMGInterfaceField.C

@@ -67,7 +67,9 @@ Foam::cyclicAMIGAMGInterfaceField::cyclicAMIGAMGInterfaceField
     GAMGInterfaceField(GAMGCp, fineInterface),
     cyclicAMIInterface_(refCast<const cyclicAMIGAMGInterface>(GAMGCp)),
     doTransform_(false),
-    rank_(0)
+    rank_(0),
+    sendRequests_(),
+    recvRequests_()
 {
     const cyclicAMILduInterfaceField& p =
         refCast<const cyclicAMILduInterfaceField>(fineInterface);
@@ -87,7 +89,9 @@ Foam::cyclicAMIGAMGInterfaceField::cyclicAMIGAMGInterfaceField
     GAMGInterfaceField(GAMGCp, doTransform, rank),
     cyclicAMIInterface_(refCast<const cyclicAMIGAMGInterface>(GAMGCp)),
     doTransform_(doTransform),
-    rank_(rank)
+    rank_(rank),
+    sendRequests_(),
+    recvRequests_()
 {}
 
 
@@ -100,7 +104,9 @@ Foam::cyclicAMIGAMGInterfaceField::cyclicAMIGAMGInterfaceField
     GAMGInterfaceField(GAMGCp, is),
     cyclicAMIInterface_(refCast<const cyclicAMIGAMGInterface>(GAMGCp)),
     doTransform_(readBool(is)),
-    rank_(readLabel(is))
+    rank_(readLabel(is)),
+    sendRequests_(),
+    recvRequests_()
 {}
 
 
@@ -114,7 +120,9 @@ Foam::cyclicAMIGAMGInterfaceField::cyclicAMIGAMGInterfaceField
     GAMGInterfaceField(GAMGCp, local),
     cyclicAMIInterface_(refCast<const cyclicAMIGAMGInterface>(GAMGCp)),
     doTransform_(false),
-    rank_(0)
+    rank_(0),
+    sendRequests_(),    // assume no requests in flight for input field
+    recvRequests_()
 {
     const auto& p = refCast<const cyclicAMILduInterfaceField>(local);
 
@@ -134,15 +142,9 @@ bool Foam::cyclicAMIGAMGInterfaceField::ready() const
             recvRequests_.start(),
             recvRequests_.size()
         )
-     && UPstream::finishedRequests
-        (
-            recvRequests1_.start(),
-            recvRequests1_.size()
-        )
     )
     {
         recvRequests_.clear();
-        recvRequests1_.clear();
 
         if
         (
@@ -151,15 +153,9 @@ bool Foam::cyclicAMIGAMGInterfaceField::ready() const
                 sendRequests_.start(),
                 sendRequests_.size()
             )
-         && UPstream::finishedRequests
-            (
-                sendRequests1_.start(),
-                sendRequests1_.size()
-            )
         )
         {
             sendRequests_.clear();
-            sendRequests1_.clear();
         }
 
         return true;
@@ -187,6 +183,7 @@ void Foam::cyclicAMIGAMGInterfaceField::initInterfaceMatrixUpdate
       ? cyclicAMIInterface_.AMI()
       : cyclicAMIInterface_.neighbPatch().AMI()
     );
+
     if (AMI.distributed() && AMI.comm() != -1)
     {
         //DebugPout<< "cyclicAMIFvPatchField::initInterfaceMatrixUpdate() :"
@@ -214,8 +211,15 @@ void Foam::cyclicAMIGAMGInterfaceField::initInterfaceMatrixUpdate
         // Transform according to the transformation tensors
         transformCoupleField(pnf, cmpt);
 
+        const auto& map =
+        (
+            cyclicAMIInterface_.owner()
+          ? AMI.tgtMap()
+          : AMI.srcMap()
+        );
+
         // Assert that all receives are known to have finished
-        if (!recvRequests_.empty() || !recvRequests1_.empty())
+        if (!recvRequests_.empty())
         {
             FatalErrorInFunction
                 << "Outstanding recv request(s) on patch "
@@ -223,63 +227,22 @@ void Foam::cyclicAMIGAMGInterfaceField::initInterfaceMatrixUpdate
                 << abort(FatalError);
         }
 
-        const auto& cache = AMI.cache();
+        // Assume that sends are also OK
+        sendRequests_.clear();
 
-        if (cache.index0() == -1 && cache.index1() == -1)
-        {
-            const auto& map =
-            (
-                cyclicAMIInterface_.owner()
-              ? AMI.tgtMap()
-              : AMI.srcMap()
-            );
-
-            // Insert send/receive requests (non-blocking). See e.g.
-            // cyclicAMIPolyPatchTemplates.C
-            const label oldWarnComm = UPstream::commWarn(AMI.comm());
-            map.send
-            (
-                pnf,
-                sendRequests_,
-                scalarSendBufs_,
-                recvRequests_,
-                scalarRecvBufs_,
-                19462+cyclicAMIInterface_.index()   // unique offset + patch index
-            );
-            UPstream::commWarn(oldWarnComm);
-        }
-        else
-        {
-            cache.setDirection(cyclicAMIInterface_.owner());
-
-            if (cache.index0() != -1)
-            {
-                const auto& map0 = cache.cTgtMapPtr0()();
-                map0.send
-                (
-                    pnf,
-                    sendRequests_,
-                    scalarSendBufs_,
-                    recvRequests_,
-                    scalarRecvBufs_,
-                    19462+cyclicAMIInterface_.index()   // unique offset + patch index
-                );
-            }
-
-            if (cache.index1() != -1)
-            {
-                const auto& map1 = cache.cTgtMapPtr1()();
-                map1.send
-                (
-                    pnf,
-                    sendRequests1_,
-                    scalarSendBufs1_,
-                    recvRequests1_,
-                    scalarRecvBufs1_,
-                    19463+cyclicAMIInterface_.index()   // unique offset + patch index
-                );
-            }
-        }
+        // Insert send/receive requests (non-blocking). See e.g.
+        // cyclicAMIPolyPatchTemplates.C
+        const label oldWarnComm = UPstream::commWarn(AMI.comm());
+        map.send
+        (
+            pnf,
+            sendRequests_,
+            scalarSendBufs_,
+            recvRequests_,
+            scalarRecvBufs_,
+            19462+cyclicAMIInterface_.index()   // unique offset + patch index
+        );
+        UPstream::commWarn(oldWarnComm);
     }
 
     this->updatedMatrix(false);
@@ -298,8 +261,6 @@ void Foam::cyclicAMIGAMGInterfaceField::updateInterfaceMatrix
     const Pstream::commsTypes commsType
 ) const
 {
-    typedef multiplyWeightedOp<scalar, plusEqOp<scalar>> opType;
-
     const labelUList& faceCells = lduAddr.patchAddr(patchId);
 
     const auto& AMI =
@@ -323,12 +284,6 @@ void Foam::cyclicAMIGAMGInterfaceField::updateInterfaceMatrix
     //    << " AMI low-weight:" << AMI.applyLowWeightCorrection()
     //    << endl;
 
-    const auto& cache = AMI.cache();
-
-    // Assume that sends are also OK
-    sendRequests_.clear();
-    sendRequests1_.clear();
-
     if (AMI.distributed() && AMI.comm() != -1)
     {
         if (commsType != UPstream::commsTypes::nonBlocking)
@@ -338,118 +293,38 @@ void Foam::cyclicAMIGAMGInterfaceField::updateInterfaceMatrix
                 << exit(FatalError);
         }
 
-        if (cache.index0() == -1 && cache.index1() == -1)
-        {
-            const auto& map =
-            (
-                cyclicAMIInterface_.owner()
-              ? AMI.tgtMap()
-              : AMI.srcMap()
-            );
+        const auto& map =
+        (
+            cyclicAMIInterface_.owner()
+          ? AMI.tgtMap()
+          : AMI.srcMap()
+        );
 
-            // Receive (= copy) data from buffers into work. TBD: receive directly
-            // into slices of work.
-            solveScalarField work;
-            map.receive
-            (
-                recvRequests_,
-                scalarRecvBufs_,
-                work,
-                19462+cyclicAMIInterface_.index()   // unique offset + patch index
-            );
+        // Receive (= copy) data from buffers into work. TBD: receive directly
+        // into slices of work.
+        solveScalarField work;
+        map.receive
+        (
+            recvRequests_,
+            scalarRecvBufs_,
+            work,
+            19462+cyclicAMIInterface_.index()   // unique offset + patch index
+        );
 
-            // Receive requests all handled by last function call
-            recvRequests_.clear();
+        // Receive requests all handled by last function call
+        recvRequests_.clear();
 
-            solveScalarField pnf(faceCells.size(), Zero);
-            AMI.weightedSum
-            (
-                cyclicAMIInterface_.owner(),
-                work,
-                pnf,                // result
-                defaultValues
-            );
+        solveScalarField pnf(faceCells.size(), Zero);
+        AMI.weightedSum
+        (
+            cyclicAMIInterface_.owner(),
+            work,
+            pnf,                // result
+            defaultValues
+        );
 
-            // Add result using coefficients
-            this->addToInternalField(result, !add, faceCells, coeffs, pnf);
-        }
-        else
-        {
-            cache.setDirection(cyclicAMIInterface_.owner());
-
-            solveScalarField pnf(faceCells.size());
-            solveScalarField work(faceCells.size());
-
-            if (cache.index0() != -1)
-            {
-                // Receive (= copy) data from buffers into work. TBD: receive directly
-                // into slices of work.
-                work = Zero;
-                cache.cTgtMapPtr0()().receive
-                (
-                    recvRequests_,
-                    scalarRecvBufs_,
-                    work,
-                    19462+cyclicAMIInterface_.index()   // unique offset + patch index
-                );
-
-                // Receive requests all handled by last function call
-                recvRequests_.clear();
-
-                pnf = Zero;
-                AMIInterpolation::weightedSum
-                (
-                    AMI.lowWeightCorrection(),
-                    cache.cSrcAddress0(),
-                    cache.cSrcWeights0(),
-                    cache.cSrcWeightsSum0(),
-                    work,
-                    opType(plusEqOp<scalar>()),
-                    pnf,
-                    defaultValues
-                );
-
-                pnf *= (1-cache.weight());
-
-                // Add result using coefficients
-                this->addToInternalField(result, !add, faceCells, coeffs, pnf);
-            }
-
-            if (cache.index1() != -1)
-            {
-                // Receive (= copy) data from buffers into work. TBD: receive directly
-                // into slices of work.
-                work = Zero;
-                cache.cTgtMapPtr1()().receive
-                (
-                    recvRequests1_,
-                    scalarRecvBufs1_,
-                    work,
-                    19463+cyclicAMIInterface_.index()   // unique offset + patch index
-                );
-
-                // Receive requests all handled by last function call
-                recvRequests1_.clear();
-
-                pnf = Zero;
-                AMIInterpolation::weightedSum
-                (
-                    AMI.lowWeightCorrection(),
-                    cache.cSrcAddress1(),
-                    cache.cSrcWeights1(),
-                    cache.cSrcWeightsSum1(),
-                    work,
-                    opType(plusEqOp<scalar>()),
-                    pnf,
-                    defaultValues
-                );
-
-                pnf *= cache.weight();
-
-                // Add result using coefficients
-                this->addToInternalField(result, !add, faceCells, coeffs, pnf);
-            }
-        }
+        // Add result using coefficients
+        this->addToInternalField(result, !add, faceCells, coeffs, pnf);
     }
     else
     {
@@ -462,68 +337,17 @@ void Foam::cyclicAMIGAMGInterfaceField::updateInterfaceMatrix
         // Transform according to the transformation tensors
         transformCoupleField(work, cmpt);
 
-        solveScalarField pnf(faceCells.size());
+        solveScalarField pnf(faceCells.size(), Zero);
+        AMI.weightedSum
+        (
+            cyclicAMIInterface_.owner(),
+            work,
+            pnf,                // result
+            defaultValues
+        );
 
-        if (cache.index0() == -1 && cache.index1() == -1)
-        {
-            pnf = Zero;
-            AMI.weightedSum
-            (
-                cyclicAMIInterface_.owner(),
-                work,
-                pnf,                // result
-                defaultValues
-            );
-
-            // Add result using coefficients
-            this->addToInternalField(result, !add, faceCells, coeffs, pnf);
-        }
-        else
-        {
-            cache.setDirection(cyclicAMIInterface_.owner());
-
-            if (cache.index0() != -1)
-            {
-                pnf = Zero;
-                AMIInterpolation::weightedSum
-                (
-                    AMI.lowWeightCorrection(),
-                    cache.cSrcAddress0(),
-                    cache.cSrcWeights0(),
-                    cache.cSrcWeightsSum0(),
-                    work,
-                    opType(plusEqOp<scalar>()),
-                    pnf,
-                    defaultValues
-                );
-
-                pnf *= (1 - cache.weight());
-
-                // Add result using coefficients
-                this->addToInternalField(result, !add, faceCells, coeffs, pnf);
-            }
-
-            if (cache.index1() != -1)
-            {
-                pnf = Zero;
-                AMIInterpolation::weightedSum
-                (
-                    AMI.lowWeightCorrection(),
-                    cache.cSrcAddress1(),
-                    cache.cSrcWeights1(),
-                    cache.cSrcWeightsSum1(),
-                    work,
-                    opType(plusEqOp<scalar>()),
-                    pnf,
-                    defaultValues
-                );
-
-                pnf *= cache.weight();
-
-                // Add result using coefficients
-                this->addToInternalField(result, !add, faceCells, coeffs, pnf);
-            }
-        }
+        // Add result using coefficients
+        this->addToInternalField(result, !add, faceCells, coeffs, pnf);
     }
 
     this->updatedMatrix(true);

src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicAMIGAMGInterfaceField/cyclicAMIGAMGInterfaceField.H

@@ -83,21 +83,6 @@ class cyclicAMIGAMGInterfaceField
             mutable PtrList<List<solveScalar>> scalarRecvBufs_;
 
 
-            // Only used for AMI caching
-
-                //- Current range of send requests (non-blocking)
-                mutable labelRange sendRequests1_;
-
-                //- Current range of recv requests (non-blocking)
-                mutable labelRange recvRequests1_;
-
-                //- Scalar send buffers
-                mutable PtrList<List<solveScalar>> scalarSendBufs1_;
-
-                //- Scalar receive buffers
-                mutable PtrList<List<solveScalar>> scalarRecvBufs1_;
-
-
     // Private Member Functions
 
         //- No copy construct

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C

@@ -381,46 +381,6 @@ void Foam::cyclicAMIPolyPatch::resetAMI(const UList<point>& points) const
         return;
     }
 
-    point refPt = point::max;
-
-    bool restoredFromCache = false;
-    if (AMIPtr_->cacheActive())
-    {
-        const label comm = boundaryMesh().mesh().comm();
-
-        if (UPstream::parRun())
-        {
-            label refProci = -1;
-            if (size() > 0)
-            {
-                refProci = UPstream::myProcNo(comm);
-            }
-            reduce(refProci, maxOp<label>(), UPstream::msgType(), comm);
-
-            if (refProci == UPstream::myProcNo(comm))
-            {
-                refPt = points[meshPoints()[0]];
-            }
-            reduce(refPt, minOp<point>(), UPstream::msgType(), comm);
-        }
-        else
-        {
-            refPt = points[meshPoints()[0]];
-        }
-
-
-        // Sets cache indices to use and time interpolation weight
-        restoredFromCache = AMIPtr_->restoreCache(refPt);
-
-        if (returnReduceOr(restoredFromCache, comm))
-        {
-            // Restored AMI weight and addressing from cache - all done
-            Info<< "AMI: weights and addresses restored from cache" << endl;
-
-            return;
-        }
-    }
-
     const cyclicAMIPolyPatch& nbr = neighbPatch();
     const pointField srcPoints(points, meshPoints());
     pointField nbrPoints(points, nbr.meshPoints());
@@ -472,7 +432,7 @@ void Foam::cyclicAMIPolyPatch::resetAMI(const UList<point>& points) const
         meshTools::writeOBJ(osN, nbrPatch0.localFaces(), nbrPoints);
 
         OFstream osO(t.path()/name() + "_ownerPatch.obj");
-        meshTools::writeOBJ(osO, this->localFaces(), srcPoints);
+        meshTools::writeOBJ(osO, this->localFaces(), localPoints());
     }
 
     // Construct/apply AMI interpolation to determine addressing and weights
@@ -486,8 +446,6 @@ void Foam::cyclicAMIPolyPatch::resetAMI(const UList<point>& points) const
     {
         AMIPtr_->checkSymmetricWeights(true);
     }
-
-    AMIPtr_->addToCache(refPt);
 }
 
 

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H

@@ -99,7 +99,7 @@ protected:
         //- Index of other half
         mutable label nbrPatchID_;
 
-        //- Particle displacement fraction across AMI
+        //- Particle displacement fraction accross AMI
         const scalar fraction_;
 
 
@@ -166,9 +166,6 @@ protected:
             //- Temporary storage for AMI face centres
             mutable vectorField faceCentres0_;
 
-        // Cache
-        bool cacheAMI_;
-
 
     // Protected Member Functions
 
@@ -523,14 +520,9 @@ public:
                 (
                     const Field<Type>& fld,
                     labelRange& sendRequests,
-                    labelRange& recvRequests,
                     PtrList<List<Type>>& sendBuffers,
-                    PtrList<List<Type>>& recvBuffers,
-
-                    labelRange& sendRequests1,
-                    labelRange& recvRequests1,
-                    PtrList<List<Type>>& sendBuffers1,
-                    PtrList<List<Type>>& recvBuffers1
+                    labelRange& recvRequests,
+                    PtrList<List<Type>>& recvBuffers
                 ) const;
 
                 template<class Type>
@@ -538,14 +530,9 @@ public:
                 (
                     const Field<Type>& fld,
                     labelRange& sendRequests,
-                    labelRange& recvRequests,
                     PtrList<List<Type>>& sendBuffers,
-                    PtrList<List<Type>>& recvBuffers,
-
-                    labelRange& sendRequests1,
-                    labelRange& recvRequests1,
-                    PtrList<List<Type>>& sendBuffers1,
-                    PtrList<List<Type>>& recvBuffers1
+                    labelRange& recvRequests,
+                    PtrList<List<Type>>& recvBuffers
                 ) const;
 
                 template<class Type>
@@ -554,8 +541,6 @@ public:
                     const Field<Type>& localFld,
                     const labelRange& requests,  // The receive requests
                     const PtrList<List<Type>>& recvBuffers,
-                    const labelRange& requests1,  // The receive requests
-                    const PtrList<List<Type>>& recvBuffers1,
                     const UList<Type>& defaultValues
                 ) const;
 

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchTemplates.C

@@ -63,7 +63,6 @@ Foam::tmp<Foam::Field<Type>> Foam::cyclicAMIPolyPatch::interpolate
 
     if constexpr (transform_supported)
     {
-        // Only creates the co-ord system if using periodic AMI
         cs.reset(cylindricalCS());
     }
 
@@ -131,7 +130,7 @@ Foam::tmp<Foam::Field<Type>> Foam::cyclicAMIPolyPatch::interpolate
         const tensorField ownT(cs().R(this->faceCentres()));
 
         Field<Type> localDeflt(defaultValues.size());
-        if (defaultValues.size() != 0 && defaultValues.size() == size())
+        if (defaultValues.size() == size())
         {
             // Transform default values into cylindrical coords (using
             // *this faceCentres)
@@ -177,77 +176,27 @@ void Foam::cyclicAMIPolyPatch::initInterpolateUntransformed
 (
     const Field<Type>& fld,
     labelRange& sendRequests,
-    labelRange& recvRequests,
     PtrList<List<Type>>& sendBuffers,
-    PtrList<List<Type>>& recvBuffers,
-
-    labelRange& sendRequests1,
-    labelRange& recvRequests1,
-    PtrList<List<Type>>& sendBuffers1,
-    PtrList<List<Type>>& recvBuffers1
+    labelRange& recvRequests,
+    PtrList<List<Type>>& recvBuffers
 ) const
 {
     const auto& AMI = (owner() ? this->AMI() : neighbPatch().AMI());
 
     if (AMI.distributed() && AMI.comm() != -1)
     {
-        const auto& cache = AMI.cache();
+        const auto& map = (owner() ? AMI.tgtMap() : AMI.srcMap());
 
-        if (cache.index0() == -1 && cache.index1() == -1)
-        {
-            // No caching
-
-            const auto& map = (owner() ? AMI.tgtMap() : AMI.srcMap());
-
-            // Insert send/receive requests (non-blocking)
-            map.send
-            (
-                fld,
-                sendRequests,
-                sendBuffers,
-                recvRequests,
-                recvBuffers,
-                3894+this->index()  // unique offset + patch index
-            );
-        }
-        else
-        {
-            // Caching is active
-
-            cache.setDirection(owner());
-
-            if (cache.index0() != -1)
-            {
-                const auto& map0 = cache.cTgtMapPtr0()();
-
-                // Insert send/receive requests (non-blocking)
-                map0.send
-                (
-                    fld,
-                    sendRequests,
-                    sendBuffers,
-                    recvRequests,
-                    recvBuffers,
-                    3894+this->index()  // unique offset + patch index
-                );
-            }
-
-            if (cache.index1() != -1)
-            {
-                const auto& map1 = cache.cTgtMapPtr1()();
-
-                // Insert send/receive requests (non-blocking)
-                map1.send
-                (
-                    fld,
-                    sendRequests1,
-                    sendBuffers1,
-                    recvRequests1,
-                    recvBuffers1,
-                    3895+this->index()  // unique offset + patch index
-                );
-            }
-        }
+        // Insert send/receive requests (non-blocking)
+        map.send
+        (
+            fld,
+            sendRequests,
+            sendBuffers,
+            recvRequests,
+            recvBuffers,
+            3894+this->index()  // unique offset + patch index
+        );
     }
 }
 
@@ -257,14 +206,9 @@ void Foam::cyclicAMIPolyPatch::initInterpolate
 (
     const Field<Type>& fld,
     labelRange& sendRequests,
-    labelRange& recvRequests,
     PtrList<List<Type>>& sendBuffers,
-    PtrList<List<Type>>& recvBuffers,
-
-    labelRange& sendRequests1,
-    labelRange& recvRequests1,
-    PtrList<List<Type>>& sendBuffers1,
-    PtrList<List<Type>>& recvBuffers1
+    labelRange& recvRequests,
+    PtrList<List<Type>>& recvBuffers
 ) const
 {
     const auto& AMI = (owner() ? this->AMI() : neighbPatch().AMI());
@@ -277,50 +221,46 @@ void Foam::cyclicAMIPolyPatch::initInterpolate
     // Can rotate fields (vector and non-spherical tensors)
     constexpr bool transform_supported = is_rotational_vectorspace_v<Type>;
 
+    [[maybe_unused]]
+    autoPtr<coordSystem::cylindrical> cs;
+
     if constexpr (transform_supported)
     {
-        // Only creates the co-ord system if using periodic AMI
-        // - convert to cylindrical coordinate system
-        auto cs = cylindricalCS();
-
-        if (cs)
-        {
-            Field<Type> localFld(fld.size());
-            const cyclicAMIPolyPatch& nbrPp = this->neighbPatch();
-            const tensorField nbrT(cs().R(nbrPp.faceCentres()));
-            Foam::invTransform(localFld, nbrT, fld);
-
-            initInterpolateUntransformed
-            (
-                localFld,
-                sendRequests,
-                recvRequests,
-                sendBuffers,
-                recvBuffers,
-
-                sendRequests1,
-                recvRequests1,
-                sendBuffers1,
-                recvBuffers1
-            );
-
-            return;
-        }
+        cs.reset(cylindricalCS());
     }
 
-    initInterpolateUntransformed
-    (
-        fld,
-        sendRequests,
-        recvRequests,
-        sendBuffers,
-        recvBuffers,
+    if (!cs)
+    {
+        initInterpolateUntransformed
+        (
+            fld,
+            sendRequests,
+            sendBuffers,
+            recvRequests,
+            recvBuffers
+        );
+    }
+    else if constexpr (transform_supported)
+    {
+        const cyclicAMIPolyPatch& nbrPp = this->neighbPatch();
 
-        sendRequests1,
-        recvRequests1,
-        sendBuffers1,
-        recvBuffers1
-    );
+        Field<Type> localFld(fld.size());
+
+        // Transform to cylindrical coords
+        {
+            const tensorField nbrT(cs().R(nbrPp.faceCentres()));
+            Foam::invTransform(localFld, nbrT, fld);
+        }
+
+        initInterpolateUntransformed
+        (
+            localFld,
+            sendRequests,
+            sendBuffers,
+            recvRequests,
+            recvBuffers
+        );
+    }
 }
 
 
@@ -330,21 +270,14 @@ Foam::tmp<Foam::Field<Type>> Foam::cyclicAMIPolyPatch::interpolate
     const Field<Type>& localFld,
     const labelRange& requests,
     const PtrList<List<Type>>& recvBuffers,
-    const labelRange& requests1,
-    const PtrList<List<Type>>& recvBuffers1,
     const UList<Type>& defaultValues
 ) const
 {
-    // Note: cannot be localFld.size() -> is set to null for distributed AMI
     auto tresult = tmp<Field<Type>>::New(this->size(), Zero);
 
     const auto& AMI = (owner() ? this->AMI() : neighbPatch().AMI());
 
-    const auto& cache = AMI.cache();
-    cache.setDirection(owner());
-
     Field<Type> work;
-    Field<Type> work1;
     if (AMI.distributed())
     {
         if (AMI.comm() == -1)
@@ -352,157 +285,71 @@ Foam::tmp<Foam::Field<Type>> Foam::cyclicAMIPolyPatch::interpolate
             return tresult;
         }
 
-        if (cache.index0() == -1 && cache.index1() == -1)
-        {
-            // No caching
-            const auto& map = (owner() ? AMI.tgtMap() : AMI.srcMap());
+        const auto& map = (owner() ? AMI.tgtMap() : AMI.srcMap());
 
-            // Receive (= copy) data from buffers into work. TBD: receive
-            // directly into slices of work.
-            map.receive
-            (
-                requests,
-                recvBuffers,
-                work,
-                3894+this->index()  // unique offset + patch index
-            );
-        }
-        else
-        {
-            // Using AMI cache
-
-            if (cache.index0() != -1)
-            {
-                cache.cTgtMapPtr0()().receive
-                (
-                    requests,
-                    recvBuffers,
-                    work,
-                    3894+this->index()  // unique offset + patch index
-                );
-            }
-
-            if (cache.index1() != -1)
-            {
-                cache.cTgtMapPtr1()().receive
-                (
-                    requests1,
-                    recvBuffers1,
-                    work1,
-                    3895+this->index()  // unique offset + patch index
-                );
-            }
-        }
+        // Receive (= copy) data from buffers into work. TBD: receive directly
+        // into slices of work.
+        map.receive
+        (
+            requests,
+            recvBuffers,
+            work,
+            3894+this->index()  // unique offset + patch index
+        );
     }
-
     const Field<Type>& fld = (AMI.distributed() ? work : localFld);
-    const Field<Type>& fld1 = (AMI.distributed() ? work1 : localFld);
 
     // Rotate fields (vector and non-spherical tensors)
     constexpr bool transform_supported = is_rotational_vectorspace_v<Type>;
 
-    // Rotate fields (vector and non-spherical tensors) for periodic AMI
-    tensorField ownTransform;
-    Field<Type> localDeflt;
+    [[maybe_unused]]
+    autoPtr<coordSystem::cylindrical> cs;
 
+    // Rotate fields (vector and non-spherical tensors)
     if constexpr (transform_supported)
     {
-        // Only creates the co-ord system if using periodic AMI
-        // - convert to cylindrical coordinate system
-        auto cs = cylindricalCS();
-
-        if (cs)
-        {
-            ownTransform = cs().R(this->faceCentres());
-            localDeflt = defaultValues;
-
-            if (defaultValues.size() == size())
-            {
-                // Transform default values into cylindrical coords (using
-                // *this faceCentres)
-                // We get in UList (why? Copied from cyclicAMI). Convert to
-                // Field so we can use transformField routines.
-                const SubField<Type> defaultSubFld(defaultValues);
-                const Field<Type>& defaultFld(defaultSubFld);
-                Foam::invTransform(localDeflt, ownTransform, defaultFld);
-            }
-        }
+        cs.reset(cylindricalCS());
     }
 
-    const auto& localDefaultValues =
-        localDeflt.size() ? localDeflt : defaultValues;
-
-    if (cache.index0() == -1 && cache.index1() == -1)
+    if (!cs)
     {
-        // No caching
         AMI.weightedSum
         (
             owner(),
             fld,
             tresult.ref(),
-            localDefaultValues
+            defaultValues
+        );
+    }
+    else if constexpr (transform_supported)
+    {
+        const tensorField ownT(cs().R(this->faceCentres()));
+
+        Field<Type> localDeflt(defaultValues.size());
+        if (defaultValues.size() == size())
+        {
+            // Transform default values into cylindrical coords (using
+            // *this faceCentres)
+            // We get in UList (why? Copied from cyclicAMI). Convert to
+            // Field so we can use transformField routines.
+            const SubField<Type> defaultSubFld(defaultValues);
+            const Field<Type>& defaultFld(defaultSubFld);
+            Foam::invTransform(localDeflt, ownT, defaultFld);
+        }
+
+        AMI.weightedSum
+        (
+            owner(),
+            fld,
+            tresult.ref(),
+            localDeflt
         );
 
         // Transform back
-        if (ownTransform.size())
-        {
-            Foam::transform(tresult.ref(), ownTransform, tresult());
-        }
-
-        return tresult;
+        Foam::transform(tresult.ref(), ownT, tresult());
     }
-    else
-    {
-        if (cache.index0() != -1)
-        {
-            AMIInterpolation::weightedSum
-            (
-                AMI.lowWeightCorrection(),
-                cache.cSrcAddress0(),
-                cache.cSrcWeights0(),
-                cache.cSrcWeightsSum0(),
-                fld,
-                multiplyWeightedOp<Type, plusEqOp<Type>>(plusEqOp<Type>()),
-                tresult.ref(),
-                localDefaultValues
-            );
 
-            if (ownTransform.size())
-            {
-                Foam::transform(tresult.ref(), ownTransform, tresult());
-            }
-
-            // Assuming cache weight is zero when index1 is inactive (==-1)
-            tresult.ref() *= (1 - cache.weight());
-        }
-
-        if (cache.index1() != -1)
-        {
-            auto tresult1 = tmp<Field<Type>>::New(this->size(), Zero);
-
-            AMIInterpolation::weightedSum
-            (
-                AMI.lowWeightCorrection(),
-                cache.cSrcAddress1(),
-                cache.cSrcWeights1(),
-                cache.cSrcWeightsSum1(),
-                fld1,
-                multiplyWeightedOp<Type, plusEqOp<Type>>(plusEqOp<Type>()),
-                tresult1.ref(),
-                localDefaultValues
-            );
-
-            if (ownTransform.size())
-            {
-                Foam::transform(tresult1.ref(), ownTransform, tresult1());
-            }
-
-            tresult1.ref() *= cache.weight();
-            tresult.ref() += tresult1();
-        }
-
-        return tresult;
-    }
+    return tresult;
 }
 
 

src/meshTools/Make/files

@@ -281,7 +281,6 @@ processorLOD/faceBox/faceBox.C
 AMI=AMIInterpolation
 $(AMI)/AMIInterpolation/AMIInterpolation.C
 $(AMI)/AMIInterpolation/AMIInterpolationNew.C
-$(AMI)/AMIInterpolation/AMICache.C
 $(AMI)/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.C
 $(AMI)/AMIInterpolation/advancingFrontAMI/advancingFrontAMIParallelOps.C
 $(AMI)/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.C

src/meshTools/algorithms/MeshWave/FaceCellWave.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2025 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -95,102 +95,6 @@ namespace Foam
                 }
             }
     };
-
-    template<class Type, class TrackingData>
-    class interpolate
-    {
-        //- Combine operator for AMIInterpolation
-
-        FaceCellWave<Type, TrackingData>& solver_;
-
-        const cyclicAMIPolyPatch& patch_;
-
-    public:
-
-            interpolate
-            (
-                FaceCellWave<Type, TrackingData>& solver,
-                const cyclicAMIPolyPatch& patch
-            )
-            :
-                solver_(solver),
-                patch_(patch)
-            {}
-
-
-            void operator()
-            (
-                Type& x,
-                const label localFacei,
-                const label f0i,
-                const Type& y0,
-                const label f1i,
-                const Type& y1,
-                const scalar weight
-            ) const
-            {
-                if (y0.valid(solver_.data()))
-                {
-                    if (y1.valid(solver_.data()))
-                    {
-                        x.interpolate
-                        (
-                            solver_.mesh(),
-                            patch_.faceCentres()[localFacei],
-                            f0i,
-                            y0,
-                            f1i,
-                            y1,
-                            weight,
-                            solver_.propagationTol(),
-                            solver_.data()
-                        );
-                    }
-                    else
-                    {
-                        // Convert patch face into mesh face
-                        label meshFacei = -1;
-                        if (patch_.owner())
-                        {
-                            meshFacei = patch_.start() + f0i;
-                        }
-                        else
-                        {
-                            meshFacei = patch_.neighbPatch().start() + f0i;
-                        }
-                        x.updateFace
-                        (
-                            solver_.mesh(),
-                            meshFacei,
-                            y0,
-                            solver_.propagationTol(),
-                            solver_.data()
-                        );
-                    }
-                }
-                else if (y1.valid(solver_.data()))
-                {
-                    // Convert patch face into mesh face
-                    label meshFacei = -1;
-                    if (patch_.owner())
-                    {
-                        meshFacei = patch_.start() + f1i;
-                    }
-                    else
-                    {
-                        meshFacei = patch_.neighbPatch().start() + f1i;
-                    }
-                    x.updateFace
-                    (
-                        solver_.mesh(),
-                        meshFacei,
-                        y1,
-                        solver_.propagationTol(),
-                        solver_.data()
-                    );
-                }
-            }
-    };
 }
 
 
@@ -880,43 +784,18 @@ void Foam::FaceCellWave<Type, TrackingData>::handleAMICyclicPatches()
                 // Transfer sendInfo to cycPatch
                 combine<Type, TrackingData> cmb(*this, cycPatch);
 
-                // Linear interpolation
-                interpolate<Type, TrackingData> interp(*this, cycPatch);
-
-                const auto& AMI =
-                (
-                    cycPatch.owner()
-                  ? cycPatch.AMI()
-                  : cycPatch.neighbPatch().AMI()
-                );
-
                 if (cycPatch.applyLowWeightCorrection())
                 {
-                    const List<Type> defVals
+                    List<Type> defVals
                     (
                         cycPatch.patchInternalList(allCellInfo_)
                     );
-                    AMI.interpolate
-                    (
-                        cycPatch.owner(),
-                        sendInfo,
-                        cmb,
-                        interp,
-                        receiveInfo,
-                        defVals
-                    );
+
+                    cycPatch.interpolate(sendInfo, cmb, receiveInfo, defVals);
                 }
                 else
                 {
-                    AMI.interpolate
-                    (
-                        cycPatch.owner(),
-                        sendInfo,
-                        cmb,
-                        interp,
-                        receiveInfo,
-                        UList<Type>::null() // no default values needed
-                    );
+                    cycPatch.interpolate(sendInfo, cmb, receiveInfo);
                 }
             }
 

src/meshTools/cellClassification/cellInfo.H

@@ -197,23 +197,6 @@ public:
             template<class TrackingData>
             inline bool equal(const cellInfo&, TrackingData& td) const;
 
-            //- Interpolate between two values (lerp). Returns true if
-            //- causes changes. Not sure if needs to be specialised between
-            //- face and cell and what index is needed...
-            template<class TrackingData>
-            inline bool interpolate
-            (
-                const polyMesh&,
-                const point& pt,
-                const label i0,
-                const cellInfo& f0,
-                const label i1,
-                const cellInfo& f1,
-                const scalar weight,
-                const scalar tol,
-                TrackingData& td
-            );
-
 
     // Member Operators
 

src/meshTools/cellClassification/cellInfoI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020,2025 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -249,31 +249,6 @@ inline bool Foam::cellInfo::equal
 }
 
 
-template<class TrackingData>
-inline bool Foam::cellInfo::interpolate
-(
-    const polyMesh&,
-    const point& pt,
-    const label i0,
-    const cellInfo& f0,
-    const label i1,
-    const cellInfo& f1,
-    const scalar weight,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    if (!f0.valid(td))
-    {
-        return update(f1, -1, -1, -1, -1, td);
-    }
-    else
-    {
-        return update(f0, -1, -1, -1, -1, td);
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::cellInfo::operator==

src/meshTools/cellDist/wallPoint/wallPoint.H

@@ -76,16 +76,6 @@ class wallPoint
 
     // Private Member Functions
 
-        // Update this with w2 if distance less
-        template<class TrackingData>
-        inline bool update
-        (
-            const wallPoint& w2,
-            const scalar dist2,
-            const scalar tol,
-            TrackingData& td
-        );
-
         //- Evaluate distance to point.
         //  Update distSqr, origin from whomever is nearer pt.
         //  \return true if w2 is closer to point, false otherwise.
@@ -216,41 +206,10 @@ public:
                 TrackingData& td
             );
 
-            //- Interpolate different values
-            template<class TrackingData>
-            wallPoint
-            (
-                const polyMesh&,
-                const scalar weight,
-                const label face0,
-                const wallPoint& info0,
-                const label face1,
-                const wallPoint& info1,
-                const scalar tol,
-                TrackingData& td
-            );
-
             //- Test for equality, with TrackingData
             template<class TrackingData>
             inline bool equal(const wallPoint&, TrackingData& td) const;
 
-            //- Interpolate between two values (lerp). Returns true if
-            //- causes changes. Not sure if needs to be specialised between
-            //- face and cell and what index is needed...
-            template<class TrackingData>
-            inline bool interpolate
-            (
-                const polyMesh&,
-                const point& pt,
-                const label i0,
-                const wallPoint& f0,
-                const label i1,
-                const wallPoint& f1,
-                const scalar weight,
-                const scalar tol,
-                TrackingData& td
-            );
-
 
     // Member Operators
 

src/meshTools/cellDist/wallPoint/wallPointI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020,2025 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,12 +35,21 @@ License
 template<class TrackingData>
 inline bool Foam::wallPoint::update
 (
+    const point& pt,
     const wallPoint& w2,
-    const scalar dist2,
     const scalar tol,
     TrackingData& td
 )
 {
+    //Already done in calling algorithm
+    //if (w2.origin() == origin_)
+    //{
+    //    // Shortcut. Same input so same distance.
+    //    return false;
+    //}
+
+    const scalar dist2 = magSqr(pt - w2.origin());
+
     if (!valid(td))
     {
         // current not yet set so use any value
@@ -74,26 +83,6 @@ inline bool Foam::wallPoint::update
 }
 
 
-// Update this with w2 if w2 nearer to pt.
-template<class TrackingData>
-inline bool Foam::wallPoint::update
-(
-    const point& pt,
-    const wallPoint& w2,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    return update
-    (
-        w2,
-        magSqr(pt - w2.origin()),
-        tol,
-        td
-    );
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 inline Foam::wallPoint::wallPoint()
@@ -271,58 +260,6 @@ inline bool Foam::wallPoint::equal
 }
 
 
-template<class TrackingData>
-inline bool Foam::wallPoint::interpolate
-(
-    const polyMesh&,
-    const point& pt,
-    const label i0,
-    const wallPoint& f0,
-    const label i1,
-    const wallPoint& f1,
-    const scalar weight,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    if (f0.valid(td))
-    {
-        if (f1.valid(td))
-        {
-            // What do we do about the origin? Should be consistent with the
-            // distance? But then interpolating between a point 'on the left'
-            // and a point 'on the right' the interpolate might just be on
-            // top of us so suddenly setting the distance to be zero...
-
-            const scalar dist2 =
-                sqr(lerp(sqrt(f0.distSqr()), sqrt(f1.distSqr()), weight));
-
-            if (f0.distSqr() < f1.distSqr())
-            {
-                // Update with interpolated distance and f0 origin
-                return update(f0, dist2, tol, td);
-            }
-            else
-            {
-                return update(f1, dist2, tol, td);
-            }
-        }
-        else
-        {
-            return update(f0, f0.distSqr(), tol, td);
-        }
-    }
-    else if (f1.valid(td))
-    {
-        return update(f1, f1.distSqr(), tol, td);
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::wallPoint::operator==

src/meshTools/meshStructure/topoDistanceData.H

@@ -196,23 +196,6 @@ public:
                 TrackingData& td
             ) const;
 
-            //- Interpolate between two values (lerp). Returns true if
-            //- causes changes. Not sure if needs to be specialised between
-            //- face and cell and what index is needed...
-            template<class TrackingData>
-            inline bool interpolate
-            (
-                const polyMesh&,
-                const point& pt,
-                const label i0,
-                const topoDistanceData<Type>& f0,
-                const label i1,
-                const topoDistanceData<Type>& f1,
-                const scalar weight,
-                const scalar tol,
-                TrackingData& td
-            );
-
 
     // Member Operators
 

src/meshTools/meshStructure/topoDistanceDataI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2025 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -199,38 +199,6 @@ inline bool Foam::topoDistanceData<Type>::equal
 }
 
 
-template<class Type>
-template<class TrackingData>
-inline bool Foam::topoDistanceData<Type>::interpolate
-(
-    const polyMesh&,
-    const point& pt,
-    const label i0,
-    const topoDistanceData<Type>& f0,
-    const label i1,
-    const topoDistanceData<Type>& f1,
-    const scalar weight,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    // Topological distance - how to interpolate?
-    if (!valid(td))
-    {
-        if (f0.valid(td))
-        {
-            this->operator=(f0);
-        }
-        else
-        {
-            this->operator=(f1);
-        }
-        return true;
-    }
-    return false;
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Type>

src/meshTools/regionSplit/minData.H

@@ -174,23 +174,6 @@ public:
             template<class TrackingData>
             inline bool equal(const minData&, TrackingData& td) const;
 
-            //- Interpolate between two values (lerp). Returns true if
-            //- causes changes. Not sure if needs to be specialised between
-            //- face and cell and what index is needed...
-            template<class TrackingData>
-            inline bool interpolate
-            (
-                const polyMesh&,
-                const point& pt,
-                const label i0,
-                const minData& f0,
-                const label i1,
-                const minData& f1,
-                const scalar weight,
-                const scalar tol,
-                TrackingData& td
-            );
-
 
     // Member Operators
 

src/meshTools/regionSplit/minDataI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2014-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2025 OpenCFD Ltd.
+    Copyright (C) 2015-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -173,35 +173,6 @@ inline bool Foam::minData::equal
 }
 
 
-template<class TrackingData>
-inline bool Foam::minData::interpolate
-(
-    const polyMesh& mesh,
-    const point& pt,
-    const label i0,
-    const minData& f0,
-    const label i1,
-    const minData& f1,
-    const scalar weight,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    if (f0.valid(td))
-    {
-        return updateFace(mesh, -1, f0, tol, td);
-    }
-    if (f1.valid(td))
-    {
-        return updateFace(mesh, -1, f1, tol, td);
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::minData::operator==

src/overset/cellCellStencil/inverseDistance/meshToMeshData.H

@@ -200,23 +200,6 @@ public:
                 TrackingData&
             ) const;
 
-            //- Interpolate between two values (lerp). Returns true if
-            //- causes changes. Not sure if needs to be specialised between
-            //- face and cell and what index is needed...
-            template<class TrackingData>
-            inline bool interpolate
-            (
-                const polyMesh&,
-                const point& pt,
-                const label i0,
-                const meshToMeshData& f0,
-                const label i1,
-                const meshToMeshData& f1,
-                const scalar weight,
-                const scalar tol,
-                TrackingData& td
-            );
-
 
     // Member Operators
 

src/overset/cellCellStencil/inverseDistance/meshToMeshDataI.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2025 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -205,35 +205,6 @@ inline bool Foam::meshToMeshData::equal
 }
 
 
-template<class TrackingData>
-inline bool Foam::meshToMeshData::interpolate
-(
-    const polyMesh& mesh,
-    const point& pt,
-    const label i0,
-    const meshToMeshData& f0,
-    const label i1,
-    const meshToMeshData& f1,
-    const scalar weight,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    if (f0.valid(td))
-    {
-        return updateFace(mesh, -1, f0, tol, td);
-    }
-    if (f1.valid(td))
-    {
-        return updateFace(mesh, -1, f1, tol, td);;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::meshToMeshData::operator==

src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2023 OpenCFD Ltd.
+    Copyright (C) 2020-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,9 +26,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "dynamicContactAngleForce.H"
-#include "addToRunTimeSelectionTable.H"
 #include "Function1.H"
 #include "distributionModel.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,6 +59,7 @@ dynamicContactAngleForce::dynamicContactAngleForce
 )
 :
     contactAngleForce(typeName, film, dict),
+    thetaPtr_(nullptr),
     U_vs_thetaPtr_
     (
         Function1<scalar>::NewIfPresent
@@ -80,46 +81,56 @@ dynamicContactAngleForce::dynamicContactAngleForce
         )
     ),
     rndGen_(label(0)),
-    distribution_
-    (
-        distributionModel::New
-        (
-            coeffDict_.subDict("distribution"),
-            rndGen_
-        )
-    )
+    distributionPtr_(nullptr)
 {
     if (U_vs_thetaPtr_ && T_vs_thetaPtr_)
     {
         FatalIOErrorInFunction(dict)
-            << "Entries Utheta and Ttheta could not be used together"
+            << "Only one of Utheta or Ttheta should be provided; "
+            << "both inputs cannot be used together."
             << abort(FatalIOError);
     }
+
+
+    thetaPtr_.emplace
+    (
+        IOobject
+        (
+            IOobject::scopedName(typeName, "theta"),
+            film.regionMesh().time().timeName(),
+            film.regionMesh().thisDb(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        film.regionMesh(),
+        dimensionedScalar(dimless, Zero)
+    );
+
+
+    if (coeffDict_.findEntry("distribution"))
+    {
+        distributionPtr_ = distributionModel::New
+        (
+            coeffDict_.subDict("distribution"),
+            rndGen_
+        );
+    }
 }
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+dynamicContactAngleForce::~dynamicContactAngleForce()
+{}  // distributionModel was forward declared
+
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 tmp<areaScalarField> dynamicContactAngleForce::theta() const
 {
-    auto ttheta = tmp<areaScalarField>::New
-    (
-        IOobject
-        (
-            IOobject::scopedName(typeName, "theta"),
-            film().regionMesh().time().timeName(),
-            film().regionMesh().thisDb(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            IOobject::NO_REGISTER
-        ),
-        film().regionMesh(),
-        dimensionedScalar(dimless, Zero)
-    );
-    areaScalarField& theta = ttheta.ref();
+    areaScalarField& theta = thetaPtr_.ref();
     scalarField& thetai = theta.ref();
 
-
     if (U_vs_thetaPtr_)
     {
         // Initialize with the function of film speed
@@ -136,13 +147,16 @@ tmp<areaScalarField> dynamicContactAngleForce::theta() const
         thetai = T_vs_thetaPtr_->value(T());
     }
 
-    // Add the stochastic perturbation
-    forAll(thetai, facei)
+    if (distributionPtr_)
     {
-        thetai[facei] += distribution_->sample();
+        // Add the stochastic perturbation
+        forAll(thetai, facei)
+        {
+            thetai[facei] += distributionPtr_->sample();
+        }
     }
 
-    return ttheta;
+    return thetaPtr_();
 }
 
 

src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -112,6 +112,9 @@ class dynamicContactAngleForce
 {
     // Private Data
 
+        //- Contact-angle field in degrees
+        mutable autoPtr<areaScalarField> thetaPtr_;
+
         //- Contact angle as a function of film speed
         autoPtr<Function1<scalar>> U_vs_thetaPtr_;
 
@@ -121,17 +124,8 @@ class dynamicContactAngleForce
         //- Random number generator
         Random rndGen_;
 
-        //- Parcel size PDF model
-        const autoPtr<distributionModel> distribution_;
-
-
-    // Private Member Functions
-
-        //- No copy construct
-        dynamicContactAngleForce(const dynamicContactAngleForce&) = delete;
-
-        //- No copy assignment
-        void operator=(const dynamicContactAngleForce&) = delete;
+        //- Stochastic perturbation model for contact angle
+        autoPtr<distributionModel> distributionPtr_;
 
 
 protected:
@@ -148,7 +142,7 @@ public:
 
     // Constructors
 
-        //- Construct from surface film model
+        //- Construct from surface film model and dictionary
         dynamicContactAngleForce
         (
             liquidFilmBase& film,
@@ -157,7 +151,7 @@ public:
 
 
     //- Destructor
-    virtual ~dynamicContactAngleForce() = default;
+    virtual ~dynamicContactAngleForce();
 };
 
 

tutorials/incompressible/pimpleFoam/laminar/mixerVesselAMI2D/mixerVesselAMI2D/system/blockMeshDict.m4

@@ -788,14 +788,9 @@ boundary
 
     AMI1
     {
-        cacheSize       360;
-
         type            cyclicAMI;
         neighbourPatch  AMI2;
         transform       noOrdering;
-        rotationCentre  (0 0 0);
-        rotationAxis    (0 0 1);
-
 /* optional
         surface
         {
@@ -820,13 +815,9 @@ boundary
 
     AMI2
     {
-        cacheSize       360;
-
         type            cyclicAMI;
         neighbourPatch  AMI1;
         transform       noOrdering;
-        rotationCentre  (0 0 0);
-        rotationAxis    (0 0 1);
 /* optional
         surface
         {

tutorials/incompressible/simpleFoam/pitzDaily_fused/Allrun

@@ -19,4 +19,12 @@ sed "s/GAUSS/fusedGauss/g" system/fvSchemes.template > system/fvSchemes
 runParallel -s fusedGauss $(getApplication)
 
 
+echo "Updating fvSolution to cache gradients"
+sed 's/GRADIENT/"grad(.*)"/g' system/fvSolution.template > system/fvSolution
+
+runParallel -s fusedGaussCached $(getApplication)
+
+# Restore
+cp system/fvSolution.template system/fvSolution
+
 #------------------------------------------------------------------------------

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/controlDict

@@ -16,6 +16,15 @@ FoamFile
 
 libs            (fusedFiniteVolume);
 
+/*
+DebugSwitches
+{
+    fusedGauss  1;
+    Gauss       1;
+    solution    1;
+}
+*/
+
 application     simpleFoam;
 
 startFrom       startTime;

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/fvSolution

@@ -55,5 +55,11 @@ relaxationFactors
     }
 }
 
+cache
+{
+    // Cache all grads
+    GRADIENT;
+}
+
 
 // ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/fvSolution.template

@@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2506                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          GAMG;
+        tolerance       1e-06;
+        relTol          0.1;
+        smoother        GaussSeidel;
+    }
+
+    "(U|k|epsilon|omega|f|v2)"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    consistent      yes;
+
+    residualControl
+    {
+        p               1e-2;
+        U               1e-3;
+        "(k|epsilon|omega|f|v2)" 1e-3;
+    }
+}
+
+relaxationFactors
+{
+    equations
+    {
+        U               0.9; // 0.9 is more stable but 0.95 more convergent
+        ".*"            0.9; // 0.9 is more stable but 0.95 more convergent
+    }
+}
+
+cache
+{
+    // Cache all grads
+    GRADIENT;
+}
+
+
+// ************************************************************************* //