ENH: AMI improvements

- Added new faceAreaWeightAMI2D AMIMethod:
  - performs intersection using a new 2D triangle class;
  - candidate face matches set using an AABBTree method (vs advancing front for
    faceAreaWeightAMI).

  - Use by setting the AMIMethod entry when specifying the AMI in the
    constant/polyMesh/boundary file, e.g.

        AMI
        {
            type            cyclicACMI;
            AMIMethod       faceAreaWeightAMI2D; // new method
            Cbb             0.1; // optional coefficient
            nFaces          1000;
            startFace       100000;
            matchTolerance  0.0001;
            transform       noOrdering;
            neighbourPatch  AMI1;
            nonOverlapPatch AMI1_non_overlap;
        }

  - The optional Cbb coeffcient controls the size of the bounding box used when
    looking for candidate pairs; the value of 0.1 is the default and worked well
    for a large range of test cases. For badly matched AMI patches this may need
    to be increased.

- Deprecated the partialFaceAreaWeightAMI class - primarily used by ACMI:
  - functionality now offered by the AMI variants.

src/OpenFOAM/primitives/Vector2D/Vector2D.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.

src/OpenFOAM/primitives/Vector2D/Vector2DI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.

src/functionObjects/field/mapFields/mapFields.H

@@ -98,7 +98,6 @@ Usage
     \verbatim
       nearestFaceAMI
       faceAreaWeightAMI
-      partialFaceAreaWeightAMI
     \endverbatim
 
     The inherited entries are elaborated in:

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

@@ -1246,7 +1246,7 @@ void Foam::AMIInterpolation::write(Ostream& os) const
 
     if (reverseTarget_)
     {
-        os.writeEntry("flipNormals", reverseTarget_);
+        os.writeEntry("reverseTarget", reverseTarget_);
     }
 
     if (lowWeightCorrection_ > 0)

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -364,6 +364,9 @@ public:
             //- Access to the requireMatch flag
             inline bool setRequireMatch(const bool flag);
 
+            //- Return true if requireMatch and lowWeightCorrectionin active
+            inline bool mustMatchFaces() const;
+
             //- Access to the reverseTarget flag
             inline bool reverseTarget() const;
 

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -73,7 +73,7 @@ inline bool Foam::AMIInterpolation::distributed() const
 
 inline bool Foam::AMIInterpolation::requireMatch() const
 {
-    return requireMatch_ && lowWeightCorrection_ < 0;
+    return requireMatch_;
 }
 
 
@@ -84,6 +84,12 @@ inline bool Foam::AMIInterpolation::setRequireMatch(const bool flag)
 }
 
 
+inline bool Foam::AMIInterpolation::mustMatchFaces() const
+{
+    return requireMatch_ && !applyLowWeightCorrection();
+}
+
+
 inline bool Foam::AMIInterpolation::reverseTarget() const
 {
     return reverseTarget_;

src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.C

@@ -55,7 +55,7 @@ void Foam::advancingFrontAMI::checkPatches() const
     }
 
 
-    if (conformal())
+    if (requireMatch_)
     {
         const scalar maxBoundsError = 0.05;
 
@@ -345,6 +345,27 @@ void Foam::advancingFrontAMI::triangulatePatch
 }
 
 
+void Foam::advancingFrontAMI::nonConformalCorrection()
+{
+    if (!requireMatch_ && distributed())
+    {
+        scalarList newTgtMagSf(std::move(tgtMagSf_));
+
+        // Assign default sizes. Override selected values with calculated
+        // values. This is to support ACMI where some of the target faces
+        // are never used (so never get sent over and hence never assigned
+        // to)
+        tgtMagSf_ = tgtPatch0().magFaceAreas();
+
+        for (const labelList& smap : this->extendedTgtMapPtr_->subMap())
+        {
+            UIndirectList<scalar>(tgtMagSf_, smap) =
+                UIndirectList<scalar>(newTgtMagSf, smap);
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::advancingFrontAMI::advancingFrontAMI
@@ -421,7 +442,8 @@ bool Foam::advancingFrontAMI::calculate
 {
     if (AMIInterpolation::calculate(srcPatch, tgtPatch, surfPtr))
     {
-        // Create a representation of the target patch that covers the source patch
+        // Create a representation of the target patch that covers the source
+        // patch
         if (distributed())
         {
             createExtendedTgtPatch();
@@ -454,10 +476,4 @@ bool Foam::advancingFrontAMI::calculate
 }
 
 
-bool Foam::advancingFrontAMI::conformal() const
-{
-    return true;
-}
-
-
 // ************************************************************************* //

src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.H

@@ -199,6 +199,9 @@ protected:
                 List<scalar>& magSf
             ) const;
 
+            //- Correction for non-conformal interpolations, e.g. for ACMI
+            virtual void nonConformalCorrection();
+
 
 public:
 
@@ -249,10 +252,6 @@ public:
         //- Labels of faces that are not overlapped by any target faces
         //  Note: this should be empty for correct functioning
         inline const labelList& srcNonOverlap() const;
-
-        //- Flag to indicate that interpolation patches are conformal, i.e.
-        //- should fully cover each other
-        virtual bool conformal() const;
 };
 
 

src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.C

@@ -115,6 +115,9 @@ void Foam::faceAreaWeightAMI::calcAddressing
     // Reset starting seed
     label startSeedi = 0;
 
+    // Should all faces be matched?
+    const bool mustMatch = mustMatchFaces();
+
     bool continueWalk = true;
     DynamicList<label> nonOverlapFaces;
     do
@@ -154,7 +157,7 @@ void Foam::faceAreaWeightAMI::calcAddressing
             mapFlag,
             seedFaces,
             visitedFaces,
-            requireMatch_ && (lowWeightCorrection_ < 0)
+            mustMatch
             // pass in nonOverlapFaces for failed tree search?
         );
     } while (continueWalk);
@@ -675,7 +678,7 @@ bool Foam::faceAreaWeightAMI::calculate
         }
 
         // Check for badly covered faces
