From a13e00b5c469e21a6b91f783f021419a8fdc8e40 Mon Sep 17 00:00:00 2001
From: Andrew Heather <>
Date: Thu, 7 Mar 2019 20:25:23 +0000
Subject: [PATCH] ENH: AMI - multiple updates

- start of work to create a 1-to-1 face mapping across AMI patches
- faces are inserted according to the AMI addressing based on Horacio's method
- removed 'updated' flag and reworked some demand driven updates
- updated to handle 'walking' through baffles
- use bitSet instead of boolList
- moved update of meshPhi to movePoints() functions at fvPatch level
- moved scaling of areas to movePoints() functions at fvPatch level
- rehomed topology change code to own file
- added warning re: geometry construction

ACMI
- split srcMask into srcMask and srcAreaMask
  - former in range 0-1, and latter has bounding or tol to (1-tol) to avoid
    sigFpe's
---
 src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C |   1 +
 .../polyPatches/polyPatch/polyPatch.C         |   1 +
 .../polyPatches/polyPatch/polyPatch.H         |  14 +
 src/dynamicFvMesh/Make/files                  |   5 +
 .../dynamicMotionSolverFvMesh.C               |   3 +
 .../dynamicMotionSolverFvMeshAMI.C            | 216 ++++++
 .../dynamicMotionSolverFvMeshAMI.H            | 121 ++++
 src/finiteVolume/Make/options                 |   3 +-
 .../cyclicACMI/cyclicACMIFvPatchField.H       |   2 +-
 src/finiteVolume/fvMesh/fvMesh.C              |  21 +-
 src/finiteVolume/fvMesh/fvMesh.H              |   2 +
 .../constraint/cyclicACMI/cyclicACMIFvPatch.C | 143 ++--
 .../constraint/cyclicACMI/cyclicACMIFvPatch.H |  24 +-
 .../constraint/cyclicAMI/cyclicAMIFvPatch.C   |  77 +-
 .../constraint/cyclicAMI/cyclicAMIFvPatch.H   |   3 +
 .../surfaceInterpolation.C                    |   1 -
 .../AMIInterpolation/AMIInterpolation.C       | 142 +++-
 .../AMIInterpolation/AMIInterpolation.H       |  35 +-
 .../AMIInterpolation/AMIInterpolationI.H      |  16 +
 .../AMIMethod/AMIMethod/AMIMethod.C           |  18 +-
 .../AMIMethod/AMIMethod/AMIMethod.H           |  12 +-
 .../AMIMethod/directAMI/directAMI.C           |   8 +-
 .../AMIMethod/directAMI/directAMI.H           |   3 +-
 .../faceAreaWeightAMI/faceAreaWeightAMI.C     | 300 ++++----
 .../faceAreaWeightAMI/faceAreaWeightAMI.H     |  21 +-
 .../AMIMethod/mapNearestAMI/mapNearestAMI.C   |   8 +-
 .../AMIMethod/mapNearestAMI/mapNearestAMI.H   |   3 +-
 .../partialFaceAreaWeightAMI.C                |  16 +-
 .../partialFaceAreaWeightAMI.H                |  30 +-
 .../AMIInterpolation/findNearestMaskedOp.H    |  62 ++
 .../faceAreaIntersect/faceAreaIntersect.C     | 135 ++--
 .../faceAreaIntersect/faceAreaIntersect.H     |  19 +-
 .../faceAreaIntersect/faceAreaIntersectI.H    |   9 +
 .../cyclicACMIPolyPatch/cyclicACMIPolyPatch.C | 379 +++++-----
 .../cyclicACMIPolyPatch/cyclicACMIPolyPatch.H |  54 +-
 .../cyclicACMIPolyPatchI.H                    |  18 +-
 .../cyclicAMIPolyPatch/cyclicAMIPolyPatch.C   | 299 +++++---
 .../cyclicAMIPolyPatch/cyclicAMIPolyPatch.H   |  95 ++-
 .../cyclicAMIPolyPatch/cyclicAMIPolyPatchI.H  |  63 ++
 .../cyclicAMIPolyPatchTopologyChange.C        | 676 ++++++++++++++++++
 src/meshTools/Make/files                      |   1 +
 src/meshTools/Make/options                    |   3 +-
 src/meshTools/polyTopoChange/polyTopoChange.H |  27 +-
 43 files changed, 2435 insertions(+), 654 deletions(-)
 create mode 100644 src/dynamicFvMesh/dynamicMotionSolverFvMeshAMI/dynamicMotionSolverFvMeshAMI.C
 create mode 100644 src/dynamicFvMesh/dynamicMotionSolverFvMeshAMI/dynamicMotionSolverFvMeshAMI.H
 create mode 100644 src/meshTools/AMIInterpolation/AMIInterpolation/findNearestMaskedOp.H
 create mode 100644 src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchTopologyChange.C

diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C b/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C
index 80232d4498..e5818171d5 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C
@@ -45,6 +45,7 @@ void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
         << "Updating addressing and (optional) pointMesh/pointFields" << endl;
 
     // Update boundaryMesh (note that patches themselves already ok)
+//    boundary_.updateMesh(mpm);
     boundary_.updateMesh();
 
     // Update zones
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
index 3e915b3e7e..6db7dc22fc 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
@@ -63,6 +63,7 @@ void Foam::polyPatch::movePoints(PstreamBuffers&, const pointField& p)
     primitivePatch::movePoints(p);
 }
 
+
 void Foam::polyPatch::updateMesh(PstreamBuffers&)
 {
     primitivePatch::clearGeom();
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
index 08bbf0c732..f7f01abb09 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
@@ -55,6 +55,7 @@ namespace Foam
 // Forward declarations
 class polyBoundaryMesh;
 class polyPatch;
+class polyTopoChange;
 class PstreamBuffers;
 
 Ostream& operator<<(Ostream&, const polyPatch&);
@@ -419,6 +420,19 @@ public:
             labelList& rotation
         ) const;
 
+        //- For dynamic mesh cases - return true if this patch will change the
+        //- topology
+        virtual bool changeTopology() const
+        {
+            return false;
+        }
+
+        //- Collect topology changes in a polyTopoChange object
+        virtual bool setTopology(polyTopoChange&)
+        {
+            return false;
+        }
+
 
     // Member operators
 
diff --git a/src/dynamicFvMesh/Make/files b/src/dynamicFvMesh/Make/files
index b8432c07c9..794af6f3d0 100644
--- a/src/dynamicFvMesh/Make/files
+++ b/src/dynamicFvMesh/Make/files
@@ -10,4 +10,9 @@ dynamicMotionSolverListFvMesh/dynamicMotionSolverListFvMesh.C
 simplifiedDynamicFvMesh/simplifiedDynamicFvMeshes.C
 simplifiedDynamicFvMesh/simplifiedDynamicFvMesh.C
 
+
+
+dynamicMotionSolverFvMeshAMI/dynamicMotionSolverFvMeshAMI.C
+
+
 LIB = $(FOAM_LIBBIN)/libdynamicFvMesh
