diff --git a/applications/test/boundSphere/Make/files b/applications/test/boundSphere/Make/files
new file mode 100644
index 0000000000..e02a3c7a27
--- /dev/null
+++ b/applications/test/boundSphere/Make/files
@@ -0,0 +1,3 @@
+Test-boundSphere.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-boundSphere
diff --git a/applications/test/boundSphere/Make/options b/applications/test/boundSphere/Make/options
new file mode 100644
index 0000000000..7ce182425d
--- /dev/null
+++ b/applications/test/boundSphere/Make/options
@@ -0,0 +1,5 @@
+EXE_INC = \
+    -I$(LIB_SRC)/fileFormats/lnInclude
+
+EXE_LIBS = \
+    -lfileFormats
diff --git a/applications/test/boundSphere/Test-boundSphere.C b/applications/test/boundSphere/Test-boundSphere.C
new file mode 100644
index 0000000000..a5cf251b21
--- /dev/null
+++ b/applications/test/boundSphere/Test-boundSphere.C
@@ -0,0 +1,224 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "argList.H"
+#include "boundSphere.H"
+#include "clock.H"
+#include "cpuTime.H"
+#include "OBJstream.H"
+#include "triFace.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::validArgs.append("number of points");
+    argList::addOption
+    (
+        "nTests",
+        "value",
+        "number of tests - default is 1"
+    );
+    argList::addOption
+    (
+        "seed",
+        "value",
+        "random number generator seed - default is clock time"
+    );
+
+    Foam::argList args(argc, argv);
+    Info<< nl;
+
+    const label nPs = args.argRead<label>(1);
+    const label nTests = args.optionLookupOrDefault<label>("nTests", 1);
+
+    // Initialise the random number generator
+    label seed;
+    if (args.optionFound("seed"))
+    {
+        seed = args.optionRead<label>("seed");
+    }
+    else
+    {
+        seed = clock::getTime();
+        Info<< "Seeding random number generator with value " << seed
+            << nl << endl;
+    }
+
+    // Do tests
+    List<point> ps;
+    Tuple2<point, scalar> sphere;
+    cpuTime time;
+    for (label testi = 0; testi < nTests; ++ testi)
+    {
+        Random rndGen(seed + testi);
+
+        // Get random points
+        ps.resize(nPs);
+        forAll(ps, i)
+        {
+            ps[i] = point::uniform(1);
+
+            while (magSqr(ps[i]) > 1)
+            {
+                ps[i] = 2*(rndGen.sample01<point>() - point::uniform(0.5));
+            }
+        }
+
+        // Determine the bounds
+        sphere = boundSphere(ps, rndGen);
+
+        // Check
+        forAll(ps, i)
+        {
+            if
+            (
+                magSqr(ps[i] - sphere.first())
+              > (1 + rootSmall)*sqr(sphere.second())
+            )
+            {
+                FatalErrorInFunction
+                    << "Sphere does not bound point " << i << '/' << nPs << " ("
+                    << magSqr(ps[i] - sphere.first())/sqr(sphere.second())
+                    << ")" << exit(FatalError);
+            }
+        }
+
+        if (nTests != 1)
+        {
+            Info<< "\rSuccessfully computed " << testi + 1 << " bound spheres";
+        }
+    }
+
+    // Report
+    if (nTests == 1)
+    {
+        Info<< "Bounding sphere with centre = " << sphere.first()
+            << ", radius = " << sphere.second() << ", calculated in "
+            << time.cpuTimeIncrement() << " s" << nl << endl;
+
+        // Write the points
+        OBJstream osPs(args[0] + "_points.obj");
+        Info<< "Writing points to " << osPs.name() << endl;
+        forAll(ps, i)
+        {
+            osPs.write(ps[i]);
+        }
+
+        // Write the sphere
+        OBJstream osSphere(args[0] + "_sphere.obj");
+        Info<< "Writing sphere to " << osSphere.name() << endl;
+
+        using namespace constant::mathematical;
+
+        static const label nRings = 20;
+        static const label nRing1Points = 4;
+
+        DynamicList<point> spherePs;
+        DynamicList<face> sphereTris;
+        auto appendP = [&](const scalar x, const scalar y, const scalar z)
+        {
+            spherePs.append(sphere.first() + sphere.second()*point(x, y, z));
+            spherePs.append(sphere.first() + sphere.second()*point(x, y, -z));
+        };
+        auto appendTri = [&](const label a, const label b, const label c)
+        {
+            sphereTris.append(triFace(2*a, 2*b, 2*c));
+            sphereTris.append(triFace(2*a + 1, 2*b + 1, 2*c + 1));
+        };
+
+        appendP(0, 0, 1);
+
+        label nPoints0 = 0, nPoints = 1;
+
+        for (label ringi = 1; ringi < nRings; ++ ringi)
+        {
+            const scalar rXY = Foam::sin(pi/2*ringi/(nRings - 1));
+
+            const label ringPN0 = max((ringi - 1)*nRing1Points, 1);
+            const label ringPN = ringi*nRing1Points;
+
+            label ringPi0 = 0, ringPi = 0;
+
+            for (label i = 0; i < nRing1Points; ++ i)
+            {
+                const scalar theta0 = 2*pi*i/nRing1Points;
+                const scalar theta1 = 2*pi*(i + 1)/nRing1Points;
+
+                for (label j = 0; j < ringi; ++ j)
+                {
+                    const scalar f = scalar(j)/ringi;
+                    const scalar theta = (1 - f)*theta0 + f*theta1;
+
+                    const scalar x = rXY*Foam::cos(theta);
+                    const scalar y = rXY*Foam::sin(theta);
+                    const scalar z = Foam::sqrt(max(1 - sqr(x) - sqr(y), 0));
+
+                    appendP(x, y, z);
+
+                    appendTri
+                    (
+                        nPoints0 + (ringPi0 % ringPN0),
+                        nPoints + ringPi,
+                        nPoints + ((ringPi + 1) % ringPN)
+                    );
+
+                    if (ringi > 1 && j < ringi - 1)
+                    {
+                        appendTri
+                        (
+                            nPoints0 + (ringPi0 % ringPN0),
+                            nPoints + ((ringPi + 1) % ringPN),
+                            nPoints0 + ((ringPi0 + 1) % ringPN0)
+                        );
+                    }
+
+                    if (j < ringi - 1) ringPi0 ++;
+                    ringPi ++;
+                }
+            }
+
+            nPoints0 = nPoints;
+            nPoints += ringPi;
+        }
+
+        osSphere.write(faceList(sphereTris), pointField(spherePs), false);
+    }
+    else
+    {
+        Info<< nl << nl << "Bounding spheres calculated in "
+            << time.cpuTimeIncrement() << " s" << endl;
+    }
+
+    Info<< nl << "End" << nl << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
+
diff --git a/applications/test/patchToPatch/Test-patchToPatch.C b/applications/test/patchToPatch/Test-patchToPatch.C
index 3708c2a910..86fdea3d85 100644
--- a/applications/test/patchToPatch/Test-patchToPatch.C
+++ b/applications/test/patchToPatch/Test-patchToPatch.C
@@ -39,18 +39,66 @@ int main(int argc, char *argv[])
     argList::validArgs.append("source");
     argList::validArgs.append("target");
     argList::validArgs.append("method");
+    argList::validArgs.append("reverse");
+
+    argList::addOption("sourceCase", "dir", "The case with the source patch");
 
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createPolyMesh.H"
 
