ENH: extend faPatchField, fvPatchField to support non-allocating snGrad()

Co-authored-by: <Mark.Olesen@keysight.com>

src/finiteArea/fields/faPatchFields/basic/basicSymmetry/basicSymmetryFaPatchField.C

@@ -82,6 +82,41 @@ Foam::basicSymmetryFaPatchField<Type>::basicSymmetryFaPatchField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+void Foam::basicSymmetryFaPatchField<Type>::snGrad(UList<Type>& result) const
+{
+    if constexpr (!is_rotational_vectorspace_v<Type>)
+    {
+        // Rotational-invariant type : treat like zero-gradient
+        result = Foam::zero{};
+    }
+    else
+    {
+        // Get patch internal field, stored temporarily in result
+        this->patchInternalField(result);
+        const auto& pif = result;
+
+        tmp<vectorField> tnHat = this->patch().edgeNormals();
+        const auto& nHat = tnHat();
+
+        const auto& dc = this->patch().deltaCoeffs();
+
+        const label len = result.size();
+
+        // (dc/2.0)*(transform(I - 2.0*sqr(nHat), iF) - iF);
+
+        for (label i = 0; i < len; ++i)
+        {
+            result[i] =
+            (
+                (0.5*dc[i])
+              * (transform(I - 2.0*sqr(nHat[i]), pif[i]) - pif[i])
+            );
+        }
+    }
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::basicSymmetryFaPatchField<Type>::snGrad() const
@@ -93,17 +128,9 @@ Foam::basicSymmetryFaPatchField<Type>::snGrad() const
     }
     else
     {
-        tmp<vectorField> nHat = this->patch().edgeNormals();
+        auto tresult = tmp<Field<Type>>::New(this->size());
-
+        this->snGrad(static_cast<UList<Type>&>(tresult.ref()));
-        const auto& dc = this->patch().deltaCoeffs();
+        return tresult;
-
-        const Field<Type> pif(this->patchInternalField());
-
-        return
-        (
-            (0.5*dc)
-          * (transform(I - 2.0*sqr(nHat), pif) - pif)
-        );
     }
 }
 

src/finiteArea/fields/faPatchFields/basic/basicSymmetry/basicSymmetryFaPatchField.H

@@ -119,21 +119,20 @@ public:
 
     // Member Functions
 
-        // Evaluation functions
+        //- Retrieve patch-normal gradient at boundary
+        virtual tmp<Field<Type>> snGrad() const;
 
-            //- Return gradient at boundary
+        //- Retrieve patch-normal gradient at boundary
-            virtual tmp<Field<Type>> snGrad() const;
+        virtual void snGrad(UList<Type>& result) const;
 
-            //- Evaluate the patch field
+        //- Evaluate the patch field
-            //  Default argument needed to allow call in constructors
+        virtual void evaluate
-            //  HJ, 30/Jun/2009
+        (
-            virtual void evaluate
+            const Pstream::commsTypes commsType = Pstream::commsTypes::buffered
-            (
+        );
-                const Pstream::commsTypes commsType = Pstream::commsTypes::buffered
-            );
 
-            //- Return face-gradient transform diagonal
+        //- Return face-gradient transform diagonal
-            virtual tmp<Field<Type>> snGradTransformDiag() const;
+        virtual tmp<Field<Type>> snGradTransformDiag() const;
 
 
     // Member Operators