-        if (restartUncoveredSourceFace_) // && requireMatch_???
+        if (restartUncoveredSourceFace_) //  && mustMatchFaces())
         {
             restartUncoveredSourceFace
             (
@@ -775,7 +778,9 @@ bool Foam::faceAreaWeightAMI::calculate
     }
 
     // Convert the weights from areas to normalised values
-    normaliseWeights(conformal(), true);
+    normaliseWeights(requireMatch_, true);
+
+    nonConformalCorrection();
 
     upToDate_ = true;
 

src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.H

@@ -30,6 +30,8 @@ Class
 Description
     Face area weighted Arbitrary Mesh Interface (AMI) method
 
+    Searching is performed using an advancing front.
+
 SourceFiles
     faceAreaWeightAMI.C
 
@@ -61,7 +63,6 @@ private:
         const bool restartUncoveredSourceFace_;
 
 
-
 protected:
 
     // Protected Member Functions
@@ -69,15 +70,6 @@ protected:
         //- No copy assignment
         void operator=(const faceAreaWeightAMI&) = delete;
 
-        //- Initialise the geometry
-        void initGeom
-        (
-            const primitivePatch& srcPatch,
-            const primitivePatch& tgtPatch,
-            const globalIndex& globalTgtFaces,
-            labelList& extendedTgtFaceIDs
-        );
-
 
         // Marching front
 

src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.C

@@ -0,0 +1,526 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "faceAreaWeightAMI2D.H"
+#include "profiling.H"
+#include "OBJstream.H"
+#include "addToRunTimeSelectionTable.H"
+#include "triangle2D.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(faceAreaWeightAMI2D, 0);
+    addToRunTimeSelectionTable(AMIInterpolation, faceAreaWeightAMI2D, dict);
+    addToRunTimeSelectionTable
+    (
+        AMIInterpolation,
+        faceAreaWeightAMI2D,
+        component
+    );
+}
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+void Foam::faceAreaWeightAMI2D::writeNoMatch
+(
+    const label srcFacei,
+    const labelList& tgtFaceCandidates,
+    const boundBox& srcFaceBb
+) const
+{
+    Info<< "NO MATCH for source face " << srcFacei << endl;
+    {
+        const auto& src = this->srcPatch();
+        const auto& tgt = this->tgtPatch(); // might be the extended patch!
+
+        OFstream os("no_match_" + Foam::name(srcFacei) + ".obj");
+
+        const pointField& srcPoints = src.points();
+        const pointField& tgtPoints = tgt.points();
+
+        label np = 0;
+
+        // Write source face
+        const face& srcF = src[srcFacei];
+        string faceStr = "f";
+        for (const label pointi : srcF)
+        {
+            const point& p = srcPoints[pointi];
+            os  << "v " << p.x() << " " << p.y() << " " << p.z() << nl;
+            ++np;
+            faceStr += " " + Foam::name(np);
+        }
+        os << faceStr.c_str() << " " << (np - srcF.size() + 1) << nl;
+
+        // Write target faces as lines
+        for (const label tgtFacei : tgtFaceCandidates)
+        {
+            const face& tgtF = tgt[tgtFacei];
+            forAll(tgtF, pointi)
+            {
+                const point& p = tgtPoints[tgtF[pointi]];
+                os  << "v " << p.x() << " " << p.y() << " " << p.z() << nl;
+                ++np;
+                if (pointi)
+                {
+                    os << "l " << np-1 << " " << np << nl;
+                }
+            }
+            os << "l " << (np - tgtF.size() + 1) << " " << np << nl;
+        }
+    }
+
+    {
+        OFstream os("no_match_" + Foam::name(srcFacei) + "_bb.obj");
+
+        const pointField points(srcFaceBb.points());
+        for (const point& p : points)
+        {
+            os  << "v " << p.x() << " " << p.y() << " " << p.z() << endl;
+        }
+        os  << "l 1 2" << nl;
+        os  << "l 2 4" << nl;
+        os  << "l 4 3" << nl;
+        os  << "l 3 1" << nl;
+        os  << "l 5 6" << nl;
+        os  << "l 6 8" << nl;
+        os  << "l 8 7" << nl;
+        os  << "l 7 5" << nl;
+        os  << "l 5 1" << nl;
+        os  << "l 6 2" << nl;
+        os  << "l 8 4" << nl;
+        os  << "l 7 3" << nl;
+    }
+}
+
+
+void Foam::faceAreaWeightAMI2D::storeInterArea
+(
+    const label srcFacei,
+    const label tgtFacei,
+    DynamicList<label>& srcAddr,
+    DynamicList<scalar>& srcWght,
+    DynamicList<vector>& srcCtr,
+    DynamicList<label>& tgtAddr,
+    DynamicList<scalar>& tgtWght
+) const
+{
+    addProfiling(ami, "faceAreaWeightAMI2D::calcInterArea");
+
+    // Quick reject if either face has zero area
+    if
+    (
+        (srcMagSf_[srcFacei] < ROOTVSMALL)
+     || (tgtMagSf_[tgtFacei] < ROOTVSMALL)
+    )
+    {
+        return;
+    }
+
+    const auto& srcPatch = this->srcPatch();
+    const auto& tgtPatch = this->tgtPatch();
+
+    const pointField& srcPoints = srcPatch.points();
+    const pointField& tgtPoints = tgtPatch.points();
+
+    const auto& srcTris = srcTris_[srcFacei];
+    const auto& tgtTris = tgtTris_[tgtFacei];
+
+    const auto& srcFaceNormals = srcPatch.faceNormals();
+
+    //triangle2D::debug = 1;
+
+    scalar area = 0;
+    vector centroid = Zero;
+
+    for (const auto& tris : srcTris)
+    {
+        const vector& origin = srcPoints[tris[0]];
+        const vector p10(srcPoints[tris[1]] - origin);
+        const vector p20(srcPoints[tris[2]] - origin);
+        const vector axis1(p10/(mag(p10) + ROOTVSMALL));
+        const vector axis2(srcFaceNormals[srcFacei]^axis1);
+
+        triangle2D s
+        (
+            vector2D(0, 0),
+            vector2D(axis1&p10, axis2&p10),
+            vector2D(axis1&p20, axis2&p20)
+        );
+
+        for (const auto& trit : tgtTris)
+        {
+            // Triangle t has opposite orientation wrt triangle s
+            triangle2D t
+            (
+                tgtPoints[trit[0]] - origin,
+                tgtPoints[trit[2]] - origin,
+                tgtPoints[trit[1]] - origin,
+                axis1,
+                axis2
+            );
+
+            scalar da = 0;
+            vector2D c(Zero);
+
+            if (t.snapClosePoints(s) == 3)
+            {
+                c = s.centre();
+                da = mag(s.area());
+            }
+            else
+            {
+                s.interArea(t, c, da);
+            }
+
+            area += da;
+            centroid += da*(origin + c.x()*axis1 + c.y()*axis2);
+        }
+    }
+
+    //triangle2D::debug = 0;
+
+    if (area/srcMagSf_[srcFacei] > triangle2D::relTol)
+    {
+        centroid /= area + ROOTVSMALL;
+
+        srcAddr.append(tgtFacei);
+        srcWght.append(area);
+        srcCtr.append(centroid);
+
+        tgtAddr.append(srcFacei);
+        tgtWght.append(area);
+    }
+}
+
+
+Foam::labelList Foam::faceAreaWeightAMI2D::overlappingTgtFaces
+(
+    const AABBTree<face>& tree,
+    const List<boundBox>& tgtFaceBbs,
+    const boundBox& srcFaceBb
+) const
+{
+    labelHashSet faceIds(16);
+
+    const auto& treeBb = tree.boundBoxes();
+    const auto& treeAddr = tree.addressing();
+
+    forAll(treeBb, boxi)
+    {
+        const auto& tbb = treeBb[boxi];
+
+        if (srcFaceBb.overlaps(tbb))
+        {
+            const auto& boxAddr = treeAddr[boxi];
+
+            for (const auto& tgtFacei : boxAddr)
+            {
+                if (srcFaceBb.overlaps(tgtFaceBbs[tgtFacei]))
+                {
+                    faceIds.insert(tgtFacei);
+                }
+            }
+        }
+    }
+
+    return faceIds.toc();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::faceAreaWeightAMI2D::faceAreaWeightAMI2D
+(
+    const dictionary& dict,
+    const bool reverseTarget
+)
+:
+    advancingFrontAMI(dict, reverseTarget),
+    Cbb_(dict.getCheckOrDefault<scalar>("Cbb", 0.1, scalarMinMax::ge(SMALL)))
+{}
+
+
+Foam::faceAreaWeightAMI2D::faceAreaWeightAMI2D
+(
+    const bool requireMatch,
+    const bool reverseTarget,
+    const scalar lowWeightCorrection,
+    const faceAreaIntersect::triangulationMode triMode,
+    const bool restartUncoveredSourceFace
+)
+:
+    advancingFrontAMI
+    (
+        requireMatch,
+        reverseTarget,
+        lowWeightCorrection,
+        triMode
+    ),
+    Cbb_(0.1)
+{}
+
+
+Foam::faceAreaWeightAMI2D::faceAreaWeightAMI2D
+(
+    const faceAreaWeightAMI2D& ami
+)
+:
+    advancingFrontAMI(ami),
+    Cbb_(ami.Cbb_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::faceAreaWeightAMI2D::calculate
+(
+    const primitivePatch& srcPatch,
+    const primitivePatch& tgtPatch,
+    const autoPtr<searchableSurface>& surfPtr
+)
+{
+    if (upToDate_)
+    {
+        return false;
+    }
+
+    addProfiling(ami, "faceAreaWeightAMI2D::calculate");
+
+    advancingFrontAMI::calculate(srcPatch, tgtPatch, surfPtr);
+
+    const auto& src = this->srcPatch();
+    const auto& tgt = this->tgtPatch(); // might be the extended patch!
+
+    bool validSize = true;
+
+    // Check that patch sizes are valid
+    if (!src.size())
+    {
+        validSize = false;
+    }
+    else if (!tgt.size())
+    {
+        WarningInFunction
+            << src.size() << " source faces but no target faces" << endl;
+
+        validSize = false;
+    }
+
+    srcCentroids_.setSize(srcAddress_.size());
+
+    // Temporary storage for addressing and weights
+    List<DynamicList<label>> tgtAddr(tgt.size());
+    List<DynamicList<scalar>> tgtWght(tgt.size());
+
+    // Find an approximate face length scale
+    const scalar lf(Cbb_*Foam::sqrt(gAverage(srcMagSf_)));
+
+    // Expansion to apply to source face bounding box
+    const vector d(lf*vector::one);
+
+    if (validSize)
+    {
+        // Create the tgt tree
+        const bool equalBinSize = true;
+        const label maxLevel = 10;
+        const label minBinSize = 4;
+        AABBTree<face> tree
+        (
+            tgt,
+            tgt.points(),
+            equalBinSize,
+            maxLevel,
+            minBinSize
+        );
+
+        const auto& tgtPoints = tgt.points();
+        List<boundBox> tgtFaceBbs(tgt.size());
+        forAll(tgt, facei)
+        {
+            tgtFaceBbs[facei] = boundBox(tgtPoints, tgt[facei], false);
+        }
+
+        DynamicList<label> nonOverlapFaces;
+
+        const auto& srcPoints = src.points();
+
+        forAll(src, srcFacei)
+        {
+            const face& srcFace = src[srcFacei];
+
+            treeBoundBox srcFaceBb(srcPoints, srcFace);
+            srcFaceBb.min() -= d;
+            srcFaceBb.max() += d;
+
+            const labelList tgtFaces
+            (
+                overlappingTgtFaces(tree, tgtFaceBbs, srcFaceBb)
+            );
+
+            DynamicList<label> srcAddr(tgtFaces.size());
+            DynamicList<scalar> srcWght(tgtFaces.size());
+            DynamicList<point> srcCtr(tgtFaces.size());
+
+            for (const label tgtFacei : tgtFaces)
+            {
+                storeInterArea
+                (
+                    srcFacei,
+                    tgtFacei,
+                    srcAddr,
+                    srcWght,
+                    srcCtr,
+                    tgtAddr[tgtFacei],
+                    tgtWght[tgtFacei]
+                );
+            }
+
+            if (mustMatchFaces() && srcAddr.empty())
+            {
+                if (debug) writeNoMatch(srcFacei, tgtFaces, srcFaceBb);
+
+                // FatalErrorInFunction
+                //    << "Unable to find match for source face " << srcFace
+                //    << exit(FatalError);
+            }
+
+            srcAddress_[srcFacei].transfer(srcAddr);
+            srcWeights_[srcFacei].transfer(srcWght);
+            srcCentroids_[srcFacei].transfer(srcCtr);
+        }
+
+        srcNonOverlap_.transfer(nonOverlapFaces);
+
+        if (debug && !srcNonOverlap_.empty())
+        {
+            Pout<< "    AMI: " << srcNonOverlap_.size()
+                << " non-overlap faces identified"
+                << endl;
+        }
+    }
+
+    // Transfer data to persistent storage
+
+    forAll(tgtAddr, i)
+    {
+        tgtAddress_[i].transfer(tgtAddr[i]);
+        tgtWeights_[i].transfer(tgtWght[i]);
+    }
+
+    if (distributed())
+    {
+        const primitivePatch& srcPatch0 = this->srcPatch0();
+        const primitivePatch& tgtPatch0 = this->tgtPatch0();
+
+        // Create global indexing for each original patch
+        globalIndex globalSrcFaces(srcPatch0.size());
+        globalIndex globalTgtFaces(tgtPatch0.size());
+
+        for (labelList& addressing : srcAddress_)
+        {
+            for (label& addr : addressing)
+            {
+                addr = extendedTgtFaceIDs_[addr];
+            }
+        }
+
+        for (labelList& addressing : tgtAddress_)
+        {
+            globalSrcFaces.inplaceToGlobal(addressing);
+        }
+
+        // Send data back to originating procs. Note that contributions
+        // from different processors get added (ListOps::appendEqOp)
+
+        mapDistributeBase::distribute
+        (
+            Pstream::commsTypes::nonBlocking,
+            List<labelPair>(),
+            tgtPatch0.size(),
+            extendedTgtMapPtr_->constructMap(),
+            false,                      // has flip
+            extendedTgtMapPtr_->subMap(),
+            false,                      // has flip
+            tgtAddress_,
+            labelList(),
+            ListOps::appendEqOp<label>(),
+            flipOp()
+        );
+
+        mapDistributeBase::distribute
+        (
+            Pstream::commsTypes::nonBlocking,
+            List<labelPair>(),
+            tgtPatch0.size(),
+            extendedTgtMapPtr_->constructMap(),
+            false,
+            extendedTgtMapPtr_->subMap(),
+            false,
+            tgtWeights_,
+            scalarList(),
+            ListOps::appendEqOp<scalar>(),
+            flipOp()
+        );
+
+        // Note: using patch face areas calculated by the AMI method
+        extendedTgtMapPtr_->reverseDistribute(tgtPatch0.size(), tgtMagSf_);
+
+        // Cache maps and reset addresses
+        List<Map<label>> cMapSrc;
+        srcMapPtr_.reset
+        (
+            new mapDistribute(globalSrcFaces, tgtAddress_, cMapSrc)
+        );
+
+        List<Map<label>> cMapTgt;
+        tgtMapPtr_.reset
+        (
+            new mapDistribute(globalTgtFaces, srcAddress_, cMapTgt)
+        );
+    }
+
+    // Convert the weights from areas to normalised values
+    normaliseWeights(requireMatch_, true);
+
+    nonConformalCorrection();
+
+    upToDate_ = true;
+
+    return upToDate_;
+}
+
+
+void Foam::faceAreaWeightAMI2D::write(Ostream& os) const
+{
+    advancingFrontAMI::write(os);
+}
+
+
+// ************************************************************************* //

src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.H

@@ -5,8 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
-    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,20 +24,22 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::partialFaceAreaWeightAMI
+    Foam::faceAreaWeightAMI2D
 
 Description
-    Partial face area weighted Arbitrary Mesh Interface (AMI) method
+    Face area weighted Arbitrary Mesh Interface (AMI) method that performs
+    the intersection of src and tgt face area in 2D.
 
 SourceFiles
-    partialFaceAreaWeightAMI.C
+    faceAreaWeightAMI2D.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef partialFaceAreaWeightAMI_H
+#ifndef faceAreaWeightAMI2D_H
-#define partialFaceAreaWeightAMI_H
+#define faceAreaWeightAMI2D_H
 
-#include "faceAreaWeightAMI.H"
+#include "advancingFrontAMI.H"
+#include "AABBTree.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,58 +47,80 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                  Class partialFaceAreaWeightAMI Declaration
+                     Class faceAreaWeightAMI2D Declaration
 \*---------------------------------------------------------------------------*/
 
-class partialFaceAreaWeightAMI
+class faceAreaWeightAMI2D
 :
-    public faceAreaWeightAMI
+    public advancingFrontAMI
 {
-
-private:
-
-    // Private Member Functions
-
-        //- No copy assignment
-        void operator=(const partialFaceAreaWeightAMI&) = delete;
-
-
 protected:
 
+    // Protected Data
+
+        //- Face bounding box factor
+        scalar Cbb_;
+
+
     // Protected Member Functions
 
-        //- Set the source and target seed faces
+        //- No copy assignment
-        virtual bool setNextFaces
+        void operator=(const faceAreaWeightAMI2D&) = delete;
+
+        //- Helper function to write non-matched source faces to the set
+        //- of candidate faces
+        void writeNoMatch
         (
-            label& startSeedi,
+            const label srcFacei,
-            label& srcFacei,
+            const labelList& tgtFaceCandidates,
-            label& tgtFacei,
+            const boundBox& srcFaceBb
-            const bitSet& mapFlag,
-            labelList& seedFaces,
-            const DynamicList<label>& visitedFaces,
-            const bool errorOnNotFound = true
         ) const;
 
 
+        // Evaluation
+
+            //- Calculate and store the area of intersection between source and
+            //- target faces
+            void storeInterArea
+            (
+                const label srcFacei,
+                const label tgtFacei,
+                DynamicList<label>& srcAddr,
+                DynamicList<scalar>& srcWght,
+                DynamicList<vector>& srcCtr,
+                DynamicList<label>& tgtAddr,
+                DynamicList<scalar>& tgtWght
+            ) const;
+
+
+            //- Return the set of tgt face IDs that overlap the src face bb
+            labelList overlappingTgtFaces
+            (
+                const AABBTree<face>& tree,
+                const List<boundBox>& tgtFaceBbs,
+                const boundBox& srcFaceBb
+            ) const;
+
+
 public:
 
     //- Runtime type information
-    TypeName("partialFaceAreaWeightAMI");
+    TypeName("faceAreaWeightAMI2D");
 
 
     // Constructors
 
         //- Construct from dictionary
-        partialFaceAreaWeightAMI
+        faceAreaWeightAMI2D
         (
             const dictionary& dict,
             const bool reverseTarget = false
         );
 
         //- Construct from components
-        partialFaceAreaWeightAMI
+        faceAreaWeightAMI2D
         (
-            const bool requireMatch = false,
+            const bool requireMatch,
             const bool reverseTarget = false,
             const scalar lowWeightCorrection = -1,
             const faceAreaIntersect::triangulationMode triMode =
@@ -106,25 +129,21 @@ public:
         );
 
         //- Construct as copy
-        partialFaceAreaWeightAMI(const partialFaceAreaWeightAMI& ami);
+        faceAreaWeightAMI2D(const faceAreaWeightAMI2D& ami);
 
         //- Construct and return a clone
         virtual autoPtr<AMIInterpolation> clone() const
         {
-            return
+            return autoPtr<AMIInterpolation>(new faceAreaWeightAMI2D(*this));
-                autoPtr<AMIInterpolation>(new partialFaceAreaWeightAMI(*this));
         }
 
 
     //- Destructor
-    virtual ~partialFaceAreaWeightAMI() = default;
+    virtual ~faceAreaWeightAMI2D() = default;
 
 
     // Member Functions
 
-        //- Flag to indicate that interpolation patches are conformal
-        virtual bool conformal() const;
-
         //- Update addressing, weights and (optional) centroids
         virtual bool calculate
         (
@@ -132,6 +151,9 @@ public:
             const primitivePatch& tgtPatch,
             const autoPtr<searchableSurface>& surfPtr = nullptr
         );
+
+        //- Write
+        virtual void write(Ostream& os) const;
 };
 
 

src/meshTools/AMIInterpolation/AMIInterpolation/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.C

@@ -1,159 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2013-2016 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/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "partialFaceAreaWeightAMI.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(partialFaceAreaWeightAMI, 0);
-    addToRunTimeSelectionTable
-    (
-        AMIInterpolation,
-        partialFaceAreaWeightAMI,
-        dict
-    );
-    addToRunTimeSelectionTable
-    (
-        AMIInterpolation,
-        partialFaceAreaWeightAMI,
-        component
-    );
-}
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-bool Foam::partialFaceAreaWeightAMI::setNextFaces
-(
-    label& startSeedi,
-    label& srcFacei,
-    label& tgtFacei,
-    const bitSet& mapFlag,
-    labelList& seedFaces,
-    const DynamicList<label>& visitedFaces,
-    const bool errorOnNotFound
-) const
-{
-    return faceAreaWeightAMI::setNextFaces
-    (
-        startSeedi,
-        srcFacei,
-        tgtFacei,
-        mapFlag,
-        seedFaces,
-        visitedFaces,
-        false // no error on not found
-    );
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::partialFaceAreaWeightAMI::partialFaceAreaWeightAMI
-(
-    const dictionary& dict,
-    const bool reverseTarget
-)
-:
-    faceAreaWeightAMI(dict, reverseTarget)
-{}
-
-
-Foam::partialFaceAreaWeightAMI::partialFaceAreaWeightAMI
-(
-    const bool requireMatch,
-    const bool reverseTarget,
-    const scalar lowWeightCorrection,
-    const faceAreaIntersect::triangulationMode triMode,
-    const bool restartUncoveredSourceFace
-)
-:
-    faceAreaWeightAMI
-    (
-        requireMatch,
-        reverseTarget,
-        lowWeightCorrection,
-        triMode,
-        restartUncoveredSourceFace
-    )
-{}
-
-
-Foam::partialFaceAreaWeightAMI::partialFaceAreaWeightAMI
-(
-    const partialFaceAreaWeightAMI& ami
-)
-:
-    faceAreaWeightAMI(ami)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::partialFaceAreaWeightAMI::conformal() const
-{
-    return false;
-}
-
-
-bool Foam::partialFaceAreaWeightAMI::calculate
-(
-    const primitivePatch& srcPatch,
-    const primitivePatch& tgtPatch,
-    const autoPtr<searchableSurface>& surfPtr
-)
-{
-    if (faceAreaWeightAMI::calculate(srcPatch, tgtPatch, surfPtr))
-    {
-        if (distributed())
-        {
-            scalarList newTgtMagSf(std::move(tgtMagSf_));
-
-            // Assign default sizes. Override selected values with calculated
-            // values. This is to support ACMI where some of the target faces
-            // are never used (so never get sent over and hence never assigned
-            // to)
-            tgtMagSf_ = tgtPatch0().magFaceAreas();
-
-            for (const labelList& smap : this->extendedTgtMapPtr_->subMap())
-            {
-                UIndirectList<scalar>(tgtMagSf_, smap) =
-                    UIndirectList<scalar>(newTgtMagSf, smap);
-            }
-        }
-
-        return true;
-    }
-
-    return false;
-}
-
-
-// ************************************************************************* //

src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -63,7 +63,6 @@ SourceFiles
 #include "cyclicAMIPolyPatch.H"
 #include "AMIPatchToPatchInterpolation.H"
 #include "polyBoundaryMesh.H"
-#include "partialFaceAreaWeightAMI.H"
 #include "PatchFunction1.H"
 #include "uniformDimensionedFields.H"
 #include "vectorList.H"
@@ -197,7 +196,7 @@ public:
             const polyBoundaryMesh& bm,
             const word& patchType,
             const transformType transform = UNKNOWN,
-            const word& defaultAMIMethod = partialFaceAreaWeightAMI::typeName
+            const word& defaultAMIMethod = faceAreaWeightAMI::typeName
         );
 
         //- Construct from dictionary
@@ -208,7 +207,7 @@ public:
             const label index,
             const polyBoundaryMesh& bm,
             const word& patchType,
-            const word& defaultAMIMethod = partialFaceAreaWeightAMI::typeName
+            const word& defaultAMIMethod = faceAreaWeightAMI::typeName
         );
 
         //- Construct as copy, resetting the boundary mesh

