ENH: Friedrich: refactor model to use cornerDetectionModel framework

- replace inline corner detection logic with cornerDetectionModel integration

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/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/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/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/Make/files

@@ -25,11 +25,20 @@ $(kinematic)/injectionModel/injectionModelList/injectionModelList.C
 $(kinematic)/injectionModel/injectionModel/injectionModel.C
 $(kinematic)/injectionModel/injectionModel/injectionModelNew.C
 
-$(kinematic)/injectionModel/filmSeparation/filmSeparation.C
-$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModel.C
-$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModelNew.C
-$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/OwenRyleyModel/OwenRyleyModel.C
-$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.C
+/* Film separation models */
+
+filmSeparation=$(kinematic)/injectionModel/filmSeparation
+filmSeparationModels=$(filmSeparation)/filmSeparationModels
+$(filmSeparation)/filmSeparation.C
+$(filmSeparationModels)/filmSeparationModel/filmSeparationModel.C
+$(filmSeparationModels)/filmSeparationModel/filmSeparationModelNew.C
+$(filmSeparationModels)/OwenRyleyModel/OwenRyleyModel.C
+$(filmSeparationModels)/FriedrichModel/FriedrichModel.C
+cornerDetectionModels=$(filmSeparation)/cornerDetectionModels
+$(cornerDetectionModels)/cornerDetectionModel/cornerDetectionModel.C
+$(cornerDetectionModels)/cornerDetectionModel/cornerDetectionModelNew.C
+$(cornerDetectionModels)/fluxBased/fluxBased.C
+$(cornerDetectionModels)/geometryBased/geometryBased.C
 
 $(kinematic)/injectionModel/BrunDrippingInjection/BrunDrippingInjection.C
 

src/regionFaModels/Make/options

@@ -11,7 +11,8 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/transportModels
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/fileFormats/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
@@ -24,4 +25,5 @@ LIB_LIBS = \
     -lthermophysicalProperties \
     -lspecie \
     -lfaOptions \