-    const polyPatch& srcPatch = mesh.boundaryMesh()[args[1]];
-    const polyPatch& tgtPatch = mesh.boundaryMesh()[args[2]];
+    // Optionally read a different mesh for the source
+    autoPtr<Time> srcRunTimePtr;
+    autoPtr<polyMesh> srcMeshPtr;
+    if (args.optionFound("sourceCase"))
+    {
+        const string tgtCase = getEnv("FOAM_CASE");
+        const string tgtCaseName = getEnv("FOAM_CASENAME");
+
+        fileName sourceCase = args["sourceCase"];
+        sourceCase.clean();
+        const fileName sourceCaseName =
+            Pstream::parRun()
+          ? fileName(sourceCase.name())/args.caseName().name()
+          : fileName(sourceCase.name());
+
+        setEnv("FOAM_CASE", sourceCase, true);
+        setEnv("FOAM_CASENAME", sourceCase.name(), true);
+
+        srcRunTimePtr.set
+        (
+            new Time(sourceCase.path(), sourceCaseName)
+        );
+
+        setEnv("FOAM_CASE", tgtCase, true);
+        setEnv("FOAM_CASENAME", tgtCaseName, true);
+
+        srcMeshPtr.set
+        (
+            new polyMesh
+            (
+                Foam::IOobject
+                (
+                    Foam::polyMesh::defaultRegion,
+                    srcRunTimePtr->timeName(),
+                    srcRunTimePtr(),
+                    Foam::IOobject::MUST_READ
+                )
+            )
+        );
+    }
+
+    const polyMesh& srcMesh = srcMeshPtr.valid() ? srcMeshPtr() : mesh;
+    const polyMesh& tgtMesh = mesh;
+
+    const polyPatch& srcPatch = srcMesh.boundaryMesh()[args[1]];
+    const polyPatch& tgtPatch = tgtMesh.boundaryMesh()[args[2]];
     const word& method = args[3];
+    const bool reverse = args.argRead<bool>(4);
 
     cpuTime time;
 
-    patchToPatch::New(method, false)->update
+    patchToPatch::New(method, reverse)->update
     (
         srcPatch,
         srcPatch.pointNormals(),
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 6706ba39ae..1ba10ca52d 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -715,16 +715,12 @@ $(interpolationWeights)/interpolationWeights/interpolationWeights.C
 $(interpolationWeights)/linearInterpolationWeights/linearInterpolationWeights.C
 $(interpolationWeights)/splineInterpolationWeights/splineInterpolationWeights.C
 
-
-
 algorithms/indexedOctree/indexedOctreeName.C
 algorithms/indexedOctree/treeDataCell.C
 algorithms/indexedOctree/volumeType.C
-algorithms/polygonTriangulate/polygonTriangulate.C
-
-
 algorithms/dynamicIndexedOctree/dynamicIndexedOctreeName.C
 algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.C
+algorithms/polygonTriangulate/polygonTriangulate.C
 
 meshes/data/data.C
 meshes/Residuals/residuals.C
diff --git a/src/OpenFOAM/algorithms/boundSphere/boundSphere.H b/src/OpenFOAM/algorithms/boundSphere/boundSphere.H
new file mode 100644
index 0000000000..30939b7821
--- /dev/null
+++ b/src/OpenFOAM/algorithms/boundSphere/boundSphere.H
@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Functions for constructing bounding spheres of lists of points
+
+SourceFiles
+    boundSphere.C
+    boundSphereTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef boundSphere_H
+#define boundSphere_H
+
+#include "point.H"
+#include "Random.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+//- Return whether or not the given sphere is valid
+inline bool isValidBoundSphere(const Tuple2<point, scalar>& sphere);
+
+//- Compute a sphere of four points or less where every point intersects the
+//  sphere's surface
+template<class PointField>
+Tuple2<point, scalar> intersectBoundSphere
+(
+    const PointField& ps,
+    const FixedList<label, 4>& pis,
+    const label nPs
+);
+
+//- Compute a bounding sphere of four points or less
+template<class PointField>
+Tuple2<point, scalar> trivialBoundSphere
+(
+    const PointField& ps,
+    const FixedList<label, 4>& pis,
+    const label nPs
+);
+
+//- Compute a bounding sphere for an arbitrary number of points recursively
+//  using Weizl's algorithm
+template<class PointField>
+Tuple2<point, scalar> weizlBoundSphere
+(
+    const PointField& ps,
+    List<label>& pis,
+    const label nPs,
+    FixedList<label, 4>& boundaryPis,
+    const label nBoundaryPs
+);
+
+//- Compute a bounding sphere for an arbitrary number of points, and given an
+//  engine with which to randomise Weizl's algorithm
+template<class PointField>
+Tuple2<point, scalar> boundSphere(const PointField& ps, Random& rndGen);
+
+//- Compute a bounding sphere for an arbitrary number of points
+template<class PointField>
+Tuple2<point, scalar> boundSphere(const PointField& ps);
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "boundSphereTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/algorithms/boundSphere/boundSphereTemplates.C b/src/OpenFOAM/algorithms/boundSphere/boundSphereTemplates.C
new file mode 100644
index 0000000000..e9b3bbc93d
--- /dev/null
+++ b/src/OpenFOAM/algorithms/boundSphere/boundSphereTemplates.C
@@ -0,0 +1,301 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "boundSphere.H"
+#include "labelPair.H"
+#include "triPointRef.H"
+#include "tetPointRef.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline bool Foam::isValidBoundSphere(const Tuple2<point, scalar>& sphere)
+{
+    return sphere.second() >= 0;
+}
+
+
+template<class PointField>
+Foam::Tuple2<Foam::point, Foam::scalar>
+Foam::intersectBoundSphere
+(
+    const PointField& ps,
+    const FixedList<label, 4>& pis,
+    const label nPs
+)
+{
+    switch (nPs)
+    {
+        case 0:
+            return Tuple2<point, scalar>(point::uniform(NaN), -vGreat);
+
+        case 1:
+            return Tuple2<point, scalar>(ps[pis[0]], 0);
+
+        case 2:
+            // Return a mid point-centred sphere
+            return Tuple2<point, scalar>
+            (
+                (ps[pis[0]] + ps[pis[1]])/2,
+                mag(ps[pis[0]] - ps[pis[1]])/2
+            );
+
+        case 3:
+        {
+            // Return the sphere that intersects the triangle
+            return
+                triPointRef
+                (
+                    ps[pis[0]],
+                    ps[pis[1]],
+                    ps[pis[2]]
+                ).circumCircle();
+        }
+
+        case 4:
+        {
+            // Return the sphere that intersects the tetrahedron
+            return
+                tetPointRef
+                (
+                    ps[pis[0]],
+                    ps[pis[1]],
+                    ps[pis[2]],
+                    ps[pis[3]]
+                ).circumSphere();
+        }
+
+        default:
+            FatalErrorInFunction
+                << "Cannot compute the intersect bounding sphere of more than "
+                << "four points" << exit(FatalError);
+            return Tuple2<point, scalar>(point::uniform(NaN), NaN);
+    }
+}
+
+
+template<class PointField>
+Foam::Tuple2<Foam::point, Foam::scalar>
+Foam::trivialBoundSphere
+(
+    const PointField& ps,
+    const FixedList<label, 4>& pis,
+    const label nPs
+)
+{
+    static const scalar tol = 1 + 2*sqrt(small*rootSmall);
+
+    // Search for spheres that intersect sub-sets of the points that bound all
+    // the other points
+    switch (nPs)
+    {
+        case 0:
+        case 1:
+        case 2:
+            break;
+
+        case 3:
+        {
+            // Test the spheres that intersect the edges
+            for (label i = 0; i < 3; ++ i)
+            {
+                const point& p0 = ps[pis[i]];
+                const point& p1 = ps[pis[(i + 1) % 3]];
+                const point& pOpp = ps[pis[(i + 2) % 3]];
+
+                const point c = (p0 + p1)/2;
+                const scalar rSqr = magSqr(p0 - p1)/4;
+
+                if (magSqr(pOpp - c) <= tol*rSqr)
+                {
+                    return Tuple2<point, scalar>(c, Foam::sqrt(rSqr));
+                }
+            }
+
+            break;
+        }
+
+        case 4:
+        {
+            // Test the spheres that intersect the edges
+            static const FixedList<labelPair, 6> p01s =
+                {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {1, 3}, {2, 3}};
+            static const FixedList<labelPair, 6> pOpp01s =
+                {{2, 3}, {1, 3}, {1, 2}, {0, 3}, {0, 2}, {0, 1}};
+            for (label i = 0; i < 6; ++ i)
+            {
+                const point& p0 = ps[pis[p01s[i].first()]];
+                const point& p1 = ps[pis[p01s[i].second()]];
+                const point& pOpp0 = ps[pis[pOpp01s[i].first()]];
+                const point& pOpp1 = ps[pis[pOpp01s[i].second()]];
+
+                const point c = (p0 + p1)/2;
+                const scalar rSqr = magSqr(p0 - p1)/4;
+
+                if
+                (
+                    magSqr(pOpp0 - c) <= tol*rSqr
+                 && magSqr(pOpp1 - c) <= tol*rSqr
+                )
+                {
+                    return Tuple2<point, scalar>(c, Foam::sqrt(rSqr));
+                }
+            }
+
+            // Test the spheres that intersect the triangles
+            point minC = point::uniform(NaN);
+            scalar minR = vGreat;
+            for (label i = 0; i < 4; ++ i)
+            {
+                const triPointRef tri
+                (
+                    ps[pis[i]],
+                    ps[pis[(i + 1) % 4]],
+                    ps[pis[(i + 2) % 4]]
+                );
+                const point& pOpp = ps[pis[(i + 3) % 4]];
+
+                const Tuple2<point, scalar> circ = tri.circumCircle();
+                const point& c = circ.first();
+                const scalar r = circ.second();
+
+                if (magSqr(pOpp - c) <= tol*sqr(r) && r < minR)
+                {
+                    minC = c;
+                    minR = r;
+                }
+            }
+            if (minR != vGreat)
+            {
+                return Tuple2<point, scalar>(minC, minR);
+            }
+
+            break;
+        }
+
+        default:
+            FatalErrorInFunction
+                << "Cannot compute the trivial bounding sphere of more than "
+                << "four points" << exit(FatalError);
+            return Tuple2<point, scalar>(point::uniform(NaN), NaN);
+    }
+
+    // Return the sphere that intersects the points
+    return intersectBoundSphere(ps, pis, nPs);
+}
+
+
+template<class PointField>
+Foam::Tuple2<Foam::point, Foam::scalar> Foam::weizlBoundSphere
+(
+    const PointField& ps,
+    List<label>& pis,
+    const label nPs,
+    FixedList<label, 4>& boundaryPis,
+    const label nBoundaryPs
+)
+{
+    static const scalar tol = 1 + 2*sqrt(small*rootSmall);
+
+    Tuple2<point, scalar> sphere;
+
+    if (nBoundaryPs != 0)
+    {
+        sphere = intersectBoundSphere(ps, boundaryPis, nBoundaryPs);
+    }
+
+    if (nBoundaryPs == 4)
+    {
+        return sphere;
+    }
+
+    for (label i = 0; i < nPs; ++ i)
+    {
+        if
+        (
+            (nBoundaryPs == 0 && i == 0)
+         || (magSqr(ps[pis[i]] - sphere.first()) > tol*sqr(sphere.second()))
+        )
+        {
+            boundaryPis[nBoundaryPs] = pis[i];
+
+            sphere = weizlBoundSphere(ps, pis, i, boundaryPis, nBoundaryPs + 1);
+
+            // Move the limiting point to the start of the list so that the
+            // sphere grows as quickly as possible in the recursive calls
+            Swap(pis[0], pis[i]);
+        }
+    }
+
+    return sphere;
+}
+
+
+template<class PointField>
+Foam::Tuple2<Foam::point, Foam::scalar>
+Foam::boundSphere(const PointField& ps, Random& rndGen)
+{
+    if (ps.size() <= 4)
+    {
+        static const FixedList<label, 4> pis({0, 1, 2, 3});
+        return trivialBoundSphere(ps, pis, ps.size());
+    }
+    else
+    {
+        labelList pis(identity(ps.size()));
+        forAll(pis, i)
+        {
+            Swap(pis[i], pis[rndGen.sampleAB<label>(i, ps.size())]);
+        }
+        FixedList<label, 4> boundaryPis({-1, -1, -1, -1});
+        return weizlBoundSphere(ps, pis, ps.size(), boundaryPis, 0);
+    }
+}
+
+
+template<class PointField>
+Foam::Tuple2<Foam::point, Foam::scalar>
+Foam::boundSphere(const PointField& ps)
+{
+    if (ps.size() <= 4)
+    {
+        static const FixedList<label, 4> pis({0, 1, 2, 3});
+        return trivialBoundSphere(ps, pis, ps.size());
+    }
+    else
+    {
+        labelList pis(identity(ps.size()));
+        Random rndGen(0);
+        forAll(pis, i)
+        {
+            Swap(pis[i], pis[rndGen.sampleAB<label>(i, ps.size())]);
+        }
+        FixedList<label, 4> boundaryPis({-1, -1, -1, -1});
+        return weizlBoundSphere(ps, pis, ps.size(), boundaryPis, 0);
+    }
+}
+
+
+// ************************************************************************* //
+
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index d4eb7119cc..fdbdb6565c 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -220,6 +220,7 @@ patchToPatchTools/patchToPatchTools.C
 patchToPatch/patchToPatch/patchToPatch.C
 patchToPatch/patchToPatch/patchToPatchParallelOps.C
 patchToPatch/matching/matchingPatchToPatch.C
+patchToPatch/nearby/nearbyPatchToPatch.C
 patchToPatch/nearest/nearestPatchToPatch.C
 patchToPatch/inverseDistance/inverseDistancePatchToPatch.C
 patchToPatch/intersection/intersectionPatchToPatch.C
diff --git a/src/meshTools/meshToMesh/meshToMesh.C b/src/meshTools/meshToMesh/meshToMesh.C
index bfef48fd81..dd2d6b142e 100644
--- a/src/meshTools/meshToMesh/meshToMesh.C
+++ b/src/meshTools/meshToMesh/meshToMesh.C
@@ -265,11 +265,11 @@ void Foam::meshToMesh::constructNoCuttingPatches
         srcToTgtPatchIDs_.transfer(srcToTgtPatchIDs);
     }
 