src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,7 +26,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "cyclicPeriodicAMIPolyPatch.H"
-#include "partialFaceAreaWeightAMI.H"
 #include "addToRunTimeSelectionTable.H"
 
 // For debugging
@@ -560,14 +559,16 @@ Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
         bm,
         patchType,
         transform,
-        partialFaceAreaWeightAMI::typeName
+        faceAreaWeightAMI::typeName
     ),
     periodicPatchName_(word::null),
     periodicPatchID_(-1),
     nTransforms_(0),
     nSectors_(0),
     maxIter_(36)
-{}
+{
+    AMIPtr_->setRequireMatch(false);
+}
 
 
 Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
@@ -586,14 +587,16 @@ Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
         index,
         bm,
         patchType,
-        partialFaceAreaWeightAMI::typeName
+        faceAreaWeightAMI::typeName
     ),
     periodicPatchName_(dict.lookup("periodicPatch")),
     periodicPatchID_(-1),
     nTransforms_(dict.getOrDefault<label>("nTransforms", 0)),
     nSectors_(dict.getOrDefault<label>("nSectors", 0)),
     maxIter_(dict.getOrDefault<label>("maxIter", 36))
-{}
+{
+    AMIPtr_->setRequireMatch(false);
+}
 
 
 Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
@@ -608,7 +611,9 @@ Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
     nTransforms_(pp.nTransforms_),
     nSectors_(pp.nSectors_),
     maxIter_(pp.maxIter_)