src/finiteArea/fields/faPatchFields/basic/coupled/coupledFaPatchField.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -107,12 +108,40 @@ Foam::coupledFaPatchField<Type>::coupledFaPatchField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+void Foam::coupledFaPatchField<Type>::snGrad
+(
+    UList<Type>& result
+) const
+{
+    // Get patch neighbour field, store temporarily in result
+    this->patchNeighbourField(result);
+    const auto& pnf = result;
+
+    // Same as patchInternalField(...), assuming edgeFaces are an indirection
+    // into internal field, but without additional storage...
+    const auto& addr = this->patch().edgeFaces();
+    const auto& iF = this->primitiveField();
+
+    const auto& deltaCoeffs = this->patch().deltaCoeffs();
+
+    // snGrad = deltaCoeffs * (patchNeighbourField - patchInternalField)
+
+    const label len = result.size();
+
+    for (label i = 0; i < len; ++i)
+    {
+        result[i] = deltaCoeffs[i]*(pnf[i] - iF[addr[i]]);
+    }
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>> Foam::coupledFaPatchField<Type>::snGrad() const
 {
-    return
+    auto tresult = tmp<Field<Type>>::New(this->size());
-        (patchNeighbourField() - this->patchInternalField())
+    this->snGrad(static_cast<UList<Type>&>(tresult.ref()));
-       *this->patch().deltaCoeffs();
+    return tresult;
 }
 
 

src/finiteArea/fields/faPatchFields/basic/coupled/coupledFaPatchField.H

@@ -148,6 +148,9 @@ public:
             //- Return patch-normal gradient
             virtual tmp<Field<Type>> snGrad() const;
 
+            //- Retrieve patch-normal gradient
+            virtual void snGrad(UList<Type>& result) const;
+
             //- Initialise the evaluation of the patch field
             virtual void initEvaluate
             (

src/finiteArea/fields/faPatchFields/basic/zeroGradient/zeroGradientFaPatchField.H

@@ -119,14 +119,22 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
-        // Evaluation functions
+        // Evaluation Functions
 
-            //- Return gradient at boundary
+            //- Return patch-normal gradient
             virtual tmp<Field<Type>> snGrad() const
             {
-                return tmp<Field<Type>>::New(this->size(), Zero);
+                // zero-gradient
+                return tmp<Field<Type>>::New(this->size(), Foam::zero{});
+            }
+
+            //- Retrieve patch-normal gradient [contiguous storage]
+            virtual void snGrad(UList<Type>& result) const
+            {
+                // zero-gradient
+                result = Foam::zero{};
             }
 
             //- Evaluate the patch field
@@ -137,25 +145,25 @@ public:
             );
 
             //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the value of this patchField with given weights
+            //- evaluation of the value of this patchField with given weights
             virtual tmp<Field<Type>> valueInternalCoeffs
             (
                 const tmp<scalarField>&
             ) const;
 
             //- Return the matrix source coefficients corresponding to the
-            //  evaluation of the value of this patchField with given weights
+            //- evaluation of the value of this patchField with given weights
             virtual tmp<Field<Type>> valueBoundaryCoeffs
             (
                 const tmp<scalarField>&
             ) const;
 
             //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
+            //- evaluation of the gradient of this patchField
             virtual tmp<Field<Type>> gradientInternalCoeffs() const;
 
             //- Return the matrix source coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
+            //- evaluation of the gradient of this patchField
             virtual tmp<Field<Type>> gradientBoundaryCoeffs() const;
 };
 

src/finiteArea/fields/faPatchFields/constraint/empty/emptyFaPatchField.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -125,7 +125,7 @@ public:
     virtual ~emptyFaPatchField() = default;
 
 
-    // Member functions
+    // Member Functions
 
         // Mapping functions
 
@@ -145,7 +145,7 @@ public:
             {}
 
 
-        // Evaluation functions
+        // Evaluation Functions
 
             //- Update the coefficients associated with the patch field
             //  This only checks to see the case is actually 1D or 2D
@@ -154,7 +154,7 @@ public:
 
 
             //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the value of this patchField with given weights
+            //- evaluation of the value of this patchField with given weights
             virtual tmp<Field<Type>> valueInternalCoeffs
             (
                 const tmp<scalarField>&
@@ -164,7 +164,7 @@ public:
             }
 
             //- Return the matrix source coefficients corresponding to the
-            //  evaluation of the value of this patchField with given weights
+            //- evaluation of the value of this patchField with given weights
             virtual tmp<Field<Type>> valueBoundaryCoeffs
             (
                 const tmp<scalarField>&
@@ -174,20 +174,31 @@ public:
             }
 
             //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
+            //- evaluation of the gradient of this patchField
             tmp<Field<Type>> gradientInternalCoeffs() const
             {
                 return tmp<Field<Type>>::New();
             }
 
             //- Return the matrix source coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
+            //- evaluation of the gradient of this patchField
             tmp<Field<Type>> gradientBoundaryCoeffs() const
             {
                 return tmp<Field<Type>>::New();
             }
 
 
+    // Contiguous storage
+
+        //- Retrieve patch-normal gradient [contiguous storage].
+        //- Placeholder value is zero (treated like zero-gradient).
+        virtual void snGrad(UList<Type>& result) const
+        {
+            // Treat like zero-gradient
+            result = Foam::zero{};
+        }
+
+
     // Member Functions
 
         //- Write without "value" entry!

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.C

@@ -305,9 +305,35 @@ void Foam::processorFaPatchField<Type>::evaluate
 }
 
 