-    // calculate cell addressing and weights
-    calculateCellToCells(methodName);
-
     // calculate patch addressing and weights
     calculatePatchToPatches(methodName);
+
+    // calculate cell addressing and weights
+    calculateCellToCells(methodName);
 }
 
 
@@ -295,12 +295,12 @@ void Foam::meshToMesh::constructFromCuttingPatches
             );
     }
 
-    // calculate cell addressing and weights
-    calculateCellToCells(methodName);
-
     // calculate patch addressing and weights
     calculatePatchToPatches(methodName);
 
+    // calculate cell addressing and weights
+    calculateCellToCells(methodName);
+
     // set IDs of cutting patches on target mesh
     tgtCuttingPatchIDs_.setSize(tgtCuttingPatches.size());
     forAll(tgtCuttingPatchIDs_, i)
diff --git a/src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.C b/src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.C
index 445f5b9804..922f3acaaf 100644
--- a/src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.C
+++ b/src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.C
@@ -542,13 +542,22 @@ void Foam::patchToPatches::intersection::initialise
 }
 
 
-Foam::labelList Foam::patchToPatches::intersection::trimLocalTgt
+Foam::labelList Foam::patchToPatches::intersection::finaliseLocal
 (
-    const primitiveOldTimePatch& localTgtPatch
+    const primitiveOldTimePatch& srcPatch,
+    const vectorField& srcPointNormals,
+    const vectorField& srcPointNormals0,
+    const primitiveOldTimePatch& tgtPatch
 )
 {
     const labelList newToOldLocalTgtFace =
-        patchToPatch::trimLocalTgt(localTgtPatch);
+        patchToPatch::finaliseLocal
+        (
+            srcPatch,
+            srcPointNormals,
+            srcPointNormals0,
+            tgtPatch
+        );
 
     tgtCouples_ = List<DynamicList<couple>>(tgtCouples_, newToOldLocalTgtFace);
 
diff --git a/src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.H b/src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.H
index 704e1d7d65..5c61eb543f 100644
--- a/src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.H
+++ b/src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.H
@@ -328,11 +328,15 @@ private:
             const primitiveOldTimePatch& tgtPatch
         );
 
-        //- Trim the local target patch so that only parts that actually
-        //  intersect the source remain
-        virtual labelList trimLocalTgt
+        //- Finalise the intersection locally. Trims the local target patch
+        //  so that only parts that actually intersect the source remain.
+        //  Returns new-to-old map for trimming target-associated data.
+        virtual labelList finaliseLocal
         (
-            const primitiveOldTimePatch& localTgtPatch
+            const primitiveOldTimePatch& srcPatch,
+            const vectorField& srcPointNormals,
+            const vectorField& srcPointNormals0,
+            const primitiveOldTimePatch& tgtPatch
         );
 
         //- Send data that resulted from an intersection between the source