diff --git a/src/dynamicFvMesh/dynamicMotionSolverFvMesh/dynamicMotionSolverFvMesh.C b/src/dynamicFvMesh/dynamicMotionSolverFvMesh/dynamicMotionSolverFvMesh.C
index ba1d636969..5b7a90cddc 100644
--- a/src/dynamicFvMesh/dynamicMotionSolverFvMesh/dynamicMotionSolverFvMesh.C
+++ b/src/dynamicFvMesh/dynamicMotionSolverFvMesh/dynamicMotionSolverFvMesh.C
@@ -92,6 +92,9 @@ const Foam::motionSolver& Foam::dynamicMotionSolverFvMesh::motion() const
 
 bool Foam::dynamicMotionSolverFvMesh::update()
 {
+    // Scan through AMI patches and update
+
+
     fvMesh::movePoints(motionPtr_->newPoints());
 
     volVectorField* Uptr = getObjectPtr<volVectorField>("U");
diff --git a/src/dynamicFvMesh/dynamicMotionSolverFvMeshAMI/dynamicMotionSolverFvMeshAMI.C b/src/dynamicFvMesh/dynamicMotionSolverFvMeshAMI/dynamicMotionSolverFvMeshAMI.C
new file mode 100644
index 0000000000..5de3713a53
--- /dev/null
+++ b/src/dynamicFvMesh/dynamicMotionSolverFvMeshAMI/dynamicMotionSolverFvMeshAMI.C
@@ -0,0 +1,216 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019-2020 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 "dynamicMotionSolverFvMeshAMI.H"
+#include "addToRunTimeSelectionTable.H"
+#include "motionSolver.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "cyclicAMIPolyPatch.H"
+#include "polyTopoChange.H"
+#include "MeshObject.H"
+#include "lduMesh.H"
+
+#include "processorFvPatch.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(dynamicMotionSolverFvMeshAMI, 0);
+    addToRunTimeSelectionTable
+    (
+        dynamicFvMesh,
+        dynamicMotionSolverFvMeshAMI,
+        IOobject
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::dynamicMotionSolverFvMeshAMI::dynamicMotionSolverFvMeshAMI
+(
+    const IOobject& io
+)
+:
+    dynamicFvMesh(io),
+    motionPtr_(motionSolver::New(*this))
+{}
+
+
+Foam::dynamicMotionSolverFvMeshAMI::dynamicMotionSolverFvMeshAMI
+(
+    const IOobject& io,
+    pointField&& points,
+    faceList&& faces,
+    labelList&& allOwner,
+    labelList&& allNeighbour,
+    const bool syncPar
+)
+:
+    dynamicFvMesh
+    (
+        io,
+        std::move(points),
+        std::move(faces),
+        std::move(allOwner),
+        std::move(allNeighbour),
+        syncPar
+    ),
+    motionPtr_(motionSolver::New(*this))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::motionSolver& Foam::dynamicMotionSolverFvMeshAMI::motion() const
+{
+    return *motionPtr_;
+}
+
+
+bool Foam::dynamicMotionSolverFvMeshAMI::update()
+{
+    // Mesh not moved/changed yet
+    moving(false);
+    topoChanging(false);
+
+    if (debug)
+    {
+        for (const fvPatch& fvp : boundary())
+        {
+            if (!isA<processorFvPatch>(fvp))
+            {
+                Info<< "1 --- patch:" << fvp.patch().name()
+                    << " area:" << gSum(fvp.magSf()) << endl;
+            }
+        }
+    }
+
+    pointField newPoints(motionPtr_->curPoints());
+
+    polyBoundaryMesh& pbm = const_cast<polyBoundaryMesh&>(boundaryMesh());
+
+    // Scan all patches and see if we want to apply a mesh topology  update
+    bool changeRequired = false;
+    for (polyPatch& pp : pbm)
+    {
+        DebugInfo
+            << "pre-topology change: patch " << pp.name()
+            << " size:" << returnReduce(pp.size(), sumOp<label>())
+            << " mag(faceAreas):" << gSum(mag(pp.faceAreas())) << endl;
+
+        //changeRequired = pp.changeTopology(newPoints) || changeRequired;
+        changeRequired = pp.changeTopology() || changeRequired;
+    }
+
+    reduce(changeRequired, orOp<bool>());
+
+    if (changeRequired)
+    {
+        polyTopoChange polyTopo(*this);
+
+        // Set new point positions in polyTopo object
+        polyTopo.movePoints(newPoints);
+
+        // Accumulate the patch-based mesh changes on the current mesh
+        // Note:
+        // - updates the AMIs using the new points
+        // - creates a topo change object that removes old added faces and
+        //   adds the new faces
+        for (polyPatch& pp : pbm)
+        {
+            pp.setTopology(polyTopo);
+        }
+
+        // Update geometry
+        // Note
+        // - changeMesh leads to polyMesh::resetPrimitives which will also
+        //   trigger polyBoundaryMesh::updateMesh (init and update) and
+        //   ::calcGeometry (with topoChanging = false)
+        // - BUT: mesh still corresponds to original (non-extended mesh) so
+        //   we want to bypass these calls...
+        // - after changes topoChanging = true
+        autoPtr<mapPolyMesh> map =
+            polyTopo.changeMesh
+            (
+                *this,
+                true       // We will be calling movePoints after this update
+            );
+
+        // Apply topology change - update fv geometry and map fields
+        // - polyMesh::updateMesh
+        //   - fires initUpdateMesh and updateMesh in AMI BCs - called before
+        //     mapFields
+        // - AMI addressing must be up-to-date - used by, e.g. FaceCellWave
+        // - will trigger (again) polyBoundaryMesh::updateMesh (init and update)
+        updateMesh(map());
+
+        // Move points and update derived properties
+        // Note:
+        // - resets face areas based on raw point locations!
+        // - polyBoundaryMesh::updateMesh (init and update)
+        // Note:
+        // - processorPolyPatches will trigger calculation of faceCentres
+        //   (and therefore cell volumes), so need to update faceAreas in
+        //   initMovePoints since proc patches will be evaluated later than
+        //   AMI patches
+        if (map().hasMotionPoints())
+        {
+            movePoints(map().preMotionPoints());
+        }
+    }
+    else
+    {
+        fvMesh::movePoints(newPoints);
+    }
+
+    volVectorField* Uptr = getObjectPtr<volVectorField>("U");
+
+    if (Uptr)
+    {
+        Uptr->correctBoundaryConditions();
+    }
+
+    if (debug)
+    {
+        for (const fvPatch& fvp : boundary())
+        {
+            if (!isA<processorFvPatch>(fvp))
+            {
+                Info<< "2 --- patch:" << fvp.patch().name()
+                    << " area:" << gSum(fvp.magSf()) << endl;
+            }
+        }
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicFvMesh/dynamicMotionSolverFvMeshAMI/dynamicMotionSolverFvMeshAMI.H b/src/dynamicFvMesh/dynamicMotionSolverFvMeshAMI/dynamicMotionSolverFvMeshAMI.H
new file mode 100644
index 0000000000..e21f5623ed
--- /dev/null
+++ b/src/dynamicFvMesh/dynamicMotionSolverFvMeshAMI/dynamicMotionSolverFvMeshAMI.H
@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2020 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::dynamicMotionSolverFvMeshAMI
+
+Description
+    The dynamicMotionSolverFvMeshAMI
+
+SourceFiles
+    dynamicMotionSolverFvMeshAMI.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dynamicMotionSolverFvMeshAMI_H
+#define dynamicMotionSolverFvMeshAMI_H
+
+#include "dynamicFvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class motionSolver;
+
+/*---------------------------------------------------------------------------*\
+               Class dynamicMotionSolverFvMeshAMI Declaration
+\*---------------------------------------------------------------------------*/
+
+class dynamicMotionSolverFvMeshAMI
+:
+    public dynamicFvMesh
+{
+    // Private data
+
+        autoPtr<motionSolver> motionPtr_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        dynamicMotionSolverFvMeshAMI
+        (
+            const dynamicMotionSolverFvMeshAMI&
+        ) = delete;
+
+        //- No copy assignment
+        void operator=(const dynamicMotionSolverFvMeshAMI&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("dynamicMotionSolverFvMeshAMI");
+
+
+    // Constructors
+
+        //- Construct from IOobject
+        dynamicMotionSolverFvMeshAMI(const IOobject& io);
+
+        //- Construct from components without boundary.
+        //  Boundary is added using addFvPatches() member function
+        dynamicMotionSolverFvMeshAMI
+        (
+            const IOobject& io,
+            pointField&& points,
+            faceList&& faces,
+            labelList&& allOwner,
+            labelList&& allNeighbour,
+            const bool syncPar = true
+        );
+
+
+    //- Destructor
+    virtual ~dynamicMotionSolverFvMeshAMI() = default;
+
+
+    // Member Functions
+
+        //- Return the motionSolver
+        const motionSolver& motion() const;
+
+        //- Update the mesh for both mesh motion and topology change
+        virtual bool update();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/Make/options b/src/finiteVolume/Make/options
index 851f94b6b3..ea11ab0d45 100644
--- a/src/finiteVolume/Make/options
+++ b/src/finiteVolume/Make/options
@@ -1,7 +1,8 @@
 EXE_INC = \
     -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 LIB_LIBS = \
     -lOpenFOAM \
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H
index eae2dd49da..db093b3940 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H
@@ -173,7 +173,7 @@ public:
 
             //- Return true if this patch field fixes a value
             //  Needed to check if a level has to be specified while solving
-            //  Poissons equations
+            //  Poisson equations
             virtual bool fixesValue() const
             {
                 const scalarField& mask =
diff --git a/src/finiteVolume/fvMesh/fvMesh.C b/src/finiteVolume/fvMesh/fvMesh.C
index 9aff27408a..e02ae63dcc 100644
--- a/src/finiteVolume/fvMesh/fvMesh.C
+++ b/src/finiteVolume/fvMesh/fvMesh.C
@@ -726,6 +726,8 @@ void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap)
 
 Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
 {
+    DebugInFunction << endl;
+
     // Grab old time volumes if the time has been incremented
     // This will update V0, V00
     if (curTimeIndex_ < time().timeIndex())
@@ -733,6 +735,14 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
         storeOldVol(V());
     }
 
+
+    // Move the polyMesh and set the mesh motion fluxes to the swept-volumes
+
+    scalar rDeltaT = 1.0/time().deltaTValue();
+
+    tmp<scalarField> tsweptVols = polyMesh::movePoints(p);
+    scalarField& sweptVols = tsweptVols.ref();
+
     if (!phiPtr_)
     {
         // Create mesh motion flux
@@ -761,14 +771,6 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
     }
 
     surfaceScalarField& phi = *phiPtr_;
-
-    // Move the polyMesh and set the mesh motion fluxes to the swept-volumes
-
-    scalar rDeltaT = 1.0/time().deltaTValue();
-
-    tmp<scalarField> tsweptVols = polyMesh::movePoints(p);
-    scalarField& sweptVols = tsweptVols.ref();
-
     phi.primitiveFieldRef() =
         scalarField::subField(sweptVols, nInternalFaces());
     phi.primitiveFieldRef() *= rDeltaT;
@@ -776,7 +778,6 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
     const fvPatchList& patches = boundary();
 
     surfaceScalarField::Boundary& phibf = phi.boundaryFieldRef();
-
     forAll(patches, patchi)
     {
         phibf[patchi] = patches[patchi].patchSlice(sweptVols);
@@ -805,6 +806,8 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
 
 void Foam::fvMesh::updateMesh(const mapPolyMesh& mpm)
 {
+    DebugInFunction << endl;
+
     // Update polyMesh. This needs to keep volume existent!
     polyMesh::updateMesh(mpm);
 
diff --git a/src/finiteVolume/fvMesh/fvMesh.H b/src/finiteVolume/fvMesh/fvMesh.H
index 53e57c7f68..0e35ef9f9d 100644
--- a/src/finiteVolume/fvMesh/fvMesh.H
+++ b/src/finiteVolume/fvMesh/fvMesh.H
@@ -91,6 +91,8 @@ class fvMesh
     public fvSolution,
     public data
 {
+protected:
+
     // Private data
 
         //- Boundary mesh
diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicACMI/cyclicACMIFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicACMI/cyclicACMIFvPatch.C
index c4427d5484..23a12796ad 100644
--- a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicACMI/cyclicACMIFvPatch.C
+++ b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicACMI/cyclicACMIFvPatch.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,9 +27,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "cyclicACMIFvPatch.H"
-#include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
 #include "transform.H"
+#include "surfaceFields.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -41,45 +43,13 @@ namespace Foam
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
-void Foam::cyclicACMIFvPatch::updateAreas() const
+void Foam::cyclicACMIFvPatch::resetPatchAreas(const fvPatch& fvp) const
 {
-    if (cyclicACMIPolyPatch_.updated())
-    {
-        if (debug)
-        {
-            Pout<< "cyclicACMIFvPatch::updateAreas() : updating fv areas for "
-                << name() << " and " << this->nonOverlapPatch().name()
-                << endl;
-        }
+    const_cast<vectorField&>(fvp.Sf()) = fvp.patch().faceAreas();
+    const_cast<scalarField&>(fvp.magSf()) = mag(fvp.patch().faceAreas());
 
-        // owner couple
-        const_cast<vectorField&>(Sf()) = patch().faceAreas();
-        const_cast<scalarField&>(magSf()) = mag(patch().faceAreas());
-
-        // owner non-overlapping
-        const fvPatch& nonOverlapPatch = this->nonOverlapPatch();
-        const_cast<vectorField&>(nonOverlapPatch.Sf()) =
-            nonOverlapPatch.patch().faceAreas();
-        const_cast<scalarField&>(nonOverlapPatch.magSf()) =
-            mag(nonOverlapPatch.patch().faceAreas());
-
-        // neighbour couple
-        const cyclicACMIFvPatch& nbrACMI = neighbPatch();
-        const_cast<vectorField&>(nbrACMI.Sf()) =
-            nbrACMI.patch().faceAreas();
-        const_cast<scalarField&>(nbrACMI.magSf()) =
-            mag(nbrACMI.patch().faceAreas());
-
-        // neighbour non-overlapping
-        const fvPatch& nbrNonOverlapPatch = nbrACMI.nonOverlapPatch();
-        const_cast<vectorField&>(nbrNonOverlapPatch.Sf()) =
-            nbrNonOverlapPatch.patch().faceAreas();
-        const_cast<scalarField&>(nbrNonOverlapPatch.magSf()) =
-            mag(nbrNonOverlapPatch.patch().faceAreas());
-
-        // set the updated flag
-        cyclicACMIPolyPatch_.setUpdated(false);
-    }
+    DebugPout
+        << fvp.patch().name() << " area:" << sum(fvp.magSf()) << endl;
 }
 
 
@@ -100,13 +70,13 @@ void Foam::cyclicACMIFvPatch::makeWeights(scalarField& w) const
             )
         );
 
-        scalar tol = cyclicACMIPolyPatch::tolerance();
+        const scalar tol = cyclicACMIPolyPatch::tolerance();
 
 
         forAll(deltas, facei)
         {
-            scalar di = deltas[facei];
-            scalar dni = nbrDeltas[facei];
+            scalar di = mag(deltas[facei]);
+            scalar dni = mag(nbrDeltas[facei]);
 
             if (dni < tol)
             {
@@ -147,7 +117,6 @@ Foam::tmp<Foam::vectorField> Foam::cyclicACMIFvPatch::delta() const
 
         vectorField nbrPatchD(interpolate(nbrPatch.coupledFvPatch::delta()));
 
-
         tmp<vectorField> tpdv(new vectorField(patchD.size()));
         vectorField& pdv = tpdv.ref();
 
@@ -201,4 +170,94 @@ Foam::tmp<Foam::labelField> Foam::cyclicACMIFvPatch::internalFieldTransfer
 }
 
 
+void Foam::cyclicACMIFvPatch::movePoints()
+{
+    if (!cyclicACMIPolyPatch_.owner())
+    {
+        return;
+    }
+
+    // Set the patch face areas to be consistent with the changes made at the
+    // polyPatch level
+
+    const fvPatch& nonOverlapPatch = this->nonOverlapPatch();
+    const cyclicACMIFvPatch& nbrACMI = neighbPatch();
+    const fvPatch& nbrNonOverlapPatch = nbrACMI.nonOverlapPatch();
+
+    resetPatchAreas(*this);
+    resetPatchAreas(nonOverlapPatch);
+    resetPatchAreas(nbrACMI);
+    resetPatchAreas(nbrNonOverlapPatch);
+
+    // Scale the mesh flux
+
+    const labelListList& newSrcAddr = AMI().srcAddress();
+    const labelListList& newTgtAddr = AMI().tgtAddress();
+
+    const fvMesh& mesh = boundaryMesh().mesh();
+    surfaceScalarField& meshPhi = const_cast<fvMesh&>(mesh).setPhi();
+    surfaceScalarField::Boundary& meshPhiBf = meshPhi.boundaryFieldRef();
+
+    // Note: phip and phiNonOverlapp will be different sizes if new faces
+    // have been added
+    scalarField& phip = meshPhiBf[cyclicACMIPolyPatch_.index()];
+    scalarField& phiNonOverlapp =
+        meshPhiBf[nonOverlapPatch.patch().index()];
+
+    forAll(phip, facei)
+    {
+        if (newSrcAddr[facei].empty())
+        {
+            // AMI patch with no connection to other coupled faces
+            phip[facei] = 0.0;
+        }
+        else
+        {
+            // Scale the mesh flux according to the area fraction
+            const face& fAMI = cyclicACMIPolyPatch_.localFaces()[facei];
+
+            // Note: using raw point locations to calculate the geometric
+            // area - faces areas are currently scaled (decoupled from
+            // mesh points)
+            const scalar geomArea = fAMI.mag(cyclicACMIPolyPatch_.localPoints());
+            phip[facei] *= magSf()[facei]/geomArea;
+        }
+    }
+
+    forAll(phiNonOverlapp, facei)
+    {
+        const scalar w = 1.0 - cyclicACMIPolyPatch_.srcMask()[facei];
+        phiNonOverlapp[facei] *= w;
+    }
+
+    scalarField& nbrPhip = meshPhiBf[nbrACMI.patch().index()];
+    scalarField& nbrPhiNonOverlapp =
+        meshPhiBf[nbrNonOverlapPatch.patch().index()];
+
+    forAll(nbrPhip, facei)
+    {
+        if (newTgtAddr[facei].empty())
+        {
+            nbrPhip[facei] = 0.0;
+        }
+        else
+        {
+            const face& fAMI = nbrACMI.patch().localFaces()[facei];
+
+            // Note: using raw point locations to calculate the geometric
+            // area - faces areas are currently scaled (decoupled from
+            // mesh points)
+            const scalar geomArea = fAMI.mag(nbrACMI.patch().localPoints());
+            nbrPhip[facei] *= nbrACMI.magSf()[facei]/geomArea;
+        }
+    }
+
+    forAll(nbrPhiNonOverlapp, facei)
+    {
+        const scalar w = 1.0 - cyclicACMIPolyPatch_.tgtMask()[facei];
+        nbrPhiNonOverlapp[facei] *= w;
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicACMI/cyclicACMIFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicACMI/cyclicACMIFvPatch.H
index 037df53fe3..c3aed02df7 100644
--- a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicACMI/cyclicACMIFvPatch.H
+++ b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicACMI/cyclicACMIFvPatch.H
@@ -65,12 +65,15 @@ protected:
 
     // Protected Member functions
 
-        //- Update the patch areas after AMI update
-        void updateAreas() const;
+        //- Helper function to reset the FV patch areas from the primitive patch
+        void resetPatchAreas(const fvPatch& fvp) const;
 
         //- Make patch weighting factors
         void makeWeights(scalarField&) const;
 
+        //- Correct patches after moving points
+        virtual void movePoints();
+
 
 public:
 
@@ -134,12 +137,7 @@ public:
             //- Return a reference to the AMI interpolator
             virtual const AMIPatchToPatchInterpolation& AMI() const
             {
-                const AMIPatchToPatchInterpolation& AMI =
-                    cyclicACMIPolyPatch_.AMI();
-
-                updateAreas();
-
-                return AMI;
+                return cyclicACMIPolyPatch_.AMI();
             }
 
             //- Are the cyclic planes parallel
@@ -169,8 +167,8 @@ public:
             }
 
             //- Return true if this patch is coupled. This is equivalent
-            //  to the coupledPolyPatch::coupled() if parallel running or
-            //  both sides present, false otherwise
+            //- to the coupledPolyPatch::coupled() if parallel running or
+            //- both sides present, false otherwise
             virtual bool coupled() const;
 
             //- Return delta (P to N) vectors across coupled patch
@@ -183,8 +181,6 @@ public:
                 const Field<Type>& fld
             ) const
             {
-                updateAreas();
-
                 return
                     cyclicACMIPolyPatch_.cyclicAMIPolyPatch::interpolate
                     (
@@ -192,7 +188,7 @@ public:
                     );
             }
 
-            //- Interpolate (make sure to have uptodate areas)
+            //- Interpolate (make sure to have up-to-date areas)
             template<class Type>
             tmp<Field<Type>> interpolate
             (
@@ -206,7 +202,7 @@ public:
         // Interface transfer functions
 
             //- Return the values of the given internal data adjacent to
-            //  the interface as a field
+            //- the interface as a field
             virtual tmp<labelField> interfaceInternalField
             (
                 const labelUList& internalData
diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C
index 356c61db10..a9c03e340c 100644
--- a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C
+++ b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C
@@ -29,6 +29,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
 #include "transform.H"
+#include "surfaceFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -82,8 +83,9 @@ void Foam::cyclicAMIFvPatch::makeWeights(scalarField& w) const
 
         forAll(deltas, facei)
         {
-            scalar di = deltas[facei];
-            scalar dni = nbrDeltas[facei];
+            // Note use of mag
+            scalar di = mag(deltas[facei]);
+            scalar dni = mag(nbrDeltas[facei]);
 
             w[facei] = dni/(di + dni);
         }
@@ -174,4 +176,75 @@ Foam::tmp<Foam::labelField> Foam::cyclicAMIFvPatch::internalFieldTransfer
 }
 
 
+void Foam::cyclicAMIFvPatch::movePoints()
+{
+    if (!owner() || !cyclicAMIPolyPatch_.createAMIFaces())
+    {
+        // Only manipulating patch face areas and mesh motion flux if the AMI
+        // creates additional faces
+        return;
+    }
+
+    // Update face data based on values set by the AMI manipulations
+    const_cast<vectorField&>(Sf()) = cyclicAMIPolyPatch_.faceAreas();
+    const_cast<vectorField&>(Cf()) = cyclicAMIPolyPatch_.faceCentres();
+    const_cast<scalarField&>(magSf()) = mag(Sf());
+
+    const cyclicAMIFvPatch& nbr = neighbPatch();
+    const_cast<vectorField&>(nbr.Sf()) = nbr.cyclicAMIPatch().faceAreas();
+    const_cast<vectorField&>(nbr.Cf()) = nbr.cyclicAMIPatch().faceCentres();
+    const_cast<scalarField&>(nbr.magSf()) = mag(nbr.Sf());
+
+
+    // Set consitent mesh motion flux
+    // TODO: currently maps src mesh flux to tgt - update to
+    // src = src + mapped(tgt) and tgt = tgt + mapped(src)?
+
+    const fvMesh& mesh = boundaryMesh().mesh();
+    surfaceScalarField& meshPhi = const_cast<fvMesh&>(mesh).setPhi();
+    surfaceScalarField::Boundary& meshPhiBf = meshPhi.boundaryFieldRef();
+
+    if (cyclicAMIPolyPatch_.owner())
+    {
+        scalarField& phip = meshPhiBf[patch().index()];
+        forAll(phip, facei)
+        {
+            const face& f = cyclicAMIPolyPatch_.localFaces()[facei];
+
+            // Note: using raw point locations to calculate the geometric
+            // area - faces areas are currently scaled by the AMI weights
+            // (decoupled from mesh points)
+            const scalar geomArea = f.mag(cyclicAMIPolyPatch_.localPoints());
+
+            const scalar scaledArea = magSf()[facei];
+            phip[facei] *= scaledArea/geomArea;
+        }
+
+        scalarField srcMeshPhi(phip);
+        if (Pstream::parRun())
+        {
+            AMI().srcMap().distribute(srcMeshPhi);
+        }
+
+        const labelListList& tgtToSrcAddr = AMI().tgtAddress();
+        scalarField& nbrPhip = meshPhiBf[nbr.index()];
+
+        forAll(tgtToSrcAddr, tgtFacei)
+        {
+            // Note: now have 1-to-1 mapping so tgtToSrcAddr[tgtFacei] is size 1
+            const label srcFacei = tgtToSrcAddr[tgtFacei][0];
+            nbrPhip[tgtFacei] = -srcMeshPhi[srcFacei];
+        }
+
+        DebugInfo
+            << "patch:" << patch().name()
+            << " sum(area):" << gSum(magSf())
+            << " min(mag(faceAreas):" << gMin(magSf())
+            << " sum(meshPhi):" << gSum(phip) << nl
+            << " sum(nbrMeshPhi):" << gSum(nbrPhip) << nl
+            << endl;
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.H
index 2b64389598..4245dd4be8 100644
--- a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.H
+++ b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.H
@@ -68,6 +68,9 @@ protected:
         //- Make patch weighting factors
         void makeWeights(scalarField&) const;
 
+        //- Correct patches after moving points
+        virtual void movePoints();
+
 
 public:
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
index 2b31e5c66b..c716e1456a 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
@@ -172,7 +172,6 @@ void Foam::surfaceInterpolation::makeWeights() const
 
     // ... and reference to the internal field of the weighting factors
     scalarField& w = weights.primitiveFieldRef();
-
     forAll(owner, facei)
     {
         // Note: mag in the dot-product.
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
index b7c2978294..090a3812a2 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
@@ -33,6 +33,8 @@ License
 #include "flipOp.H"
 #include "profiling.H"
 
+#define DEBUGAMI(msg){Pout<< "[" << __FILE__ << ":" << __LINE__ << "]: " << #msg << "=" << msg << endl;}
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 template<class SourcePatch, class TargetPatch>
@@ -49,10 +51,12 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolationMethodNames_
     { interpolationMethod::imPartialFaceAreaWeight, "partialFaceAreaWeightAMI" }
 });
 
+
 template<class SourcePatch, class TargetPatch>
 bool Foam::AMIInterpolation<SourcePatch, TargetPatch>::cacheIntersections_ =
     false;
 
+
 template<class SourcePatch, class TargetPatch>
 template<class Patch>
 Foam::tmp<Foam::scalarField>
@@ -264,7 +268,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
         }
     }
 
-
     // Agglomerate weights and indices
     if (targetMapPtr.valid())
     {
@@ -298,7 +301,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
         //           the slots are equal to face indices.
         // A mapDistribute has:
         // - a subMap : these are face indices
-        // - a constructMap : these are from 'transferred-date' to slots
+        // - a constructMap : these are from 'transferred-data' to slots
 
         labelListList tgtSubMap(Pstream::nProcs());
 
@@ -616,10 +619,12 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
     srcAddress_(),
     srcWeights_(),
     srcWeightsSum_(),
+    srcCentroids_(),
     tgtMagSf_(),
     tgtAddress_(),
     tgtWeights_(),
     tgtWeightsSum_(),
+    tgtCentroids_(),
     triMode_(triMode),
     srcMapPtr_(nullptr),
     tgtMapPtr_(nullptr)
@@ -649,10 +654,12 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
     srcAddress_(),
     srcWeights_(),
     srcWeightsSum_(),
+    srcCentroids_(),
     tgtMagSf_(),
     tgtAddress_(),
     tgtWeights_(),
     tgtWeightsSum_(),
+    tgtCentroids_(),
     triMode_(triMode),
     srcMapPtr_(nullptr),
     tgtMapPtr_(nullptr)
@@ -683,10 +690,12 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
     srcAddress_(),
     srcWeights_(),
     srcWeightsSum_(),
+    srcCentroids_(),
     tgtMagSf_(),
     tgtAddress_(),
     tgtWeights_(),
     tgtWeightsSum_(),
+    tgtCentroids_(),
     triMode_(triMode),
     srcMapPtr_(nullptr),
     tgtMapPtr_(nullptr)
@@ -717,10 +726,12 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
     srcAddress_(),
     srcWeights_(),
     srcWeightsSum_(),
+    srcCentroids_(),
     tgtMagSf_(),
     tgtAddress_(),
     tgtWeights_(),
     tgtWeightsSum_(),
+    tgtCentroids_(),
     triMode_(triMode),
     srcMapPtr_(nullptr),
     tgtMapPtr_(nullptr)
@@ -834,13 +845,6 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-template<class SourcePatch, class TargetPatch>
-Foam::AMIInterpolation<SourcePatch, TargetPatch>::~AMIInterpolation()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class SourcePatch, class TargetPatch>
@@ -932,6 +936,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
         (
             srcAddress_,
             srcWeights_,
+            srcCentroids_,
             tgtAddress_,
             tgtWeights_
         );
@@ -1009,9 +1014,16 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
         AMIPtr->normaliseWeights(true, *this);
 
         // Cache maps and reset addresses
-        List<Map<label>> cMap;
-        srcMapPtr_.reset(new mapDistribute(globalSrcFaces, tgtAddress_, cMap));
-        tgtMapPtr_.reset(new mapDistribute(globalTgtFaces, srcAddress_, cMap));
+        List<Map<label>> cMapSrc;
+        srcMapPtr_.reset
+        (
+            new mapDistribute(globalSrcFaces, tgtAddress_, cMapSrc)
+        );
+        List<Map<label>> cMapTgt;
+        tgtMapPtr_.reset
+        (
+            new mapDistribute(globalTgtFaces, srcAddress_, cMapTgt)
+        );
     }
     else
     {
@@ -1033,6 +1045,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
         (
             srcAddress_,
             srcWeights_,
+            srcCentroids_,
             tgtAddress_,
             tgtWeights_
         );
@@ -1056,6 +1069,42 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
 }
 
 
+template<class SourcePatch, class TargetPatch>
+void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
+(
+    autoPtr<mapDistribute>&& srcToTgtMap,
+    autoPtr<mapDistribute>&& tgtToSrcMap,
+    labelListList&& srcAddress,
+    scalarListList&& srcWeights,
+    labelListList&& tgtAddress,
+    scalarListList&& tgtWeights
+)
+{
+    DebugInFunction<< endl;
+
+    srcAddress_.transfer(srcAddress);
+    srcWeights_.transfer(srcWeights);
+    tgtAddress_.transfer(tgtAddress);
+    tgtWeights_.transfer(tgtWeights);
+
+    // Reset the sums of the weights
+    srcWeightsSum_.setSize(srcWeights_.size());
+    forAll(srcWeights_, facei)
+    {
+        srcWeightsSum_[facei] = sum(srcWeights_[facei]);
+    }
+
+    tgtWeightsSum_.setSize(tgtWeights_.size());
+    forAll(tgtWeights_, facei)
+    {
+        tgtWeightsSum_[facei] = sum(tgtWeights_[facei]);
+    }
+
+    srcMapPtr_ = srcToTgtMap;
+    tgtMapPtr_ = tgtToSrcMap;
+}
+
+
 template<class SourcePatch, class TargetPatch>
 void Foam::AMIInterpolation<SourcePatch, TargetPatch>::append
 (
@@ -1397,7 +1446,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToSource
             FatalErrorInFunction
                 << "Employing default values when sum of weights falls below "
                 << lowWeightCorrection_
-                << " but supplied default field size is not equal to target "
+                << " but supplied default field size is not equal to source "
                 << "patch size" << nl
                 << "    default values = " << defaultValues.size() << nl
                 << "    source patch   = " << srcAddress_.size() << nl
@@ -1693,6 +1742,73 @@ const
 }
 
 
+template<class SourcePatch, class TargetPatch>
+bool Foam::AMIInterpolation<SourcePatch, TargetPatch>::checkSymmetricWeights
+(
+    const bool log
+) const
+{
+    if (Pstream::parRun() && (singlePatchProc_ == -1))
+    {
+        Log << "Checks only valid for serial running (currently)" << endl;
+
+        return true;
+    }
+
+    bool symmetricSrc = true;
+
+    Log << "    Checking for missing src face in tgt lists" << nl;
+
+    forAll(srcAddress_, srcFacei)
+    {
+        const labelList& tgtIds = srcAddress_[srcFacei];
+        for (const label tgtFacei : tgtIds)
+        {
+            if (!tgtAddress_[tgtFacei].found(srcFacei))
+            {
+                symmetricSrc = false;
+
+                Log << "       srcFacei:" << srcFacei
+                    << " not found in tgtToSrc list for tgtFacei:"
+                    << tgtFacei << nl;
+            }
+        }
+    }
+
+    if (symmetricSrc)
+    {
+        Log << "    - symmetric" << endl;
+    }
+
+    bool symmetricTgt = true;
+
+    Log << "    Checking for missing tgt face in src lists" << nl;
+
+    forAll(tgtAddress_, tgtFacei)
+    {
+        const labelList& srcIds = tgtAddress_[tgtFacei];
+        for (const label srcFacei : srcIds)
+        {
+            if (!srcAddress_[srcFacei].found(tgtFacei))
+            {
+                symmetricTgt = false;
+
+                Log << "       tgtFacei:" << tgtFacei
+                    << " not found in srcToTgt list for srcFacei:"
+                    << srcFacei << nl;
+            }
+        }
+    }
+
+    if (symmetricTgt)
+    {
+        Log << "    - symmetric" << endl;
+    }
+
+    return symmetricSrc && symmetricTgt;
+}
+
+
 template<class SourcePatch, class TargetPatch>
 void Foam::AMIInterpolation<SourcePatch, TargetPatch>::writeFaceConnectivity
 (
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
index fcf418ff6b..2e57a3c7dd 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
@@ -61,6 +61,7 @@ SourceFiles
 #include "globalIndex.H"
 #include "ops.H"
 #include "Enum.H"
+#include "pointList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -146,6 +147,9 @@ private:
             //- Sum of weights of target faces per source face
             scalarField srcWeightsSum_;
 
+            //- Centroid of target faces per source face
+            pointListList srcCentroids_;
+
 
         // Target patch
 
@@ -161,6 +165,9 @@ private:
             //- Sum of weights of source faces per target face
             scalarField tgtWeightsSum_;
 
+            //- Centroid of source faces per target face
+            pointListList tgtCentroids_;
+
 
         //- Face triangulation mode
         const faceAreaIntersect::triangulationMode triMode_;
@@ -347,7 +354,7 @@ public:
 
 
     //- Destructor
-    ~AMIInterpolation();
+    ~AMIInterpolation() = default;
 
     // Typedef to SourcePatch type this AMIInterpolation is instantiated on
     typedef SourcePatch sourcePatchType;
@@ -399,6 +406,12 @@ public:
                 //- patch weights (i.e. the sum before normalisation)
                 inline scalarField& srcWeightsSum();
 
+                //- Return const access to source patch face centroids
+                inline const pointListList& srcCentroids() const;
+
+                //- Return access to source patch face centroids
+                inline pointListList& srcCentroids();
+
                 //- Source map pointer - valid only if singlePatchProc = -1
                 //- This gets source data into a form to be consumed by
                 //- tgtAddress, tgtWeights
@@ -448,6 +461,22 @@ public:
                 const TargetPatch& tgtPatch
             );
 
+            void update
+            (
+                autoPtr<mapDistribute>&& srcToTgtMap,
+                autoPtr<mapDistribute>&& tgtToSrcMap,
+                labelListList&& srcAddress,
+                scalarListList&& srcWeights,
+                labelListList&& tgtAddress,
+                scalarListList&& tgtWeights
+            );
+
+            void setAreas(const scalarList& srcMagSf, const scalarList& tgtMagSf)
+            {
+                srcMagSf_ = srcMagSf;
+                tgtMagSf_ = tgtMagSf;
+            }
+
             //- Append additional addressing and weights
             void append
             (
@@ -582,6 +611,10 @@ public:
 
         // Checks
 
+            //- Check if src addresses are present in tgt addresses and
+            //- viceversa
+            bool checkSymmetricWeights(const bool log) const;
+
             //- Write face connectivity as OBJ file
             void writeFaceConnectivity
             (
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H
index 9cd5e51ff1..027a504521 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H
@@ -115,6 +115,22 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcWeightsSum()
 }
 
 
+template<class SourcePatch, class TargetPatch>
+inline const Foam::pointListList&
+Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcCentroids() const
+{
+    return srcCentroids_;
+}
+
+
+template<class SourcePatch, class TargetPatch>
+inline Foam::pointListList&
+Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcCentroids()
+{
+    return srcCentroids_;
+}
+
+
 template<class SourcePatch, class TargetPatch>
 inline const Foam::mapDistribute&
 Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcMap() const
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C
index 3b1f06d405..5a8f11ec74 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C
@@ -31,6 +31,8 @@ License
 #include "mapDistribute.H"
 #include "unitConversion.H"
 
+#include "findNearestMaskedOp.H"
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class SourcePatch, class TargetPatch>
@@ -231,17 +233,25 @@ void Foam::AMIMethod<SourcePatch, TargetPatch>::resetTree()
 template<class SourcePatch, class TargetPatch>
 Foam::label Foam::AMIMethod<SourcePatch, TargetPatch>::findTargetFace
 (
-    const label srcFacei
+    const label srcFacei,
+    const UList<label>& excludeFaces,
+    const label srcFacePti
 ) const
 {
     label targetFacei = -1;
 
     const pointField& srcPts = srcPatch_.points();
     const face& srcFace = srcPatch_[srcFacei];
-    const point srcPt = srcFace.centre(srcPts);
-    const scalar srcFaceArea = srcMagSf_[srcFacei];
 
-    pointIndexHit sample = treePtr_->findNearest(srcPt, 10.0*srcFaceArea);
+    findNearestMaskedOp<TargetPatch> fnOp(*treePtr_, excludeFaces);
+
+    boundBox bb(srcPts, srcFace, false);
+
+    const point srcPt =
+        srcFacePti == -1 ? bb.centre() : srcPts[srcFace[srcFacePti]];
+
+    const pointIndexHit sample =
+        treePtr_->findNearest(srcPt, magSqr(bb.max() - bb.centre()), fnOp);
 
     if (sample.hit())
     {
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.H
index bf569aad56..e019c384ce 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.H
@@ -45,6 +45,7 @@ SourceFiles
 #include "treeDataPrimitivePatch.H"
 #include "treeBoundBoxList.H"
 #include "runTimeSelectionTables.H"
+#include "pointList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -101,7 +102,7 @@ protected:
         List<scalar> tgtMagSf_;
 
         //- Labels of faces that are not overlapped by any target faces
-        //- (should be empty for correct functioning)
+        //- (should be empty for correct functioning for fully covered AMIs)
         labelList srcNonOverlap_;
 
         //- Octree used to find face seeds
@@ -145,8 +146,12 @@ protected:
             //- Reset the octree for the target patch face search
             void resetTree();
 
-            //- Find face on target patch that overlaps source face
-            label findTargetFace(const label srcFacei) const;
+            label findTargetFace
+            (
+                const label srcFacei,
+                const UList<label>& excludeFaces = UList<label>::null(),
+                const label srcFacePti = -1
+            ) const;
 
             //- Add faces neighbouring facei to the ID list
             void appendNbrFaces
@@ -251,6 +256,7 @@ public:
             (
                 labelListList& srcAddress,
                 scalarListList& srcWeights,
+                pointListList& srcCentroids,
                 labelListList& tgtAddress,
                 scalarListList& tgtWeights,
                 label srcFacei = -1,
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C
index 894b33fe0f..14989fdfbe 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C
@@ -205,13 +205,6 @@ Foam::directAMI<SourcePatch, TargetPatch>::directAMI
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-template<class SourcePatch, class TargetPatch>
-Foam::directAMI<SourcePatch, TargetPatch>::~directAMI()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class SourcePatch, class TargetPatch>
@@ -219,6 +212,7 @@ void Foam::directAMI<SourcePatch, TargetPatch>::calculate
 (
     labelListList& srcAddress,
     scalarListList& srcWeights,
+    pointListList& srcCentroids,
     labelListList& tgtAddress,
     scalarListList& tgtWeights,
     label srcFacei,
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H
index a451cad7da..a07aa4f617 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H
@@ -118,7 +118,7 @@ public:
 
 
     //- Destructor
-    virtual ~directAMI();
+    virtual ~directAMI() = default;
 
 
     // Member Functions
@@ -130,6 +130,7 @@ public:
             (
                 labelListList& srcAddress,
                 scalarListList& srcWeights,
+                pointListList& srcCentroids,
                 labelListList& tgtAddress,
                 scalarListList& tgtWeights,
                 label srcFacei = -1,
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.C
index 3cbe9cd65f..f770dfbb58 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.C
@@ -28,6 +28,7 @@ License
 
 #include "faceAreaWeightAMI.H"
 #include "profiling.H"
+#include "OBJstream.H"
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
@@ -36,6 +37,7 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calcAddressing
 (
     List<DynamicList<label>>& srcAddr,
     List<DynamicList<scalar>>& srcWght,
+    List<DynamicList<point>>& srcCtr,
     List<DynamicList<label>>& tgtAddr,
     List<DynamicList<scalar>>& tgtWght,
     label srcFacei,
@@ -60,7 +62,7 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calcAddressing
     seedFaces[srcFacei] = tgtFacei;
 
     // list to keep track of whether src face can be mapped
-    boolList mapFlag(nFacesRemaining, true);
+    bitSet mapFlag(nFacesRemaining, true);
 
     // reset starting seed
     label startSeedi = 0;
@@ -68,6 +70,9 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calcAddressing
     DynamicList<label> nonOverlapFaces;
     do
     {
+        nbrFaces.clear();
+        visitedFaces.clear();
+
         // Do advancing front starting from srcFacei,tgtFacei
         bool faceProcessed = processSourceFace
         (
@@ -79,11 +84,12 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calcAddressing
 
             srcAddr,
             srcWght,
+            srcCtr,
             tgtAddr,
             tgtWght
         );
 
-        mapFlag[srcFacei] = false;
+        mapFlag.unset(srcFacei);
 
         nFacesRemaining--;
 
@@ -122,9 +128,10 @@ bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
     // list of faces currently visited for srcFacei to avoid multiple hits
     DynamicList<label>& visitedFaces,
 
-    // temporary storage for addressing and weights
+    // temporary storage for addressing, weights and centroid
     List<DynamicList<label>>& srcAddr,
     List<DynamicList<scalar>>& srcWght,
+    List<DynamicList<point>>& srcCtr,
     List<DynamicList<label>>& tgtAddr,
     List<DynamicList<scalar>>& tgtWght
 )
@@ -136,9 +143,6 @@ bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
         return false;
     }
 
-    nbrFaces.clear();
-    visitedFaces.clear();
-
     // append initial target face and neighbours
     nbrFaces.append(tgtStartFacei);
     this->appendNbrFaces
@@ -158,16 +162,24 @@ bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
         // process new target face
         label tgtFacei = nbrFaces.remove();
         visitedFaces.append(tgtFacei);
-        scalar area = interArea(srcFacei, tgtFacei);
+
+        scalar interArea = 0;
+        vector interCentroid(Zero);
+        calcInterArea(srcFacei, tgtFacei, interArea, interCentroid);
 
         // store when intersection fractional area > tolerance
-        if (area/this->srcMagSf_[srcFacei] > faceAreaIntersect::tolerance())
+        if
+        (
+            interArea/this->srcMagSf_[srcFacei]
+          > faceAreaIntersect::tolerance()
+        )
         {
             srcAddr[srcFacei].append(tgtFacei);
-            srcWght[srcFacei].append(area);
+            srcWght[srcFacei].append(interArea);
+            srcCtr[srcFacei].append(interCentroid);
 
             tgtAddr[tgtFacei].append(srcFacei);
-            tgtWght[tgtFacei].append(area);
+            tgtWght[tgtFacei].append(interArea);
 
             this->appendNbrFaces
             (
@@ -199,10 +211,10 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
     label& startSeedi,
     label& srcFacei,
     label& tgtFacei,
-    const boolList& mapFlag,
+    const bitSet& mapFlag,
     labelList& seedFaces,
     const DynamicList<label>& visitedFaces,
-    bool errorOnNotFound
+    const bool errorOnNotFound
 ) const
 {
     addProfiling(ami, "faceAreaWeightAMI::setNextFaces");
@@ -216,15 +228,14 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
     bool valuesSet = false;
     for (label faceS: srcNbrFaces)
     {
-        if (mapFlag[faceS] && seedFaces[faceS] == -1)
+        if (mapFlag.test(faceS) && seedFaces[faceS] == -1)
         {
             for (label faceT : visitedFaces)
             {
                 const scalar threshold =
                     this->srcMagSf_[faceS]*faceAreaIntersect::tolerance();
 
-                // store when intersection fractional area > tolerance
-                // Check that faces have enough overlap for robust walking
+                // Store when intersection fractional area > threshold
                 if (overlaps(faceS, faceT, threshold))
                 {
                     seedFaces[faceS] = faceT;
@@ -248,25 +259,31 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
     else
     {
         // try to use existing seed
-        bool foundNextSeed = false;
-        for (label facei = startSeedi; facei < mapFlag.size(); ++facei)
+        label facei = startSeedi;
+        if (!mapFlag.test(startSeedi))
         {
-            if (mapFlag[facei])
+            facei = mapFlag.find_next(facei);
+        }
+        const label startSeedi0 = facei;
+
+        bool foundNextSeed = false;
+        while (facei != -1)
+        {
+            if (!foundNextSeed)
             {
-                if (!foundNextSeed)
-                {
-                    startSeedi = facei;
-                    foundNextSeed = true;
-                }
-
-                if (seedFaces[facei] != -1)
-                {
-                    srcFacei = facei;
-                    tgtFacei = seedFaces[facei];
-
-                    return;
-                }
+                startSeedi = facei;
+                foundNextSeed = true;
             }
+
+            if (seedFaces[facei] != -1)
+            {
+                srcFacei = facei;
+                tgtFacei = seedFaces[facei];
+
+                return;
+            }
+
+            facei = mapFlag.find_next(facei);
         }
 
         // perform new search to find match
@@ -276,25 +293,18 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
                 << "target face" << endl;
         }
 
-        foundNextSeed = false;
-        for (label facei = startSeedi; facei < mapFlag.size(); ++facei)
+        facei = startSeedi0;
+        while (facei != -1)
         {
-            if (mapFlag[facei])
+            srcFacei = facei;
+            tgtFacei = this->findTargetFace(srcFacei, visitedFaces);
+
+            if (tgtFacei >= 0)
             {
-                if (!foundNextSeed)
-                {
-                    startSeedi = facei + 1;
-                    foundNextSeed = true;
-                }
-
-                srcFacei = facei;
-                tgtFacei = this->findTargetFace(srcFacei);
-
-                if (tgtFacei >= 0)
-                {
-                    return;
-                }
+                return;
             }
+
+            facei = mapFlag.find_next(facei);
         }
 
         if (errorOnNotFound)
@@ -307,16 +317,16 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
 
 
 template<class SourcePatch, class TargetPatch>
-Foam::scalar Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
+void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calcInterArea
 (
     const label srcFacei,
-    const label tgtFacei
+    const label tgtFacei,
+    scalar& area,
+    vector& centroid
 ) const
 {
     addProfiling(ami, "faceAreaWeightAMI::interArea");
 
-    scalar area = 0;
-
     const pointField& srcPoints = this->srcPatch_.points();
     const pointField& tgtPoints = this->tgtPatch_.points();
 
@@ -331,7 +341,7 @@ Foam::scalar Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
      || (this->tgtMagSf_[tgtFacei] < ROOTVSMALL)
     )
     {
-        return area;
+        return;
     }
 
     // create intersection object
@@ -359,7 +369,7 @@ Foam::scalar Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
 
     if (magN > ROOTVSMALL)
     {
-        area = inter.calc(src, tgt, n/magN);
+        inter.calc(src, tgt, n/magN, area, centroid);
     }
     else
     {
@@ -371,13 +381,24 @@ Foam::scalar Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
             << endl;
     }
 
+    if
+    (
+        AMIInterpolation<SourcePatch, TargetPatch>::cacheIntersections_
+     && debug
+    )
+    {
+        static OBJstream tris("intersectionTris.obj");
+        const auto& triPts = inter.triangles();
+        for (const auto& tp : triPts)
+        {
+            tris.write(triPointRef(tp[0], tp[1], tp[2]), false);
+        }
+    }
 
     if ((debug > 1) && (area > 0))
     {
         this->writeIntersectionOBJ(area, src, tgt, srcPoints, tgtPoints);
     }
-
-    return area;
 }
 
 
@@ -452,95 +473,70 @@ restartUncoveredSourceFace
 (
     List<DynamicList<label>>& srcAddr,
     List<DynamicList<scalar>>& srcWght,
+    List<DynamicList<point>>& srcCtr,
     List<DynamicList<label>>& tgtAddr,
     List<DynamicList<scalar>>& tgtWght
 )
 {
     addProfiling(ami, "faceAreaWeightAMI::restartUncoveredSourceFace");
 
-    // Collect all src faces with a low weight
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Note: exclude faces in srcNonOverlap_ for ACMI?
+
+    label nBelowMinWeight = 0;
+    const scalar minWeight = 0.95;
+
+    // list of tgt face neighbour faces
+    DynamicList<label> nbrFaces(10);
+
+    // list of faces currently visited for srcFacei to avoid multiple hits
+    DynamicList<label> visitedFaces(10);
 
-    labelHashSet lowWeightFaces(100);
     forAll(srcWght, srcFacei)
     {
-        scalar s = sum(srcWght[srcFacei]);
-        scalar t = s/this->srcMagSf_[srcFacei];
+        const scalar s = sum(srcWght[srcFacei]);
+        const scalar t = s/this->srcMagSf_[srcFacei];
 
-        if (t < 0.5)
+        if (t < minWeight)
         {
-            lowWeightFaces.insert(srcFacei);
-        }
-    }
+            ++nBelowMinWeight;
 
-    if (debug)
-    {
-        Pout<< "faceAreaWeightAMI: restarting search on "
-            << lowWeightFaces.size() << " faces since sum of weights < 0.5"
-            << endl;
-    }
+            const face& f = this->srcPatch_[srcFacei];
 
-    if (lowWeightFaces.size() > 0)
-    {
-        // Erase all the lowWeight source faces from the target
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        DynamicList<label> okSrcFaces(10);
-        DynamicList<scalar> okSrcWeights(10);
-        forAll(tgtAddr, tgtFacei)
-        {
-            okSrcFaces.clear();
-            okSrcWeights.clear();
-            DynamicList<label>& srcFaces = tgtAddr[tgtFacei];
-            DynamicList<scalar>& srcWeights = tgtWght[tgtFacei];
-            forAll(srcFaces, i)
+            forAll(f, fpi)
             {
-                if (!lowWeightFaces.found(srcFaces[i]))
+                const label tgtFacei =
+                    this->findTargetFace(srcFacei, srcAddr[srcFacei], fpi);
+
+                if (tgtFacei != -1)
                 {
-                    okSrcFaces.append(srcFaces[i]);
-                    okSrcWeights.append(srcWeights[i]);
+                    nbrFaces.clear();
+                    visitedFaces = srcAddr[srcFacei];
+
+                    (void)processSourceFace
+                    (
+                        srcFacei,
+                        tgtFacei,
+
+                        nbrFaces,
+                        visitedFaces,
+
+                        srcAddr,
+                        srcWght,
+                        srcCtr,
+                        tgtAddr,
+                        tgtWght
+                    );
                 }
             }
-            if (okSrcFaces.size() < srcFaces.size())
-            {
-                srcFaces.transfer(okSrcFaces);
-                srcWeights.transfer(okSrcWeights);
-            }
         }
+    }
 
-
-
-        // Restart search from best hit
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        // list of tgt face neighbour faces
-        DynamicList<label> nbrFaces(10);
-
-        // list of faces currently visited for srcFacei to avoid multiple hits
-        DynamicList<label> visitedFaces(10);
-
-        for (const label srcFacei : lowWeightFaces)
-        {
-            const label tgtFacei = this->findTargetFace(srcFacei);
-            if (tgtFacei != -1)
-            {
-                //bool faceProcessed =
-                processSourceFace
-                (
-                    srcFacei,
-                    tgtFacei,
-
-                    nbrFaces,
-                    visitedFaces,
-
-                    srcAddr,
-                    srcWght,
-                    tgtAddr,
-                    tgtWght
-                );
-                // ? Check faceProcessed to see if restarting has worked.
-            }
-        }
+    if (debug && nBelowMinWeight)
+    {
+        WarningInFunction
+            << "Restarted search on " << nBelowMinWeight
+            << " faces since sum of weights < " << minWeight
+            << endl;
     }
 }
 
@@ -572,16 +568,45 @@ Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::faceAreaWeightAMI
 {
     this->triangulatePatch(srcPatch, srcTris_, this->srcMagSf_);
     this->triangulatePatch(tgtPatch, tgtTris_, this->tgtMagSf_);
+
+    if (debug)
+    {
+        static label nAMI = 0;
+
+        // Write out triangulated surfaces as OBJ files
+        OBJstream srcTriObj("srcTris_" + Foam::name(nAMI) + ".obj");
+        const pointField& srcPts = srcPatch.points();
+        forAll(srcTris_, facei)
+        {
+            const DynamicList<face>& faces = srcTris_[facei];
+            for (const face& f : faces)
+            {
+                srcTriObj.write
+                (
+                    triPointRef(srcPts[f[0]], srcPts[f[1]], srcPts[f[2]])
+                );
+            }
+        }
+
+        OBJstream tgtTriObj("tgtTris_" + Foam::name(nAMI) + ".obj");
+        const pointField& tgtPts = tgtPatch.points();
+        forAll(tgtTris_, facei)
+        {
+            const DynamicList<face>& faces = tgtTris_[facei];
+            for (const face& f : faces)
+            {
+                tgtTriObj.write
+                (
+                    triPointRef(tgtPts[f[0]], tgtPts[f[1]], tgtPts[f[2]])
+                );
+            }
+        }
+
+        ++nAMI;
+    }
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-template<class SourcePatch, class TargetPatch>
-Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::~faceAreaWeightAMI()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class SourcePatch, class TargetPatch>
@@ -589,6 +614,7 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
 (
     labelListList& srcAddress,
     scalarListList& srcWeights,
+    pointListList& srcCentroids,
     labelListList& tgtAddress,
     scalarListList& tgtWeights,
     label srcFacei,
@@ -608,6 +634,9 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
             tgtFacei
         );
 
+    srcCentroids.setSize(srcAddress.size());
+
+
     if (!ok)
     {
         return;
@@ -616,6 +645,7 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
     // temporary storage for addressing and weights
     List<DynamicList<label>> srcAddr(this->srcPatch_.size());
     List<DynamicList<scalar>> srcWght(srcAddr.size());
+    List<DynamicList<point>> srcCtr(srcAddr.size());
     List<DynamicList<label>> tgtAddr(this->tgtPatch_.size());
     List<DynamicList<scalar>> tgtWght(tgtAddr.size());
 
@@ -623,6 +653,7 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
     (
         srcAddr,
         srcWght,
+        srcCtr,
         tgtAddr,
         tgtWght,
         srcFacei,
@@ -644,17 +675,18 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
         (
             srcAddr,
             srcWght,
+            srcCtr,
             tgtAddr,
             tgtWght
         );
     }
 
-
-    // transfer data to persistent storage
+    // Transfer data to persistent storage
     forAll(srcAddr, i)
     {
         srcAddress[i].transfer(srcAddr[i]);
         srcWeights[i].transfer(srcWght[i]);
+        srcCentroids[i].transfer(srcCtr[i]);
     }
     forAll(tgtAddr, i)
     {
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.H
index 0c6dec88d6..10622d727a 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.H
@@ -80,11 +80,13 @@ protected:
 
         // Marching front
 
-            //- Calculate addressing and weights using temporary storage
+            //- Calculate addressing, weights and centroids using temporary
+            //- storage
             virtual void calcAddressing
             (
                 List<DynamicList<label>>& srcAddress,
                 List<DynamicList<scalar>>& srcWeights,
+                List<DynamicList<point>>& srcCentroids,
                 List<DynamicList<label>>& tgtAddress,
                 List<DynamicList<scalar>>& tgtWeights,
                 label srcFacei,
@@ -100,6 +102,7 @@ protected:
                 DynamicList<label>& visitedFaces,
                 List<DynamicList<label>>& srcAddr,
                 List<DynamicList<scalar>>& srcWght,
+                List<DynamicList<point>>& srcCtr,
                 List<DynamicList<label>>& tgtAddr,
                 List<DynamicList<scalar>>& tgtWght
             );
@@ -109,6 +112,7 @@ protected:
             (
                 List<DynamicList<label>>& srcAddr,
                 List<DynamicList<scalar>>& srcWght,
+                List<DynamicList<point>>& srcCtr,
                 List<DynamicList<label>>& tgtAddr,
                 List<DynamicList<scalar>>& tgtWght
             );
@@ -119,20 +123,22 @@ protected:
                 label& startSeedi,
                 label& srcFacei,
                 label& tgtFacei,
-                const boolList& mapFlag,
+                const bitSet& mapFlag,
                 labelList& seedFaces,
                 const DynamicList<label>& visitedFaces,
-                bool errorOnNotFound = true
+                const bool errorOnNotFound = true
             ) const;
 
 
         // Evaluation
 
             //- Area of intersection between source and target faces
-            virtual scalar interArea
+            virtual void calcInterArea
             (
                 const label srcFacei,
-                const label tgtFacei
+                const label tgtFacei,
+                scalar& area,
+                vector& centroid
             ) const;
 
             //- Return true if faces overlap
@@ -165,18 +171,19 @@ public:
 
 
     //- Destructor
-    virtual ~faceAreaWeightAMI();
+    virtual ~faceAreaWeightAMI() = default;
 
 
     // Member Functions
 
         // Manipulation
 
-            //- Update addressing and weights
+            //- Update addressing, weights and centroids
             virtual void calculate
             (
                 labelListList& srcAddress,
                 scalarListList& srcWeights,
+                pointListList& srcCentroids,
                 labelListList& tgtAddress,
                 scalarListList& tgtWeights,
                 label srcFacei = -1,
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C
index 5bfaa347ca..a159fae7a0 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C
@@ -190,13 +190,6 @@ Foam::mapNearestAMI<SourcePatch, TargetPatch>::mapNearestAMI
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-template<class SourcePatch, class TargetPatch>
-Foam::mapNearestAMI<SourcePatch, TargetPatch>::~mapNearestAMI()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class SourcePatch, class TargetPatch>
@@ -204,6 +197,7 @@ void Foam::mapNearestAMI<SourcePatch, TargetPatch>::calculate
 (
     labelListList& srcAddress,
     scalarListList& srcWeights,
+    pointListList& srcCentroids,
     labelListList& tgtAddress,
     scalarListList& tgtWeights,
     label srcFacei,
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.H
index b59e0682d3..f347fcc767 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.H
@@ -123,7 +123,7 @@ public:
 
 
     //- Destructor
-    virtual ~mapNearestAMI();
+    virtual ~mapNearestAMI() = default;
 
 
     // Member Functions
@@ -135,6 +135,7 @@ public:
             (
                 labelListList& srcAddress,
                 scalarListList& srcWeights,
+                pointListList& srcCentroids,
                 labelListList& tgtAddress,
                 scalarListList& tgtWeights,
                 label srcFacei = -1,
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.C
index 78275e0111..5692a656ba 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.C
@@ -35,7 +35,7 @@ void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
     label& startSeedi,
     label& srcFacei,
     label& tgtFacei,
-    const boolList& mapFlag,
+    const bitSet& mapFlag,
     labelList& seedFaces,
     const DynamicList<label>& visitedFaces,
     const bool errorOnNotFound
@@ -78,14 +78,6 @@ partialFaceAreaWeightAMI
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-template<class SourcePatch, class TargetPatch>
-Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::
-~partialFaceAreaWeightAMI()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class SourcePatch, class TargetPatch>
@@ -100,6 +92,7 @@ void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
 (
     labelListList& srcAddress,
     scalarListList& srcWeights,
+    pointListList& srcCentroids,
     labelListList& tgtAddress,
     scalarListList& tgtWeights,
     label srcFacei,
@@ -122,9 +115,12 @@ void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
         return;
     }
 
+    srcCentroids.setSize(srcAddress.size());
+
     // temporary storage for addressing and weights
     List<DynamicList<label>> srcAddr(this->srcPatch_.size());
     List<DynamicList<scalar>> srcWght(srcAddr.size());
+    List<DynamicList<point>> srcCtr(srcAddr.size());
     List<DynamicList<label>> tgtAddr(this->tgtPatch_.size());
     List<DynamicList<scalar>> tgtWght(tgtAddr.size());
 
@@ -132,6 +128,7 @@ void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
     (
         srcAddr,
         srcWght,
+        srcCtr,
         tgtAddr,
         tgtWght,
         srcFacei,
@@ -143,6 +140,7 @@ void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
     {
         srcAddress[i].transfer(srcAddr[i]);
         srcWeights[i].transfer(srcWght[i]);
+        srcCentroids[i].transfer(srcCtr[i]);
     }
     forAll(tgtAddr, i)
     {
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.H
index becb1c69fb..2d0255d6dc 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.H
@@ -65,19 +65,22 @@ private:
         //- No copy assignment
         void operator=(const partialFaceAreaWeightAMI&) = delete;
 
-        // Marching front
 
-            //- Set the source and target seed faces
-            virtual void setNextFaces
-            (
-                label& startSeedi,
-                label& srcFacei,
-                label& tgtFacei,
-                const boolList& mapFlag,
-                labelList& seedFaces,
-                const DynamicList<label>& visitedFaces,
-                bool errorOnNotFound = true
-            ) const;
+protected:
+
+    // Protected Member Functions
+
+        //- Set the source and target seed faces
+        virtual void setNextFaces
+        (
+            label& startSeedi,
+            label& srcFacei,
+            label& tgtFacei,
+            const bitSet& mapFlag,
+            labelList& seedFaces,
+            const DynamicList<label>& visitedFaces,
+            const bool errorOnNotFound = true
+        ) const;
 
 
 public:
@@ -100,7 +103,7 @@ public:
 
 
     //- Destructor
-    virtual ~partialFaceAreaWeightAMI();
+    virtual ~partialFaceAreaWeightAMI() = default;
 
 
     // Member Functions
@@ -118,6 +121,7 @@ public:
             (
                 labelListList& srcAddress,
                 scalarListList& srcWeights,
+                pointListList& srcCentroids,
                 labelListList& tgtAddress,
                 scalarListList& tgtWeights,
                 label srcFacei = -1,
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/findNearestMaskedOp.H b/src/meshTools/AMIInterpolation/AMIInterpolation/findNearestMaskedOp.H
new file mode 100644
index 0000000000..4bc2507cc9
--- /dev/null
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/findNearestMaskedOp.H
@@ -0,0 +1,62 @@
+#include "indexedOctree.H"
+#include "treeDataPrimitivePatch.H"
+
+namespace Foam
+{
+
+template<class PatchType>
+class findNearestMaskedOp
+{
+    const indexedOctree<treeDataPrimitivePatch<PatchType>>& tree_;
+    const labelUList& excludeIndices_;
+
+public:
+
+    findNearestMaskedOp
+    (
+        const indexedOctree<treeDataPrimitivePatch<PatchType>>& tree,
+        const labelUList& excludeIndices
+    )
+    :
+        tree_(tree),
+        excludeIndices_(excludeIndices)
+    {}
+
+    void operator()
+    (
+        const labelUList& indices,
+        const point& sample,
+
+        scalar& nearestDistSqr,
+        label& minIndex,
+        point& nearestPoint
+    ) const
+    {
+        const treeDataPrimitivePatch<PatchType>& shape = tree_.shapes();
+        const PatchType& patch = shape.patch();
+
+        const pointField& points = patch.points();
+
+        forAll(indices, i)
+        {
+            const label index = indices[i];
+
+            if (!excludeIndices_.found(index))
+            {
+                const typename PatchType::FaceType& f = patch[index];
+
+                pointHit nearHit = f.nearestPoint(sample, points);
+                scalar distSqr = sqr(nearHit.distance());
+
+                if (distSqr < nearestDistSqr)
+                {
+                    nearestDistSqr = distSqr;
+                    minIndex = index;
+                    nearestPoint = nearHit.rawPoint();
+                }
+            }
+        }
+    }
+};
+
+} // End namespace Foam
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
index 97428ec0e3..8976e57b58 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
@@ -219,13 +219,15 @@ void Foam::faceAreaIntersect::triSliceWithPlane
 }
 
 
-Foam::scalar Foam::faceAreaIntersect::triangleIntersect
+void Foam::faceAreaIntersect::triangleIntersect
 (
     const triPoints& src,
     const point& tgt0,
     const point& tgt1,
     const point& tgt2,
-    const vector& n
+    const vector& n,
+    scalar& area,
+    vector& centroid
 ) const
 {
     // Work storage
@@ -241,7 +243,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
     const scalar srcArea(triArea(src));
     if (srcArea < ROOTVSMALL)
     {
-        return 0.0;
+        return;
     }
 
     // Typical length scale
@@ -255,7 +257,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
         scalar s = mag(tgt1 - tgt0);
         if (s < ROOTVSMALL)
         {
-            return 0.0;
+            return;
         }
 
         // Note: outer product with n pre-scaled with edge length. This is
@@ -268,7 +270,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
             // Triangle either zero edge length (s == 0) or
             // perpendicular to face normal n. In either case zero
             // overlap area
-            return 0.0;
+            return;
         }
         else
         {
@@ -279,7 +281,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
 
     if (nWorkTris1 == 0)
     {
-        return 0.0;
+        return;
     }
 
     // Edge 1
@@ -290,7 +292,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
         scalar s = mag(tgt2 - tgt1);
         if (s < ROOTVSMALL)
         {
-            return 0.0;
+            return;
         }
         const vector n1((tgt1 - tgt2)^(-s*n));
         const scalar magSqrN1(magSqr(n1));
@@ -300,7 +302,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
             // Triangle either zero edge length (s == 0) or
             // perpendicular to face normal n. In either case zero
             // overlap area
-            return 0.0;
+            return;
         }
         else
         {
@@ -315,7 +317,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
 
             if (nWorkTris2 == 0)
             {
-                return 0.0;
+                return;
             }
         }
     }
@@ -328,7 +330,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
         scalar s = mag(tgt2 - tgt0);
         if (s < ROOTVSMALL)
         {
-            return 0.0;
+            return;
         }
         const vector n2((tgt2 - tgt0)^(-s*n));
         const scalar magSqrN2(magSqr(n2));
@@ -338,7 +340,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
             // Triangle either zero edge length (s == 0) or
             // perpendicular to face normal n. In either case zero
             // overlap area
-            return 0.0;
+            return;
         }
         else
         {
@@ -353,23 +355,25 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
 
             if (nWorkTris1 == 0)
             {
-                return 0.0;
+                return;
             }
             else
             {
                 // Calculate area of sub-triangles
-                scalar area = 0.0;
                 for (label i = 0; i < nWorkTris1; ++i)
                 {
-                    area += triArea(workTris1[i]);
+                    // Area of intersection
+                    const scalar currArea = triArea(workTris1[i]);
+                    area += currArea;
+
+                    // Area-weighted centroid of intersection
+                    centroid += currArea*triCentroid(workTris1[i]);
 
                     if (cacheTriangulation_)
                     {
                         triangles_.append(workTris1[i]);
                     }
                 }
-
-                return area;
             }
         }
     }
@@ -443,12 +447,13 @@ void Foam::faceAreaIntersect::triangulate
 }
 
 
-
-Foam::scalar Foam::faceAreaIntersect::calc
+void Foam::faceAreaIntersect::calc
 (
     const face& faceA,
     const face& faceB,
-    const vector& n
+    const vector& n,
+    scalar& area,
+    vector& centroid
 ) const
 {
     if (cacheTriangulation_)
@@ -456,8 +461,10 @@ Foam::scalar Foam::faceAreaIntersect::calc
         triangles_.clear();
     }
 
+    area = 0.0;
+    centroid = vector::zero;
+
     // Intersect triangles
-    scalar totalArea = 0.0;
     for (const face& triA : trisA_)
     {
         triPoints tpA = getTriPoints(pointsA_, triA, false);
@@ -466,32 +473,38 @@ Foam::scalar Foam::faceAreaIntersect::calc
         {
             if (reverseB_)
             {
-                totalArea +=
-                    triangleIntersect
-                    (
-                        tpA,
-                        pointsB_[triB[0]],
-                        pointsB_[triB[1]],
-                        pointsB_[triB[2]],
-                        n
-                    );
+                triangleIntersect
+                (
+                    tpA,
+                    pointsB_[triB[0]],
+                    pointsB_[triB[1]],
+                    pointsB_[triB[2]],
+                    n,
+                    area,
+                    centroid
+                );
             }
             else
             {
-                totalArea +=
-                    triangleIntersect
-                    (
-                        tpA,
-                        pointsB_[triB[2]],
-                        pointsB_[triB[1]],
-                        pointsB_[triB[0]],
-                        n
-                    );
+                triangleIntersect
+                (
+                    tpA,
+                    pointsB_[triB[2]],
+                    pointsB_[triB[1]],
+                    pointsB_[triB[0]],
+                    n,
+                    area,
+                    centroid
+                );
             }
         }
     }
 
-    return totalArea;
+    // Area weighed centroid
+    if (area > 0)
+    {
+        centroid /= area;
+    }
 }
 
 
@@ -503,8 +516,10 @@ bool Foam::faceAreaIntersect::overlaps
     const scalar threshold
 ) const
 {
+    scalar area = 0.0;
+    vector centroid(Zero);
+
     // Intersect triangles
-    scalar totalArea = 0.0;
     for (const face& triA : trisA_)
     {
         const triPoints tpA = getTriPoints(pointsA_, triA, false);
@@ -513,30 +528,32 @@ bool Foam::faceAreaIntersect::overlaps
         {
             if (reverseB_)
             {
-                totalArea +=
-                    triangleIntersect
-                    (
-                        tpA,
-                        pointsB_[triB[0]],
-                        pointsB_[triB[1]],
-                        pointsB_[triB[2]],
-                        n
-                    );
+                triangleIntersect
+                (
+                    tpA,
+                    pointsB_[triB[0]],
+                    pointsB_[triB[1]],
+                    pointsB_[triB[2]],
+                    n,
+                    area,
+                    centroid
+                );
             }
             else
             {
-                totalArea +=
-                    triangleIntersect
-                    (
-                        tpA,
-                        pointsB_[triB[2]],
-                        pointsB_[triB[1]],
-                        pointsB_[triB[0]],
-                        n
-                    );
+                triangleIntersect
+                (
+                    tpA,
+                    pointsB_[triB[2]],
+                    pointsB_[triB[1]],
+                    pointsB_[triB[0]],
+                    n,
+                    area,
+                    centroid
+                );
             }
 
-            if (totalArea > threshold)
+            if (area > threshold)
             {
                 return true;
             }
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
index f0ba9b3e92..916d40e36f 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
@@ -97,6 +97,7 @@ private:
 
     // Static data members
 
+        //- Tolerance
         static scalar tol;
 
 
@@ -132,6 +133,9 @@ private:
         //- Return triangle area
         inline scalar triArea(const triPoints& t) const;
 
+        //- Return triangle centre
+        inline vector triCentroid(const triPoints& t) const;
+
         //- Slice triangle with plane and generate new cut sub-triangles
         void triSliceWithPlane
         (
@@ -143,13 +147,15 @@ private:
         ) const;
 
         //- Return area of intersection of triangles src and tgt
-        scalar triangleIntersect
+        void triangleIntersect
         (
             const triPoints& src,
             const point& tgt0,
             const point& tgt1,
             const point& tgt2,
-            const vector& n
+            const vector& n,
+            scalar& area,
+            vector& centroid
         ) const;
 
 
@@ -199,12 +205,15 @@ public:
             DynamicList<face>& faces
         );
 
-        //- Return area of intersection of faceA with faceB
-        scalar calc
+        //- Return area of intersection of faceA with faceB and effective
+        //- face centre
+        void calc
         (
             const face& faceA,
             const face& faceB,
-            const vector& n
+            const vector& n,
+            scalar& area,
+            vector& centroid
         ) const;
 
         //- Return area of intersection of faceA with faceB
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
index 12113a594d..836b87d762 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
@@ -113,6 +113,15 @@ inline Foam::scalar Foam::faceAreaIntersect::triArea(const triPoints& t) const
 }
 
 
+inline Foam::vector Foam::faceAreaIntersect::triCentroid
+(
+    const triPoints& t
+) const
+{
+    return (t[0] + t[1] + t[2])/3;
+}
+
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::scalar& Foam::faceAreaIntersect::tolerance()
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C
index 3cd26ef2b8..cbba50fb34 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C
+++ b/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2018 OpenCFD Ltd.
+    Copyright (C) 2017-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,195 +41,225 @@ namespace Foam
     addToRunTimeSelectionTable(polyPatch, cyclicACMIPolyPatch, dictionary);
 }
 
-const Foam::scalar Foam::cyclicACMIPolyPatch::tolerance_ = 1e-10;
-
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+void Foam::cyclicACMIPolyPatch::reportCoverage
+(
+    const word& name,
+    const scalarField& weightSum
+) const
+{
+    label nUncovered = 0;
+    label nCovered = 0;
+    for (const scalar sum : weightSum)
+    {
+        if (sum < tolerance_)
+        {
+            ++nUncovered;
+        }
+        else if (sum > scalar(1) - tolerance_)
+        {
+            ++nCovered;
+        }
+    }
+    reduce(nUncovered, sumOp<label>());
+    reduce(nCovered, sumOp<label>());
+    label nTotal = returnReduce(weightSum.size(), sumOp<label>());
+
+    Info<< "ACMI: Patch " << name << " uncovered/blended/covered = "
+        << nUncovered << ", " << nTotal-nUncovered-nCovered
+        << ", " << nCovered << endl;
+}
+
+
 void Foam::cyclicACMIPolyPatch::resetAMI
 (
-    const AMIPatchToPatchInterpolation::interpolationMethod&
+    const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod
 ) const
 {
-    if (owner())
-    {
-        const polyPatch& nonOverlapPatch = this->nonOverlapPatch();
+    resetAMI(boundaryMesh().mesh().points(), AMIMethod);
+}
 
+
+void Foam::cyclicACMIPolyPatch::resetAMI
+(
+    const UList<point>& points,
+    const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod
+) const
+{
+    if (!owner())
+    {
+        return;
+    }
+
+    const polyPatch& nonOverlapPatch = this->nonOverlapPatch();
+
+    if (debug)
+    {
+        Pout<< "cyclicACMIPolyPatch::resetAMI : recalculating weights"
+            << " for " << name() << " and " << nonOverlapPatch.name()
+            << endl;
+    }
+
+
+WarningInFunction<< "DEACTIVATED clearGeom()" << endl;
+//    if (boundaryMesh().mesh().hasCellCentres())
+    if (0 && boundaryMesh().mesh().hasCellCentres())
+    {
         if (debug)
         {
-            Pout<< "cyclicACMIPolyPatch::resetAMI : recalculating weights"
+            Pout<< "cyclicACMIPolyPatch::resetAMI : clearing cellCentres"
                 << " for " << name() << " and " << nonOverlapPatch.name()
                 << endl;
         }
 
-        if (boundaryMesh().mesh().hasCellCentres())
-        {
-            if (debug)
-            {
-                Pout<< "cyclicACMIPolyPatch::resetAMI : clearing cellCentres"
-                    << " for " << name() << " and " << nonOverlapPatch.name()
-                    << endl;
-            }
-
-            //WarningInFunction
-            //    << "The mesh already has cellCentres calculated when"
-            //    << " resetting ACMI " << name() << "." << endl
-            //    << "This is a problem since ACMI adapts the face areas"
-            //    << " (to close cells) so this has" << endl
-            //    << "to be done before cell centre calculation." << endl
-            //    << "This can happen if e.g. the cyclicACMI is after"
-            //    << " any processor patches in the boundary." << endl;
-            const_cast<polyMesh&>
-            (
-                boundaryMesh().mesh()
-            ).primitiveMesh::clearGeom();
-        }
-
-
-        // Trigger re-building of faceAreas
-        (void)boundaryMesh().mesh().faceAreas();
-
-
-        // Calculate the AMI using partial face-area-weighted. This leaves
-        // the weights as fractions of local areas (sum(weights) = 1 means
-        // face is fully covered)
-        cyclicAMIPolyPatch::resetAMI
+        //WarningInFunction
+        //    << "The mesh already has cellCentres calculated when"
+        //    << " resetting ACMI " << name() << "." << endl
+        //    << "This is a problem since ACMI adapts the face areas"
+        //    << " (to close cells) so this has" << endl
+        //    << "to be done before cell centre calculation." << endl
+        //    << "This can happen if e.g. the cyclicACMI is after"
+        //    << " any processor patches in the boundary." << endl;
+        const_cast<polyMesh&>
         (
-            AMIPatchToPatchInterpolation::imPartialFaceAreaWeight
-        );
+            boundaryMesh().mesh()
+        ).primitiveMesh::clearGeom();
+    }
 
-        AMIPatchToPatchInterpolation& AMI =
-            const_cast<AMIPatchToPatchInterpolation&>(this->AMI());
 
-        // Output some stats. AMIInterpolation will have already output the
-        // average weights ("sum(weights) min:1 max:1 average:1")
+    // Trigger re-building of faceAreas
+    (void)boundaryMesh().mesh().faceAreas();
+
+
+    // Calculate the AMI using partial face-area-weighted. This leaves
+    // the weights as fractions of local areas (sum(weights) = 1 means
+    // face is fully covered)
+    cyclicAMIPolyPatch::resetAMI
+    (
+        points,
+        AMIPatchToPatchInterpolation::imPartialFaceAreaWeight
+    );
+
+    const AMIPatchToPatchInterpolation& AMI = this->AMI();
+
+    // Output some statistics
+    reportCoverage("source", AMI.srcWeightsSum());
+    reportCoverage("target", AMI.tgtWeightsSum());
+
+    // Set the mask fields
+    // Note:
+    // - assumes that the non-overlap patches are decomposed using the same
+    //   decomposition as the coupled patches (per side)
+    srcMask_ = min(scalar(1), max(scalar(0), AMI.srcWeightsSum()));
+    tgtMask_ = min(scalar(1), max(scalar(0), AMI.tgtWeightsSum()));
+
+    if (debug)
+    {
+        Pout<< "resetAMI" << endl;
         {
-            const scalarField& wghtsSum = AMI.srcWeightsSum();
-
-            label nUncovered = 0;
-            label nCovered = 0;
-            forAll(wghtsSum, facei)
-            {
-                scalar sum = wghtsSum[facei];
-                if (sum < tolerance_)
-                {
-                    nUncovered++;
-                }
-                else if (sum > scalar(1)-tolerance_)
-                {
-                    nCovered++;
-                }
-            }
-            reduce(nUncovered, sumOp<label>());
-            reduce(nCovered, sumOp<label>());
-            label nTotal = returnReduce(wghtsSum.size(), sumOp<label>());
-
-            Info<< "ACMI: Patch source uncovered/blended/covered = "
-                << nUncovered << ", " << nTotal-nUncovered-nCovered
-                << ", " << nCovered << endl;
+            const cyclicACMIPolyPatch& patch = *this;
+            Pout<< "patch:" << patch.name() << " size:" << patch.size()
+                << " non-overlap patch: " << patch.nonOverlapPatch().name()
+                << " size:" << patch.nonOverlapPatch().size()
+                << " mask size:" << patch.srcMask().size() << endl;
         }
         {
-            const scalarField& wghtsSum = AMI.tgtWeightsSum();
-
-            label nUncovered = 0;
-            label nCovered = 0;
-            forAll(wghtsSum, facei)
-            {
-                scalar sum = wghtsSum[facei];
-                if (sum < tolerance_)
-                {
-                    nUncovered++;
-                }
-                else if (sum > scalar(1)-tolerance_)
-                {
-                    nCovered++;
-                }
-            }
-            reduce(nUncovered, sumOp<label>());
-            reduce(nCovered, sumOp<label>());
-            label nTotal = returnReduce(wghtsSum.size(), sumOp<label>());
-
-            Info<< "ACMI: Patch target uncovered/blended/covered = "
-                << nUncovered << ", " << nTotal-nUncovered-nCovered
-                << ", " << nCovered << endl;
+            const cyclicACMIPolyPatch& patch = this->neighbPatch();
+            Pout<< "patch:" << patch.name() << " size:" << patch.size()
+                << " non-overlap patch: " << patch.nonOverlapPatch().name()
+                << " size:" << patch.nonOverlapPatch().size()
+                << " mask size:" << patch.neighbPatch().tgtMask().size()
+                << endl;
         }
-
-        srcMask_ =
-            min(scalar(1) - tolerance_, max(tolerance_, AMI.srcWeightsSum()));
-
-        tgtMask_ =
-            min(scalar(1) - tolerance_, max(tolerance_, AMI.tgtWeightsSum()));
+    }
+}
 
 
-        // Adapt owner side areas. Note that in uncoupled situations (e.g.
-        // decomposePar) srcMask, tgtMask can be zero size.
-        if (srcMask_.size())
+void Foam::cyclicACMIPolyPatch::scalePatchFaceAreas()
+{
+    if (!owner() || !canResetAMI())
+    {
+        return;
+    }
+
+    scalePatchFaceAreas(*this);
+    scalePatchFaceAreas(this->neighbPatch());
+}
+
+
+void Foam::cyclicACMIPolyPatch::scalePatchFaceAreas
+(
+    const cyclicACMIPolyPatch& acmipp
+)
+{
+    // Primitive patch face areas have been cleared/reset based on the raw
+    // points - need to reset to avoid double-accounting of face areas
+
+    DebugPout
+        << "rescaling non-overlap patch areas" << endl;
+
+    const scalar maxTol = scalar(1) - tolerance_;
+    const scalarField& mask = acmipp.mask();
+
+    const polyPatch& nonOverlapPatch = acmipp.nonOverlapPatch();
+    vectorField::subField noSf = nonOverlapPatch.faceAreas();
+
+    if (mask.size() != noSf.size())
+    {
+        WarningInFunction
+            << "Inconsistent sizes for patch: " << acmipp.name()
+            << " - not manipulating patches" << nl
+            << " - size: " << size() << nl
+            << " - non-overlap patch size: " << noSf.size() << nl
+            << " - mask size: " << mask.size() << nl
+            << "This is OK for decomposition but should be considered fatal "
+            << "at run-time" << endl;
+
+        return;
+    }
+
+    forAll(noSf, facei)
+    {
+        const scalar w = min(maxTol, max(tolerance_, mask[facei]));
+        noSf[facei] *= scalar(1) - w;
+    }
+
+    if (!createAMIFaces_)
+    {
+        // Note: for topological update (createAMIFaces_ = true)
+        // AMI coupled patch face areas are updated as part of the topological
+        // updates, e.g. by the calls to cyclicAMIPolyPatch's setTopology and
+        // initMovePoints
+        DebugPout
+            << "scaling coupled patch areas" << endl;
+
+        // Scale the coupled patch face areas
+        vectorField::subField Sf = acmipp.faceAreas();
+
+        forAll(Sf, facei)
         {
-            vectorField::subField Sf = faceAreas();
-            vectorField::subField noSf = nonOverlapPatch.faceAreas();
-
-            forAll(Sf, facei)
-            {
-                Sf[facei] *= srcMask_[facei];
-                noSf[facei] *= 1.0 - srcMask_[facei];
-            }
-        }
-        // Adapt slave side areas
-        if (tgtMask_.size())
-        {
-            const cyclicACMIPolyPatch& cp =
-                refCast<const cyclicACMIPolyPatch>(this->neighbPatch());
-            const polyPatch& pp = cp.nonOverlapPatch();
-
-            vectorField::subField Sf = cp.faceAreas();
-            vectorField::subField noSf = pp.faceAreas();
-
-            forAll(Sf, facei)
-            {
-                Sf[facei] *= tgtMask_[facei];
-                noSf[facei] *= 1.0 - tgtMask_[facei];
-            }
+            Sf[facei] *= max(tolerance_, mask[facei]);
         }
 
         // Re-normalise the weights since the effect of overlap is already
-        // accounted for in the area.
+        // accounted for in the area
+        auto& weights = const_cast<scalarListList&>(acmipp.weights());
+        auto& weightsSum = const_cast<scalarField&>(acmipp.weightsSum());
+        forAll(weights, i)
         {
-            scalarListList& srcWeights = AMI.srcWeights();
-            scalarField& srcWeightsSum = AMI.srcWeightsSum();
-            forAll(srcWeights, i)
+            scalarList& wghts = weights[i];
+            if (wghts.size())
             {
-                scalarList& wghts = srcWeights[i];
-                if (wghts.size())
-                {
-                    scalar& sum = srcWeightsSum[i];
+                scalar& sum = weightsSum[i];
 
-                    forAll(wghts, j)
-                    {
-                        wghts[j] /= sum;
-                    }
-                    sum = 1.0;
+                forAll(wghts, j)
+                {
+                    wghts[j] /= sum;
                 }
+                sum = 1.0;
             }
         }
-        {
-            scalarListList& tgtWeights = AMI.tgtWeights();
-            scalarField& tgtWeightsSum = AMI.tgtWeightsSum();
-            forAll(tgtWeights, i)
-            {
-                scalarList& wghts = tgtWeights[i];
-                if (wghts.size())
-                {
-                    scalar& sum = tgtWeightsSum[i];
-                    forAll(wghts, j)
-                    {
-                        wghts[j] /= sum;
-                    }
-                    sum = 1.0;
-                }
-            }
-        }
-
-        // Set the updated flag
-        updated_ = true;
     }
 }
 
@@ -244,9 +274,9 @@ void Foam::cyclicACMIPolyPatch::initGeometry(PstreamBuffers& pBufs)
     // Note: calculates transformation and triggers face centre calculation
     cyclicAMIPolyPatch::initGeometry(pBufs);
 
-    // Initialise the AMI early to make sure we adapt the face areas before the
-    // cell centre calculation gets triggered.
-    resetAMI();
+    // On start-up there are no topological updates so scale the face areas
+    // - Note: resetAMI called by cyclicAMIPolyPatch::initGeometry
+    scalePatchFaceAreas();
 }
 
 
@@ -272,10 +302,10 @@ void Foam::cyclicACMIPolyPatch::initMovePoints
     }
 
     // Note: calculates transformation and triggers face centre calculation
+    // - Note: resetAMI called by cyclicAMIPolyPatch::initMovePoints
     cyclicAMIPolyPatch::initMovePoints(pBufs, p);
 
-    // Initialise the AMI early. See initGeometry.
-    resetAMI();
+    scalePatchFaceAreas();
 }
 
 
@@ -289,6 +319,8 @@ void Foam::cyclicACMIPolyPatch::movePoints
     {
         Pout<< "cyclicACMIPolyPatch::movePoints : " << name() << endl;
     }
+
+    // When topology is changing, this will scale the duplicate AMI faces
     cyclicAMIPolyPatch::movePoints(pBufs, p);
 }
 
@@ -352,8 +384,7 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     nonOverlapPatchName_(word::null),
     nonOverlapPatchID_(-1),
     srcMask_(),
-    tgtMask_(),
-    updated_(false)
+    tgtMask_()
 {
     AMIRequireMatch_ = false;
 
@@ -375,8 +406,7 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     nonOverlapPatchName_(dict.lookup("nonOverlapPatch")),
     nonOverlapPatchID_(-1),
     srcMask_(),
-    tgtMask_(),
-    updated_(false)
+    tgtMask_()
 {
     AMIRequireMatch_ = false;
 
@@ -403,8 +433,7 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     nonOverlapPatchName_(pp.nonOverlapPatchName_),
     nonOverlapPatchID_(-1),
     srcMask_(),
-    tgtMask_(),
-    updated_(false)
+    tgtMask_()
 {
     AMIRequireMatch_ = false;
 
@@ -428,8 +457,7 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     nonOverlapPatchName_(nonOverlapPatchName),
     nonOverlapPatchID_(-1),
     srcMask_(),
-    tgtMask_(),
-    updated_(false)
+    tgtMask_()
 {
     AMIRequireMatch_ = false;
 
@@ -459,19 +487,12 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     nonOverlapPatchName_(pp.nonOverlapPatchName_),
     nonOverlapPatchID_(-1),
     srcMask_(),
-    tgtMask_(),
-    updated_(false)
+    tgtMask_()
 {
     AMIRequireMatch_ = false;
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::cyclicACMIPolyPatch::~cyclicACMIPolyPatch()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 const Foam::cyclicACMIPolyPatch& Foam::cyclicACMIPolyPatch::neighbPatch() const
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.H b/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.H
index 5f2bd9dfad..c858c14e04 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.H
+++ b/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -60,9 +60,6 @@ private:
 
     // Private data
 
-        //- Fraction of face area below which face is considered disconnected
-        static const scalar tolerance_;
-
         //- Name of non-overlapping patch
         const word nonOverlapPatchName_;
 
@@ -75,21 +72,40 @@ private:
         //- Mask/weighting for target patch
         mutable scalarField tgtMask_;
 
-        //- Flag to indicate that AMI has been updated
-        mutable bool updated_;
-
 
 protected:
 
     // Protected Member Functions
 
-        //- Reset the AMI interpolator
+        //- Report coverage statics, e.g. number of uncovered/blended/covered
+        //- faces
+        void reportCoverage
+        (
+            const word& name,
+            const scalarField& weightSum
+        ) const;
+
+        //- Reset the AMI interpolator, supply patch points
+        virtual void resetAMI
+        (
+            const UList<point>& points,
+            const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod =
+                AMIPatchToPatchInterpolation::imFaceAreaWeight
+        ) const;
+
+        //-  Reset the AMI interpolator, use current patch points
         virtual void resetAMI
         (
             const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod =
                 AMIPatchToPatchInterpolation::imFaceAreaWeight
         ) const;
 
+        //- Scale patch face areas to maintain physical area
+        virtual void scalePatchFaceAreas();
+
+        //- Scale patch face areas to maintain physical area
+        virtual void scalePatchFaceAreas(const cyclicACMIPolyPatch& acmipp);
+
         //- Initialise the calculation of the patch geometry
         virtual void initGeometry(PstreamBuffers&);
 
@@ -111,12 +127,6 @@ protected:
         //- Clear geometry
         virtual void clearGeom();
 
-        //- Return the mask/weighting for the source patch
-        virtual const scalarField& srcMask() const;
-
-        //- Return the mask/weighting for the target patch
-        virtual const scalarField& tgtMask() const;
-
 
 public:
 
@@ -235,19 +245,13 @@ public:
 
 
     //- Destructor
-    virtual ~cyclicACMIPolyPatch();
+    virtual ~cyclicACMIPolyPatch() = default;
 
 
     // Member Functions
 
         // Access
 
-            //- Reset the updated flag
-            inline void setUpdated(bool flag) const;
-
-            //- Return access to the updated flag
-            inline bool updated() const;
-
             //- Return a reference to the neighbour patch
             virtual const cyclicACMIPolyPatch& neighbPatch() const;
 
@@ -263,9 +267,15 @@ public:
             //- Return a reference to the non-overlapping patch
             inline polyPatch& nonOverlapPatch();
 
-            //- Mask field where 1 = overlap, 0 = no-overlap
+            //- Mask field where 1 = overlap(coupled), 0 = no-overlap
             inline const scalarField& mask() const;
 
+            //- Return the mask/weighting for the source patch
+            virtual const scalarField& srcMask() const;
+
+            //- Return the mask/weighting for the target patch
+            virtual const scalarField& tgtMask() const;
+
             //- Overlap tolerance
             inline static scalar tolerance();
 
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchI.H b/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchI.H
index 42b4119667..50f4022ed6 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchI.H
+++ b/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchI.H
@@ -27,18 +27,6 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-inline void Foam::cyclicACMIPolyPatch::setUpdated(const bool flag) const
-{
-    updated_ = flag;
-}
-
-
-inline bool Foam::cyclicACMIPolyPatch::updated() const
-{
-    return updated_;
-}
-
-
 inline const Foam::word& Foam::cyclicACMIPolyPatch::nonOverlapPatchName() const
 {
     return nonOverlapPatchName_;
@@ -69,10 +57,8 @@ inline const Foam::scalarField& Foam::cyclicACMIPolyPatch::mask() const
     {
         return srcMask_;
     }
-    else
-    {
-        return neighbPatch().tgtMask();
-    }
+
+    return neighbPatch().tgtMask();
 }
 
 
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
index 51f14ccb5e..04cb5b6bb3 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
+++ b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
@@ -27,13 +27,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "cyclicAMIPolyPatch.H"
-#include "transformField.H"
 #include "SubField.H"
-#include "polyMesh.H"
 #include "Time.H"
-#include "addToRunTimeSelectionTable.H"
 #include "faceAreaIntersect.H"
-#include "ops.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,6 +42,7 @@ namespace Foam
     addToRunTimeSelectionTable(polyPatch, cyclicAMIPolyPatch, dictionary);
 }
 
+const Foam::scalar Foam::cyclicAMIPolyPatch::tolerance_ = 1e-10;
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
@@ -295,77 +293,114 @@ void Foam::cyclicAMIPolyPatch::resetAMI
     const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod
 ) const
 {
-    if (owner())
+    resetAMI(boundaryMesh().mesh().points(), AMIMethod);
+}
+
+
+void Foam::cyclicAMIPolyPatch::resetAMI
+(
+    const UList<point>& points,
+    const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod
+) const
+{
+    DebugInFunction << endl;
+
+    if (!owner())
     {
-        AMIPtr_.clear();
+        return;
+    }
 
-        const polyPatch& nbr = neighbPatch();
-        pointField nbrPoints
-        (
-            neighbPatch().boundaryMesh().mesh().points(),
-            neighbPatch().meshPoints()
-        );
+    AMIPtr_.clear();
 
-        if (debug)
+    const cyclicAMIPolyPatch& nbr = neighbPatch();
+    pointField srcPoints(points, meshPoints());
+    pointField nbrPoints(points, neighbPatch().meshPoints());
+
+    if (debug)
+    {
+        const Time& t = boundaryMesh().mesh().time();
+        OFstream os(t.path()/name() + "_neighbourPatch-org.obj");
+        meshTools::writeOBJ(os, neighbPatch().localFaces(), nbrPoints);
+    }
+
+    label patchSize0 = size();
+    label nbrPatchSize0 = nbr.size();
+
+    if (createAMIFaces_)
+    {
+        // AMI is created based on the original patch faces (non-extended patch)
+        if (srcFaceIDs_.size())
         {
-            const Time& t = boundaryMesh().mesh().time();
-            OFstream os(t.path()/name() + "_neighbourPatch-org.obj");
-            meshTools::writeOBJ(os, neighbPatch().localFaces(), nbrPoints);
+            patchSize0 = srcFaceIDs_.size();
         }
-
-        // Transform neighbour patch to local system
-        transformPosition(nbrPoints);
-        primitivePatch nbrPatch0
-        (
-            SubList<face>
-            (
-                nbr.localFaces(),
-                nbr.size()
-            ),
-            nbrPoints
-        );
-
-        if (debug)
+        if (tgtFaceIDs_.size())
         {
-            const Time& t = boundaryMesh().mesh().time();
-            OFstream osN(t.path()/name() + "_neighbourPatch-trans.obj");
-            meshTools::writeOBJ(osN, nbrPatch0.localFaces(), nbrPoints);
-
-            OFstream osO(t.path()/name() + "_ownerPatch.obj");
-            meshTools::writeOBJ(osO, this->localFaces(), localPoints());
+            nbrPatchSize0 = tgtFaceIDs_.size();
         }
+    }
 
-        // Construct/apply AMI interpolation to determine addressing and weights
-        AMIPtr_.reset
+    // Transform neighbour patch to local system
+    transformPosition(nbrPoints);
+    primitivePatch nbrPatch0
+    (
+        SubList<face>
         (
-            new AMIPatchToPatchInterpolation
-            (
-                *this,
-                nbrPatch0,
-                surfPtr(),
-                faceAreaIntersect::tmMesh,
-                AMIRequireMatch_,
-                AMIMethod,
-                AMILowWeightCorrection_,
-                AMIReverse_
-            )
-        );
+            nbr.localFaces(),
+            nbrPatchSize0
+        ),
+        nbrPoints
+    );
+    primitivePatch patch0
+    (
+        SubList<face>
+        (
+            localFaces(),
+            patchSize0
+        ),
+        srcPoints
+    );
 
-        if (debug)
-        {
-            Pout<< "cyclicAMIPolyPatch : " << name()
-                << " constructed AMI with " << nl
-                << "    " << "srcAddress:" << AMIPtr_().srcAddress().size()
-                << nl
-                << "    " << "tgAddress :" << AMIPtr_().tgtAddress().size()
-                << nl << endl;
-        }
+
+    if (debug)
+    {
+        const Time& t = boundaryMesh().mesh().time();
+        OFstream osN(t.path()/name() + "_neighbourPatch-trans.obj");
+        meshTools::writeOBJ(osN, nbrPatch0.localFaces(), nbrPoints);
+
+        OFstream osO(t.path()/name() + "_ownerPatch.obj");
+        meshTools::writeOBJ(osO, this->localFaces(), localPoints());
+    }
+
+    // Construct/apply AMI interpolation to determine addressing and weights
+    AMIPtr_.reset
+    (
+        new AMIPatchToPatchInterpolation
+        (
+            patch0, // *this,
+            nbrPatch0,
+            surfPtr(),
+            faceAreaIntersect::tmMesh,
+            AMIRequireMatch_,
+            AMIMethod,
+            AMILowWeightCorrection_,
+            AMIReverse_
+        )
+    );
+
+    // Set the updated flag
+    updated_ = true;
+
+    if (debug)
+    {
+        AMIPtr_->checkSymmetricWeights(true);
     }
 }
 
 
 void Foam::cyclicAMIPolyPatch::calcTransforms()
 {
+    DebugInFunction << endl;
+
     const cyclicAMIPolyPatch& half0 = *this;
     vectorField half0Areas(half0.size());
     forAll(half0, facei)
@@ -402,22 +437,29 @@ void Foam::cyclicAMIPolyPatch::calcTransforms()
 
 void Foam::cyclicAMIPolyPatch::initGeometry(PstreamBuffers& pBufs)
 {
-    // The AMI is no longer valid. Leave it up to demand-driven calculation
-    AMIPtr_.clear();
+    DebugInFunction << endl;
+
+    // Only resetting the AMI if not creating additional AMI faces
+    // - When creating additional faces the AMI is updated externally by the
+    //   dynamic mesh via the setTopology function.
+    if (!createAMIFaces_ && canResetAMI())
+    {
+        resetAMI(AMIMethod_);
+    }
 
     polyPatch::initGeometry(pBufs);
 
     // Early calculation of transforms so e.g. cyclicACMI can use them.
     // Note: also triggers primitiveMesh face centre. Note that cell
     // centres should -not- be calculated
-    // since e.g. cyclicACMI override face areas
+    // since e.g. cyclicACMI overrides face areas
     calcTransforms();
 }
 
 
 void Foam::cyclicAMIPolyPatch::calcGeometry(PstreamBuffers& pBufs)
 {
-    // All geometry done inside initGeometry
+    DebugInFunction << endl;
 }
 
 
@@ -427,14 +469,31 @@ void Foam::cyclicAMIPolyPatch::initMovePoints
     const pointField& p
 )
 {
-    // The AMI is no longer valid. Leave it up to demand-driven calculation
-    AMIPtr_.clear();
-
-    polyPatch::initMovePoints(pBufs, p);
+    DebugInFunction << endl;
 
     // See below. Clear out any local geometry
     primitivePatch::movePoints(p);
 
+    // Note: processorPolyPatch::initMovePoints calls
+    // processorPolyPatch::initGeometry which will trigger calculation of
+    // patch faceCentres() and cell volumes...
+
+    if (createAMIFaces_)
+    {
+        // Note: AMI should have been updated in setTopology
+
+        // faceAreas() and faceCentres() have been reset and will be
+        // recalculated on-demand using the mesh points and no longer
+        // correspond to the scaled areas!
+        restoreScaledGeometry();
+
+        // deltas need to be recalculated to use new face centres!
+    }
+    else
+    {
+        resetAMI(p, AMIMethod_);
+    }
+
     // Early calculation of transforms. See above.
     calcTransforms();
 }
@@ -446,30 +505,55 @@ void Foam::cyclicAMIPolyPatch::movePoints
     const pointField& p
 )
 {
-    polyPatch::movePoints(pBufs, p);
+    DebugInFunction << endl;
 
-    // All transformation tensors already done in initMovePoints
+//    Note: not calling movePoints since this will undo our manipulations!
+//    polyPatch::movePoints(pBufs, p);
+/*
+    polyPatch::movePoints
+     -> primitivePatch::movePoints
+        -> primitivePatch::clearGeom:
+    deleteDemandDrivenData(localPointsPtr_);
+    deleteDemandDrivenData(faceCentresPtr_);
+    deleteDemandDrivenData(faceAreasPtr_);
+    deleteDemandDrivenData(magFaceAreasPtr_);
+    deleteDemandDrivenData(faceNormalsPtr_);
+    deleteDemandDrivenData(pointNormalsPtr_);
+*/
 }
 
 
 void Foam::cyclicAMIPolyPatch::initUpdateMesh(PstreamBuffers& pBufs)
 {
-    // The AMI is no longer valid. Leave it up to demand-driven calculation
-    AMIPtr_.clear();
+    DebugInFunction << endl;
 
     polyPatch::initUpdateMesh(pBufs);
+
+    if (createAMIFaces_ && boundaryMesh().mesh().topoChanging() && owner())
+    {
+        setAMIFaces();
+    }
 }
 
 
 void Foam::cyclicAMIPolyPatch::updateMesh(PstreamBuffers& pBufs)
 {
+    DebugInFunction << endl;
+
+    // Note: this clears out cellCentres(), faceCentres() and faceAreas()
     polyPatch::updateMesh(pBufs);
 }
 
 
 void Foam::cyclicAMIPolyPatch::clearGeom()
 {
-    AMIPtr_.clear();
+    DebugInFunction << endl;
+
+    if (!updatingAMI_)
+    {
+//DEBUG("*** CLEARING AMI ***");
+        AMIPtr_.clear();
+    }
     polyPatch::clearGeom();
 }
 
@@ -488,6 +572,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
 )
 :
     coupledPolyPatch(name, size, start, index, bm, patchType, transform),
+    updated_(false),
     nbrPatchName_(word::null),
     nbrPatchID_(-1),
     rotationAxis_(Zero),
@@ -501,7 +586,13 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     AMIRequireMatch_(true),
     AMILowWeightCorrection_(-1.0),
     surfPtr_(nullptr),
-    surfDict_(fileName("surface"))
+    surfDict_(fileName("surface")),
+    createAMIFaces_(false),
+    updatingAMI_(true),
+    srcFaceIDs_(),
+    tgtFaceIDs_(),
+    faceAreas0_(),
+    faceCentres0_()
 {
     // Neighbour patch might not be valid yet so no transformation
     // calculation possible
@@ -518,6 +609,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
 )
 :
     coupledPolyPatch(name, dict, index, bm, patchType),
+    updated_(false),
     nbrPatchName_(dict.getOrDefault<word>("neighbourPatch", "")),
     coupleGroup_(dict),
     nbrPatchID_(-1),
@@ -545,7 +637,13 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     AMIRequireMatch_(true),
     AMILowWeightCorrection_(dict.getOrDefault("lowWeightCorrection", -1.0)),
     surfPtr_(nullptr),
-    surfDict_(dict.subOrEmptyDict("surface"))
+    surfDict_(dict.subOrEmptyDict("surface")),
+    createAMIFaces_(false),
+    updatingAMI_(true),
+    srcFaceIDs_(),
+    tgtFaceIDs_(),
+    faceAreas0_(),
+    faceCentres0_()
 {
     if (nbrPatchName_ == word::null && !coupleGroup_.valid())
     {
@@ -605,6 +703,20 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
 
     // Neighbour patch might not be valid yet so no transformation
     // calculation possible
+
+    // If topology change, recover the sizes of the original patches
+    label srcSize0 = 0;
+    if (dict.readIfPresent("srcSize", srcSize0))
+    {
+        srcFaceIDs_.setSize(srcSize0);
+        createAMIFaces_ = true;
+    }
+    label tgtSize0 = 0;
+    if (dict.readIfPresent("tgtSize", tgtSize0))
+    {
+        tgtFaceIDs_.setSize(tgtSize0);
+        createAMIFaces_ = true;
+    }
 }
 
 
@@ -615,6 +727,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
 )
 :
     coupledPolyPatch(pp, bm),
+    updated_(false),
     nbrPatchName_(pp.nbrPatchName_),
     coupleGroup_(pp.coupleGroup_),
     nbrPatchID_(-1),
@@ -629,7 +742,13 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     AMIRequireMatch_(pp.AMIRequireMatch_),
     AMILowWeightCorrection_(pp.AMILowWeightCorrection_),
     surfPtr_(nullptr),
-    surfDict_(pp.surfDict_)
+    surfDict_(pp.surfDict_),
+    createAMIFaces_(false),
+    updatingAMI_(true),
+    srcFaceIDs_(),
+    tgtFaceIDs_(),
+    faceAreas0_(),
+    faceCentres0_()
 {
     // Neighbour patch might not be valid yet so no transformation
     // calculation possible
@@ -647,6 +766,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
 )
 :
     coupledPolyPatch(pp, bm, index, newSize, newStart),
+    updated_(false),
     nbrPatchName_(nbrPatchName),
     coupleGroup_(pp.coupleGroup_),
     nbrPatchID_(-1),
@@ -661,7 +781,13 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     AMIRequireMatch_(pp.AMIRequireMatch_),
     AMILowWeightCorrection_(pp.AMILowWeightCorrection_),
     surfPtr_(nullptr),
-    surfDict_(pp.surfDict_)
+    surfDict_(pp.surfDict_),
+    createAMIFaces_(false),
+    updatingAMI_(true),
+    srcFaceIDs_(),
+    tgtFaceIDs_(),
+    faceAreas0_(),
+    faceCentres0_()
 {
     if (nbrPatchName_ == name())
     {
@@ -686,6 +812,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
 )
 :
     coupledPolyPatch(pp, bm, index, mapAddressing, newStart),
+    updated_(false),
     nbrPatchName_(pp.nbrPatchName_),
     coupleGroup_(pp.coupleGroup_),
     nbrPatchID_(-1),
@@ -700,13 +827,13 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     AMIRequireMatch_(pp.AMIRequireMatch_),
     AMILowWeightCorrection_(pp.AMILowWeightCorrection_),
     surfPtr_(nullptr),
-    surfDict_(pp.surfDict_)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::cyclicAMIPolyPatch::~cyclicAMIPolyPatch()
+    surfDict_(pp.surfDict_),
+    createAMIFaces_(false),
+    updatingAMI_(true),
+    srcFaceIDs_(),
+    tgtFaceIDs_(),
+    faceAreas0_(),
+    faceCentres0_()
 {}
 
 
@@ -1105,6 +1232,12 @@ void Foam::cyclicAMIPolyPatch::write(Ostream& os) const
     {
         surfDict_.writeEntry(surfDict_.dictName(), os);
     }
+
+    if (createAMIFaces_)
+    {
+        os.writeEntry("srcSize", srcFaceIDs_.size());
+        os.writeEntry("tgtSize", tgtFaceIDs_.size());
+    }
 }
 
 
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H
index 25259ffacd..84fab30b99 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H
+++ b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H
@@ -74,6 +74,9 @@ protected:
 
     // Protected data
 
+        //- Flag to indicate that AMI has been updated
+        mutable bool updated_;
+
         //- Name of other half
         mutable word nbrPatchName_;
 
@@ -129,9 +132,47 @@ protected:
         const dictionary surfDict_;
 
 
+        // Change of topology as AMI is updated
+
+            //- Flag to indicate that new AMI faces will created
+            //  Set by the call to changeTopology
+            mutable bool createAMIFaces_;
+
+            mutable bool updatingAMI_;
+
+            labelListList srcFaceIDs_;
+
+            labelListList tgtFaceIDs_;
+
+            //- Temporary storage for AMI face areas
+            mutable vectorField faceAreas0_;
+
+            //- Temporary storage for AMI face centres
+            mutable vectorField faceCentres0_;
+
+
     // Protected Member Functions
 
-        //- Reset the AMI interpolator
+        // Topology change
+
+            virtual bool removeAMIFaces(polyTopoChange& topoChange);
+
+            virtual bool addAMIFaces(polyTopoChange& topoChange);
+
+            virtual void setAMIFaces();
+
+            virtual void restoreScaledGeometry();
+
+
+        //- Reset the AMI interpolator, supply patch points
+        virtual void resetAMI
+        (
+            const UList<point>& points,
+            const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod =
+                AMIPatchToPatchInterpolation::imFaceAreaWeight
+        ) const;
+
+        //-  Reset the AMI interpolator, use current patch points
         virtual void resetAMI
         (
             const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod =
@@ -274,15 +315,45 @@ public:
 
 
     //- Destructor
-    virtual ~cyclicAMIPolyPatch();
+    virtual ~cyclicAMIPolyPatch() = default;
 
 
     // Member Functions
 
         // Access
 
+            //- Tolerance used e.g. for area calculations/limits
+            static const scalar tolerance_;
+
+            //- Flag to indicate whether the AMI can be reset
+            inline bool canResetAMI() const;
+
+             //- Reset the updated flag
+            inline void setUpdated(bool flag) const;
+
+            //- Return access to the updated flag
+            inline bool updated() const;
+
+            //- Return access to the createAMIFaces flag
+            inline bool createAMIFaces() const;
+
+            //- Return access to the updated flag
+            inline bool updatingAMI() const;
+
+void clearAMI() const
+{
+    AMIPtr_.clear();
+}
+
+            //- Return true if this patch changes the mesh topology
+            //  True when createAMIFaces is true
+            virtual bool changeTopology() const;
+
+            //- Set topology changes in the polyTopoChange object
+            virtual bool setTopology(polyTopoChange& topoChange);
+
             //- Is patch 'coupled'. Note that on AMI the geometry is not
-            //  coupled but the fields are!
+            //- coupled but the fields are!
             virtual bool coupled() const
             {
                 return false;
@@ -306,9 +377,23 @@ public:
             //- Return a reference to the AMI interpolator
             const AMIPatchToPatchInterpolation& AMI() const;
 
+            //- Helper function to return the weights
+            inline const scalarListList& weights() const;
+
+            //- Helper function to return the weights sum
+            inline const scalarField& weightsSum() const;
+
             //- Return true if applying the low weight correction
             bool applyLowWeightCorrection() const;
 
+            //- Return access to the initial face areas
+            //  Used for topology change
+            inline vectorField& faceAreas0() const;
+
+            //- Return access to the initial face centres
+            //  Used for topology change
+            inline vectorField& faceCentres0() const;
+
 
             // Transformations
 
@@ -388,7 +473,7 @@ public:
         );
 
         //- Initialize ordering for primitivePatch. Does not
-        //  refer to *this (except for name() and type() etc.)
+        //- refer to *this (except for name() and type() etc.)
         virtual void initOrder
         (
             PstreamBuffers&,
@@ -409,7 +494,7 @@ public:
         ) const;
 
         //- Return face index on neighbour patch which shares point p
-        //  following trajectory vector n
+        //- following trajectory vector n
         label pointFace
         (
             const label facei,
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchI.H b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchI.H
index ccc36fc300..6d55bf1841 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchI.H
+++ b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchI.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,6 +28,36 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+inline bool Foam::cyclicAMIPolyPatch::canResetAMI() const
+{
+    return Pstream::parRun() || !boundaryMesh().mesh().time().processorCase();
+}
+
+
+inline void Foam::cyclicAMIPolyPatch::setUpdated(const bool flag) const
+{
+    updated_ = flag;
+}
+
+
+inline bool Foam::cyclicAMIPolyPatch::updated() const
+{
+    return updated_;
+}
+
+
+inline bool Foam::cyclicAMIPolyPatch::createAMIFaces() const
+{
+    return createAMIFaces_;
+}
+
+
+inline bool Foam::cyclicAMIPolyPatch::updatingAMI() const
+{
+    return updatingAMI_;
+}
+
+
 inline const Foam::word& Foam::cyclicAMIPolyPatch::neighbPatchName() const
 {
     if (nbrPatchName_.empty())
@@ -39,6 +70,38 @@ inline const Foam::word& Foam::cyclicAMIPolyPatch::neighbPatchName() const
     return nbrPatchName_;
 }
 
+inline const Foam::scalarListList& Foam::cyclicAMIPolyPatch::weights() const
+{
+    if (owner())
+    {
+        return AMI().srcWeights();
+    }
+
+    return neighbPatch().AMI().tgtWeights();
+}
+
+
+inline const Foam::scalarField& Foam::cyclicAMIPolyPatch::weightsSum() const
+{
+    if (owner())
+    {
+        return AMI().srcWeightsSum();
+    }
+
+    return neighbPatch().AMI().tgtWeightsSum();
+}
+
+
+inline Foam::vectorField& Foam::cyclicAMIPolyPatch::faceAreas0() const
+{
+    return faceAreas0_;
+}
+
+
+inline Foam::vectorField& Foam::cyclicAMIPolyPatch::faceCentres0() const
+{
+    return faceCentres0_;
+}
 
 inline const Foam::vector& Foam::cyclicAMIPolyPatch::rotationAxis() const
 {
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchTopologyChange.C b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchTopologyChange.C
new file mode 100644
index 0000000000..f7e33910da
--- /dev/null
+++ b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchTopologyChange.C
@@ -0,0 +1,676 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     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 "cyclicAMIPolyPatch.H"
+#include "SubField.H"
+#include "vectorList.H"
+#include "polyTopoChange.H"
+
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
+
+void Foam::cyclicAMIPolyPatch::restoreScaledGeometry()
+{
+    DebugInFunction << endl;
+
+    // Note: only used for topology update (createAMIFaces_ flag)
+    if (!createAMIFaces_)
+    {
+        FatalErrorInFunction
+            << "Attempted to perform topology update when createAMIFaces_ "
+            << "flag is set to false"
+            << abort(FatalError);
+    }
+
+    if (boundaryMesh().mesh().hasCellVolumes())
+    {
+        WarningInFunction
+            << "Mesh already has volumes set!"
+            << endl;
+    }
+
+    vectorField::subField faceAreas = this->faceAreas();
+    vectorField::subField faceCentres = this->faceCentres();
+
+    if (debug)
+    {
+        Info<< "Patch:" << name() << " before: sum(mag(faceAreas)):"
+            << gSum(mag(faceAreas)) << endl;
+        Info<< "Patch:" << name() << " before: sum(mag(faceAreas0)):"
+            << gSum(mag(faceAreas0_)) << endl;
+    }
+
+    faceAreas = faceAreas0_;
+    Info<< "WARNING: ==============================" << endl;
+    Info<< "WARNING: Not updating face cell centres" << endl;
+    Info<< "WARNING: ==============================" << endl;
+    //faceCentres = faceCentres0_;
+    faceAreas0_.clear();
+    faceCentres0_.clear();
+
+    if (debug)
+    {
+        Info<< "Patch:" << name() << " after: sum(mag(faceAreas)):"
+            << gSum(mag(faceAreas)) << endl;
+        Info<< "Patch:" << name() << " after: sum(mag(faceAreas0)):"
+            << gSum(mag(faceAreas0_)) << endl;
+    }
+}
+
+
+bool Foam::cyclicAMIPolyPatch::removeAMIFaces(polyTopoChange& topoChange)
+{
+    DebugInFunction << endl;
+
+    // Note: only used for topology update (createAMIFaces_ flag)
+    if (!createAMIFaces_)
+    {
+        FatalErrorInFunction
+            << "Attempted to perform topology update when createAMIFaces_ "
+            << "flag is set to false"
+            << abort(FatalError);
+    }
+
+    if (!owner())
+    {
+        return false;
+    }
+
+    bool changeRequired = false;
+
+    // Remove any faces that we inserted to create the 1-to-1 match...
+
+    const cyclicAMIPolyPatch& nbr = neighbPatch();
+
+    const label newSrcFaceStart = srcFaceIDs_.size();
+
+    if (newSrcFaceStart != 0)
+    {
+        for (label facei = newSrcFaceStart; facei < size(); ++facei)
+        {
+            changeRequired = true;
+            label meshFacei = start() + facei;
+            topoChange.removeFace(meshFacei, -1);
+        }
+    }
+
+    const label newTgtFaceStart = tgtFaceIDs_.size();
+
+    if (newTgtFaceStart != 0)
+    {
+        for (label facei = newTgtFaceStart; facei < nbr.size(); ++facei)
+        {
+            changeRequired = true;
+            label meshFacei = nbr.start() + facei;
+            topoChange.removeFace(meshFacei, -1);
+        }
+    }
+
+    srcFaceIDs_.clear();
+    tgtFaceIDs_.clear();
+
+    return changeRequired;
+}
+
+
+bool Foam::cyclicAMIPolyPatch::addAMIFaces(polyTopoChange& topoChange)
+{
+    DebugInFunction << endl;
+
+    // Note: only used for topology update (createAMIFaces_ flag = true)
+    if (!createAMIFaces_)
+    {
+        FatalErrorInFunction
+            << "Attempted to perform topology update when createAMIFaces_ "
+            << "flag is set to false"
+            << abort(FatalError);
+    }
+
+    bool changedFaces = false;
+    const cyclicAMIPolyPatch& nbr = neighbPatch();
+
+    polyMesh& mesh = const_cast<polyMesh&>(boundaryMesh().mesh());
+    const faceZoneMesh& faceZones = mesh.faceZones();
+
+    // First face address and weight are used to manipulate the
+    // original face - all other addresses and weights are used to
+    // create additional faces
+    const labelListList& srcToTgtAddr = AMI().srcAddress();
+    const labelListList& tgtToSrcAddr = AMI().tgtAddress();
+
+    const label nSrcFace = srcToTgtAddr.size();
+    const label nTgtFace = tgtToSrcAddr.size();
+
+    srcFaceIDs_.setSize(nSrcFace);
+    tgtFaceIDs_.setSize(nTgtFace);
+
+    label nNewSrcFaces = 0;
+    forAll(srcToTgtAddr, srcFacei)
+    {
+        const labelList& tgtAddr = srcToTgtAddr[srcFacei];
+
+        // No tgt faces linked to srcFacei (ACMI)
+        if (tgtAddr.empty()) continue;
+
+        srcFaceIDs_[srcFacei].setSize(tgtAddr.size());
+        srcFaceIDs_[srcFacei][0] = srcFacei;
+
+        const label meshFacei = start() + srcFacei;
+        for (label addri = 1; addri < tgtAddr.size(); ++addri)
+        {
+            changedFaces = true;
+
+            // Note: new faces reuse originating face points
+            // - but areas are scaled by the weights (later)
+
+            // New source face for each target face address
+            srcFaceIDs_[srcFacei][addri] = nNewSrcFaces + nSrcFace;
+            ++nNewSrcFaces;
+            (void)topoChange.addFace
+            (
+                mesh.faces()[meshFacei],        // modified face
+                mesh.faceOwner()[meshFacei],    // owner
+                -1,                             // neighbour
+                -1,                             // master point
+                -1,                             // master edge
+                meshFacei,                      // master face
+                false,                          // face flip
+                index(),                        // patch for face
+                faceZones.whichZone(meshFacei), // zone for original face
+                false                           // face flip in zone
+            );
+        }
+    }
+
+    label nNewTgtFaces = 0;
+    forAll(tgtToSrcAddr, tgtFacei)
+    {
+        const labelList& srcAddr = tgtToSrcAddr[tgtFacei];
+
+        // No src faces linked to tgtFacei (ACMI)
+        if (srcAddr.empty()) continue;
+
+        tgtFaceIDs_[tgtFacei].setSize(srcAddr.size());
+        tgtFaceIDs_[tgtFacei][0] = tgtFacei;
+
+        const label meshFacei = nbr.start() + tgtFacei;
+        for (label addri = 1; addri < srcAddr.size(); ++addri)
+        {
+            changedFaces = true;
+
+            // Note: new faces reuse originating face points
+            // - but areas are scaled by the weights (later)
+
+            // New target face for each source face address
+            tgtFaceIDs_[tgtFacei][addri] = nNewTgtFaces + nTgtFace;
+            ++nNewTgtFaces;
+
+            (void)topoChange.addFace
+            (
+                mesh.faces()[meshFacei],        // modified face
+                mesh.faceOwner()[meshFacei],    // owner
+                -1,                             // neighbour
+                -1,                             // master point
+                -1,                             // master edge
+                meshFacei,                      // master face
+                false,                          // face flip
+                nbr.index(),                    // patch for face
+                faceZones.whichZone(meshFacei), // zone for original face
+                false                           // face flip in zone
+            );
+        }
+    }
+
+    Info<< "AMI: Patch " << name() << " additional faces: "
+        << returnReduce(nNewSrcFaces, sumOp<label>()) << nl
+        << "AMI: Patch " << nbr.name() << " additional faces: "
+        << returnReduce(nNewTgtFaces, sumOp<label>())
+        << endl;
+
+    if (debug)
+    {
+        Pout<< "New faces - " << name() << ": " << nNewSrcFaces
+            << " "  << nbr.name() << ": " << nNewTgtFaces << endl;
+    }
+
+    return returnReduce(changedFaces, orOp<bool>());
+}
+
+
+void Foam::cyclicAMIPolyPatch::setAMIFaces()
+{
+    // Create new mesh faces so that there is a 1-to-1 correspondence
+    // between faces on each side of the AMI
+
+    // Note: only used for topology update (createAMIFaces_ flag)
+    if (!createAMIFaces_)
+    {
+        FatalErrorInFunction
+            << "Attempted to perform topology update when createAMIFaces_ "
+            << "flag is set to false"
+            << abort(FatalError);
+    }
+
+
+    DebugInFunction << endl;
+
+    if (!owner())
+    {
+        return;
+    }
+
+    const cyclicAMIPolyPatch& nbr = neighbPatch();
+
+    vectorField& nbrFaceAreas0 = nbr.faceAreas0();
+    vectorField& nbrFaceCentres0 = nbr.faceCentres0();
+
+    // Scale the new face areas and set the centroids
+    // Note:
+    // - storing local copies so that they can be re-applied after the call to
+    //   movePoints that will reset any changes to the areas and centroids
+    //
+    // - For AMI, src and tgt patches should be the same
+    // - For ACMI they are likely to be different!
+    faceAreas0_ = faceAreas();
+    faceCentres0_ = faceCentres();
+    nbrFaceAreas0 = nbr.faceAreas();
+    nbrFaceCentres0 = nbr.faceCentres();
+
+    // Original AMI info (based on the mesh state when the AMI was evaluated)
+    const labelListList& srcToTgtAddr0 = AMIPtr_->srcAddress();
+    const labelListList& tgtToSrcAddr0 = AMIPtr_->tgtAddress();
+    const pointListList& srcCtr0 = AMIPtr_->srcCentroids();
+    const scalarListList& srcToTgtWght0 = AMIPtr_->srcWeights();
+
+    // New addressing on new mesh (extended by polyTopoChange)
+    labelListList srcToTgtAddr1(size(), labelList());
+    labelListList tgtToSrcAddr1(nbr.size(), labelList());
+
+    // Need to calc new parallel maps (mesh has changed since AMI was computed)
+    autoPtr<mapDistribute> srcToTgtMap1;
+    autoPtr<mapDistribute> tgtToSrcMap1;
+
+    if (AMIPtr_->singlePatchProc() == -1)
+    {
+       // Parallel running
+
+        // Global index based on old patch sizes (when AMI was computed)
+        globalIndex globalSrcFaces0(srcToTgtAddr0.size());
+        globalIndex globalTgtFaces0(tgtToSrcAddr0.size());
+
+        // Global index based on new patch sizes
+        globalIndex globalSrcFaces1(size());
+        globalIndex globalTgtFaces1(nbr.size());
+
+
+        // Gather source side info
+        // =======================
+
+        // Note: using new global index for addressing, and distributed using
+        // the old AMI map
+        labelListList newTgtGlobalFaces(tgtFaceIDs_);
+        forAll(newTgtGlobalFaces, tgtFacei)
+        {
+            globalTgtFaces1.inplaceToGlobal(newTgtGlobalFaces[tgtFacei]);
+        }
+        AMIPtr_->tgtMap().distribute(newTgtGlobalFaces);
+
+        // Now have new tgt face indices for each src face
+
+        labelList globalSrcFaceIDs(identity(srcToTgtAddr0.size()));
+        globalSrcFaces0.inplaceToGlobal(globalSrcFaceIDs);
+        AMIPtr_->srcMap().distribute(globalSrcFaceIDs);
+        // globalSrcFaceIDs now has remote data for each srcFacei0 known to the
+        // tgt patch
+
+        List<List<point>> globalSrcCtrs0(srcCtr0);
+        AMIPtr_->srcMap().distribute(globalSrcCtrs0);
+
+        labelList globalTgtFaceIDs(identity(tgtToSrcAddr0.size()));
+        globalTgtFaces0.inplaceToGlobal(globalTgtFaceIDs);
+        AMIPtr_->tgtMap().distribute(globalTgtFaceIDs);
+        // globalTgtFaceIDs now has remote data for each tgtFacei0 known to the
+        // src patch
+
+        // For debug - send tgt face centres and compare against mapped src
+        // face centres
+        //List<List<point>> globalTgtCtrs0(tgtCtr0);
+        //AMIPtr_->tgtMap().distribute(globalTgtCtrs0);
+
+        labelListList globalTgtToSrcAddr(tgtToSrcAddr0);
+        forAll(tgtToSrcAddr0, tgtFacei0)
+        {
+            forAll(tgtToSrcAddr0[tgtFacei0], addri)
+            {
+                const label globalSrcFacei =
+                    globalSrcFaceIDs[tgtToSrcAddr0[tgtFacei0][addri]];
+                globalTgtToSrcAddr[tgtFacei0][addri] = globalSrcFacei;
+            }
+        }
+        AMIPtr_->tgtMap().distribute(globalTgtToSrcAddr);
+
+        labelListList globalSrcToTgtAddr(srcToTgtAddr0);
+        forAll(srcToTgtAddr0, srcFacei0)
+        {
+            forAll(srcToTgtAddr0[srcFacei0], addri)
+            {
+                const label globalTgtFacei =
+                    globalTgtFaceIDs[srcToTgtAddr0[srcFacei0][addri]];
+                globalSrcToTgtAddr[srcFacei0][addri] = globalTgtFacei;
+            }
+        }
+        AMIPtr_->srcMap().distribute(globalSrcToTgtAddr);
+
+        label nError = 0;
+        forAll(srcToTgtAddr0, srcFacei0)
+        {
+            const labelList& newSrcFaces = srcFaceIDs_[srcFacei0];
+
+            forAll(newSrcFaces, i)
+            {
+                const label srcFacei1 = newSrcFaces[i];
+
+                // What index did srcFacei0 appear in tgtToSrc0 list?
+                // - if first index, all ok
+                // - else tgt face has been moved to according to tgtFaceIDs_
+                const label tgtFacei0 = srcToTgtAddr0[srcFacei0][i];
+                const label addri =
+                    globalTgtToSrcAddr[tgtFacei0].find
+                    (
+                        globalSrcFaceIDs[srcFacei0]
+                    );
+
+                if (addri == -1)
+                {
+                    ++nError;
+                    continue;
+
+                    if (debug)
+                    {
+                        Pout<< "Unable to find global source face "
+                            << globalSrcFaceIDs[srcFacei0]
+                            << " in globalTgtToSrcAddr[" << tgtFacei0 << "]: "
+                            << globalTgtToSrcAddr[tgtFacei0]
+                            << endl;
+                    }
+                }
+
+                const label tgtFacei1 = newTgtGlobalFaces[tgtFacei0][addri];
+
+                // Sanity check to see that we've picked the correct face
+                // point tgtCtr0(globalTgtCtrs0[tgtFacei0][addri]);
+                // Pout<< "srcCtr:" << srcCtr0[srcFacei0][i]
+                //     << " tgtCtr:" << tgtCtr0 << endl;
+
+                srcToTgtAddr1[srcFacei1] = labelList(1, tgtFacei1);
+                faceAreas0_[srcFacei1] *= srcToTgtWght0[srcFacei0][i];
+                faceCentres0_[srcFacei1] = srcCtr0[srcFacei0][i];
+            }
+        }
+
+        if (nError)
+        {
+            FatalErrorInFunction
+                << "Unable to find " << nError << " global source faces"
+                << abort(FatalError);
+        }
+
+
+        // Gather Target side info
+        // =======================
+
+        labelListList newSrcGlobalFaces(srcFaceIDs_);
+        forAll(newSrcGlobalFaces, srcFacei)
+        {
+            globalSrcFaces1.inplaceToGlobal(newSrcGlobalFaces[srcFacei]);
+        }
+
+        AMIPtr_->srcMap().distribute(newSrcGlobalFaces);
+
+        // Now have new src face indices for each tgt face
+        forAll(tgtToSrcAddr0, tgtFacei0)
+        {
+            const labelList& newTgtFaces = tgtFaceIDs_[tgtFacei0];
+            forAll(newTgtFaces, i)
+            {
+                const label srcFacei0 = tgtToSrcAddr0[tgtFacei0][i];
+
+                const label addri =
+                    globalSrcToTgtAddr[srcFacei0].find
+                    (
+                        globalTgtFaceIDs[tgtFacei0]
+                    );
+
+                if (addri == -1)
+                {
+                    ++nError;
+                    continue;
+
+                    if (debug)
+                    {
+                        Pout<< "Unable to find global target face "
+                            << globalTgtFaceIDs[tgtFacei0]
+                            << " in globalSrcToTgtAddr[" << srcFacei0 << "]: "
+                            << globalSrcToTgtAddr[srcFacei0]
+                            << endl;
+                    }
+                }
+
+                const label srcFacei1 = newSrcGlobalFaces[srcFacei0][addri];
+
+                // Sanity check to see that we've picked the correct face
+                point srcCtr0(globalSrcCtrs0[srcFacei0][addri]);
+                reverseTransformPosition(srcCtr0, srcFacei0);
+
+                const label tgtFacei1 = newTgtFaces[i];
+                tgtToSrcAddr1[tgtFacei1] = labelList(1, srcFacei1);
+                nbrFaceCentres0[tgtFacei1] = srcCtr0;
+            }
+        }
+
+        if (nError)
+        {
+            FatalErrorInFunction
+                << "Unable to find " << nError << " global target faces"
+                << abort(FatalError);
+        }
+
+        // Update the maps
+        {
+            List<Map<label>> cMap;
+            srcToTgtMap1.reset
+            (
+                new mapDistribute(globalSrcFaces1, tgtToSrcAddr1, cMap)
+            );
+        }
+        {
+            List<Map<label>> cMap;
+            tgtToSrcMap1.reset
+            (
+                new mapDistribute(globalTgtFaces1, srcToTgtAddr1, cMap)
+            );
+        }
+
+        // Reset tgt patch areas using the new map
+        vectorList newSrcGlobalFaceAreas(faceAreas0_);
+
+        srcToTgtMap1->distribute(newSrcGlobalFaceAreas);
+        forAll(nbrFaceAreas0, tgtFacei)
+        {
+            if (!tgtToSrcAddr1[tgtFacei].empty())
+            {
+                const label srcFacei = tgtToSrcAddr1[tgtFacei][0];
+                nbrFaceAreas0[tgtFacei] = -newSrcGlobalFaceAreas[srcFacei];
+            }
+        }
+    }
+    else
+    {
+        label nError = 0;
+        forAll(srcToTgtAddr0, srcFacei0)
+        {
+            const labelList& srcFaceTgtAddr0 = srcToTgtAddr0[srcFacei0];
+            const scalarList& srcFaceTgtWght0 = srcToTgtWght0[srcFacei0];
+            const pointList& srcFaceTgtCtr0 = srcCtr0[srcFacei0];
+            forAll(srcFaceTgtAddr0, addri)
+            {
+                const label srcFacei1 = srcFaceIDs_[srcFacei0][addri];
+
+                // Find which slot srcFacei0 appears in tgt->src addressing
+                const label tgtFacei0 = srcFaceTgtAddr0[addri];
+                const label tgtAddri0 =
+                    tgtToSrcAddr0[tgtFacei0].find(srcFacei0);
+
+                if (tgtAddri0 == -1)
+                {
+                    ++nError;
+                    continue;
+
+                    if (debug)
+                    {
+                        Pout<< "Unable to find source face " << srcFacei0
+                            << " in tgtToSrcAddr0[" << tgtFacei0 << "]: "
+                            << tgtToSrcAddr0[tgtFacei0]
+                            << endl;
+                    }
+                }
+
+                const label tgtFacei1 = tgtFaceIDs_[tgtFacei0][tgtAddri0];
+
+                faceAreas0_[srcFacei1] *= srcFaceTgtWght0[addri];
+                nbrFaceAreas0[tgtFacei1] = -faceAreas0_[srcFacei1];
+
+                point pt(srcFaceTgtCtr0[addri]);
+                faceCentres0_[srcFacei1] = pt;
+                reverseTransformPosition(pt, srcFacei0);
+                nbrFaceCentres0[tgtFacei1] = pt;
+
+                // SANITY CHECK
+                // Info<< "srcPt:" << srcFaceCentres[srcFacei1]
+                //     << " tgtPt:" << tgtFaceCentres[tgtFacei1] << endl;
+
+                srcToTgtAddr1[srcFacei1] = labelList(1, tgtFacei1);
+                tgtToSrcAddr1[tgtFacei1] = labelList(1, srcFacei1);
+            }
+        }
+
+        if (nError)
+        {
+            FatalErrorInFunction
+                << "Unable to find " << nError
+                << " source faces in tgtToSrcAddr0"
+                << abort(FatalError);
+        }
+    }
+
+    scalarListList newSrcToTgtWeights(srcToTgtAddr1.size());
+    forAll(srcToTgtAddr1, facei)
+    {
+        if (srcToTgtAddr1[facei].size())
+        {
+            newSrcToTgtWeights[facei] = scalarList(1, scalar(1));
+        }
+        else
+        {
+            // No connection - effect of face removed by setting area to a
+            // a small value
+            faceAreas0_[facei] *= tolerance_;
+        }
+    }
+
+    scalarListList newTgtToSrcWeights(tgtToSrcAddr1.size());
+    forAll(tgtToSrcAddr1, facei)
+    {
+        if (tgtToSrcAddr1[facei].size())
+        {
+            newTgtToSrcWeights[facei] = scalarList(1, scalar(1));
+        }
+        else
+        {
+            // No connection - effect of face removed by setting area to a
+            // a small value
+            nbrFaceAreas0[facei] *= tolerance_;
+        }
+    }
+
+    // Update the AMI addressing and weights to reflect the new 1-to-1
+    // correspondence
+    AMIPtr_->update
+    (
+        std::move(srcToTgtMap1),
+        std::move(tgtToSrcMap1),
+        std::move(srcToTgtAddr1),
+        std::move(newSrcToTgtWeights),
+        std::move(tgtToSrcAddr1),
+        std::move(newTgtToSrcWeights)
+    );
+
+    AMIPtr_->setAreas(mag(faceAreas0_), mag(nbrFaceAreas0));
+
+    if (debug)
+    {
+        Pout<< "cyclicAMIPolyPatch : " << name()
+            << " constructed AMI with " << nl
+            << "    " << "srcAddress:" << AMIPtr_().srcAddress().size()
+            << nl
+            << "    " << "tgAddress :" << AMIPtr_().tgtAddress().size()
+            << nl << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::cyclicAMIPolyPatch::changeTopology() const
+{
+    DebugInFunction << endl;
+
+    createAMIFaces_ = true;
+
+    return true;
+}
+
+
+bool Foam::cyclicAMIPolyPatch::setTopology(polyTopoChange& topoChange)
+{
+    DebugInFunction << endl;
+
+    if (createAMIFaces_ && owner())
+    {
+        // Calculate the AMI using the new points
+        // Note: mesh still has old points
+        resetAMI(topoChange.points(), AMIMethod_);
+
+        removeAMIFaces(topoChange);
+
+        addAMIFaces(topoChange);
+
+        return true;
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index 22b12ee488..c89a1afb64 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -271,6 +271,7 @@ AMICycPatches=$(AMI)/patches/cyclicAMI
 $(AMICycPatches)/cyclicAMILduInterfaceField/cyclicAMILduInterface.C
 $(AMICycPatches)/cyclicAMILduInterfaceField/cyclicAMILduInterfaceField.C
 $(AMICycPatches)/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
+$(AMICycPatches)/cyclicAMIPolyPatch/cyclicAMIPolyPatchTopologyChange.C
 $(AMICycPatches)/cyclicAMIPointPatch/cyclicAMIPointPatch.C
 $(AMICycPatches)/cyclicAMIPointPatchField/cyclicAMIPointPatchFields.C
 
diff --git a/src/meshTools/Make/options b/src/meshTools/Make/options
index 9ee5884e59..5b4f9ec841 100644
--- a/src/meshTools/Make/options
+++ b/src/meshTools/Make/options
@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/fileFormats/lnInclude \
-    -I$(LIB_SRC)/surfMesh/lnInclude
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 LIB_LIBS = \
     -lfileFormats \
diff --git a/src/meshTools/polyTopoChange/polyTopoChange.H b/src/meshTools/polyTopoChange/polyTopoChange.H
index d861a2d7da..613f6d398f 100644
--- a/src/meshTools/polyTopoChange/polyTopoChange.H
+++ b/src/meshTools/polyTopoChange/polyTopoChange.H
@@ -94,7 +94,7 @@ class IOobject;
 template<class T, class Container> class CompactListList;
 
 /*---------------------------------------------------------------------------*\
-                           Class polyTopoChange Declaration
+                       Class polyTopoChange Declaration
 \*---------------------------------------------------------------------------*/
 
 class polyTopoChange
@@ -121,7 +121,7 @@ class polyTopoChange
             DynamicList<label> pointMap_;
 
             //- For all original and added points contains new point label.
-            //  (used to map return value of addPoint to new mesh point)
+            //- (used to map return value of addPoint to new mesh point)
             DynamicList<label> reversePointMap_;
 
             //- Zone of point
@@ -324,10 +324,10 @@ class polyTopoChange
         );
 
         //- Remove all unused/removed points/faces/cells and update
-        //  face ordering (always), cell ordering (bandcompression,
-        //  orderCells=true),
-        //  point ordering (sorted into internal and boundary points,
-        //  orderPoints=true)
+        //- face ordering (always), cell ordering (bandcompression,
+        //- orderCells=true),
+        //- point ordering (sorted into internal and boundary points,
+        //- orderPoints=true)
         void compact
         (
             const bool orderCells,
@@ -433,7 +433,7 @@ public:
     // Constructors
 
         //- Construct without mesh. Either specify nPatches or use
-        //  setNumPatches before trying to make a mesh (makeMesh, changeMesh)
+        //- setNumPatches before trying to make a mesh (makeMesh, changeMesh)
         polyTopoChange(const label nPatches, const bool strict = true);
 
         //- Construct from mesh. Adds all points/face/cells from mesh.
@@ -471,15 +471,15 @@ public:
             }
 
             //- Is point removed?
-            //  Considered removed if point is GREAT.
+            //- Considered removed if point is GREAT.
             inline bool pointRemoved(const label pointi) const;
 
             //- Is face removed?
-            //  Considered removed if face is empty
+            //- Considered removed if face is empty
             inline bool faceRemoved(const label facei) const;
 
             //- Is cell removed?
-            //  Considered removed if the cellMap is -2
+            //- Considered removed if the cellMap is -2
             inline bool cellRemoved(const label celli) const;
 
 
@@ -489,7 +489,7 @@ public:
             void clear();
 
             //- Add all points/faces/cells of mesh. Additional offset for patch
-            //  or zone ids.
+            //- or zone ids.
             void addMesh
             (
                 const polyMesh& mesh,
@@ -500,7 +500,7 @@ public:
             );
 
             //- Explicitly pre-size the dynamic storage for expected mesh
-            //  size for if construct-without-mesh
+            //- size for if construct-without-mesh
             void setCapacity
             (
                 const label nPoints,
@@ -590,7 +590,7 @@ public:
             void removeCell(const label celli, const label mergeCelli);
 
             //- Explicitly set the number of patches if construct-without-mesh
-            //  used.
+            //- used.
             inline void setNumPatches(const label nPatches);
 
         // Other
@@ -628,7 +628,6 @@ public:
                 const bool orderCells = false,
                 const bool orderPoints = false
             );
-
 };