+template<class Type>
+void Foam::processorFaPatchField<Type>::snGrad
+(
+    UList<Type>& result
+) const
+{
+    // Get patch internal field, store temporarily in result
+    this->patchInternalField(result);
+    const auto& pif = result;
+
+    const Field<Type>& pnf = *this;
+    const auto& deltaCoeffs = this->patch().deltaCoeffs();
+
+    const label len = result.size();
+
+    for (label i = 0; i < len; ++i)
+    {
+        result[i] = deltaCoeffs[i]*(pnf[i] - pif[i]);
+    }
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>> Foam::processorFaPatchField<Type>::snGrad() const
 {
+    // OR
+    // auto tfld = tmp<Field<Type>>::New(this->size());
+    // this->snGrad(static_cast<UList<Type>&>(tfld.ref()));
+    // return tfld;
     return this->patch().deltaCoeffs()*(*this - this->patchInternalField());
 }
 

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.H

@@ -210,6 +210,9 @@ public:
         //- Return patch-normal gradient
         virtual tmp<Field<Type>> snGrad() const;
 
+        //- Retrieve patch-normal gradient
+        virtual void snGrad(UList<Type>& result) const;
+
 
         // Coupled interface functionality
 

src/finiteArea/fields/faPatchFields/faPatchField/faPatchField.C

@@ -202,10 +202,31 @@ void Foam::faPatchField<Type>::check(const faPatchField<Type>& rhs) const
 }
 
 
+template<class Type>
+void Foam::faPatchField<Type>::snGrad(UList<Type>& result) const
+{
+    // Get patch internal field, store temporarily in result
+    this->patchInternalField(result);
+    const auto& pif = result;
+
+    const Field<Type>& pfld = *this;
+    const auto& dc = patch().deltaCoeffs();
+
+    const label len = result.size();
+
+    for (label i = 0; i < len; ++i)
+    {
+        result[i] = dc[i]*(pfld[i] - pif[i]);
+    }
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>> Foam::faPatchField<Type>::snGrad() const
 {
-    return (*this - patchInternalField())*patch().deltaCoeffs();
+    auto tfld = tmp<Field<Type>>::New(this->size());
+    this->snGrad(static_cast<UList<Type>&>(tfld.ref()));
+    return tfld;
 }
 
 

src/finiteArea/fields/faPatchFields/faPatchField/faPatchField.H

@@ -582,6 +582,31 @@ public:
             //- Return patch-normal gradient
             virtual tmp<Field<Type>> snGrad() const;
 
+            //- Retrieve patch-normal gradient [contiguous storage]
+            virtual void snGrad(UList<Type>& result) const;
+
+            //- Return patch-normal gradient for coupled-patches
+            //- using the deltaCoeffs provided
+            virtual tmp<Field<Type>> snGrad
+            (
+                const scalarField& deltaCoeffs
+            ) const
+            {
+                NotImplemented;
+                return *this;
+            }
+
+            //- Retrieve patch-normal gradient for coupled-patches
+            //- using the deltaCoeffs provided [contiguous storage]
+            virtual void snGrad
+            (
+                const scalarField& deltaCoeffs,
+                UList<Type>&
+            ) const
+            {
+                NotImplemented;
+            }
+
             //- Return internal field next to patch
             virtual tmp<Field<Type>> patchInternalField() const;
 
@@ -677,6 +702,69 @@ public:
             }
 
 
+        // Contiguous storage
+
+            //- Retrieve the matrix diagonal coefficients corresponding to the
+            //- evaluation of the value of this patchField with given weights
+            virtual void valueInternalCoeffs
+            (
+                const tmp<Field<scalar>>&,
+                UList<Type>&
+            ) const
+            {
+                NotImplemented;
+            }
+
+            //- Retrieve the matrix source coefficients corresponding to the
+            //- evaluation of the value of this patchField with given weights
+            virtual void valueBoundaryCoeffs
+            (
+                const tmp<Field<scalar>>&,
+                UList<Type>&
+            ) const
+            {
+                NotImplemented;
+            }
+
+            //- Retrieve the matrix diagonal coefficients corresponding to the
+            //- evaluation of the gradient of this patchField
+            virtual void gradientInternalCoeffs(UList<Type>&) const
+            {
+                NotImplemented;
+            }
+
+            //- Retrieve the matrix diagonal coefficients corresponding to the
+            //- evaluation of the gradient of this coupled patchField
+            //- using the deltaCoeffs provided
+            virtual void gradientInternalCoeffs
+            (
+                const scalarField& deltaCoeffs,
+                UList<Type>&
+            ) const
+            {
+                NotImplemented;
+            }
+
+            //- Retrieve the matrix source coefficients corresponding to the
+            //- evaluation of the gradient of this patchField
+            virtual void gradientBoundaryCoeffs(UList<Type>&) const
+            {
+                NotImplemented;
+            }
+
+            //- Retrieve the matrix source coefficients corresponding to the
+            //- evaluation of the gradient of this coupled patchField
+            //- using the deltaCoeffs provided
+            virtual void gradientBoundaryCoeffs
+            (
+                const scalarField& deltaCoeffs,
+                UList<Type>&
+            ) const
+            {
+                NotImplemented;
+            }
+
+
     // Other
 
         //- Write