diff --git a/src/meshTools/patchToPatch/inverseDistance/inverseDistancePatchToPatch.C b/src/meshTools/patchToPatch/inverseDistance/inverseDistancePatchToPatch.C
index 3b31271341..b65c188617 100644
--- a/src/meshTools/patchToPatch/inverseDistance/inverseDistancePatchToPatch.C
+++ b/src/meshTools/patchToPatch/inverseDistance/inverseDistancePatchToPatch.C
@@ -39,31 +39,15 @@ namespace patchToPatches
 }
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * Private Static Member Functions * * * * * * * * * * //
 
-Foam::treeBoundBox Foam::patchToPatches::inverseDistance::srcBox
+bool Foam::patchToPatches::inverseDistance::rayHitsFace
 (
-    const face& srcFace,
-    const pointField& srcPoints,
-    const vectorField& srcPointNormals
-) const
-{
-    const treeBoundBox bb(srcPoints, srcFace);
-
-    const point c = bb.midpoint();
-    const scalar l = bb.maxDim();
-
-    return treeBoundBox(c - l*vector::one, c + l*vector::one);
-}
-
-
-bool Foam::patchToPatches::inverseDistance::inside
-(
-    const face& f,
-    const pointField& ps,
     const point& p,
-    const vector& r
-) const
+    const vector& r,
+    const face& f,
+    const pointField& ps
+)
 {
     using namespace constant::mathematical;
 
@@ -75,100 +59,23 @@ bool Foam::patchToPatches::inverseDistance::inside
     {
         const vector& a = T & (ps[f[i]] - p);
         const vector& b = T & (ps[f[f.fcIndex(i)]] - p);
-        const scalar magAB = sqrt(magSqr(a)*magSqr(b));
-        angle +=
-            (reverse() ? +1 : -1)
-           *sign(r & (a ^ b))
-           *acos(min(max(-1, (a & b)/magAB), +1));
+
+        const scalar meanMagSqrAB = (magSqr(a) + magSqr(b))/2;
+        const scalar geometricMeanMagSqrAB = sqrt(magSqr(a)*magSqr(b));
+
+        // This indicates that we have hit a point to within round off error
+        if (geometricMeanMagSqrAB < small*meanMagSqrAB) return true;
+
+        angle -=
+            sign(r & (a ^ b))
+           *acos(min(max(-1, (a & b)/geometricMeanMagSqrAB), +1));
     }
 
     return pi < angle && angle < 3*pi;
 }
 
 
-bool Foam::patchToPatches::inverseDistance::intersectFaces
-(
-    const primitivePatch& patch,
-    const primitivePatch& otherPatch,
-    const label facei,
-    const label otherFacei,
-    DynamicList<label>& faceOtherFaces,
-    DynamicList<scalar>& faceWeights
-) const
-{
-    const face& f = otherPatch[otherFacei];
-    const pointField& ps = otherPatch.points();
-
-    const point& p = patch.faceCentres()[facei];
-    const vector& r = patch.faceNormals()[facei];
-
-    bool intersectsOther = inside(f, ps, p, r);
-
-    if (!intersectsOther)
-    {
-        forAll(otherPatch.faceFaces()[otherFacei], otherFaceFacei)
-        {
-            const label otherFacej =
-                otherPatch.faceFaces()[otherFacei][otherFaceFacei];
-
-            const face& g = otherPatch[otherFacej];
-
-            if (inside(g, ps, p, r))
-            {
-                intersectsOther = true;
-                break;
-            }
-        }
-    }
-
-    if (intersectsOther)
-    {
-        faceOtherFaces.append(otherFacei);
-        faceWeights.append
-        (
-            1/max(mag(p - otherPatch.faceCentres()[otherFacei]), vSmall)
-        );
-    }
-
-    return intersectsOther;
-}
-
-
-bool Foam::patchToPatches::inverseDistance::intersectFaces
-(
-    const primitiveOldTimePatch& srcPatch,
-    const vectorField& srcPointNormals,
-    const vectorField& srcPointNormals0,
-    const primitiveOldTimePatch& tgtPatch,
-    const label srcFacei,
-    const label tgtFacei
-)
-{
-    const bool srcCouples =
-        intersectFaces
-        (
-            srcPatch,
-            tgtPatch,
-            srcFacei,
-            tgtFacei,
-            srcLocalTgtFaces_[srcFacei],
-            srcWeights_[srcFacei]
-        );
-
-    const bool tgtCouples =
-        intersectFaces
-        (
-            tgtPatch,
-            srcPatch,
-            tgtFacei,
-            srcFacei,
-            tgtLocalSrcFaces_[tgtFacei],
-            tgtWeights_[tgtFacei]
-        );
-
-    return srcCouples || tgtCouples;
-}
-
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::patchToPatches::inverseDistance::initialise
 (
@@ -178,7 +85,7 @@ void Foam::patchToPatches::inverseDistance::initialise
     const primitiveOldTimePatch& tgtPatch
 )
 {
-    patchToPatch::initialise
+    nearby::initialise
     (
         srcPatch,
         srcPointNormals,
@@ -187,17 +94,122 @@ void Foam::patchToPatches::inverseDistance::initialise
     );
 
     srcWeights_.resize(srcPatch.size());
+    forAll(srcWeights_, i)
+    {
+        srcWeights_[i].clear();
+    }
+
     tgtWeights_.resize(tgtPatch.size());
+    forAll(tgtWeights_, i)
+    {
+        tgtWeights_[i].clear();
+    }
 }
 
 