-{}
+{
+    AMIPtr_->setRequireMatch(false);
+}
 
 
 Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
@@ -627,7 +632,9 @@ Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
     nTransforms_(pp.nTransforms_),
     nSectors_(pp.nSectors_),
     maxIter_(pp.maxIter_)
-{}
+{
+    AMIPtr_->setRequireMatch(false);
+}
 
 
 Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
@@ -645,7 +652,9 @@ Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
     nTransforms_(pp.nTransforms_),
     nSectors_(pp.nSectors_),
     maxIter_(pp.maxIter_)
-{}
+{
+    AMIPtr_->setRequireMatch(false);
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

src/meshTools/AMIInterpolation/triangle2D/triangle2D.C

@@ -0,0 +1,303 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2021 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 "triangle2D.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+int Foam::triangle2D::debug = 0;
+
+Foam::scalar Foam::triangle2D::relTol = 1e-8;
+
+Foam::scalar Foam::triangle2D::absTol = 1e-10;
+
+Foam::FixedList<Foam::vector2D, 7> Foam::triangle2D::work_
+(
+    vector2D::zero
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::triangle2D::triangle2D
+(
+    const vector2D& a,
+    const vector2D& b,
+    const vector2D& c,
+    const bool orient
+)
+:
+    FixedList<vector2D, 3>({a, b, c}),
+    area_(area(a, b, c))
+{
+    if (orient && area_ < 0)
+    {
+        // Inverted triangle
+        triangle2D& tri = *this;
+        vector2D tmp = tri[0];
+        tri[0] = tri[2];
+        tri[2] = tmp;
+
+        area_ = mag(area_);
+    }
+}
+
+
+Foam::triangle2D::triangle2D
+(
+    const vector& a3d,
+    const vector& b3d,
+    const vector& c3d,
+    const vector& axis1,
+    const vector& axis2,
+    const bool orient
+)
+:
+    triangle2D
+    (
+        vector2D(axis1 & a3d, axis2 & a3d),
+        vector2D(axis1 & b3d, axis2 & b3d),
+        vector2D(axis1 & c3d, axis2 & c3d),
+        orient
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::label Foam::triangle2D::snapClosePoints(const triangle2D& triB)
+{
+    label nSnapped = 0;
+    triangle2D& triA = *this;
+
+    FixedList<bool, 3> match(true);
+
+    for (auto& a : triA)
+    {
+        forAll(triB, tb)
+        {
+            if (match[tb] && a.isClose(triB[tb]))
+            {
+                a = triB[tb];
+                match[tb] = false;
+                ++nSnapped;
+                break;
+            }
+        }
+    }
+
+    return nSnapped;
+}
+
+
+void Foam::triangle2D::interArea
+(
+    const triangle2D& triB,
+    vector2D& centre,
+    scalar& area
+) const
+{
+    const triangle2D& triA = *this;
+
+    // Potential short-cut if the triangles are the same-ish.  Happens rarely
+    // for moving mesh cases.
+    // if (nClosePoints(triA, triB) == 3)
+    // {
+    //     centre = triA.centre();
+    //     area = triA.area();
+    //     return;
+    // }
+
+    if (debug)
+    {
+        static int nInter = 0;
+        OFstream os("intersection-tris-"+Foam::name(nInter++)+".obj");
+        writeOBJ(os, triA, 0);
+        writeOBJ(os, triB, 3);
+        Info<< "written " << os.name() << endl;
+    }
+
+
+    // Use work_ to store intersections
+
+    // Current number of intersections
+    int nPoint = 0;
+
+    static FixedList<vector2D, 7> work2;
+
+    // Clipped triangle starts as triA
+    work2[0] = triA[0];
+    work2[1] = triA[1];
+    work2[2] = triA[2];
+
+    int nPoint2 = 3;
+
+
+    vector2D intersection(Zero);
+
+    // Cut triA using triB's edges as clipping planes
+    for (int i0 = 0; i0 <= 2; ++i0)
+    {
+        if (debug)
+        {
+            Info<< "\nstarting points:";
+            for (int i = 0; i < nPoint2; ++i)
+            {
+                Info<< work2[i];
+            }
+            Info<< endl;
+        }
+
+        // Clipping plane points
+        const label i1 = (i0 + 1) % 3;
+        const vector2D& c0 = triB[i0];
+        const vector2D& c1 = triB[i1];
+
+        nPoint = 0;
+
+        // Test against all intersection poly points
+        for (int j = 0; j < nPoint2; ++j)
+        {
+            const vector2D& p0 = work2[j];
+            const vector2D& p1 = work2[(j+1) % nPoint2];
+
+            if (triangle2D(c0, c1, p0).order() == 1)
+            {
+                if (triangle2D(c0, c1, p1).order() == 1)
+                {
+                    work_[nPoint++] = p1;
+                }
+                else
+                {
+                    if (lineIntersectionPoint(p0, p1, c0, c1, intersection))
+                    {
+                        work_[nPoint++] = intersection;
+                    }
+                }
+            }
+            else
+            {
+                if (triangle2D(c0, c1, p1).order() == 1)
+                {
+                    if (lineIntersectionPoint(p0, p1, c0, c1, intersection))
+                    {
+                        work_[nPoint++] = intersection;
+                    }
+                    work_[nPoint++] = p1;
+                }
+            }
+        }
+
+        work2 = work_;
+        nPoint2 = nPoint;
+    }
+
+    if (debug)
+    {
+        static int n = 0;
+        OFstream os("intersection-poly-"+Foam::name(n++)+".obj");
+        for (int i = 0; i < nPoint; ++i)
+        {
+            os << "v " << work_[i].x() << " " << work_[i].y() << " 0" << nl;
+        }
+        os  << "l";
+        for (int i = 0; i < nPoint; ++i)
+        {
+            os << " " << (i + 1);
+        }
+        os  << " 1" << endl;
+        Info<< "written " << os.name() << endl;
+
+        Info<< "Intersection points:" << endl;
+        for (int i = 0; i < nPoint; ++i)
+        {
+            Info<< "    " << work_[i] << endl;
+        }
+    }
+
+    // Calculate the area of the clipped triangle
+    scalar twoArea = 0;
+    centre = vector2D::zero;
+    if (nPoint)
+    {
+        for (int i = 0; i < nPoint - 1; ++i)
+        {
+            twoArea += work_[i].x()*work_[i+1].y();
+            twoArea -= work_[i].y()*work_[i+1].x();
+            centre += work_[i];
+        }
+        twoArea += work_[nPoint-1].x()*work_[0].y();
+        twoArea -= work_[nPoint-1].y()*work_[0].x();
+
+        centre += work_[nPoint - 1];
+        centre /= scalar(nPoint);
+    }
+
+    area = 0.5*twoArea;
+}
+
+
+Foam::scalar Foam::triangle2D::interArea(const triangle2D& triB) const
+{
+    vector2D dummyCentre(Zero);
+    scalar area;
+
+    interArea(triB, dummyCentre, area);
+
+    return area;
+}
+
+
+bool Foam::triangle2D::overlaps(const triangle2D& triB) const
+{
+    const triangle2D& triA = *this;
+
+    // True if any of triB's edges intersect a triA edge
+    for (int i = 0; i < 3; ++i)
+    {
+        int i1 = (i + 1) % 3;
+
+        for (int j = 0; j < 3; ++j)
+        {
+            int j1 = (j + 1) % 3;
+
+            if (lineIntersects(triA[i], triA[i1], triB[j], triB[j1]))
+            {
+                return true;
+            }
+        }
+    }
+
+    return
+        (nClosePoints(triA, triB) == 3) // same tri
+     || triA.contains(triB)             // triA contains triB
+     || triB.contains(triA);            // triB contains triA
+}
+
+
+// ************************************************************************* //

src/meshTools/AMIInterpolation/triangle2D/triangle2D.H

@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2021 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::triangle2D
+
+Description
+    2-D triangle and queries
+
+SourceFiles
+    triangle2DI.H
+    triangle2D.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef triangle2D_H
+#define triangle2D_H
+
+#include "vector.H"
+#include "vector2D.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class triangle2D Declaration
+\*---------------------------------------------------------------------------*/
+
+class triangle2D
+:
+    public FixedList<vector2D, 3>
+{
+    // Private Data
+
+        //- Area
+        scalar area_;
+
+        //- Helper for intersections
+        static FixedList<vector2D, 7> work_;
+
+
+public:
+
+    static int debug;
+
+    //- Relative tolerance
+    static scalar relTol;
+
+    //- Absolute tolerance
+    static scalar absTol;
+
+    //- Construct from 3 2-D points
+    triangle2D
+    (
+        const vector2D& a,
+        const vector2D& b,
+        const vector2D& c,
+        const bool orient = false
+    );
+
+    //- Construct from 3 3-D points and axes
+    triangle2D
+    (
+        const vector& a3d,
+        const vector& b3d,
+        const vector& c3d,
+        const vector& axis1,
+        const vector& axis2,
+        const bool orient = false
+    );
+
+
+    // Member Functions
+
+        //- Returns:
+        //-   1 if points are ordered in anti-clockwise direction
+        //-  -1 if points are ordered in clockwise direction
+        //-   0 if the triangle  has collapsed to a line
+        inline label order() const;
+
+        //- Write the triangle in OBJ format
+        inline static void writeOBJ
+        (
+            Ostream& os,
+            const triangle2D& tri,
+            label offset
+        );
+
+        //- Return the number of similar points between two triangles
+        inline static label nClosePoints
+        (
+            const triangle2D& triA,
+            const triangle2D& triB
+        );
+
+        //- Return the signed area
+        inline static scalar area
+        (
+            const vector2D& a,
+            const vector2D& b,
+            const vector2D& c
+        );
+
+        //- Set the intersection between a line and segment
+        //- Return true if lines intersect
+        inline static bool lineSegmentIntersectionPoint
+        (
+            const vector2D& lp1,
+            const vector2D& lp2,
+            const vector2D& sp1,
+            const vector2D& sp2,
+            vector2D& intersection
+        );
+
+        //- Set the intersection between two lines
+        //- Return true if lines intersect
+        inline static bool lineIntersectionPoint
+        (
+            const vector2D& a,
+            const vector2D& b,
+            const vector2D& c,
+            const vector2D& d,
+            vector2D& intersection
+        );
+
+        //- Return true if lines ab and cd intersect
+        inline static bool lineIntersects
+        (
+            const vector2D& a,
+            const vector2D& b,
+            const vector2D& c,
+            const vector2D& d
+        );
+
+        //- Snap [this] triangle's points to those of triB if they are within
+        //- absTol
+        //  Returns the number of snapped points
+        label snapClosePoints(const triangle2D& triB);
+
+        //- Return the intersection centre and area
+        void interArea
+        (
+            const triangle2D& triB,
+            vector2D& centre,
+            scalar& area
+        ) const;
+
+        //- Return the intersection area
+        scalar interArea(const triangle2D& triB) const;
+
+        //- Return true if triB overlaps
+        bool overlaps(const triangle2D& triB) const;
+
+        //- Return the triangle area
+        inline scalar area() const noexcept;
+
+        //- Return the triangle centre
+        inline vector2D centre() const;
+
+        //- Return true if tri is within this triangle
+        inline bool contains(const triangle2D& tri) const;
+
+        //- Return true if triB is the same as this triangle
+        inline bool isSame(const triangle2D& triB) const;
+
+        //- Return true if t point p is inside this triangle
+        inline bool pointInside(const vector2D& p) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "triangle2DI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/AMIInterpolation/triangle2D/triangle2DI.H

@@ -0,0 +1,241 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2021 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 "OFstream.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline void Foam::triangle2D::writeOBJ
+(
+    Ostream& os,
+    const triangle2D& tri,
+    label offset
+)
+{
+    os  << "v " << tri[0].x() << " " << tri[0].y() << " 0" << nl
+        << "v " << tri[1].x() << " " << tri[1].y() << " 0" << nl
+        << "v " << tri[2].x() << " " << tri[2].y() << " 0" << nl
+        << "l"
+        << " " << 1 + offset
+        << " " << 2 + offset
+        << " " << 3 + offset
+        << " " << 1 + offset
+        << endl;
+}
+
+
+inline Foam::label Foam::triangle2D::nClosePoints
+(
+    const triangle2D& triA,
+    const triangle2D& triB
+)
+{
+    label nClose = 0;
+
+    FixedList<bool, 3> match(true);
+
+    for (auto& a : triA)
+    {
+        forAll(triB, tb)
+        {
+            if (match[tb] && a.isClose(triB[tb]))
+            {
+                match[tb] = false;
+                ++nClose;
+                break;
+            }
+        }
+    }
+
+    return nClose;
+}
+
+
+inline Foam::scalar Foam::triangle2D::area
+(
+    const vector2D& a,
+    const vector2D& b,
+    const vector2D& c
+)
+{
+    const vector2D e1(a - c);
+    const vector2D e2(b - c);
+
+    return 0.5*e1.perp(e2);
+}
+
+
+inline Foam::vector2D Foam::triangle2D::centre() const
+{
+    const triangle2D& tri = *this;
+
+    return (1/3.0)*(tri[0] + tri[1] + tri[2]);
+}
+
+
+inline bool Foam::triangle2D::lineSegmentIntersectionPoint
+(
+    const vector2D& lp1,
+    const vector2D& lp2,
+    const vector2D& sp1,
+    const vector2D& sp2,
+    vector2D& intersection
+)
+{
+    const vector2D r(lp2 - lp1);
+    const vector2D s(sp2 - sp1);
+
+    const scalar rcs = r.perp(s);
+
+    if (mag(rcs) > ROOTVSMALL)
+    {
+        const scalar u = (sp1 - lp1).perp(r)/rcs;
+
+        if (u >= -relTol && u <= 1+relTol)
+        {
+            intersection = sp1 + u*s;
+            return true;
+        }
+    }
+
+/*
+        // Collapsed to line
+        if (mag(triangle2D(lp1, lp2, sp1).area()) < absTol)
+        {
+            intersection = sp1;
+            return true;
+        }
+
+        if (mag(triangle2D(lp1, lp2, sp2).area()) < absTol)
+        {
+            intersection = sp2;
+            return true;
+        }
+*/
+
+    if (debug)
+    {
+        OFstream os("bad-intersection.obj");
+        os  << "v " << lp1.x() << " " << lp1.y() << " 0" << nl
+            << "v " << lp2.x() << " " << lp2.y() << " 0" << nl
+            << "v " << sp1.x() << " " << sp1.y() << " 0" << nl
+            << "v " << sp2.x() << " " << sp2.y() << " 0" << nl
+            << "l 1 2"
+            << "l 3 4"
+            << endl;
+    }
+
+    return false;
+}
+
+
+inline bool Foam::triangle2D::lineIntersectionPoint
+(
+    const vector2D& a,
+    const vector2D& b,
+    const vector2D& c,
+    const vector2D& d,
+    vector2D& intersection
+)
+{
+    return lineSegmentIntersectionPoint(c, d, a, b, intersection);
+}
+
+
+inline bool Foam::triangle2D::lineIntersects
+(
+    const vector2D& a,
+    const vector2D& b,
+    const vector2D& c,
+    const vector2D& d
+)
+{
+    if
+    (
+        (triangle2D(a, c, d).order() == triangle2D(b, c, d).order())
+     || (triangle2D(a, b, c).order() == triangle2D(a, b, d).order())
+    )
+    {
+        return false;
+    }
+
+
+    DebugInfo<< "line " << a << b << " intersects " << c << d << endl;
+
+    return true;
+}
+
+
+inline Foam::scalar Foam::triangle2D::area() const noexcept
+{
+    return area_;
+}
+
+
+inline Foam::label Foam::triangle2D::order() const
+{
+    return mag(area_) < SMALL ? 0 : sign(area_);
+}
+
+
+inline bool Foam::triangle2D::contains(const triangle2D& tri) const
+{
+    return
+        pointInside(tri[0])
+     && pointInside(tri[1])
+     && pointInside(tri[2]);
+}
+
+
+inline bool Foam::triangle2D::isSame(const triangle2D& triB) const
+{
+    const triangle2D& triA = *this;
+
+    return
+        triA[0].isClose(triB[0])
+     && triA[1].isClose(triB[1])
+     && triA[2].isClose(triB[2]);
+}
+
+
+inline bool Foam::triangle2D::pointInside(const vector2D& p) const
+{
+    const triangle2D& triA = *this;
+
+    for (int i = 0; i < 3; ++i)
+    {
+        if ((triA[(i + 1) % 3] - triA[i]).perp(p - triA[i]) < 0)
+        {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/meshTools/Make/files

@@ -269,7 +269,6 @@ $(AMI)/AMIInterpolation/AMIInterpolationNew.C
 $(AMI)/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.C
 $(AMI)/AMIInterpolation/advancingFrontAMI/advancingFrontAMIParallelOps.C
 $(AMI)/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.C
-$(AMI)/AMIInterpolation/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.C
 $(AMI)/AMIInterpolation/nearestFaceAMI/nearestFaceAMI.C
 $(AMI)/faceAreaIntersect/faceAreaIntersect.C
 $(AMI)/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.C
@@ -277,6 +276,10 @@ $(AMI)/GAMG/interfaceFields/cyclicAMIGAMGInterfaceField/cyclicAMIGAMGInterfaceFi
 $(AMI)/GAMG/interfaces/cyclicACMIGAMGInterface/cyclicACMIGAMGInterface.C
 $(AMI)/GAMG/interfaceFields/cyclicACMIGAMGInterfaceField/cyclicACMIGAMGInterfaceField.C
 
+$(AMI)/triangle2D/triangle2D.C
+$(AMI)/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.C
+
+
 AMICycPatches=$(AMI)/patches/cyclicAMI
 $(AMICycPatches)/cyclicAMILduInterfaceField/cyclicAMILduInterface.C
 $(AMICycPatches)/cyclicAMILduInterfaceField/cyclicAMILduInterfaceField.C