src/finiteVolume/fields/fvPatchFields/basic/basicSymmetry/basicSymmetryFvPatchField.C

@@ -82,6 +82,41 @@ Foam::basicSymmetryFvPatchField<Type>::basicSymmetryFvPatchField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+void Foam::basicSymmetryFvPatchField<Type>::snGrad(UList<Type>& result) const
+{
+    if constexpr (!is_rotational_vectorspace_v<Type>)
+    {
+        // Rotational-invariant type : treat like zero-gradient
+        result = Foam::zero{};
+    }
+    else
+    {
+        // Get patch internal field, stored temporarily in result
+        this->patchInternalField(result);
+        const auto& pif = result;
+
+        tmp<vectorField> tnHat = this->patch().nf();
+        const auto& nHat = tnHat();
+
+        const auto& dc = this->patch().deltaCoeffs();
+
+        const label len = result.size();
+
+        // (dc/2.0)*(transform(I - 2.0*sqr(nHat), iF) - iF);
+
+        for (label i = 0; i < len; ++i)
+        {
+            result[i] =
+            (
+                (0.5*dc[i])
+              * (transform(I - 2.0*sqr(nHat[i]), pif[i]) - pif[i])
+            );
+        }
+    }
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::basicSymmetryFvPatchField<Type>::snGrad() const
@@ -93,17 +128,9 @@ Foam::basicSymmetryFvPatchField<Type>::snGrad() const
     }
     else
     {
-        tmp<vectorField> nHat = this->patch().nf();
+        auto tresult = tmp<Field<Type>>::New(this->size());
-
+        this->snGrad(static_cast<UList<Type>&>(tresult.ref()));
-        const auto& dc = this->patch().deltaCoeffs();
+        return tresult;
-
-        const Field<Type> pif(this->patchInternalField());
-
-        return
-        (
-            (0.5*dc)
-          * (transform(I - 2.0*sqr(nHat), pif) - pif)
-        );
     }
 }
 

src/finiteVolume/fields/fvPatchFields/basic/basicSymmetry/basicSymmetryFvPatchField.H

@@ -119,9 +119,12 @@ public:
 
     // Member Functions
 
-        //- Return gradient at boundary
+        //- Retrieve patch-normal gradient at boundary
         virtual tmp<Field<Type>> snGrad() const;
 