-Foam::labelList Foam::patchToPatches::inverseDistance::trimLocalTgt
+void Foam::patchToPatches::inverseDistance::generateWeights
 (
-    const primitiveOldTimePatch& localTgtPatch
+    const primitiveOldTimePatch& srcPatch,
+    const primitiveOldTimePatch& tgtPatch
 )
 {
+    auto generate = []
+    (
+        const primitiveOldTimePatch& patch,
+        const primitiveOldTimePatch& otherPatch,
+        const bool reverse,
+        List<DynamicList<label>>& otherFaces,
+        List<DynamicList<scalar>>& weights
+    )
+    {
+        forAll(otherFaces, facei)
+        {
+            if (otherFaces[facei].empty()) continue;
+
+            label otherFacei = -1;
+
+            // Find the other face that "contains" this face's centre
+            forAll(otherFaces[facei], i)
+            {
+                if
+                (
+                    rayHitsFace
+                    (
+                        patch.faceCentres()[facei],
+                        (reverse ? -1 : +1)*patch.faceNormals()[facei],
+                        otherPatch[otherFaces[facei][i]],
+                        otherPatch.points()
+                    )
+                )
+                {
+                    otherFacei = otherFaces[facei][i];
+                    break;
+                }
+            }
+
+            const point& c = patch.faceCentres()[facei];
+
+            // If the above failed, find the closest
+            if (otherFacei == -1)
+            {
+                scalar minDistSqr = vGreat;
+
+                forAll(otherFaces[facei], i)
+                {
+                    const point& otherC =
+                        otherPatch.faceCentres()[otherFaces[facei][i]];
+                    const scalar distSqr = magSqr(c - otherC);
+                    if (distSqr < minDistSqr)
+                    {
+                        minDistSqr = distSqr;
+                        otherFacei = otherFaces[facei][i];
+                    }
+                }
+            }
+
+            // Remove all faces
+            otherFaces[facei].clear();
+
+            // Add the found face and all its neighbours
+            const point& otherC = otherPatch.faceCentres()[otherFacei];
+            otherFaces[facei].append(otherFacei);
+            weights[facei].append(1/(mag(c - otherC) + rootVSmall));
+
+            forAll(otherPatch.faceFaces()[otherFacei], i)
+            {
+                const label otherFacej = otherPatch.faceFaces()[otherFacei][i];
+
+                const point& otherC = otherPatch.faceCentres()[otherFacej];
+                otherFaces[facei].append(otherFacej);
+                weights[facei].append(1/(mag(c - otherC) + rootVSmall));
+            }
+        }
+    };
+
+    generate(srcPatch, tgtPatch, reverse_, srcLocalTgtFaces_, srcWeights_);
+    generate(tgtPatch, srcPatch, reverse_, tgtLocalSrcFaces_, tgtWeights_);
+}
+
+
+Foam::labelList Foam::patchToPatches::inverseDistance::finaliseLocal
+(
+    const primitiveOldTimePatch& srcPatch,
+    const vectorField& srcPointNormals,
+    const vectorField& srcPointNormals0,
+    const primitiveOldTimePatch& tgtPatch
+)
+{
+    // Transfer weight addressing into actual addressing
+    generateWeights(srcPatch, tgtPatch);
+
     const labelList newToOldLocalTgtFace =
-        patchToPatch::trimLocalTgt(localTgtPatch);
+        nearby::finaliseLocal
+        (
+            srcPatch,
+            srcPointNormals,
+            srcPointNormals0,
+            tgtPatch
+        );
 
     tgtWeights_ = List<DynamicList<scalar>>(tgtWeights_, newToOldLocalTgtFace);
 
@@ -210,8 +222,10 @@ void Foam::patchToPatches::inverseDistance::rDistributeTgt
     const primitiveOldTimePatch& tgtPatch
 )
 {
-    patchToPatch::rDistributeTgt(tgtPatch);
+    // Let the base class reverse distribute the addressing
+    nearby::rDistributeTgt(tgtPatch);
 
+    // Reverse distribute the weights
     patchToPatchTools::rDistributeListList
     (
         tgtPatch.size(),
@@ -240,20 +254,25 @@ Foam::label Foam::patchToPatches::inverseDistance::finalise
             tgtToSrc
         );
 
+    // Transfer weight addressing into actual addressing (if not done in the
+    // finaliseLocal method above)
+    if (isSingleProcess())
+    {
+        generateWeights(srcPatch, tgtPatch);
+    }
+
+    // Normalise the weights
     forAll(srcWeights_, srcFacei)
     {
         const scalar w = sum(srcWeights_[srcFacei]);
-
         forAll(srcWeights_[srcFacei], i)
         {
             srcWeights_[srcFacei][i] /= max(w, vSmall);
         }
     }
-
     forAll(tgtWeights_, tgtFacei)
     {
         const scalar w = sum(tgtWeights_[tgtFacei]);
-
         forAll(tgtWeights_[tgtFacei], i)
         {
             tgtWeights_[tgtFacei][i] /= max(w, vSmall);
@@ -279,8 +298,10 @@ Foam::label Foam::patchToPatches::inverseDistance::finalise
             return result;
         };
 
-        Info<< "Number of source faces by number of target connections = "
+        Info<< indent
+            << "Number of source faces by number of target connections = "
             << histogram(srcLocalTgtFaces_) << nl
+            << indent
             << "Number of target faces by number of source connections = "
             << histogram(tgtLocalSrcFaces_) << endl;
     }
@@ -307,7 +328,7 @@ Foam::patchToPatches::inverseDistance::tgtWeights() const
 
 Foam::patchToPatches::inverseDistance::inverseDistance(const bool reverse)
 :
-    patchToPatch(reverse),
+    nearby(reverse),
     srcWeights_(),
     tgtWeights_()
 {}
diff --git a/src/meshTools/patchToPatch/inverseDistance/inverseDistancePatchToPatch.H b/src/meshTools/patchToPatch/inverseDistance/inverseDistancePatchToPatch.H
index 685a89ba9a..5ce25a8483 100644
--- a/src/meshTools/patchToPatch/inverseDistance/inverseDistancePatchToPatch.H
+++ b/src/meshTools/patchToPatch/inverseDistance/inverseDistancePatchToPatch.H
@@ -38,7 +38,7 @@ SourceFiles
 #ifndef inverseDistancePatchToPatch_H
 #define inverseDistancePatchToPatch_H
 
-#include "patchToPatch.H"
+#include "nearbyPatchToPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,7 +53,7 @@ namespace patchToPatches
 
 class inverseDistance
 :
-    public patchToPatch
+    public nearby
 {
     // Private Member Data
 
@@ -64,47 +64,21 @@ class inverseDistance
         List<DynamicList<scalar>> tgtWeights_;
 
 
-    // Private Member Functions
+    // Private Static Member Functions
 
-        //- Get the bound box for a source face
-        virtual treeBoundBox srcBox
+        //- Return whether a ray starting at point `p` with direction `r` hits
+        //  a face `f` with points `ps`
+        static bool rayHitsFace
         (
-            const face& srcFace,
-            const pointField& srcPoints,
-            const vectorField& srcPointNormals
-        ) const;
-
-        //- Return whether or not the face contains a point
-        bool inside
-        (
-            const face& f,
-            const pointField& ps,
             const point& p,
-            const vector& r
-        ) const;
-
-        //- Intersect two faces
-        virtual bool intersectFaces
-        (
-            const primitivePatch& patch,
-            const primitivePatch& otherPatch,
-            const label facei,
-            const label otherFacei,
-            DynamicList<label>& faceOtherFaces,
-            DynamicList<scalar>& faceWeights
-        ) const;
-
-        //- Intersect two faces
-        virtual bool intersectFaces
-        (
-            const primitiveOldTimePatch& srcPatch,
-            const vectorField& srcPointNormals,
-            const vectorField& srcPointNormals0,
-            const primitiveOldTimePatch& tgtPatch,
-            const label srcFacei,
-            const label tgtFacei
+            const vector& r,
+            const face& f,
+            const pointField& ps
         );
 
+
+    // Private Member Functions
+
         //- Initialisation
         virtual void initialise
         (
@@ -114,11 +88,23 @@ class inverseDistance
             const primitiveOldTimePatch& tgtPatch
         );
 
-        //- Trim the local target patch so that only parts that actually
-        //  intersect the source remain
-        virtual labelList trimLocalTgt
+        //- Convert the addressing generated by the base class into weights
+        //  (and addressing) appropriate for inverse distance interpolation
+        void generateWeights
         (
-            const primitiveOldTimePatch& localTgtPatch
+            const primitiveOldTimePatch& srcPatch,
+            const primitiveOldTimePatch& tgtPatch
+        );
+
+        //- Finalise the intersection locally. Trims the local target patch
+        //  so that only parts that actually intersect the source remain.
+        //  Returns new-to-old map for trimming target-associated data.
+        virtual labelList finaliseLocal
+        (
+            const primitiveOldTimePatch& srcPatch,
+            const vectorField& srcPointNormals,
+            const vectorField& srcPointNormals0,
+            const primitiveOldTimePatch& tgtPatch
         );
 
         //- Send data that resulted from an intersection between the source
diff --git a/src/meshTools/patchToPatch/nearby/nearbyPatchToPatch.C b/src/meshTools/patchToPatch/nearby/nearbyPatchToPatch.C
new file mode 100644
index 0000000000..4afa25d293
--- /dev/null
+++ b/src/meshTools/patchToPatch/nearby/nearbyPatchToPatch.C
@@ -0,0 +1,145 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "nearbyPatchToPatch.H"
+#include "boundSphere.H"
+#include "OFstream.H"
+#include "OBJstream.H"
+#include "vtkWritePolyData.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace patchToPatches
+{
+    defineTypeNameAndDebug(nearby, 0);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::treeBoundBox Foam::patchToPatches::nearby::srcBox
+(
+    const face& srcFace,
+    const pointField& srcPoints,
+    const vectorField& srcPointNormals
+) const
+{
+    const treeBoundBox bb(srcPoints, srcFace);
+
+    const point c = bb.midpoint();
+    const scalar l = bb.maxDim();
+
+    return treeBoundBox(c - l*vector::one, c + l*vector::one);
+}
+
+
+bool Foam::patchToPatches::nearby::intersectFaces
+(
+    const primitiveOldTimePatch& srcPatch,
+    const vectorField& srcPointNormals,
+    const vectorField& srcPointNormals0,
+    const primitiveOldTimePatch& tgtPatch,
+    const label srcFacei,
+    const label tgtFacei
+)
+{
+    const point& srcC = srcSpheres_[srcFacei].first();
+    const scalar srcR = srcSpheres_[srcFacei].second();
+    const point& tgtC = tgtSpheres_[tgtFacei].first();
+    const scalar tgtR = tgtSpheres_[tgtFacei].second();
+
+    if (magSqr(srcC - tgtC) < sqr(srcR + tgtR))
+    {
+        srcLocalTgtFaces_[srcFacei].append(tgtFacei);
+        tgtLocalSrcFaces_[tgtFacei].append(srcFacei);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void Foam::patchToPatches::nearby::initialise
+(
+    const primitiveOldTimePatch& srcPatch,
+    const vectorField& srcPointNormals,
+    const vectorField& srcPointNormals0,
+    const primitiveOldTimePatch& tgtPatch
+)
+{
+    patchToPatch::initialise
+    (
+        srcPatch,
+        srcPointNormals,
+        srcPointNormals0,
+        tgtPatch
+    );
+
+    srcSpheres_.resize(srcPatch.size());
+    forAll(srcPatch, srcFacei)
+    {
+        srcSpheres_[srcFacei] =
+            boundSphere
+            (
+                UIndirectList<point>(srcPatch.points(), srcPatch[srcFacei])
+            );
+    }
+
+    tgtSpheres_.resize(tgtPatch.size());
+    forAll(tgtPatch, tgtFacei)
+    {
+        tgtSpheres_[tgtFacei] =
+            boundSphere
+            (
+                UIndirectList<point>(tgtPatch.points(), tgtPatch[tgtFacei])
+            );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::patchToPatches::nearby::nearby(const bool reverse)
+:
+    patchToPatch(reverse),
+    srcSpheres_(),
+    tgtSpheres_()
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::patchToPatches::nearby::~nearby()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/patchToPatch/nearby/nearbyPatchToPatch.H b/src/meshTools/patchToPatch/nearby/nearbyPatchToPatch.H
new file mode 100644
index 0000000000..3c21c3d99f
--- /dev/null
+++ b/src/meshTools/patchToPatch/nearby/nearbyPatchToPatch.H
@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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::patchToPatches::nearby
+
+Description
+    Class to generate patchToPatch coupling geometry. Couples a face to the
+    single nearby opposite face only.
+
+SourceFiles
+    nearby.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef nearbyPatchToPatch_H
+#define nearbyPatchToPatch_H
+
+#include "patchToPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace patchToPatches
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class nearby Declaration
+\*---------------------------------------------------------------------------*/
+
+class nearby
+:
+    public patchToPatch
+{
+protected:
+
+    // Private Member Data
+
+        //- For each source face, the bounding sphere
+        List<Tuple2<point, scalar>> srcSpheres_;
+
+        //- For each target face, the bounding sphere
+        List<Tuple2<point, scalar>> tgtSpheres_;
+
+
+    // Private Member Functions
+
+        //- Get the bound box for a source face
+        virtual treeBoundBox srcBox
+        (
+            const face& srcFace,
+            const pointField& srcPoints,
+            const vectorField& srcPointNormals
+        ) const;
+
+        //- Intersect two faces
+        virtual bool intersectFaces
+        (
+            const primitiveOldTimePatch& srcPatch,
+            const vectorField& srcPointNormals,
+            const vectorField& srcPointNormals0,
+            const primitiveOldTimePatch& tgtPatch,
+            const label srcFacei,
+            const label tgtFacei
+        );
+
+        //- Initialisation
+        virtual void initialise
+        (
+            const primitiveOldTimePatch& srcPatch,
+            const vectorField& srcPointNormals,
+            const vectorField& srcPointNormals0,
+            const primitiveOldTimePatch& tgtPatch
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("nearby");
+
+
+    // Constructors
+
+        //- Construct from components
+        nearby(const bool reverse);
+
+
+    //- Destructor
+    ~nearby();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace patchToPatches
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/patchToPatch/nearest/nearestPatchToPatch.C b/src/meshTools/patchToPatch/nearest/nearestPatchToPatch.C
index e4f933ca1b..d26d77ced1 100644
--- a/src/meshTools/patchToPatch/nearest/nearestPatchToPatch.C
+++ b/src/meshTools/patchToPatch/nearest/nearestPatchToPatch.C
@@ -25,6 +25,11 @@ License
 
 #include "nearestPatchToPatch.H"
 #include "addToRunTimeSelectionTable.H"
+#include "boundSphere.H"
+#include "OFstream.H"
+#include "OBJstream.H"
+#include "vtkWritePolyData.H"
+#include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -40,105 +45,6 @@ namespace patchToPatches
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::treeBoundBox Foam::patchToPatches::nearest::srcBox
-(
-    const face& srcFace,
-    const pointField& srcPoints,
-    const vectorField& srcPointNormals
-) const
-{
-    const treeBoundBox bb(srcPoints, srcFace);
-
-    const point c = bb.midpoint();
-    const scalar l = bb.maxDim();
-
-    return treeBoundBox(c - l*vector::one, c + l*vector::one);
-}
-
-
-bool Foam::patchToPatches::nearest::intersectFaces
-(
-    const primitivePatch& patch,
-    const primitivePatch& otherPatch,
-    const label facei,
-    const label otherFacei,
-    DynamicList<label>& faceOtherFaces,
-    scalar& faceDistance
-) const
-{
-    auto closest = [&patch,&otherPatch]
-    (
-        const label facei,
-        const label otherFacei
-    )
-    {
-        const point& c = patch.faceCentres()[facei];
-        const point& otherC = otherPatch.faceCentres()[otherFacei];
-        const scalar distSqr = magSqr(c - otherC);
-
-        forAll(otherPatch.faceEdges()[otherFacei], otherFaceEdgei)
-        {
-            const label otherEdgei =
-                otherPatch.faceEdges()[otherFacei][otherFaceEdgei];
-
-            point otherNbrC;
-
-            if (otherPatch.edgeFaces()[otherEdgei].size() == 2)
-            {
-                const label facej =
-                    otherPatch.edgeFaces()[otherEdgei]
-                    [otherPatch.edgeFaces()[otherEdgei][0] == otherFacei];
-
-                otherNbrC = otherPatch.faceCentres()[facej];
-            }
-            else
-            {
-                const edge& e = otherPatch.edges()[otherEdgei];
-                const point& p = otherPatch.localPoints()[e[0]];
-                const vector dp = e.vec(otherPatch.localPoints());
-                const vector n = otherPatch.faceNormals()[otherFacei] ^ dp;
-
-                otherNbrC = p + ((tensor::I - 2*sqr(n)) & (otherC - p));
-            }
-
-            if (magSqr(c - otherNbrC) < distSqr)
-            {
-                return false;
-            }
-        }
-
-        return true;
-    };
-
-    if (closest(facei, otherFacei))
-    {
-        const point& c = patch.faceCentres()[facei];
-        const point& otherC = otherPatch.faceCentres()[otherFacei];
-        const scalar distSqr = magSqr(c - otherC);
-
-        if (faceOtherFaces.empty() || faceDistance > distSqr)
-        {
-            faceOtherFaces.clear();
-            faceOtherFaces.append(otherFacei);
-            faceDistance = distSqr;
-        }
-
-        return true;
-    }
-
-    const labelList& otherFaceFaces = otherPatch.faceFaces()[otherFacei];
-    forAll(otherFaceFaces, otherFaceFacei)
-    {
-        if (closest(facei, otherFaceFaces[otherFaceFacei]))
-        {
-            return true;
-        }
-    }
-
-    return false;
-}
-
-
 bool Foam::patchToPatches::nearest::intersectFaces
 (
     const primitiveOldTimePatch& srcPatch,
@@ -149,30 +55,52 @@ bool Foam::patchToPatches::nearest::intersectFaces
     const label tgtFacei
 )
 {
-    const bool srcCouples =
-        intersectFaces
+    if
+    (
+        nearby::intersectFaces
         (
             srcPatch,
+            srcPointNormals,
+            srcPointNormals0,
             tgtPatch,
             srcFacei,
-            tgtFacei,
-            srcLocalTgtFaces_[srcFacei],
-            srcDistances_[srcFacei]
-        );
+            tgtFacei
+        )
+    )
+    {
+        const scalar dSqr =
+            magSqr
+            (
+                srcPatch.faceCentres()[srcFacei]
+              - tgtPatch.faceCentres()[tgtFacei]
+            );
 
-    const bool tgtCouples =
-        intersectFaces
-        (
-            tgtPatch,
-            srcPatch,
-            tgtFacei,
-            srcFacei,
-            tgtLocalSrcFaces_[tgtFacei],
-            tgtDistances_[tgtFacei]
-        );
+        if (dSqr < srcDistances_[srcFacei])
+        {
+            srcDistances_[srcFacei] = dSqr;
+            Swap
+            (
+                srcLocalTgtFaces_[srcFacei].first(),
+                srcLocalTgtFaces_[srcFacei].last()
+            );
+        }
 
-    return srcCouples || tgtCouples;
+        if (dSqr < tgtDistances_[tgtFacei])
+        {
+            tgtDistances_[tgtFacei] = dSqr;
+            Swap
+            (
+                tgtLocalSrcFaces_[tgtFacei].first(),
+                tgtLocalSrcFaces_[tgtFacei].last()
+            );
+        }
 
+        return true;
+    }
+    else
+    {
+        return false;
+    }
 }
 
 
@@ -184,7 +112,7 @@ void Foam::patchToPatches::nearest::initialise
     const primitiveOldTimePatch& tgtPatch
 )
 {
-    patchToPatch::initialise
+    nearby::initialise
     (
         srcPatch,
         srcPointNormals,
@@ -200,13 +128,22 @@ void Foam::patchToPatches::nearest::initialise
 }
 
 
-Foam::labelList Foam::patchToPatches::nearest::trimLocalTgt
+Foam::labelList Foam::patchToPatches::nearest::finaliseLocal
 (
-    const primitiveOldTimePatch& localTgtPatch
+    const primitiveOldTimePatch& srcPatch,
+    const vectorField& srcPointNormals,
+    const vectorField& srcPointNormals0,
+    const primitiveOldTimePatch& tgtPatch
 )
 {
     const labelList newToOldLocalTgtFace =
-        patchToPatch::trimLocalTgt(localTgtPatch);
+        nearby::finaliseLocal
+        (
+            srcPatch,
+            srcPointNormals,
+            srcPointNormals0,
+            tgtPatch
+        );
 
     tgtDistances_ = List<scalar>(tgtDistances_, newToOldLocalTgtFace);
 
@@ -219,6 +156,22 @@ void Foam::patchToPatches::nearest::rDistributeTgt
     const primitiveOldTimePatch& tgtPatch
 )
 {
+    // Keep only the closest opposing face
+    forAll(srcLocalTgtFaces_, srcFacei)
+    {
+        srcLocalTgtFaces_[srcFacei].resize
+        (
+            min(srcLocalTgtFaces_[srcFacei].size(), 1)
+        );
+    }
+    forAll(tgtLocalSrcFaces_, tgtFacei)
+    {
+        tgtLocalSrcFaces_[tgtFacei].resize
+        (
+            min(tgtLocalSrcFaces_[tgtFacei].size(), 1)
+        );
+    }
+
     // Create a list-list of distances to match the addressing
     List<List<scalar>> tgtDistances(tgtLocalSrcFaces_.size());
     forAll(tgtLocalSrcFaces_, tgtFacei)
@@ -230,7 +183,7 @@ void Foam::patchToPatches::nearest::rDistributeTgt
     }
 
     // Let the base class reverse distribute the addressing
-    patchToPatch::rDistributeTgt(tgtPatch);
+    nearby::rDistributeTgt(tgtPatch);
 
     // Reverse distribute the distances
     patchToPatchTools::rDistributeListList
@@ -259,6 +212,132 @@ void Foam::patchToPatches::nearest::rDistributeTgt
 }
 
 
+Foam::label Foam::patchToPatches::nearest::finalise
+(
+    const primitiveOldTimePatch& srcPatch,
+    const vectorField& srcPointNormals,
+    const vectorField& srcPointNormals0,
+    const primitiveOldTimePatch& tgtPatch,
+    const transformer& tgtToSrc
+)
+{
+    // Keep only the closest opposing face
+    forAll(srcLocalTgtFaces_, srcFacei)
+    {
+        srcLocalTgtFaces_[srcFacei].resize
+        (
+            min(srcLocalTgtFaces_[srcFacei].size(), 1)
+        );
+    }
+    forAll(tgtLocalSrcFaces_, tgtFacei)
+    {
+        tgtLocalSrcFaces_[tgtFacei].resize
+        (
+            min(tgtLocalSrcFaces_[tgtFacei].size(), 1)
+        );
+    }
+
+    const label nCouples =
+        nearby::finalise
+        (
+            srcPatch,
+            srcPointNormals,
+            srcPointNormals0,
+            tgtPatch,
+            tgtToSrc
+        );
+
+    if (debug)
+    {
+        auto countNonEmpty = [](const List<DynamicList<label>>& ll)
+        {
+            label result = 0;
+
+            forAll(ll, i)
+            {
+                result += !ll[i].empty();
+            }
+
+            return returnReduce(result, sumOp<label>());
+        };
+
+        Info<< indent
+            << "Coupled " << countNonEmpty(srcLocalTgtFaces_)
+            << "/" << returnReduce(srcLocalTgtFaces_.size(), sumOp<label>())
+            << " source faces and " << countNonEmpty(tgtLocalSrcFaces_)
+            << "/" << returnReduce(tgtLocalSrcFaces_.size(), sumOp<label>())
+            << " target faces" << endl;
+    }
+
+    if (debug && !Pstream::parRun())
+    {
+        auto writeConnections = []
+        (
+            const primitivePatch& patch,
+            const primitivePatch& otherPatch,
+            const bool isSrc,
+            const List<DynamicList<label>>& faceLocalOtherFaces
+        )
+        {
+            const word name =
+                typeName + "_" + (isSrc ? "src" : "tgt") + "Connections";
+
+            OBJstream obj(name + ".obj");
+
+            forAll(faceLocalOtherFaces, facei)
+            {
+                const point& p = patch.faceCentres()[facei];
+                forAll(faceLocalOtherFaces[facei], i)
+                {
+                    const label otherFacei = faceLocalOtherFaces[facei][i];
+                    const point& q = otherPatch.faceCentres()[otherFacei];
+                    obj.write(linePointRef(p, q));
+                }
+            }
+        };
+
+        writeConnections(srcPatch, tgtPatch, true, srcLocalTgtFaces_);
+        writeConnections(tgtPatch, srcPatch, false, tgtLocalSrcFaces_);
+
+        auto writeNotConnected = []
+        (
+            const primitivePatch& patch,
+            const List<DynamicList<label>>& faceLocalOtherFaces,
+            const bool isSrc
+        )
+        {
+            DynamicList<label> unconnected;
+            forAll(faceLocalOtherFaces, facei)
+            {
+                if (faceLocalOtherFaces[facei].empty())
+                {
+                    unconnected.append(facei);
+                }
+            }
+
+            const word name =
+                typeName + "_" + (isSrc ? "src" : "tgt") + "NotConnected";
+
+            vtkWritePolyData::write
+            (
+                name + ".vtk",
+                name,
+                false,
+                patch.localPoints(),
+                labelList(),
+                labelListList(),
+                UIndirectList<face>(patch.localFaces(), unconnected)
+            );
+        };
+
+        writeNotConnected(srcPatch, srcLocalTgtFaces_, true);
+        writeNotConnected(tgtPatch, tgtLocalSrcFaces_, false);
+    }
+
+    return nCouples;
+}
+
+
 Foam::tmpNrc<Foam::List<Foam::DynamicList<Foam::scalar>>>
 Foam::patchToPatches::nearest::srcWeights() const
 {
@@ -307,7 +386,9 @@ Foam::patchToPatches::nearest::tgtWeights() const
 
 Foam::patchToPatches::nearest::nearest(const bool reverse)
 :
-    patchToPatch(reverse)
+    nearby(reverse),
+    srcDistances_(),
+    tgtDistances_()
 {}
 
 
diff --git a/src/meshTools/patchToPatch/nearest/nearestPatchToPatch.H b/src/meshTools/patchToPatch/nearest/nearestPatchToPatch.H
index 3e5d185171..67bdf7039d 100644
--- a/src/meshTools/patchToPatch/nearest/nearestPatchToPatch.H
+++ b/src/meshTools/patchToPatch/nearest/nearestPatchToPatch.H
@@ -36,7 +36,7 @@ SourceFiles
 #ifndef nearestPatchToPatch_H
 #define nearestPatchToPatch_H
 
-#include "patchToPatch.H"
+#include "nearbyPatchToPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,7 +51,7 @@ namespace patchToPatches
 
 class nearest
 :
-    public patchToPatch
+    public nearby
 {
 protected:
 
@@ -66,25 +66,6 @@ protected:
 
     // Private Member Functions
 
-        //- Get the bound box for a source face
-        virtual treeBoundBox srcBox
-        (
-            const face& srcFace,
-            const pointField& srcPoints,
-            const vectorField& srcPointNormals
-        ) const;
-
-        //- Intersect two faces
-        bool intersectFaces
-        (
-            const primitivePatch& patch,
-            const primitivePatch& otherPatch,
-            const label facei,
-            const label otherFacei,
-            DynamicList<label>& faceOtherFaces,
-            scalar& faceDistance
-        ) const;
-
         //- Intersect two faces
         virtual bool intersectFaces
         (
@@ -105,11 +86,15 @@ protected:
             const primitiveOldTimePatch& tgtPatch
         );
 
-        //- Trim the local target patch so that only parts that actually
-        //  intersect the source remain
-        virtual labelList trimLocalTgt
+        //- Finalise the intersection locally. Trims the local target patch
+        //  so that only parts that actually intersect the source remain.
+        //  Returns new-to-old map for trimming target-associated data.
+        virtual labelList finaliseLocal
         (
-            const primitiveOldTimePatch& localTgtPatch
+            const primitiveOldTimePatch& srcPatch,
+            const vectorField& srcPointNormals,
+            const vectorField& srcPointNormals0,
+            const primitiveOldTimePatch& tgtPatch
         );
 
         //- Send data that resulted from an intersection between the source
@@ -117,6 +102,16 @@ protected:
         //  original target processes
         virtual void rDistributeTgt(const primitiveOldTimePatch& tgtPatch);
 
+        //- Finalise the intersection
+        virtual label finalise
+        (
+            const primitiveOldTimePatch& srcPatch,
+            const vectorField& srcPointNormals,
+            const vectorField& srcPointNormals0,
+            const primitiveOldTimePatch& tgtPatch,
+            const transformer& tgtToSrc
+        );
+
         //- For each source face, the coupled target weights
         virtual tmpNrc<List<DynamicList<scalar>>> srcWeights() const;
 
diff --git a/src/meshTools/patchToPatch/patchToPatch/patchToPatch.C b/src/meshTools/patchToPatch/patchToPatch/patchToPatch.C
index 8be111c88f..e6150d7614 100644
--- a/src/meshTools/patchToPatch/patchToPatch/patchToPatch.C
+++ b/src/meshTools/patchToPatch/patchToPatch/patchToPatch.C
@@ -657,13 +657,16 @@ void Foam::patchToPatch::initialise
 }
 
 
-Foam::labelList Foam::patchToPatch::trimLocalTgt
+Foam::labelList Foam::patchToPatch::finaliseLocal
 (
-    const primitiveOldTimePatch& localTgtPatch
+    const primitiveOldTimePatch& srcPatch,
+    const vectorField& srcPointNormals,
+    const vectorField& srcPointNormals0,
+    const primitiveOldTimePatch& tgtPatch
 )
 {
     // Determine which local target faces are actually used
-    boolList oldLocalTgtFaceIsUsed(localTgtPatch.size(), false);
+    boolList oldLocalTgtFaceIsUsed(tgtPatch.size(), false);
     forAll(srcLocalTgtFaces_, srcFacei)
     {
         forAll(srcLocalTgtFaces_[srcFacei], i)
@@ -961,7 +964,13 @@ void Foam::patchToPatch::update
         }
 
         // Trim the local target patch
-        trimLocalTgt(localTTgtPatch);
+        finaliseLocal
+        (
+            srcPatch,
+            srcPointNormals,
+            srcPointNormals,
+            localTTgtPatch
+        );
 
         // Distribute the source patch
         srcMapPtr_ =
diff --git a/src/meshTools/patchToPatch/patchToPatch/patchToPatch.H b/src/meshTools/patchToPatch/patchToPatch/patchToPatch.H
index b251d293b2..681f761638 100644
--- a/src/meshTools/patchToPatch/patchToPatch/patchToPatch.H
+++ b/src/meshTools/patchToPatch/patchToPatch/patchToPatch.H
@@ -215,11 +215,15 @@ protected:
                 const primitiveOldTimePatch& tgtPatch
             );
 
-            //- Trim the local target patch so that only parts that actually
-            //  intersect the source remain
-            virtual labelList trimLocalTgt
+            //- Finalise the intersection locally. Trims the local target patch
+            //  so that only parts that actually intersect the source remain.
+            //  Returns new-to-old map for trimming target-associated data.
+            virtual labelList finaliseLocal
             (
-                const primitiveOldTimePatch& localTgtPatch
+                const primitiveOldTimePatch& srcPatch,
+                const vectorField& srcPointNormals,
+                const vectorField& srcPointNormals0,
+                const primitiveOldTimePatch& tgtPatch
             );
 
             //- Distribute the source patch so that any information needed by
diff --git a/src/meshTools/patchToPatch/rays/raysPatchToPatch.C b/src/meshTools/patchToPatch/rays/raysPatchToPatch.C
index 67cce63806..70a6c9cc7c 100644
--- a/src/meshTools/patchToPatch/rays/raysPatchToPatch.C
+++ b/src/meshTools/patchToPatch/rays/raysPatchToPatch.C
@@ -86,21 +86,30 @@ bool Foam::patchToPatches::rays::intersectFaces
 }
 
 
-Foam::labelList Foam::patchToPatches::rays::trimLocalTgt
+Foam::labelList Foam::patchToPatches::rays::finaliseLocal
 (
-    const primitiveOldTimePatch& localTgtPatch
+    const primitiveOldTimePatch& srcPatch,
+    const vectorField& srcPointNormals,
+    const vectorField& srcPointNormals0,
+    const primitiveOldTimePatch& tgtPatch
 )
 {
     const labelList newToOldLocalTgtFace =
-        patchToPatch::trimLocalTgt(localTgtPatch);
+        patchToPatch::finaliseLocal
+        (
+            srcPatch,
+            srcPointNormals,
+            srcPointNormals0,
+            tgtPatch
+        );
 
     localTgtPatchPtr_.reset
     (
         new PrimitiveOldTimePatch<faceList, pointField>
         (
-            faceList(localTgtPatch, newToOldLocalTgtFace),
-            localTgtPatch.points(),
-            localTgtPatch.points0()
+            faceList(tgtPatch, newToOldLocalTgtFace),
+            tgtPatch.points(),
+            tgtPatch.points0()
         )
     );
 
diff --git a/src/meshTools/patchToPatch/rays/raysPatchToPatch.H b/src/meshTools/patchToPatch/rays/raysPatchToPatch.H
index 9c4ae21aa5..5a01b18962 100644
--- a/src/meshTools/patchToPatch/rays/raysPatchToPatch.H
+++ b/src/meshTools/patchToPatch/rays/raysPatchToPatch.H
@@ -92,11 +92,15 @@ class rays
             const label tgtFacei
         );
 
-        //- Subset the local target patch so that only parts that actually
-        //  intersect the source remain
-        virtual labelList trimLocalTgt
+        //- Finalise the intersection locally. Trims the local target patch
+        //  so that only parts that actually intersect the source remain.
+        //  Returns new-to-old map for trimming target-associated data.
+        virtual labelList finaliseLocal
         (
-            const primitiveOldTimePatch& localTgtPatch
+            const primitiveOldTimePatch& srcPatch,
+            const vectorField& srcPointNormals,
+            const vectorField& srcPointNormals0,
+            const primitiveOldTimePatch& tgtPatch
         );
 
         //- Distribute the source patch so that any information needed by