-    -ldistributionModels
+    -ldistributionModels \
+    -lfileFormats

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

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/cornerDetectionModels/cornerDetectionModel/cornerDetectionModel.C

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "cornerDetectionModel.H"
+#include "faMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(cornerDetectionModel, 0);
+    defineRunTimeSelectionTable(cornerDetectionModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::cornerDetectionModel::dihedralAngle
+(
+    const vector& n0,
+    const vector& n1
+) const
+{
+    if (mag(n0) <= VSMALL || mag(n1) <= VSMALL)
+    {
+#ifdef FULL_DEBUG
+        WarningInFunction
+            << "Degenerate face normal magnitude (|n| ~ 0). "
+            << "Returning 0 for dihedral angle." << nl;
+#endif
+        return 0;
+    }
+
+    const scalar a = mag(n1 - n0);
+    const scalar b = mag(n1 + n0);
+
+    // The dihedral angle is calculated as 2*atan2(|n1 - n0|, |n1 + n0|),
+    // which gives the angle between the two normals n0 and n1.
+    scalar phi = scalar(2)*std::atan2(a, b);
+
+    // Clamp to [0, pi]
+    phi = max(0, min(constant::mathematical::pi, phi));
+
+    if (!std::isfinite(phi))
+    {
+#ifdef FULL_DEBUG
+        WarningInFunction
+            << "Non-finite dihedral angle computed. Returning 0." << nl;
+#endif
+        return 0;
+    }
+
+    return phi;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cornerDetectionModel::cornerDetectionModel
+(
+    const faMesh& mesh,
+    const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
+    const dictionary& dict
+)
+:
+    mesh_(mesh),
+    film_(film),
+    dict_(dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::cornerDetectionModel::~cornerDetectionModel()
+{}  // faMesh was forward declared
+
+
+// ************************************************************************* //

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/cornerDetectionModels/cornerDetectionModel/cornerDetectionModel.H

@@ -0,0 +1,211 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Namespace
+    Foam::cornerDetectionModels
+
+Description
+    A namespace for various \c cornerDetection model implementations.
+
+Class
+    Foam::cornerDetectionModel
+
+Description
+    A base class for \c cornerDetection models.
+
+SourceFiles
+    cornerDetectionModel.C
+    cornerDetectionModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_cornerDetectionModel_H
+#define Foam_cornerDetectionModel_H
+
+#include "liquidFilmBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class faMesh;
+class dictionary;
+
+/*---------------------------------------------------------------------------*\
+                    Class cornerDetectionModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class cornerDetectionModel
+{
+    // Private Data
+
+        //- Const reference to the finite-area mesh
+        const faMesh& mesh_;
+
+        //- Const reference to the film
+        const regionModels::areaSurfaceFilmModels::liquidFilmBase& film_;
+
+        //- Const reference to the dictionary
+        const dictionary& dict_;
+
+        //- Bitset for corner faces, true for corner faces
+        bitSet cornerFaces_;
+
+        //- List of corner angles [rad]
+        scalarList cornerAngles_;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Return the dihedral angle between two neighbour faces
+        //  Robust 2*atan form (handles parallel normals better than acos)
+        //  \param n0 First normal vector
+        //  \param n1 Second normal vector
+        //  \return The dihedral angle [rad]
+        scalar dihedralAngle(const vector& n0, const vector& n1) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("cornerDetectionModel");
+
+
+    // Declare runtime constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            cornerDetectionModel,
+            dictionary,
+            (
+                const faMesh& mesh,
+                const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
+                const dictionary& dict
+            ),
+            (mesh, film, dict)
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected cornerDetection model
+        static autoPtr<cornerDetectionModel> New
+        (
+            const faMesh& mesh,
+            const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
+            const dictionary& dict
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        cornerDetectionModel
+        (
+            const faMesh& mesh,
+            const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
+            const dictionary& dict
+        );
+
+
+    // Generated Methods
+
+        //- No copy construct
+        cornerDetectionModel(const cornerDetectionModel&) = delete;
+
+        //- No copy assignment
+        void operator=(const cornerDetectionModel&) = delete;
+
+
+    //- Destructor
+    virtual ~cornerDetectionModel();
+
+
+    // Member Functions
+
+    // Getters
+
+        //- Return const reference to the finite-area mesh
+        const faMesh& mesh() const noexcept { return mesh_; }
+
+        //- Return const reference to the film
+        const regionModels::areaSurfaceFilmModels::liquidFilmBase&
+            film() const noexcept { return film_; }
+
+        //- Return const reference to the dictionary
+        const dictionary& dict() const noexcept { return dict_; }
+
+        //- Return const reference to the corner faces bitset
+        const bitSet& getCornerFaces() const noexcept { return cornerFaces_; }
+
+        //- Return const reference to the corner angles list [rad]
+        const scalarList& getCornerAngles() const noexcept
+        {
+            return cornerAngles_;
+        }
+
+
+    // Setters
+
+        //- Set the corner faces bitset
+        void setCornerFaces(const bitSet& cornerFaces)
+        {
+            cornerFaces_ = cornerFaces;
+        }
+
+        //- Set the corner angles list [rad]
+        void setCornerAngles(const scalarList& cornerAngles)
+        {
+            cornerAngles_ = cornerAngles;
+        }
+
+
+    // Evaluation
+
+        //- Detect and store the corner faces
+        virtual bool detectCorners() = 0;
+
+
+    // I-O
+
+        //- Read the model dictionary
+        virtual bool read(const dictionary&) { return true; }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/cornerDetectionModels/cornerDetectionModel/cornerDetectionModelNew.C

@@ -25,38 +25,40 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "cornerDetectionModel.H"
+#include "faMesh.H"
+#include "liquidFilmBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class Type, class EvalFunction>
-bool Foam::AMICache::apply(List<Type>& result, const EvalFunction& eval) const
+Foam::autoPtr<Foam::cornerDetectionModel> Foam::cornerDetectionModel::New
+(
+    const faMesh& mesh,
+    const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
+    const dictionary& dict
+)
 {
-    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_);
+    const word modelType
+    (
+        dict.getOrDefault<word>("cornerDetectionModel", "fluxBased")
+    );
 
-        //result = (r1 - r0)*interpWeight_ + r0;
-        forAll(result, i)
-        {
-            result[i] = lerp(r0[i], r1[i], interpWeight_);
-        }
-        return true;
+    Info<< "    Selecting corner-detection model: " << modelType << nl << endl;
+
+    auto* ctorPtr = dictionaryConstructorTable(modelType);
+
+    if (!ctorPtr)
+    {
+        FatalIOErrorInLookup
+        (
+            dict,
+            "cornerDetectionModel",
+            modelType,
+            *dictionaryConstructorTablePtr_
+        ) << exit(FatalIOError);
     }
 
-    return false;
+    return autoPtr<cornerDetectionModel>(ctorPtr(mesh, film, dict));
 }
 
 

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/cornerDetectionModels/fluxBased/fluxBased.C

@@ -0,0 +1,391 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "fluxBased.H"
+#include "OBJstream.H"
+#include "processorFaPatch.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace cornerDetectionModels
+{
+    defineTypeNameAndDebug(fluxBased, 0);
+    addToRunTimeSelectionTable(cornerDetectionModel, fluxBased, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::bitSet Foam::cornerDetectionModels::fluxBased::identifyCornerEdges() const
+{
+    // Return true if face normals converge, i.e. sharp edge
+    // Face-normal vectors diverge: no separation, converge: separation (maybe)
+    const auto isCornerEdgeSharp =
+    [](
+        const vector& fcO,  // face-centre owner
+        const vector& fcN,  // face-centre neigh
+        const vector& fnO,  // face-normal owner
+        const vector& fnN   // face-normal neigh
+    ) noexcept -> bool
+    {
+        // Threshold for sharpness detection
+        constexpr scalar sharpEdgeThreshold = -1e-8;
+
+        // Relative centre and normal of the two faces sharing the edge
+        const vector relativePosition(fcN - fcO);
+        const vector relativeNormal(fnN - fnO);
+
+        // Sharp if normals converge along the centre-to-centre direction
+        return ((relativeNormal & relativePosition) < sharpEdgeThreshold);
+    };
+
+
+    // Cache the operand references
+    const areaVectorField& fc = mesh().areaCentres();
+    const areaVectorField& fn = mesh().faceAreaNormals();
+    const labelUList& own = mesh().edgeOwner();
+    const labelUList& nei = mesh().edgeNeighbour();
+
+
+    // Allocate the resource for the return object
+    bitSet cornerEdges(mesh().nEdges(), false);
+
+    // Internal edges (owner <-> neighbour)
+    forAll(nei, edgei)
+    {
+        const label faceO = own[edgei];
+        const label faceN = nei[edgei];
+
+        cornerEdges[edgei] = isCornerEdgeSharp
+        (
+            fc[faceO],
+            fc[faceN],
+            fn[faceO],
+            fn[faceN]
+        );
+    }
+
+
+    // Skip the rest of the routine if the simulation is a serial run
+    if (!Pstream::parRun()) return cornerEdges;
+
+    // Check if processor face-normal vectors diverge (no separation)
+    // or converge (separation may occur)
+    const faBoundaryMesh& patches = mesh().boundary();
+
+    for (const faPatch& fap : patches)
+    {
+        if (!isA<processorFaPatch>(fap)) continue;
+
+        const label patchi = fap.index();
+        const auto& edgeFaces = fap.edgeFaces();
+        const label internalEdge0 = fap.start();
+
+        const auto& fcp = fc.boundaryField()[patchi];
+        const auto& fnp = fn.boundaryField()[patchi];
+
+        // Processor edges (owner <-| none)
+        forAll(fnp, bndEdgei)
+        {
+            const label faceO = edgeFaces[bndEdgei];
+            const label meshEdgei = internalEdge0 + bndEdgei;
+
+            cornerEdges[meshEdgei] = isCornerEdgeSharp
+            (
+                fc[faceO],
+                fcp[bndEdgei],
+                fn[faceO],
+                fnp[bndEdgei]
+            );
+        }
+    }
+
+    return cornerEdges;
+}
+
+
+Foam::bitSet Foam::cornerDetectionModels::fluxBased::identifyCornerFaces
+(
+    const bitSet& cornerEdges
+) const
+{
+    // Marks the separating face based on edge flux sign
+    const auto markSeparation =
+    [](
+        bitSet& cornerFaces,
+        const scalar phiEdge,
+        const label faceO,
+        const label faceN = -1  /* = -1 for processor edges */
+    ) noexcept -> void
+    {
+        constexpr scalar tol = 1e-8;
+
+        // Assuming no sources/sinks at the edge
+        if (phiEdge > tol)  // From owner to neighbour
+        {
+            cornerFaces[faceO] = true;
+        }
+        else if ((phiEdge < -tol) && (faceN != -1))  // From nei to own
+        {
+            cornerFaces[faceN] = true;
+        }
+    };
+
+
+    // Cache the operand references
+    const edgeScalarField& phis = film().phi2s();
+    const labelUList& own = mesh().edgeOwner();
+    const labelUList& nei = mesh().edgeNeighbour();
+
+    // Allocate the resource for the return object
+    bitSet cornerFaces(mesh().faces().size(), false);
+
+    // Internal faces (owner <-> neighbour)
+    forAll(nei, edgei)
+    {
+        if (!cornerEdges[edgei]) continue;
+
+        markSeparation
+        (
+            cornerFaces,
+            phis[edgei],
+            own[edgei],  // faceO
+            nei[edgei]   // faceN
+        );
+    }
+
+
+    // Skip the rest of the routine if the simulation is a serial run
+    if (!Pstream::parRun()) return cornerFaces;
+
+    const faBoundaryMesh& patches = mesh().boundary();
+
+    for (const faPatch& fap : patches)
+    {
+        if (!isA<processorFaPatch>(fap)) continue;
+
+        const label patchi = fap.index();
+        const auto& edgeFaces = fap.edgeFaces();
+        const label internalEdge0 = fap.start();
+
+        const auto& phisp = phis.boundaryField()[patchi];
+
+        // Processor faces (owner <-| none)
+        forAll(phisp, bndEdgei)
+        {
+            const label faceO = edgeFaces[bndEdgei];
+            const label meshEdgei = internalEdge0 + bndEdgei;
+
+            if (!cornerEdges[meshEdgei]) continue;
+
+            markSeparation
+            (
+                cornerFaces,
+                phisp[bndEdgei],
+                faceO
+                /*faceN = -1*/
+            );
+        }
+    }
+
+    return cornerFaces;
+}
+
+
+Foam::scalarList Foam::cornerDetectionModels::fluxBased::calcCornerAngles
+(
+    const bitSet& faces,
+    const bitSet& edges
+) const
+{
+    // Cache the operand references
+    const areaVectorField& fn = mesh().faceAreaNormals();
+    const labelUList& own = mesh().edgeOwner();
+    const labelUList& nei = mesh().edgeNeighbour();
+
+    scalarList cornerFaceAngles(mesh().faces().size(), Zero);
+
+    // Internal edges (owner <-> neighbour)
+    forAll(nei, edgei)
+    {
+        if (!edges[edgei]) continue;
+
+        const label faceO = own[edgei];
+        const label faceN = nei[edgei];
+
+        // If neither adjacent face is flagged as a corner, skip the atan2 work
+        if (!faces[faceO] && !faces[faceN]) continue;
+
+        const scalar ang = this->dihedralAngle(fn[faceO], fn[faceN]);
+
+        if (faces[faceO]) cornerFaceAngles[faceO] = ang;
+        if (faces[faceN]) cornerFaceAngles[faceN] = ang;
+    }
+
+
+    // Skip the rest of the routine if the simulation is a serial run
+    if (!Pstream::parRun()) return cornerFaceAngles;
+
+    const faBoundaryMesh& patches = mesh().boundary();
+
+    for (const faPatch& fap : patches)
+    {
+        if (!isA<processorFaPatch>(fap)) continue;
+        const label patchi = fap.index();
+        const auto& edgeFaces = fap.edgeFaces();
+        const label internalEdge0 = fap.start();
+
+        const auto& fnp = fn.boundaryField()[patchi];
+
+        // Processor edges (owner <-| none)
+        forAll(fnp, bndEdgei)
+        {
+            const label faceO = edgeFaces[bndEdgei];
+            const label meshEdgei = internalEdge0 + bndEdgei;
+
+                // Only if the mesh edge and owner face are both corners
+                if (!edges[meshEdgei] || !faces[faceO]) continue;
+
+                cornerFaceAngles[faceO] =
+                    this->dihedralAngle(fn[faceO], fnp[bndEdgei]);
+        }
+    }
+
+    return cornerFaceAngles;
+}
+
+
+void Foam::cornerDetectionModels::fluxBased::writeEdgesAndFaces
+(
+    const word& prefix
+) const
+{
+    const pointField& pts = mesh().points();
+
+    const word timeName(Foam::name(mesh().time().value()));
+    const word nameEdges("fluxBased-edges-" + timeName + ".obj");
+    const word nameFaces("fluxBased-faces-" + timeName + ".obj");
+
+
+    // Write OBJ of edge faces to file
+    OBJstream osEdges(mesh().time().path()/nameEdges);
+
+    const auto& edges = mesh().edges();
+    forAll(cornerEdges_, ei)
+    {
+        if (cornerEdges_[ei])
+        {
+            const edge& e = edges[ei];
+            osEdges.write(e, pts);
+        }
+    }
+
+    // Write OBJ of corner faces to file
+    OBJstream osFaces(mesh().time().path()/nameFaces);
+
+    const bitSet& cornerFaces = this->getCornerFaces();
+    const auto& faces = mesh().faces();
+    forAll(cornerFaces, fi)
+    {
+        if (cornerFaces[fi])
+        {
+            const face& f = faces[fi];
+            osFaces.write(f, pts);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cornerDetectionModels::fluxBased::fluxBased
+(
+    const faMesh& mesh,
+    const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
+    const dictionary& dict
+)
+:
+    cornerDetectionModel(mesh, film, dict),
+    init_(false)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::cornerDetectionModels::fluxBased::detectCorners()
+{
+    if (!init_ || mesh().moving())
+    {
+        // Identify and store corner edges based on face normals
+        cornerEdges_ = identifyCornerEdges();
+        init_ = true;
+    }
+
+
+    // Identify and store corner faces based on edge flux sign
+    this->setCornerFaces(identifyCornerFaces(cornerEdges_));
+
+
+    // Calculate and store corner face angles
+    const bitSet& cornerFaces = this->getCornerFaces();
+    this->setCornerAngles
+    (
+        calcCornerAngles(cornerFaces, cornerEdges_)
+    );
+
+
+    // Write edges and faces as OBJ sets for debug purposes, if need be
+    if (debug && mesh().time().writeTime())
+    {
+        writeEdgesAndFaces();
+    }
+
+    return true;
+}
+
+
+bool Foam::cornerDetectionModels::fluxBased::read(const dictionary& dict)
+{
+    if (!cornerDetectionModel::read(dict))
+    {
+        return false;
+    }
+
+    // Force the re-identification of corner edges/faces
+    init_ = false;
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/cornerDetectionModels/fluxBased/fluxBased.H

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::cornerDetectionModels::fluxBased
+
+Description
+    Flux-based corner detection model. Marks faces at which the liquid film can
+    separate based on the flux.
+
+    The model identifies sharp edges based on the face normals of the two
+    faces sharing the edge. Then, if the edge is sharp, the flux direction is
+    evaluated to mark the face through which the flux leaves the liquid film.
+    If the edge is sharp and the flux leaves through one of the two faces
+    sharing the edge, the face is marked as a corner face, where the film can
+    separate.
+
+Usage
+    Minimal example in boundary-condition files:
+    \verbatim
+    filmSeparationCoeffs
+    {
+        // Inherited entries
+        ...
+
+        // Optional entries
+        cornerDetectionModel  fluxBased;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                     | Type | Reqd | Deflt
+      cornerDetectionModel | Corner detector model   | word | no   | fluxBased
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link cornerDetectionModel.H \endlink
+
+SourceFiles
+    fluxBased.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_cornerDetectionModels_fluxBased_H
+#define Foam_cornerDetectionModels_fluxBased_H
+
+#include "cornerDetectionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace cornerDetectionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class fluxBased Declaration
+\*---------------------------------------------------------------------------*/
+
+class fluxBased
+:
+    public cornerDetectionModel
+{
+    // Private Data
+
+        //- Flag to deduce if the object is initialised
+        bool init_;
+
+        //- Identified corner edges
+        bitSet cornerEdges_;
+
+
+    // Private Member Functions
+
+        //- Return Boolean list of identified corner edges
+        bitSet identifyCornerEdges() const;
+
+        //- Return Boolean list of identified corner faces
+        bitSet identifyCornerFaces(const bitSet& cornerEdges) const;
+
+        //- Return the list of corner angles for each edge [rad]
+        scalarList calcCornerAngles
+        (
+            const bitSet& faces,
+            const bitSet& edges
+        ) const;
+
+        // Write edges and faces as OBJ sets for debug purposes
+        void writeEdgesAndFaces(const word& prefix = "geomCorners") const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fluxBased");
+
+
+    // Constructors
+
+        //- Construct from components
+        fluxBased
+        (
+            const faMesh& mesh,
+            const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+    virtual ~fluxBased() = default;
+
+
+    // Member Functions
+
+    // Evaluation
+
+        //- Detect and store the corner faces
+        virtual bool detectCorners();
+
+
+    // I-O
+
+        //- Read the model dictionary
+        virtual bool read(const dictionary&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace cornerDetectionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/cornerDetectionModels/geometryBased/geometryBased.C

@@ -0,0 +1,586 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "geometryBased.H"
+#include "processorFaPatch.H"
+#include "unitConversion.H"
+#include "syncTools.H"
+#include "OBJstream.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace cornerDetectionModels
+{
+    defineTypeNameAndDebug(geometryBased, 0);
+    addToRunTimeSelectionTable(cornerDetectionModel, geometryBased, dictionary);
+}
+}
+
+const Foam::Enum
+<
+    Foam::cornerDetectionModels::geometryBased::cornerCurveType
+>
+Foam::cornerDetectionModels::geometryBased::cornerCurveTypeNames
+({
+    { cornerCurveType::ANY, "any" },
+    { cornerCurveType::CONCAVE , "concave" },
+    { cornerCurveType::CONVEX , "convex" }
+});
+
+const Foam::Enum
+<
+    Foam::cornerDetectionModels::geometryBased::cornerType
+>
+Foam::cornerDetectionModels::geometryBased::cornerTypeNames
+({
+    { cornerType::ALL, "sharpOrRound" },
+    { cornerType::SHARP , "sharp" },
+    { cornerType::ROUND , "round" }
+});
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::cornerDetectionModels::geometryBased::curvatureSign
+(
+    const vector& t,
+    const vector& n0,
+    const vector& n1
+) const
+{
+    // t: unit edge tangent
+    // n0: owner face unit normal
+    // n1: neighbour face unit normal
+    scalar curvature = (t & (n0 ^ n1));
+
+    // Orientation: sign of triple product t . (n0 x n1)
+    // Positive => one sense (together with outward normals, treat as "convex");
+    // mapping to convex/concave is finally gated by 'cornerCurveType_'.
+    return sign(curvature);
+}
+
+
+void Foam::cornerDetectionModels::geometryBased::classifyEdges
+(
+    bitSet& sharpEdges,
+    bitSet& roundEdges
+) const
+{
+    // Cache the operand references
+    const areaVectorField& nf = mesh().faceAreaNormals();
+    const edgeList& edges = mesh().edges();
+    const labelUList& own = mesh().edgeOwner();      // own.sz = nEdges
+    const labelUList& nei = mesh().edgeNeighbour();  // nei.sz = nInternalEdges
+    const pointField& pts = mesh().points();
+
+
+    // Convert input-angle parameters from degrees to radians
+    const scalar angSharp = degToRad(angleSharpDeg_);
+    const scalar angRoundMin = degToRad(angleRoundMinDeg_);
+    const scalar angRoundMax = degToRad(angleRoundMaxDeg_);
+
+
+    // Limit to subset of patches if requested
+    bitSet allowedFaces(mesh().nFaces(), true);
+
+
+    // Allocate the resource for the return objects
+    sharpEdges.resize(mesh().nEdges());  // internal + boundary edges
+    sharpEdges.reset();
+
+    roundEdges.resize(mesh().nEdges());
+    roundEdges.reset();
+
+
+    // Internal edges
+    const label nInternalEdges = mesh().nInternalEdges();
+    for (label ei = 0; ei < nInternalEdges; ++ei)
+    {
+        // Do not allow processing of edges shorter than 'minEdgeLength'
+        const edge& e = edges[ei];
+        const scalar le = e.mag(pts);
+        if (le <= max(minEdgeLength_, VSMALL)) continue;
+
+
+        // Do not allow processing of excluded faces
+        const label f0 = own[ei];
+        const label f1 = nei[ei];
+        if (!allowedFaces.test(f0) && !allowedFaces.test(f1)) continue;
+
+
+        // Calculate the dihedral angle and curvature per edge
+        const vector& n0 = nf[f0];
+        const vector& n1 = nf[f1];
+
+        const scalar phi = this->dihedralAngle(n0, n1); // [rad]
+        const scalar kappa = 2.0*Foam::sin(0.5*phi)/max(le, VSMALL); // [1/m]
+        const scalar R = (kappa > VSMALL ? scalar(1)/kappa : GREAT);
+
+        const vector tangent(e.unitVec(pts));
+
+        const scalar sgn = curvatureSign(tangent, n0, n1);
+
+        const bool curvatureType =
+            (cornerCurveType_ == cornerCurveType::ANY)
+         || (cornerCurveType_ == cornerCurveType::CONVEX  && sgn > 0)
+         || (cornerCurveType_ == cornerCurveType::CONCAVE && sgn < 0);
+
+        // Do not allow processing of excluded curvature-type faces
+        if (!curvatureType) continue;
+
+
+        // Sharp: dihedral above threshold
+        if (phi >= angSharp && cornerType_ != cornerType::ROUND)
+        {
+            sharpEdges.set(ei);
+            continue; // do not double-classify as round
+        }
+
+
+        // Round: small-to-moderate angle but small radius (tight fillet)
+        if
+        (
+            phi >= angRoundMin && phi <= angRoundMax
+         && R <= maxRoundRadius_
+         && cornerType_ != cornerType::SHARP
+        )
+        {
+            roundEdges.set(ei);
+        }
+    }
+
+
+    // Optional binary smoothing (edge-neighbour OR)
+
+
+    // Skip the rest of the routine if the simulation is a serial run
+    if (!Pstream::parRun()) return;
+
+    // Boundary edges
+    const faBoundaryMesh& patches = mesh().boundary();
+    for (const faPatch& fap : patches)
+    {
+        const label patchi = fap.index();
+        const label boundaryEdge0 = fap.start();
+
+        const auto& nfp = nf.boundaryField()[patchi];
+
+        if (isA<processorFaPatch>(fap))
+        {
+            forAll(nfp, bEdgei)
+            {
+                const label meshEdgei = boundaryEdge0 + bEdgei;
+
+                // Do not allow processing of edges shorter than 'minEdgeLength'
+                const edge& e = edges[meshEdgei];
+                const scalar le = e.mag(pts);
+                if (le <= max(minEdgeLength_, VSMALL)) continue;
+
+
+                // Do not allow processing of excluded faces
+                const label faceO = own[meshEdgei];
+                if (!allowedFaces.test(faceO)) continue;
+
+
+                // Fetch normal vector of owner and neigh faces
+                const vector& n0 = nf[faceO];
+                const vector& n1 = nfp[bEdgei];
+
+
+                // Calculate the dihedral angle and curvature per edge
+                const scalar phi = this->dihedralAngle(n0, n1); // [rad]
+                const scalar kappa = 2.0*Foam::sin(0.5*phi)/max(le, VSMALL);
+                const scalar R = (kappa > VSMALL ? scalar(1)/kappa : GREAT);
+
+                const vector tangent(e.unitVec(pts));
+
+                const scalar sgn = curvatureSign(tangent, n0, n1);
+
+                const bool curvatureType =
+                    (cornerCurveType_ == cornerCurveType::ANY)
+                 || (cornerCurveType_ == cornerCurveType::CONVEX  && sgn > 0)
+                 || (cornerCurveType_ == cornerCurveType::CONCAVE && sgn < 0);
+
+                // Do not allow processing of excluded curvature-type faces
+                if (!curvatureType) continue;
+
+
+                // Sharp: dihedral above threshold
+                if (phi >= angSharp && cornerType_ != cornerType::ROUND)
+                {
+                    sharpEdges.set(meshEdgei);
+                    continue; // do not double-classify as round
+                }
+
+
+                // Round: small-to-moderate angle but small radius
+                if
+                (
+                    phi >= angRoundMin && phi <= angRoundMax
+                 && R <= maxRoundRadius_
+                 && cornerType_ != cornerType::SHARP
+                )
+                {
+                    roundEdges.set(meshEdgei);
+                }
+            }
+        }
+        else
+        {
+            forAll(nfp, bEdgei)
+            {
+                const label meshEdgei = boundaryEdge0 + bEdgei;
+                const label faceO = own[meshEdgei];
+
+                if (sharpBoundaryEdges_ && allowedFaces.test(faceO))
+                {
+                    sharpEdges.set(meshEdgei);
+                }
+                // Do not allow round edges on physical boundaries
+            }
+        }
+    }
+}
+
+
+void Foam::cornerDetectionModels::geometryBased::edgesToFaces
+(
+    const bitSet& edgeMask,
+    bitSet& faceMask
+) const
+{
+    // Cache the operand references
+    const labelUList& own = mesh().edgeOwner();
+    const labelUList& nei = mesh().edgeNeighbour();
+
+
+    // Allocate the resource for the return objects
+    faceMask.resize(mesh().nFaces());
+    faceMask.reset();
+
+
+    // Internal edges
+    const label nInternalEdges = mesh().nInternalEdges();
+    for (label ei = 0; ei < nInternalEdges; ++ei)
+    {
+        if (edgeMask.test(ei))
+        {
+            // pick the intersecting owner and neighbour faces at the edge
+            faceMask.set(nei[ei]);
+        }
+    }
+
+
+    // Skip the rest of the routine if the simulation is a serial run
+    if (!Pstream::parRun()) return;
+
+
+    // Boundary edges
+    const faBoundaryMesh& patches = mesh().boundary();
+    for (const faPatch& fap : patches)
+    {
+        const label bEdge0 = fap.start();
+        const label nbEdges = fap.size();
+
+        for (label bEdgei = 0; bEdgei < nbEdges; ++bEdgei)
+        {
+            const label meshEdgei = bEdge0 + bEdgei;
+
+            if (edgeMask.test(meshEdgei))
+            {
+                faceMask.set(own[meshEdgei]);
+            }
+        }
+    }
+}
+
+
+Foam::scalarList Foam::cornerDetectionModels::geometryBased::calcCornerAngles
+(
+    const bitSet& faces,
+    const bitSet& edges
+) const
+{
+    // Cache the operand references
+    const areaVectorField& nf = mesh().faceAreaNormals();
+    const labelUList& own = mesh().edgeOwner();
+    const labelUList& nei = mesh().edgeNeighbour();
+
+
+    // Allocate the resource for the return object
+    scalarList cornerFaceAngles(mesh().faces().size(), Zero);
+
+
+    // Internal edges
+    const label nInternalEdges = mesh().nInternalEdges();
+    for (label ei = 0; ei < nInternalEdges; ++ei)
+    {
+        if (!edges[ei]) continue;
+
+        const label faceO = own[ei];
+        const label faceN = nei[ei];
+
+        // If neither adjacent face is flagged as a corner, skip the atan2 work
+        if (!faces[faceO] && !faces[faceN]) continue;
+
+        const scalar ang = this->dihedralAngle(nf[faceO], nf[faceN]);
+
+        if (faces[faceO]) cornerFaceAngles[faceO] = ang;
+        if (faces[faceN]) cornerFaceAngles[faceN] = ang;
+    }
+
+
+    // Skip the rest of the routine if the simulation is a serial run
+    if (!Pstream::parRun()) return cornerFaceAngles;
+
+
+    // Boundary edges
+    const faBoundaryMesh& patches = mesh().boundary();
+    for (const faPatch& fap : patches)
+    {
+        if (!isA<processorFaPatch>(fap)) continue;
+
+        const label patchi = fap.index();
+        const label bEdge0 = fap.start();
+
+        const auto& nfp = nf.boundaryField()[patchi];
+
+        forAll(nfp, bEdgei)
+        {
+            const label meshEdgei = bEdge0 + bEdgei;
+            const label faceO = own[meshEdgei];
+
+            // Only if the mesh edge is a corner and the owner face is a corner
+            if (!edges[meshEdgei] || !faces[faceO]) continue;
+
+            cornerFaceAngles[faceO] =
+                this->dihedralAngle(nf[faceO], nfp[bEdgei]);
+        }
+    }
+
+    return cornerFaceAngles;
+}
+
+
+void Foam::cornerDetectionModels::geometryBased::writeEdgesAndFaces() const
+{
+    // Cache the operand references
+    const auto& edges = mesh().edges();
+    const auto& faces = mesh().faces();
+    const pointField& pts = mesh().points();
+
+    // Generic writer for masked primitives (edge/face)
+    auto writeMasked =
+    [&](
+        const auto& geom,
+        const auto& mask,
+        const char* file
+    )
+    {
+        OBJstream os(mesh().time().path()/file);
+        forAll(mask, i) if (mask[i]) os.write(geom[i], pts);
+    };
+
+    const bool writeSharp =
+        (cornerType_ == cornerType::ALL || cornerType_ == cornerType::SHARP);
+    const bool writeRound =
+        (cornerType_ == cornerType::ALL || cornerType_ == cornerType::ROUND);
+
+    if (writeSharp)
+    {
+        writeMasked(edges, sharpEdges_, "sharp-edges.obj");
+        writeMasked(faces, sharpFaces_, "sharp-faces.obj");
+    }
+
+    if (writeRound)
+    {
+        writeMasked(edges, roundEdges_, "round-edges.obj");
+        writeMasked(faces, roundFaces_, "round-faces.obj");
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cornerDetectionModels::geometryBased::geometryBased
+(
+    const faMesh& mesh,
+    const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
+    const dictionary& dict
+)
+:
+    cornerDetectionModel(mesh, film, dict),
+    init_(false)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::cornerDetectionModels::geometryBased::detectCorners()
+{
+    if (!init_ || mesh().moving())
+    {
+        // Identify and store sharp/round edges/faces
+        classifyEdges(sharpEdges_, roundEdges_);
+        edgesToFaces(sharpEdges_, sharpFaces_);
+        edgesToFaces(roundEdges_, roundFaces_);
+
+
+        // Collect all operand edges/faces
+        cornerEdges_ = (sharpEdges_ | roundEdges_);
+        cornerFaces_ = (sharpFaces_ | roundFaces_);
+
+
+        // Pass the operand edges/faces to the film-separation model
+        this->setCornerFaces(cornerFaces_);
+        this->setCornerAngles
+        (
+            calcCornerAngles(cornerFaces_, cornerEdges_)
+        );
+
+
+        init_ = true;
+    }
+
+
+    // Write edges and faces as OBJ sets for debug purposes, if need be
+    if (debug && mesh().time().writeTime())
+    {
+        writeEdgesAndFaces();
+    }
+
+    return true;
+}
+
+
+bool Foam::cornerDetectionModels::geometryBased::read(const dictionary& dict)
+{
+    if (!cornerDetectionModel::read(dict))
+    {
+        return false;
+    }
+
+
+    cornerCurveType_ = cornerCurveTypeNames.getOrDefault
+    (
+        "cornerCurveType",
+        dict,
+        cornerCurveType::ANY
+    );
+
+    cornerType_ = cornerTypeNames.getOrDefault
+    (
+        "cornerType",
+        dict,
+        cornerType::ALL
+    );
+
+    angleSharpDeg_ = dict.getOrDefault<scalar>("angleSharp", 45);
+    angleRoundMinDeg_ = dict.getOrDefault<scalar>("angleRoundMin", 5);
+    angleRoundMaxDeg_ = dict.getOrDefault<scalar>("angleRoundMax", 45);
+    maxRoundRadius_ = dict.getOrDefault<scalar>("maxRoundRadius", 2e-3);
+
+    minEdgeLength_ = dict.getOrDefault<scalar>("minEdgeLength", 0);
+    nSmooth_ = dict.getOrDefault<label>("nSmooth", 0);
+
+    sharpBoundaryEdges_ = dict.getOrDefault<bool>("sharpBoundaryEdges", false);
+
+
+    // Validate the input parameters
+    if (angleSharpDeg_ <= 0 || angleSharpDeg_ >= 180)
+    {
+        FatalIOErrorInFunction(dict)
+            << "angleSharp (" << angleSharpDeg_
+            << " deg) must be in (0, 180)."
+            << exit(FatalIOError);
+    }
+
+    if
+    (
+        angleRoundMinDeg_ < 0 || angleRoundMaxDeg_ > 180
+     || angleRoundMinDeg_ > angleRoundMaxDeg_
+    )
+    {
+        FatalIOErrorInFunction(dict)
+            << "Inconsistent round-angle range: angleRoundMin="
+            << angleRoundMinDeg_ << " deg, angleRoundMax=" << angleRoundMaxDeg_
+            << " deg. Require 0 <= min <= max <= 180."
+            << exit(FatalIOError);
+    }
+
+    if (angleSharpDeg_ <= angleRoundMaxDeg_)
+    {
+        WarningInFunction
+            << "angleSharp (" << angleSharpDeg_
+            << " deg) <= angleRoundMax (" << angleRoundMaxDeg_
+            << " deg): sharp vs round thresholds overlap; "
+            << "classification may be ambiguous."
+            << nl;
+    }
+
+    if (maxRoundRadius_ < 0)
+    {
+        FatalIOErrorInFunction(dict)
+            << "maxRoundRadius must be non-negative."
+            << exit(FatalIOError);
+    }
+
+    if (minEdgeLength_ < 0)
+    {
+        FatalIOErrorInFunction(dict)
+            << "minEdgeLength must be non-negative."
+            << exit(FatalIOError);
+    }
+
+    if (nSmooth_ < 0)
+    {
+        FatalIOErrorInFunction(dict)
+            << "nSmooth must be non-negative."
+            << exit(FatalIOError);
+    }
+
+
+    sharpEdges_.clear();
+    roundEdges_.clear();
+    cornerEdges_.clear();
+    sharpFaces_.clear();
+    roundFaces_.clear();
+    cornerFaces_.clear();
+
+
+    // Force the re-identification of corner edges/faces
+    init_ = false;
+
+    return true;
+}
+
+// ************************************************************************* //

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/cornerDetectionModels/geometryBased/geometryBased.H

@@ -0,0 +1,278 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::cornerDetectionModels::geometryBased
+
+Description
+    Geometry-based corner detection model. Marks faces at which the liquid film
+    can separate based on the geometry.
+
+    The model identifies sharp edges based on the face normals of the two
+    faces sharing the edge. Then, if the edge is sharp, the curvature is
+    evaluated to mark the face through which the flux leaves the liquid film.
+    If the edge is sharp and the curvature is of the specified type, the face
+    is marked as a corner face, where the film can separate.
+
+Usage
+    Minimal example in boundary-condition files:
+    \verbatim
+    filmSeparationCoeffs
+    {
+        // Inherited entries
+        ...
+
+        // Optional entries
+        cornerDetectionModel  geometryBased;
+        cornerCurveType       <word>;
+        cornerType            <word>;
+        angleSharp            <scalar>;   // [deg]
+        angleRoundMin         <scalar>;   // [deg]
+        angleRoundMax         <scalar>;   // [deg]
+        maxRoundRadius        <scalar>;   // [m]
+        minEdgeLength         <scalar>;   // [m]
+        nSmooth               <label>;    // [no. of passes]
+        sharpBoundaryEdges    <bool>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                   | Type | Reqd | Deflt
+      cornerDetectionModel | Corner detector model | word | no   | fluxBased
+      cornerCurveType | Corner-curvature type      | word | no   | any
+      cornerType      | Corner type                | word | no   | sharpOrRound
+      angleSharp      | Sharp-angle limit [deg]    | scalar | no | 45
+      angleRoundMin   | Minimum round-angle limit [deg] | scalar | no | 5
+      angleRoundMax   | Maximum round-angle limit [deg] | scalar | no | 45
+      maxRoundRadius  | Maximum round-radius limit [m] | scalar | no | 2e-3
+      minEdgeLength   | Minimum edge length [m]     | scalar | no | 0
+      nSmooth         | No. of smoothing passes     | label | no | 0
+      sharpBoundaryEdges | Treat boundary edges as sharp | bool | no | false
+    \endtable
+
+    Options for the \c cornerCurve entry:
+    \verbatim
+        any          | Convex or concave corners
+        convex       | Convex corners
+        concave      | Concave corners
+    \endverbatim
+
+    Options for the \c cornerType entry:
+    \verbatim
+        sharp        | Sharp corners
+        round        | Round corners
+        sharpOrRound | Sharp or round corners
+    \endverbatim
+
+    The inherited entries are elaborated in:
+      - \link cornerDetectionModel.H \endlink
+
+SourceFiles
+    geometryBased.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_cornerDetectionModels_geometryBased_H
+#define Foam_cornerDetectionModels_geometryBased_H
+
+#include "cornerDetectionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace cornerDetectionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class geometryBased Declaration
+\*---------------------------------------------------------------------------*/
+
+class geometryBased
+:
+    public cornerDetectionModel
+{
+    // Private Enumerations
+
+        //- Options for the corner-curvature type
+        enum cornerCurveType : char
+        {
+            ANY = 0,    //!< "Convex or concave corners"
+            CONVEX,     //!< "Convex corners"
+            CONCAVE     //!< "Concave corners"
+        };
+
+        //- Names for cornerCurveType
+        static const Enum<cornerCurveType> cornerCurveTypeNames;
+
+
+        //- Options for the corner type
+        enum cornerType : char
+        {
+            ALL = 0,    //!< "Sharp or round corners"
+            SHARP,      //!< "Sharp corners"
+            ROUND       //!< "Round corners"
+        };
+
+        //- Names for cornerType
+        static const Enum<cornerType> cornerTypeNames;
+
+
+    // Private Data
+
+        //- Corner-curvature type
+        enum cornerCurveType cornerCurveType_;
+
+        //- Corner type
+        enum cornerType cornerType_;
+
+        //- Flag to deduce if the object is initialised
+        bool init_;
+
+        //- Bitset of edges identified as sharp
+        bitSet sharpEdges_;
+
+        //- Bitset of edges identified as round
+        bitSet roundEdges_;
+
+        //- Bitset of edges identified as a combination of sharp and round
+        bitSet cornerEdges_;
+
+        //- Bitset of faces identified as sharp
+        bitSet sharpFaces_;
+
+        //- Bitset of faces identified as round
+        bitSet roundFaces_;
+
+        //- Bitset of faces identified as a combination of sharp and round
+        bitSet cornerFaces_;
+
+        //- Sharp-angle limit
+        scalar angleSharpDeg_;
+
+        //- Minimum round-angle limit
+        scalar angleRoundMinDeg_;
+
+        //- Maximum round-angle limit
+        scalar angleRoundMaxDeg_;
+
+        //- Maximum round-radius limit
+        scalar maxRoundRadius_;
+
+        //- Minimum edge length; ignore edges shorter than this
+        scalar minEdgeLength_;
+
+        //- Number of smoothing passes on the binary edge mask
+        label nSmooth_;
+
+        //- Flag to treat one-face boundary edges as sharp
+        bool sharpBoundaryEdges_;
+
+
+    // Private Member Functions
+
+        //- Return the signed bending sense, sign(+1/-1) of curvature, across
+        //- an edge with respect to the edge tangent
+        scalar curvatureSign
+        (
+            const vector& t,
+            const vector& n0,
+            const vector& n1
+        ) const;
+
+        //- Classify edges into sharp/round according to dihedral angle and
+        //- inferred radius
+        void classifyEdges
+        (
+            bitSet& sharpEdges,
+            bitSet& roundEdges
+        ) const;
+
+        //- Convert an edge mask to a face mask (face is set if any of its
+        //- edges are set)
+        void edgesToFaces
+        (
+            const bitSet& edgeMask,
+            bitSet& faceMask
+        ) const;
+
+        //- Return the list of corner angles [rad] for each edge
+        scalarList calcCornerAngles
+        (
+            const bitSet& faces,
+            const bitSet& edges
+        ) const;
+
+        // Write edges and faces as OBJ sets for debug purposes
+        void writeEdgesAndFaces() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("geometryBased");
+
+
+    // Constructors
+
+        //- Construct from components
+        geometryBased
+        (
+            const faMesh& mesh,
+            const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+    virtual ~geometryBased() = default;
+
+
+    // Member Functions
+
+    // Evaluation
+
+        //- Detect and store the corner faces
+        virtual bool detectCorners();
+
+
+    // I-O
+
+        //- Read the model dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace cornerDetectionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.C

@@ -26,6 +26,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "FriedrichModel.H"
+#include "cornerDetectionModel.H"
 #include "processorFaPatch.H"
 #include "unitConversion.H"
 #include "addToRunTimeSelectionTable.H"
@@ -53,321 +54,6 @@ FriedrichModel::separationTypeNames
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-bitSet FriedrichModel::calcCornerEdges() const
-{
-    bitSet cornerEdges(mesh().nEdges(), false);
-
-    const areaVectorField& faceCentres = mesh().areaCentres();
-    const areaVectorField& faceNormals = mesh().faceAreaNormals();
-
-    const labelUList& own = mesh().edgeOwner();
-    const labelUList& nbr = mesh().edgeNeighbour();
-
-    // Check if internal face-normal vectors diverge (no separation)
-    // or converge (separation may occur)
-    forAll(nbr, edgei)
-    {
-        const label faceO = own[edgei];
-        const label faceN = nbr[edgei];
-
-        cornerEdges[edgei] = isCornerEdgeSharp
-        (
-            faceCentres[faceO],
-            faceCentres[faceN],
-            faceNormals[faceO],
-            faceNormals[faceN]
-        );
-    }
-
-
-    // Skip the rest of the routine if the simulation is a serial run
-    if (!Pstream::parRun()) return cornerEdges;
-
-    // Check if processor face-normal vectors diverge (no separation)
-    // or converge (separation may occur)
-    const faBoundaryMesh& patches = mesh().boundary();
-
-    for (const faPatch& fap : patches)
-    {
-        if (isA<processorFaPatch>(fap))
-        {
-            const label patchi = fap.index();
-            const auto& edgeFaces = fap.edgeFaces();
-            const label internalEdgei = fap.start();
-
-            const auto& faceCentresp = faceCentres.boundaryField()[patchi];
-            const auto& faceNormalsp = faceNormals.boundaryField()[patchi];
-
-            forAll(faceNormalsp, bndEdgei)
-            {
-                const label faceO = edgeFaces[bndEdgei];
-                const label meshEdgei = internalEdgei + bndEdgei;
-
-                cornerEdges[meshEdgei] = isCornerEdgeSharp
-                (
-                    faceCentres[faceO],
-                    faceCentresp[bndEdgei],
-                    faceNormals[faceO],
-                    faceNormalsp[bndEdgei]
-                );
-            }
-        }
-    }
-
-    return cornerEdges;
-}
-
-
-bool FriedrichModel::isCornerEdgeSharp
-(
-    const vector& faceCentreO,
-    const vector& faceCentreN,
-    const vector& faceNormalO,
-    const vector& faceNormalN
-) const
-{
-    // Calculate the relative position of centres of faces sharing an edge
-    const vector relativePosition(faceCentreN - faceCentreO);
-
-    // Calculate the relative normal of faces sharing an edge
-    const vector relativeNormal(faceNormalN - faceNormalO);
-
-    // Return true if the face normals converge, meaning that the edge is sharp
-    return ((relativeNormal & relativePosition) < -1e-8);
-}
-
-
-scalarList FriedrichModel::calcCornerAngles() const
-{
-    scalarList cornerAngles(mesh().nEdges(), Zero);
-
-    const areaVectorField& faceNormals = mesh().faceAreaNormals();
-
-    const labelUList& own = mesh().edgeOwner();
-    const labelUList& nbr = mesh().edgeNeighbour();
-
-    // Process internal edges
-    forAll(nbr, edgei)
-    {
-        if (!cornerEdges_[edgei]) continue;
-
-        const label faceO = own[edgei];
-        const label faceN = nbr[edgei];
-
-        cornerAngles[edgei] = calcCornerAngle
-        (
-            faceNormals[faceO],
-            faceNormals[faceN]
-        );
-    }
-
-
-    // Skip the rest of the routine if the simulation is a serial run
-    if (!Pstream::parRun()) return cornerAngles;
-
-    // Process processor edges
-    const faBoundaryMesh& patches = mesh().boundary();
-
-    for (const faPatch& fap : patches)
-    {
-        if (isA<processorFaPatch>(fap))
-        {
-            const label patchi = fap.index();
-            const auto& edgeFaces = fap.edgeFaces();
-            const label internalEdgei = fap.start();
-
-            const auto& faceNormalsp = faceNormals.boundaryField()[patchi];
-
-            forAll(faceNormalsp, bndEdgei)
-            {
-                const label faceO = edgeFaces[bndEdgei];
-                const label meshEdgei = internalEdgei + bndEdgei;
-
-                if (!cornerEdges_[meshEdgei]) continue;
-
-                cornerAngles[meshEdgei] = calcCornerAngle
-                (
-                    faceNormals[faceO],
-                    faceNormalsp[bndEdgei]
-                );
-            }
-        }
-    }
-
-    return cornerAngles;
-}
-
-
-scalar FriedrichModel::calcCornerAngle
-(
-    const vector& faceNormalO,
-    const vector& faceNormalN
-) const
-{
-    const scalar magFaceNormal = mag(faceNormalO)*mag(faceNormalN);
-
-    // Avoid any potential exceptions during the cosine calculations
-    if (magFaceNormal < SMALL) return 0;
-
-    scalar cosAngle = (faceNormalO & faceNormalN)/magFaceNormal;
-    cosAngle = clamp(cosAngle, -1, 1);
-
-    return std::acos(cosAngle);
-}
-
-
-bitSet FriedrichModel::calcSeparationFaces() const
-{
-    bitSet separationFaces(mesh().faces().size(), false);
-
-    const edgeScalarField& phis = film().phi2s();
-
-    const labelUList& own = mesh().edgeOwner();
-    const labelUList& nbr = mesh().edgeNeighbour();
-
-    // Process internal faces
-    forAll(nbr, edgei)
-    {
-        if (!cornerEdges_[edgei]) continue;
-
-        const label faceO = own[edgei];
-        const label faceN = nbr[edgei];
-
-        isSeparationFace
-        (
-            separationFaces,
-            phis[edgei],
-            faceO,
-            faceN
-        );
-    }
-
-
-    // Skip the rest of the routine if the simulation is a serial run
-    if (!Pstream::parRun()) return separationFaces;
-
-    // Process processor faces
-    const faBoundaryMesh& patches = mesh().boundary();
-
-    for (const faPatch& fap : patches)
-    {
-        if (isA<processorFaPatch>(fap))
-        {
-            const label patchi = fap.index();
-            const auto& edgeFaces = fap.edgeFaces();
-            const label internalEdgei = fap.start();
-
-            const auto& phisp = phis.boundaryField()[patchi];
-
-            forAll(phisp, bndEdgei)
-            {
-                const label faceO = edgeFaces[bndEdgei];
-                const label meshEdgei(internalEdgei + bndEdgei);
-
-                if (!cornerEdges_[meshEdgei]) continue;
-
-                isSeparationFace
-                (
-                    separationFaces,
-                    phisp[bndEdgei],
-                    faceO
-                );
-            }
-        }
-    }
-
-    return separationFaces;
-}
-
-
-void FriedrichModel::isSeparationFace
-(
-    bitSet& separationFaces,
-    const scalar phiEdge,
-    const label faceO,
-    const label faceN
-) const
-{
-    const scalar tol = 1e-8;
-
-    // Assuming there are no sources/sinks at the edge
-    if (phiEdge > tol)  // From owner to neighbour
-    {
-        separationFaces[faceO] = true;
-    }
-    else if ((phiEdge < -tol) && (faceN != -1))  // From neighbour to owner
-    {
-        separationFaces[faceN] = true;
-    }
-}
-
-
-scalarList FriedrichModel::calcSeparationAngles
-(
-    const bitSet& separationFaces
-) const
-{
-    scalarList separationAngles(mesh().faces().size(), Zero);
-
-    const labelUList& own = mesh().edgeOwner();
-    const labelUList& nbr = mesh().edgeNeighbour();
-
-    // Process internal faces
-    forAll(nbr, edgei)
-    {
-        if (!cornerEdges_[edgei]) continue;
-
-        const label faceO = own[edgei];
-        const label faceN = nbr[edgei];
-
-        if (separationFaces[faceO])
-        {
-            separationAngles[faceO] = cornerAngles_[edgei];
-        }
-
-        if (separationFaces[faceN])
-        {
-            separationAngles[faceN] = cornerAngles_[edgei];
-        }
-    }
-
-
-    // Skip the rest of the routine if the simulation is a serial run
-    if (!Pstream::parRun()) return separationAngles;
-
-    // Process processor faces
-    const edgeScalarField& phis = film().phi2s();
-    const faBoundaryMesh& patches = mesh().boundary();
-
-    for (const faPatch& fap : patches)
-    {
-        if (isA<processorFaPatch>(fap))
-        {
-            const label patchi = fap.index();
-            const auto& edgeFaces = fap.edgeFaces();
-            const label internalEdgei = fap.start();
-
-            const auto& phisp = phis.boundaryField()[patchi];
-
-            forAll(phisp, bndEdgei)
-            {
-                const label faceO = edgeFaces[bndEdgei];
-                const label meshEdgei(internalEdgei + bndEdgei);
-
-                if (!cornerEdges_[meshEdgei]) continue;
-
-                if (separationFaces[faceO])
-                {
-                    separationAngles[faceO] = cornerAngles_[meshEdgei];
-                }
-            }
-        }
-    }
-
-    return separationAngles;
-}
-
-
 tmp<scalarField> FriedrichModel::Fratio() const
 {
     const areaVectorField Up(film().Up());
@@ -378,10 +64,11 @@ tmp<scalarField> FriedrichModel::Fratio() const
     const areaScalarField& sigma = film().sigma();
 
     // Identify the faces where separation may occur
-    const bitSet separationFaces(calcSeparationFaces());
+    const bitSet& separationFaces = cornerDetectorPtr_->getCornerFaces();
 
     // Calculate the corner angles corresponding to the separation faces
-    const scalarList separationAngles(calcSeparationAngles(separationFaces));
+    const scalarList& separationAngles = cornerDetectorPtr_->getCornerAngles();
+
 
     // Initialize the force ratio
     auto tFratio = tmp<scalarField>::New(mesh().faces().size(), Zero);
@@ -431,7 +118,7 @@ tmp<scalarField> FriedrichModel::Fratio() const
         if (isA<processorFaPatch>(fap))
         {
             const label patchi = fap.index();
-            const label internalEdgei = fap.start();
+            const auto& edgeFaces = fap.edgeFaces();
 
             const auto& hp = h.boundaryField()[patchi];
             const auto& Ufp = Uf.boundaryField()[patchi];
@@ -445,18 +132,18 @@ tmp<scalarField> FriedrichModel::Fratio() const
                 // Skip the routine if the face is not a candidate for separation
                 if (!separationFaces[i]) continue;
 
-                const label meshEdgei = internalEdgei + i;
+                const label faceO = edgeFaces[i];
 
                 // Calculate the corner-angle trigonometric values
-                const scalar sinAngle = std::sin(cornerAngles_[meshEdgei]);
-                const scalar cosAngle = std::cos(cornerAngles_[meshEdgei]);
+                const scalar sinAngle = std::sin(separationAngles[faceO]);
+                const scalar cosAngle = std::cos(separationAngles[faceO]);
 
                 // Reynolds number (FLW:Eq. 16)
                 const scalar Re = hp[i]*mag(Ufp[i])*rhop[i]/mup[i];
 
                 // Weber number (FLW:Eq. 17)
-                const vector Urelp(Upp[i] - Ufp[i]);
-                const scalar We = hp[i]*rhop_*sqr(mag(Urelp))/(2.0*sigmap[i]);
+                const vector Urel(Upp[i] - Ufp[i]);
+                const scalar We = hp[i]*rhop_*sqr(mag(Urel))/(2.0*sigmap[i]);
 
                 // Characteristic breakup length (FLW:Eq. 15)
                 const scalar Lb =
@@ -499,13 +186,12 @@ FriedrichModel::FriedrichModel
             separationType::FULL
         )
     ),
+    cornerDetectorPtr_(cornerDetectionModel::New(mesh(), film, dict)),
     rhop_(dict.getScalar("rhop")),
     magG_(mag(film.g().value())),
     C0_(dict.getOrDefault<scalar>("C0", 0.882)),
     C1_(dict.getOrDefault<scalar>("C1", -1.908)),
-    C2_(dict.getOrDefault<scalar>("C2", 1.264)),
-    cornerEdges_(calcCornerEdges()),
-    cornerAngles_(calcCornerAngles())
+    C2_(dict.getOrDefault<scalar>("C2", 1.264))
 {
     if (rhop_ < VSMALL)
     {
@@ -523,10 +209,18 @@ FriedrichModel::FriedrichModel
 }
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+FriedrichModel::~FriedrichModel()
+{}  // cornerDetectionModel was forward declared
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 tmp<scalarField> FriedrichModel::separatedMassRatio() const
 {
+    cornerDetectorPtr_->detectCorners();
+
     tmp<scalarField> tFratio = Fratio();
     const auto& Fratio = tFratio.cref();
 
@@ -576,6 +270,25 @@ tmp<scalarField> FriedrichModel::separatedMassRatio() const
             areaFratio.primitiveFieldRef() = Fratio;
             areaFratio.write();
         }
+
+        {
+            areaScalarField cornerAngles
+            (
+                mesh().newIOobject("cornerAngles"),
+                mesh(),
+                dimensionedScalar(dimless, Zero)
+            );
+
+            const bitSet& cornerFaces = cornerDetectorPtr_->getCornerFaces();
+            const scalarList& angles = cornerDetectorPtr_->getCornerAngles();
+
+            forAll(cornerFaces, i)
+            {
+                if (!cornerFaces[i]) continue;
+                cornerAngles[i] = radToDeg(angles[i]);
+            }
+            cornerAngles.write();
+        }
     }
 
 
@@ -583,6 +296,16 @@ tmp<scalarField> FriedrichModel::separatedMassRatio() const
 }
 
 
+/*
+bool FriedrichModel::read(const dictionary& dict) const
+{
+    // Add the base-class reading later
+    // Read the film separation model dictionary
+
+    return true;
+}
+*/
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace filmSeparationModels

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2024 OpenCFD Ltd.
+    Copyright (C) 2024-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -116,11 +116,14 @@ Usage
     filmSeparationCoeffs
     {
         // Mandatory entries
-        model               Friedrich;
-        rhop                <scalar>;
+        model                 Friedrich;
+        rhop                  <scalar>;
 
         // Optional entries
-        separationType      <word>;
+        separationType        <word>;
+
+        // Inherited entries
+        cornerDetectionModel  <word>;
     }
     \endverbatim
 
@@ -130,14 +133,23 @@ Usage
       model        | Model name: Friedrich              | word | yes  | -
       rhop         | Primary-phase density            | scalar | yes  | -
       separationType | Film separation type             | word | no   | full
+      cornerDetectionModel | Corner detector model      | word | no   | flux
     \endtable
 
+    The inherited entries are elaborated in:
+      - \link cornerDetectionModel.H \endlink
+
     Options for the \c separationType entry:
     \verbatim
       full          | Full film separation (Friedrich et al., 2008)
       partial       | Partial film separation (Zhang et al., 2018)
     \endverbatim
 
+    Options for the \c cornerDetectionModel entry:
+    \verbatim
+      flux          | Flux-based corner detection algorithm
+    \endverbatim
+
 SourceFiles
     FriedrichModel.C
 
@@ -152,6 +164,10 @@ SourceFiles
 
 namespace Foam
 {
+
+// Forward Declarations
+class cornerDetectionModel;
+
 namespace filmSeparationModels
 {
 
@@ -181,6 +197,9 @@ class FriedrichModel
         //- Film separation type
         enum separationType separation_;
 
+        //- Corner-detection model
+        mutable autoPtr<cornerDetectionModel> cornerDetectorPtr_;
+
         //- Approximate uniform mass density of primary phase
         scalar rhop_;
 
@@ -196,53 +215,9 @@ class FriedrichModel
         //- Empirical constant for the partial separation model
         scalar C2_;
 
-        //- List of flags identifying sharp-corner edges where separation
-        //- may occur
-        bitSet cornerEdges_;
-
-        //- Corner angles of sharp-corner edges where separation may occur
-        scalarList cornerAngles_;
-
 
     // Private Member Functions
 
-        //- Return Boolean list of identified sharp-corner edges
-        bitSet calcCornerEdges() const;
-
-        //- Return true if the given edge is identified as sharp
-        bool isCornerEdgeSharp
-        (
-            const vector& faceCentreO,
-            const vector& faceCentreN,
-            const vector& faceNormalO,
-            const vector& faceNormalN
-        ) const;
-
-        //- Return the list of sharp-corner angles for each edge
-        scalarList calcCornerAngles() const;
-
-        //- Return the sharp-corner angle for a given edge
-        scalar calcCornerAngle
-        (
-            const vector& faceNormalO,
-            const vector& faceNormalN
-        ) const;
-
-        //- Return Boolean list of identified separation faces
-        bitSet calcSeparationFaces() const;
-
-        //- Return true if the given face is identified as a separation face
-        void isSeparationFace
-        (
-            bitSet& separationFaces,
-            const scalar phiEdge,
-            const label faceO,
-            const label faceN = -1
-        ) const;
-
-        //- Return the list of sharp-corner angles for each face
-        scalarList calcSeparationAngles(const bitSet& separationFaces) const;
-
         //- Return the film-separation force ratio per face
         tmp<scalarField> Fratio() const;
 
@@ -264,7 +239,7 @@ public:
 
 
     // Destructor
-    virtual ~FriedrichModel() = default;
+    virtual ~FriedrichModel();
 
 
     // Member Functions

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
         {