+        //- Retrieve patch-normal gradient at boundary
+        virtual void snGrad(UList<Type>& result) const;
+
         //- Evaluate the patch field
         virtual void evaluate
         (

src/finiteVolume/fields/fvPatchFields/basic/coupled/coupledFvPatchField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2023 OpenCFD Ltd.
+    Copyright (C) 2023-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -108,15 +108,42 @@ Foam::coupledFvPatchField<Type>::coupledFvPatchField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+void Foam::coupledFvPatchField<Type>::snGrad
+(
+    const scalarField& deltaCoeffs,
+    UList<Type>& result
+) const
+{
+    // Get patch neighbour field, store temporarily in result
+    this->patchNeighbourField(result);
+    const auto& pnf = result;
+
+    // Same as patchInternalField(...), assuming faceCells are an indirection
+    // into internal field, but without additional storage...
+    const auto& addr = this->patch().faceCells();
+    const auto& iF = this->primitiveField();
+
+    // snGrad = deltaCoeffs * (patchNeighbourField - patchInternalField)
+
+    const label len = result.size();
+
+    for (label i = 0; i < len; ++i)
+    {
+        result[i] = deltaCoeffs[i]*(pnf[i] - iF[addr[i]]);
+    }
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>> Foam::coupledFvPatchField<Type>::snGrad
 (
     const scalarField& deltaCoeffs
 ) const
 {
-    return
+    auto tresult = tmp<Field<Type>>::New(this->size());
-        deltaCoeffs
+    this->snGrad(deltaCoeffs, tresult.ref());
-       *(this->patchNeighbourField() - this->patchInternalField());
+    return tresult;
 }
 
 
@@ -152,6 +179,44 @@ void Foam::coupledFvPatchField<Type>::evaluate(const Pstream::commsTypes)
 }
 
 
+template<class Type>
+void Foam::coupledFvPatchField<Type>::valueInternalCoeffs
+(
+    const tmp<scalarField>& tweights,
+    UList<Type>& result
+) const
+{
+    const auto& w = tweights();
+
+    const label len = result.size();
+
+    for (label i = 0; i < len; ++i)
+    {
+        result[i] = Type(pTraits<Type>::one)*w[i];
+    }
+    tweights.clear();
+}
+
+
+template<class Type>
+void Foam::coupledFvPatchField<Type>::valueBoundaryCoeffs
+(
+    const tmp<scalarField>& tweights,
+    UList<Type>& result
+) const
+{
+    const auto& w = tweights();
+
+    const label len = result.size();
+
+    for (label i = 0; i < len; ++i)
+    {
+        result[i] = Type(pTraits<Type>::one)*(1.0 - w[i]);
+    }
+    tweights.clear();
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::coupledFvPatchField<Type>::valueInternalCoeffs
@@ -174,6 +239,38 @@ Foam::coupledFvPatchField<Type>::valueBoundaryCoeffs
 }
 
 
+template<class Type>
+void Foam::coupledFvPatchField<Type>::gradientInternalCoeffs
+(
+    const scalarField& deltaCoeffs,
+    UList<Type>& result
+) const
+{
+    const label len = result.size();
+
+    for (label i = 0; i < len; ++i)
+    {
+        result[i] = -Type(pTraits<Type>::one)*deltaCoeffs[i];
+    }
+}
+
+
+template<class Type>
+void Foam::coupledFvPatchField<Type>::gradientBoundaryCoeffs
+(
+    const scalarField& deltaCoeffs,
+    UList<Type>& result
+) const
+{
+    const label len = result.size();
+
+    for (label i = 0; i < len; ++i)
+    {
+        result[i] = Type(pTraits<Type>::one)*deltaCoeffs[i];
+    }
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::coupledFvPatchField<Type>::gradientInternalCoeffs
@@ -181,7 +278,9 @@ Foam::coupledFvPatchField<Type>::gradientInternalCoeffs
     const scalarField& deltaCoeffs
 ) const
 {
-    return -Type(pTraits<Type>::one)*deltaCoeffs;
+    auto tresult = tmp<Field<Type>>::New(deltaCoeffs.size());
+    this->gradientInternalCoeffs(deltaCoeffs, tresult.ref());
+    return tresult;
 }
 
 
@@ -194,6 +293,17 @@ Foam::coupledFvPatchField<Type>::gradientInternalCoeffs() const
 }
 
 
+template<class Type>
+void Foam::coupledFvPatchField<Type>::gradientInternalCoeffs
+(
+    UList<Type>& result
+) const
+{
+    NotImplemented;
+    this->gradientInternalCoeffs(this->patch().deltaCoeffs(), result);
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::coupledFvPatchField<Type>::gradientBoundaryCoeffs
@@ -214,6 +324,17 @@ Foam::coupledFvPatchField<Type>::gradientBoundaryCoeffs() const
 }
 
 
+template<class Type>
+void Foam::coupledFvPatchField<Type>::gradientBoundaryCoeffs
+(
+    UList<Type>& result
+) const
+{
+    NotImplemented;
+    this->gradientBoundaryCoeffs(this->patch().deltaCoeffs(), result);
+}
+
+
 template<class Type>
 void Foam::coupledFvPatchField<Type>::write(Ostream& os) const
 {

src/finiteVolume/fields/fvPatchFields/basic/coupled/coupledFvPatchField.H

@@ -208,32 +208,32 @@ public:
             }
 
             //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the value of this patchField with given weights
+            //- evaluation of the value of this patchField with given weights
             virtual tmp<Field<Type>> valueInternalCoeffs
             (
                 const tmp<scalarField>&
             ) const;
 
             //- Return the matrix source coefficients corresponding to the
-            //  evaluation of the value of this patchField with given weights
+            //- evaluation of the value of this patchField with given weights
             virtual tmp<Field<Type>> valueBoundaryCoeffs
             (
                 const tmp<scalarField>&
             ) const;
 
             //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
+            //- evaluation of the gradient of this patchField
             virtual tmp<Field<Type>> gradientInternalCoeffs
             (
                 const scalarField& deltaCoeffs
             ) const;
 
             //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
+            //- evaluation of the gradient of this patchField
             virtual tmp<Field<Type>> gradientInternalCoeffs() const;
 
             //- Return the matrix source coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
+            //- evaluation of the gradient of this patchField
             virtual tmp<Field<Type>> gradientBoundaryCoeffs
             (
                 const scalarField& deltaCoeffs
@@ -244,6 +244,65 @@ public:
             virtual tmp<Field<Type>> gradientBoundaryCoeffs() const;
 
 
+      // Contiguous storage
+
+            //- Retrieve patch-normal gradient.
+            //- Assumes faceCells is an indirection into internal field
+            virtual void snGrad
+            (
+                const scalarField& deltaCoeffs,
+                UList<Type>&
+            ) const;
+
+            //- Retrieve patch-normal gradient
+            virtual void snGrad(UList<Type>&) const
+            {
+                NotImplemented;
+            }
+
+            //- Return the matrix diagonal coefficients corresponding to the
+            //  evaluation of the value of this patchField with given weights
+            virtual void valueInternalCoeffs
+            (
+                const tmp<scalarField>&,
+                UList<Type>&
+            ) const;
+
+            //- Return the matrix source coefficients corresponding to the
+            //  evaluation of the value of this patchField with given weights
+            virtual void valueBoundaryCoeffs
+            (
+                const tmp<scalarField>&,
+                UList<Type>&
+            ) const;
+
+            //- Return the matrix diagonal coefficients corresponding to the
+            //  evaluation of the gradient of this patchField
+            virtual void gradientInternalCoeffs(UList<Type>&) const;
+
+            //- Return the matrix diagonal coefficients corresponding to the
+            //  evaluation of the gradient of this coupled patchField
+            //  using the deltaCoeffs provided
+            virtual void gradientInternalCoeffs
+            (
+                const scalarField& deltaCoeffs,
+                UList<Type>&
+            ) const;
+
+            //- Return the matrix source coefficients corresponding to the
+            //  evaluation of the gradient of this patchField
+            virtual void gradientBoundaryCoeffs(UList<Type>&) const;
+
+            //- Return the matrix source coefficients corresponding to the
+            //  evaluation of the gradient of this coupled patchField
+            //  using the deltaCoeffs provided
+            virtual void gradientBoundaryCoeffs
+            (
+                const scalarField& deltaCoeffs,
+                UList<Type>&
+            ) const;
+
+
         // Coupled interface functionality
 
             //- Update result field based on interface functionality

src/finiteVolume/fields/fvPatchFields/basic/zeroGradient/zeroGradientFvPatchField.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -128,14 +129,22 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
-        // Evaluation functions
+        // Evaluation Functions
 
-            //- Return gradient at boundary
+            //- Return patch-normal gradient
             virtual tmp<Field<Type>> snGrad() const
             {
-                return tmp<Field<Type>>::New(this->size(), Zero);
+                // zero-gradient
+                return tmp<Field<Type>>::New(this->size(), Foam::zero{});
+            }
+
+            //- Retrieve patch-normal gradient [contiguous storage]
+            virtual void snGrad(UList<Type>& result) const
+            {
+                // zero-gradient
+                result = Foam::zero{};
             }
 
             //- Evaluate the patch field
@@ -146,25 +155,25 @@ public:
             );
 
             //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the value of this patchField with given weights
+            //- evaluation of the value of this patchField with given weights
             virtual tmp<Field<Type>> valueInternalCoeffs
             (
                 const tmp<scalarField>&
             ) const;
 
             //- Return the matrix source coefficients corresponding to the
-            //  evaluation of the value of this patchField with given weights
+            //- evaluation of the value of this patchField with given weights
             virtual tmp<Field<Type>> valueBoundaryCoeffs
             (
                 const tmp<scalarField>&
             ) const;
 
             //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
+            //- evaluation of the gradient of this patchField
             virtual tmp<Field<Type>> gradientInternalCoeffs() const;
 
             //- Return the matrix source coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
+            //- evaluation of the gradient of this patchField
             virtual tmp<Field<Type>> gradientBoundaryCoeffs() const;
 };
 

src/finiteVolume/fields/fvPatchFields/constraint/empty/emptyFvPatchField.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -133,7 +134,7 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         // Mapping functions
 
@@ -153,7 +154,7 @@ public:
             {}
 
 
-        // Evaluation functions
+        // Evaluation Functions
 
             //- Update the coefficients associated with the patch field
             //  This only checks to see the case is actually 1D or 2D
@@ -182,20 +183,31 @@ public:
             }
 
             //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
+            //- evaluation of the gradient of this patchField
             tmp<Field<Type>> gradientInternalCoeffs() const
             {
                 return tmp<Field<Type>>::New();
             }
 
             //- Return the matrix source coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
+            //- evaluation of the gradient of this patchField
             tmp<Field<Type>> gradientBoundaryCoeffs() const
             {
                 return tmp<Field<Type>>::New();
             }
 
 
+    // Contiguous storage
+
+        //- Retrieve patch-normal gradient [contiguous storage].
+        //- Placeholder value is zero (treated like zero-gradient).
+        virtual void snGrad(UList<Type>& result) const
+        {
+            // Treat like zero-gradient
+            result = Foam::zero{};
+        }
+
+
     // Member Functions
 
         //- Write without "value" entry!

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.C

@@ -314,6 +314,28 @@ void Foam::processorFvPatchField<Type>::evaluate
 }
 
 
+template<class Type>
+void Foam::processorFvPatchField<Type>::snGrad
+(
+    const scalarField& deltaCoeffs,
+    UList<Type>& result
+) const
+{
+    // Get patch internal field, store temporarily in result
+    this->patchInternalField(result);
+    const auto& pif = result;
+
+    const auto& pnf = *this;
+
+    const label len = result.size();
+
+    for (label i = 0; i < len; ++i)
+    {
+        result[i] = deltaCoeffs[i]*(pnf[i] - pif[i]);
+    }
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::processorFvPatchField<Type>::snGrad
@@ -321,6 +343,10 @@ Foam::processorFvPatchField<Type>::snGrad
     const scalarField& deltaCoeffs
 ) const
 {
+    // OR
+    // auto tfld = tmp<Field<Type>>::New(this->size());
+    // this->snGrad(deltaCoeffs, tfld.ref());
+    // return tfld;
     return deltaCoeffs*(*this - this->patchInternalField());
 }
 

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H

@@ -221,6 +221,13 @@ public:
             const scalarField& deltaCoeffs
         ) const;
 
+        //- Retrieve patch-normal gradient
+        virtual void snGrad
+        (
+            const scalarField& deltaCoeffs,
+            UList<Type>& result
+        ) const;
+
 
         // Coupled interface functionality
 

src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C

@@ -221,10 +221,31 @@ void Foam::fvPatchField<Type>::check(const fvPatchField<Type>& rhs) const
 }
 
 
+template<class Type>
+void Foam::fvPatchField<Type>::snGrad(UList<Type>& result) const
+{
+    // Get patch internal field, store temporarily in result
+    this->patchInternalField(result);
+    const auto& pif = result;
+
+    const Field<Type>& pfld = *this;
+    const auto& dc = patch().deltaCoeffs();
+
+    const label len = result.size();
+
+    for (label i = 0; i < len; ++i)
+    {
+        result[i] = dc[i]*(pfld[i] - pif[i]);
+    }
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>> Foam::fvPatchField<Type>::snGrad() const
 {
-    return patch().deltaCoeffs()*(*this - patchInternalField());
+    auto tfld = tmp<Field<Type>>::New(this->size());
+    this->snGrad(static_cast<UList<Type>&>(tfld.ref()));
+    return tfld;
 }
 
 

src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H

@@ -622,8 +622,11 @@ public:
             //- Return patch-normal gradient
             virtual tmp<Field<Type>> snGrad() const;
 
+            //- Retrieve patch-normal gradient [contiguous storage]
+            virtual void snGrad(UList<Type>& result) const;
+
             //- Return patch-normal gradient for coupled-patches
-            //  using the deltaCoeffs provided
+            //- using the deltaCoeffs provided
             virtual tmp<Field<Type>> snGrad
             (
                 const scalarField& deltaCoeffs
@@ -633,15 +636,16 @@ public:
                 return *this;
             }
 
-            //- Update the coefficients associated with the patch field
+            //- Retrieve patch-normal gradient for coupled-patches
-            //  Sets Updated to true
+            //- using the deltaCoeffs provided [contiguous storage]
-            virtual void updateCoeffs();
+            virtual void snGrad
-
+            (
-            //- Update the coefficients associated with the patch field
+                const scalarField& deltaCoeffs,
-            //  with a weight field (0..1). This weight field is usually
+                UList<Type>&
-            //  provided as the amount of geometric overlap for 'duplicate'
+            ) const
-            //  patches. Sets Updated to true
+            {
-            virtual void updateWeightedCoeffs(const scalarField& weights);
+                NotImplemented;
+            }
 
             //- Return internal field next to patch
             virtual tmp<Field<Type>> patchInternalField() const;
@@ -657,6 +661,22 @@ public:
                 return *this;
             }
 
+            //- Retrieve patchField on the opposite patch of a coupled patch
+            virtual void patchNeighbourField(UList<Type>&) const
+            {
+                NotImplemented;
+            }
+
+            //- Update the coefficients associated with the patch field
+            //  Sets Updated to true
+            virtual void updateCoeffs();
+
+            //- Update the coefficients associated with the patch field
+            //  with a weight field (0..1). This weight field is usually
+            //  provided as the amount of geometric overlap for 'duplicate'
+            //  patches. Sets Updated to true
+            virtual void updateWeightedCoeffs(const scalarField& weights);
+
             //- Initialise the evaluation of the patch field
             virtual void initEvaluate
             (
@@ -690,7 +710,7 @@ public:
             {}
 
             //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the value of this patchField with given weights
+            //- evaluation of the value of this patchField with given weights
             virtual tmp<Field<Type>> valueInternalCoeffs
             (
                 const tmp<scalarField>&
@@ -701,7 +721,7 @@ public:
             }
 
             //- Return the matrix source coefficients corresponding to the
-            //  evaluation of the value of this patchField with given weights
+            //- evaluation of the value of this patchField with given weights
             virtual tmp<Field<Type>> valueBoundaryCoeffs
             (
                 const tmp<scalarField>&
@@ -712,7 +732,7 @@ public:
             }
 
             //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
+            //- evaluation of the gradient of this patchField
             virtual tmp<Field<Type>> gradientInternalCoeffs() const
             {
                 NotImplemented;
@@ -720,8 +740,8 @@ public:
             }
 
             //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the gradient of this coupled patchField
+            //- evaluation of the gradient of this coupled patchField
-            //  using the deltaCoeffs provided
+            //- using the deltaCoeffs provided
             virtual tmp<Field<Type>> gradientInternalCoeffs
             (
                 const scalarField& deltaCoeffs
@@ -732,7 +752,7 @@ public:
             }
 
             //- Return the matrix source coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
+            //- evaluation of the gradient of this patchField
             virtual tmp<Field<Type>> gradientBoundaryCoeffs() const
             {
                 NotImplemented;
@@ -740,8 +760,8 @@ public:
             }
 
             //- Return the matrix source coefficients corresponding to the
-            //  evaluation of the gradient of this coupled patchField
+            //- evaluation of the gradient of this coupled patchField
-            //  using the deltaCoeffs provided
+            //- using the deltaCoeffs provided
             virtual tmp<Field<Type>> gradientBoundaryCoeffs
             (
                 const scalarField& deltaCoeffs
@@ -751,7 +771,6 @@ public:
                 return *this;
             }
 
-
             //- Manipulate matrix
             virtual void manipulateMatrix(fvMatrix<Type>& matrix);
 
@@ -771,6 +790,69 @@ public:
             );
 
 
+        // Contiguous storage
+
+            //- Retrieve the matrix diagonal coefficients corresponding to the
+            //- evaluation of the value of this patchField with given weights
+            virtual void valueInternalCoeffs
+            (
+                const tmp<Field<scalar>>&,
+                UList<Type>&
+            ) const
+            {
+                NotImplemented;
+            }
+
+            //- Retrieve the matrix source coefficients corresponding to the
+            //- evaluation of the value of this patchField with given weights
+            virtual void valueBoundaryCoeffs
+            (
+                const tmp<Field<scalar>>&,
+                UList<Type>&
+            ) const
+            {
+                NotImplemented;
+            }
+
+            //- Retrieve the matrix diagonal coefficients corresponding to the
+            //- evaluation of the gradient of this patchField
+            virtual void gradientInternalCoeffs(UList<Type>&) const
+            {
+                NotImplemented;
+            }
+
+            //- Retrieve the matrix diagonal coefficients corresponding to the
+            //- evaluation of the gradient of this coupled patchField
+            //- using the deltaCoeffs provided
+            virtual void gradientInternalCoeffs
+            (
+                const scalarField& deltaCoeffs,
+                UList<Type>&
+            ) const
+            {
+                NotImplemented;
+            }
+
+            //- Retrieve the matrix source coefficients corresponding to the
+            //- evaluation of the gradient of this patchField
+            virtual void gradientBoundaryCoeffs(UList<Type>&) const
+            {
+                NotImplemented;
+            }
+
+            //- Retrieve the matrix source coefficients corresponding to the
+            //- evaluation of the gradient of this coupled patchField
+            //- using the deltaCoeffs provided
+            virtual void gradientBoundaryCoeffs
+            (
+                const scalarField& deltaCoeffs,
+                UList<Type>&
+            ) const
+            {
+                NotImplemented;
+            }
+
+
     // Other
 
         //- Write