ENH: AMIInterpolation - updated to perform face triangulation once only; code clean-up

applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C

@@ -233,10 +233,10 @@ int main(int argc, char *argv[])
         meshToMesh::interpolationMethod method =
             meshToMesh::interpolationMethodNames_[mapMethod];
 
-        patchMapMethod = AMIPatchToPatchInterpolation::interpolationMethodToWord
+        patchMapMethod = AMIPatchToPatchInterpolation::interpolationMethodNames_
-        (
+        [
             meshToMesh::interpolationMethodAMI(method)
-        );
+        ];
     }
 
     // Optionally override

src/functionObjects/field/mapFields/mapFields.C

@@ -89,10 +89,10 @@ void Foam::functionObjects::mapFields::createInterpolation
 
     // Lookup corresponding AMI method
     word patchMapMethodName =
-        AMIPatchToPatchInterpolation::interpolationMethodToWord
+        AMIPatchToPatchInterpolation::interpolationMethodNames_
-        (
+        [
             meshToMesh::interpolationMethodAMI(mapMethod)
-        );
+        ];
 
     // Optionally override
     if (dict.readIfPresent("patchMapMethod", patchMapMethodName))

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,26 +46,8 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolationMethodNames_
 };
 
 template<class SourcePatch, class TargetPatch>
-Foam::word
+bool Foam::AMIInterpolation<SourcePatch, TargetPatch>::cacheIntersections_ =
-Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolationMethodToWord
+    false;
-(
-    const interpolationMethod& im
-)
-{
-    return interpolationMethodNames_[im];
-}
-
-
-template<class SourcePatch, class TargetPatch>
-typename Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolationMethod
-Foam::AMIInterpolation<SourcePatch, TargetPatch>::wordTointerpolationMethod
-(
-    const word& im
-)
-{
-    return interpolationMethodNames_[im];
-}
-
 
 template<class SourcePatch, class TargetPatch>
 template<class Patch>
@@ -156,7 +138,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::projectPointsToSurface
 template<class SourcePatch, class TargetPatch>
 void Foam::AMIInterpolation<SourcePatch, TargetPatch>::normaliseWeights
 (
-    const scalarField& patchAreas,
+    const scalarList& patchAreas,
     const word& patchName,
     const labelListList& addr,
     scalarListList& wght,
@@ -237,14 +219,14 @@ template<class SourcePatch, class TargetPatch>
 void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
 (
     const autoPtr<mapDistribute>& targetMapPtr,
-    const scalarField& fineSrcMagSf,
+    const scalarList& fineSrcMagSf,
     const labelListList& fineSrcAddress,
     const scalarListList& fineSrcWeights,
 
     const labelList& sourceRestrictAddressing,
     const labelList& targetRestrictAddressing,
 
-    scalarField& srcMagSf,
+    scalarList& srcMagSf,
     labelListList& srcAddress,
     scalarListList& srcWeights,
     scalarField& srcWeightsSum,
@@ -339,9 +321,9 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
                 labelList oldToNew(targetCoarseSize, -1);
                 label newi = 0;
 
-                forAll(elems, i)
+                for (const label elemi : elems)
                 {
-                    label fineElem = elemsMap[elems[i]];
+                    label fineElem = elemsMap[elemi];
                     label coarseElem = allRestrict[fineElem];
                     if (oldToNew[coarseElem] == -1)
                     {
@@ -375,9 +357,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
             // reachable elements (using the local constructMap)
 
             const labelList& elemsMap = map.constructMap()[Pstream::myProcNo()];
-            forAll(elemsMap, i)
+            for (const label fineElem : elemsMap)
             {
-                label fineElem = elemsMap[i];
                 label coarseElem = allRestrict[fineElem];
                 tgtCompactMap[fineElem] = coarseElem;
             }
@@ -403,12 +384,12 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
                     // Get the maximum target coarse size for this set of
                     // received data.
                     label remoteTargetCoarseSize = labelMin;
-                    forAll(elems, i)
+                    for (const label elemi : elems)
                     {
                         remoteTargetCoarseSize = max
                         (
                             remoteTargetCoarseSize,
-                            allRestrict[elems[i]]
+                            allRestrict[elemi]
                         );
                     }
                     remoteTargetCoarseSize += 1;
@@ -417,9 +398,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
                     labelList oldToNew(remoteTargetCoarseSize, -1);
                     label newi = 0;
 
-                    forAll(elems, i)
+                    for (const label fineElem : elems)
                     {
-                        label fineElem = elems[i];
                         // fineElem now points to section from proci
                         label coarseElem = allRestrict[fineElem];
                         if (oldToNew[coarseElem] == -1)
@@ -563,9 +543,9 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::constructFromSurface
         if (debug)
         {
             OFstream os("amiSrcPoints.obj");
-            forAll(srcPoints, i)
+            for (const point& pt : srcPoints)
             {
-                meshTools::writeOBJ(os, srcPoints[i]);
+                meshTools::writeOBJ(os, pt);
             }
         }
 
@@ -584,9 +564,9 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::constructFromSurface
         if (debug)
         {
             OFstream os("amiTgtPoints.obj");
-            forAll(tgtPoints, i)
+            for (const point& pt : tgtPoints)
             {
-                meshTools::writeOBJ(os, tgtPoints[i]);
+                meshTools::writeOBJ(os, pt);
             }
         }
 
@@ -620,7 +600,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
     const bool reverseTarget
 )
 :
-    methodName_(interpolationMethodToWord(method)),
+    methodName_(interpolationMethodNames_[method]),
     reverseTarget_(reverseTarget),
     requireMatch_(requireMatch),
     singlePatchProc_(-999),
@@ -683,7 +663,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
     const bool reverseTarget
 )
 :
-    methodName_(interpolationMethodToWord(method)),
+    methodName_(interpolationMethodNames_[method]),
     reverseTarget_(reverseTarget),
     requireMatch_(requireMatch),
     singlePatchProc_(-999),
@@ -870,10 +850,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
         << tgtTotalSize << " target faces"
         << endl;
 
-    // Calculate face areas
-    srcMagSf_ = patchMagSf(srcPatch, triMode_);
-    tgtMagSf_ = patchMagSf(tgtPatch, triMode_);
-
     // Calculate if patches present on multiple processors
     singlePatchProc_ = calcDistribution(srcPatch, tgtPatch);
 
@@ -919,13 +895,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
                 newTgtPoints
             );
 
-        scalarField newTgtMagSf(newTgtPatch.size());
-        forAll(newTgtPatch, facei)
-        {
-            newTgtMagSf[facei] = newTgtPatch[facei].mag(newTgtPatch.points());
-        }
-
-
         // Calculate AMI interpolation
         autoPtr<AMIMethod<SourcePatch, TargetPatch>> AMIPtr
         (
@@ -934,8 +903,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
                 methodName_,
                 srcPatch,
                 newTgtPatch,
-                srcMagSf_,
-                newTgtMagSf,
                 triMode_,
                 reverseTarget_,
                 requireMatch_ && (lowWeightCorrection_ < 0)
@@ -950,6 +917,12 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
             tgtWeights_
         );
 
+        // Note: using patch face areas calculated by the AMI method
+        // - TODO: should move into the calculate method
+        srcMagSf_.transfer(AMIPtr->srcMagSf());
+        tgtMagSf_.transfer(AMIPtr->tgtMagSf());
+        map.reverseDistribute(tgtPatch.size(), tgtMagSf_);
+
         // Now
         // ~~~
         //  srcAddress_ :   per srcPatch face a list of the newTgtPatch (not
@@ -967,21 +940,19 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 
-        forAll(srcAddress_, i)
+        for (labelList& addressing : srcAddress_)
         {
-            labelList& addressing = srcAddress_[i];
+            for (label& addr : addressing)
-            forAll(addressing, addri)
             {
-                addressing[addri] = tgtFaceIDs[addressing[addri]];
+                addr = tgtFaceIDs[addr];
             }
         }
 
-        forAll(tgtAddress_, i)
+        for (labelList& addressing : tgtAddress_)
         {
-            labelList& addressing = tgtAddress_[i];
+            for (label& addr : addressing)
-            forAll(addressing, addri)
             {
-                addressing[addri] = globalSrcFaces.toGlobal(addressing[addri]);
+                addr = globalSrcFaces.toGlobal(addr);
             }
         }
 
@@ -1036,8 +1007,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
                 methodName_,
                 srcPatch,
                 tgtPatch,
-                srcMagSf_,
-                tgtMagSf_,
                 triMode_,
                 reverseTarget_,
                 requireMatch_ && (lowWeightCorrection_ < 0)
@@ -1052,6 +1021,9 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
             tgtWeights_
         );
 
+        srcMagSf_.transfer(AMIPtr->srcMagSf());
+        tgtMagSf_.transfer(AMIPtr->tgtMagSf());
+
         AMIPtr->normaliseWeights(true, *this);
     }
 
@@ -1623,9 +1595,8 @@ const
     nearest.setDistance(GREAT);
     label nearestFacei = -1;
 
-    forAll(addr, i)
+    for (const label srcFacei : addr)
     {
-        const label srcFacei = addr[i];
         const face& f = srcPatch[srcFacei];
 
         pointHit ray = f.ray(tgtPoint, n, srcPoints);
@@ -1672,9 +1643,8 @@ const
     // Target face addresses that intersect source face srcFacei
     const labelList& addr = srcAddress_[srcFacei];
 
-    forAll(addr, i)
+    for (const label tgtFacei : addr)
     {
-        const label tgtFacei = addr[i];
         const face& f = tgtPatch[tgtFacei];
 
         pointHit ray = f.ray(srcPoint, n, tgtPoints);
@@ -1719,9 +1689,8 @@ const
         const labelList& addr = srcAddress[i];
         const point& srcPt = srcPatch.faceCentres()[i];
 
-        forAll(addr, j)
+        for (const label tgtPti : addr)
         {
-            label tgtPti = addr[j];
             const point& tgtPt = tgtPatch.faceCentres()[tgtPti];
 
             meshTools::writeOBJ(os, srcPt);

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -95,17 +95,7 @@ public:
 
         static const Enum<interpolationMethod> interpolationMethodNames_;
 
-        //- Convert interpolationMethod to word representation
+        static bool cacheIntersections_;
-        static word interpolationMethodToWord
-        (
-            const interpolationMethod& method
-        );
-
-        //- Convert word to interpolationMethod
-        static interpolationMethod wordTointerpolationMethod
-        (
-            const word& method
-        );
 
         //- Calculate the patch face magnitudes for the given tri-mode
         template<class Patch>
@@ -142,7 +132,7 @@ private:
         // Source patch
 
             //- Source face areas
-            scalarField srcMagSf_;
+            scalarList srcMagSf_;
 
             //- Addresses of target faces per source face
             labelListList srcAddress_;
@@ -157,7 +147,7 @@ private:
         // Target patch
 
             //- Target face areas
-            scalarField tgtMagSf_;
+            scalarList tgtMagSf_;
 
             //- Addresses of source faces per target face
             labelListList tgtAddress_;
@@ -250,7 +240,7 @@ private:
             //     numerical error!
             static void normaliseWeights
             (
-                const scalarField& patchAreas,
+                const scalarList& patchAreas,
                 const word& patchName,
                 const labelListList& addr,
                 scalarListList& wght,
@@ -266,14 +256,14 @@ private:
             static void agglomerate
             (
                 const autoPtr<mapDistribute>& targetMap,
-                const scalarField& fineSrcMagSf,
+                const scalarList& fineSrcMagSf,
                 const labelListList& fineSrcAddress,
                 const scalarListList& fineSrcWeights,
 
                 const labelList& sourceRestrictAddressing,
                 const labelList& targetRestrictAddressing,
 
-                scalarField& srcMagSf,
+                scalarList& srcMagSf,
                 labelListList& srcAddress,
                 scalarListList& srcWeights,
                 scalarField& srcWeightsSum,
@@ -311,7 +301,7 @@ public:
             const faceAreaIntersect::triangulationMode& triMode,
             const bool requireMatch = true,
             const word& methodName =
-                interpolationMethodToWord(imFaceAreaWeight),
+                interpolationMethodNames_[imFaceAreaWeight],
             const scalar lowWeightCorrection = -1,
             const bool reverseTarget = false
         );
@@ -338,7 +328,7 @@ public:
             const faceAreaIntersect::triangulationMode& triMode,
             const bool requireMatch = true,
             const word& methodName =
-                interpolationMethodToWord(imFaceAreaWeight),
+                interpolationMethodNames_[imFaceAreaWeight],
             const scalar lowWeightCorrection = -1,
             const bool reverseTarget = false
         );
@@ -381,10 +371,10 @@ public:
             // Source patch
 
                 //- Return const access to source patch face areas
-                inline const scalarField& srcMagSf() const;
+                inline const List<scalar>& srcMagSf() const;
 
                 //- Return access to source patch face areas
-                inline scalarField& srcMagSf();
+                inline List<scalar>& srcMagSf();
 
                 //- Return const access to source patch addressing
                 inline const labelListList& srcAddress() const;
@@ -415,10 +405,10 @@ public:
             // Target patch
 
                 //- Return const access to target patch face areas
-                inline const scalarField& tgtMagSf() const;
+                inline const List<scalar>& tgtMagSf() const;
 
                 //- Return access to target patch face areas
-                inline scalarField& tgtMagSf();
+                inline List<scalar>& tgtMagSf();
 
                 //- Return const access to target patch addressing
                 inline const labelListList& tgtAddress() const;

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H

@@ -49,7 +49,7 @@ applyLowWeightCorrection() const
 
 
 template<class SourcePatch, class TargetPatch>
-inline const Foam::scalarField&
+inline const Foam::List<Foam::scalar>&
 Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcMagSf() const
 {
     return srcMagSf_;
@@ -57,7 +57,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcMagSf() const
 
 
 template<class SourcePatch, class TargetPatch>
-inline Foam::scalarField&
+inline Foam::List<Foam::scalar>&
 Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcMagSf()
 {
     return srcMagSf_;
@@ -121,7 +121,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcMap() const
 
 
 template<class SourcePatch, class TargetPatch>
-inline const Foam::scalarField&
+inline const Foam::List<Foam::scalar>&
 Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtMagSf() const
 {
     return tgtMagSf_;
@@ -129,7 +129,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtMagSf() const
 
 
 template<class SourcePatch, class TargetPatch>
-inline Foam::scalarField&
+inline Foam::List<Foam::scalar>&
 Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtMagSf()
 {
     return tgtMagSf_;

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationParallelOps.C

@@ -100,9 +100,9 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcOverlappingProcs
     {
         const treeBoundBoxList& bbp = procBb[proci];
 
-        forAll(bbp, bbi)
+        for (const treeBoundBox& tbb: bbp)
         {
-            if (bbp[bbi].overlaps(bb))
+            if (tbb.overlaps(bb))
             {
                 overlaps[proci] = true;
                 nOverlaps++;
@@ -128,7 +128,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::distributePatches
 {
     PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
 
-    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    for (label domain = 0; domain < Pstream::nProcs(); ++domain)
     {
         const labelList& sendElems = map.subMap()[domain];
 
@@ -195,7 +195,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::distributePatches
     }
 
     // Consume
-    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    for (label domain = 0; domain < Pstream::nProcs(); ++domain)
     {
         const labelList& recvElems = map.constructMap()[domain];
 
@@ -251,9 +251,8 @@ distributeAndMergePatches
         SubList<label>(tgtFaceIDs, faceIDs.size()) = faceIDs;
 
         const faceList& fcs = allFaces[Pstream::myProcNo()];
-        forAll(fcs, i)
+        for (const face& f : fcs)
         {
-            const face& f = fcs[i];
             face& newF = tgtFaces[nFaces++];
             newF.setSize(f.size());
             forAll(f, fp)
@@ -263,9 +262,9 @@ distributeAndMergePatches
         }
 
         const pointField& pts = allPoints[Pstream::myProcNo()];
-        forAll(pts, i)
+        for (const point& pt: pts)
         {
-            tgtPoints[nPoints++] = pts[i];
+            tgtPoints[nPoints++] = pt;
         }
     }
 
@@ -279,9 +278,8 @@ distributeAndMergePatches
             SubList<label>(tgtFaceIDs, faceIDs.size(), nFaces) = faceIDs;
 
             const faceList& fcs = allFaces[proci];
-            forAll(fcs, i)
+            for (const face& f : fcs)
             {
-                const face& f = fcs[i];
                 face& newF = tgtFaces[nFaces++];
                 newF.setSize(f.size());
                 forAll(f, fp)
@@ -291,9 +289,9 @@ distributeAndMergePatches
             }
 
             const pointField& pts = allPoints[proci];
-            forAll(pts, i)
+            for (const point& pt: pts)
             {
-                tgtPoints[nPoints++] = pts[i];
+                tgtPoints[nPoints++] = pt;
             }
         }
     }
@@ -319,9 +317,9 @@ distributeAndMergePatches
         }
 
         tgtPoints.transfer(newTgtPoints);
-        forAll(tgtFaces, i)
+        for (face& f : tgtFaces)
         {
-            inplaceRenumber(oldToNew, tgtFaces[i]);
+            inplaceRenumber(oldToNew, f);
         }
     }
 }
@@ -447,7 +445,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcProcMap
             label nRecv = sendSizes[proci][Pstream::myProcNo()];
             constructMap[proci].setSize(nRecv);
 
-            for (label i = 0; i < nRecv; i++)
+            for (label i = 0; i < nRecv; ++i)
             {
                 constructMap[proci][i] = segmentI++;
             }

src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -167,9 +167,9 @@ void Foam::AMIMethod<SourcePatch, TargetPatch>::writeIntersectionOBJ
 
     OFstream os("areas" + name(count) + ".obj");
 
-    forAll(f1pts, i)
+    for (const point& pt : f1pts)
     {
-        meshTools::writeOBJ(os, f1pts[i]);
+        meshTools::writeOBJ(os, pt);
     }
     os<< "l";
     forAll(f1pts, i)
@@ -178,19 +178,19 @@ void Foam::AMIMethod<SourcePatch, TargetPatch>::writeIntersectionOBJ
     }
     os<< " 1" << endl;
 
-
+    for (const point& pt : f2pts)
-    forAll(f2pts, i)
     {
-        meshTools::writeOBJ(os, f2pts[i]);
+        meshTools::writeOBJ(os, pt);
     }
     os<< "l";
+    const label n = f1pts.size();
     forAll(f2pts, i)
     {
-        os<< " " << f1pts.size() + i + 1;
+        os<< " " << n + i + 1;
     }
-    os<< " " << f1pts.size() + 1 << endl;
+    os<< " " << n + 1 << endl;
 
-    count++;
+    ++count;
 }
 
 
@@ -265,16 +265,17 @@ void Foam::AMIMethod<SourcePatch, TargetPatch>::appendNbrFaces
     DynamicList<label>& faceIDs
 ) const
 {
+    static const scalar thetaCos = Foam::cos(degToRad(89.0));
+
     const labelList& nbrFaces = patch.faceFaces()[facei];
 
     // filter out faces already visited from face neighbours
-    forAll(nbrFaces, i)
+    for (const label nbrFacei : nbrFaces)
     {
-        label nbrFacei = nbrFaces[i];
         bool valid = true;
-        forAll(visitedFaces, j)
+        for (const label visitedFacei : visitedFaces)
         {
-            if (nbrFacei == visitedFaces[j])
+            if (nbrFacei == visitedFacei)
             {
                 valid = false;
                 break;
@@ -283,9 +284,9 @@ void Foam::AMIMethod<SourcePatch, TargetPatch>::appendNbrFaces
 
         if (valid)
         {
-            forAll(faceIDs, j)
+            for (const label testFacei : faceIDs)
             {
-                if (nbrFacei == faceIDs[j])
+                if (nbrFacei == testFacei)
                 {
                     valid = false;
                     break;
@@ -301,7 +302,7 @@ void Foam::AMIMethod<SourcePatch, TargetPatch>::appendNbrFaces
 
             const scalar cosI = n1 & n2;
 
-            if (cosI > Foam::cos(degToRad(89.0)))
+            if (cosI > thetaCos)
             {
                 faceIDs.append(nbrFacei);
             }
@@ -310,6 +311,52 @@ void Foam::AMIMethod<SourcePatch, TargetPatch>::appendNbrFaces
 }
 
 
+template<class SourcePatch, class TargetPatch>
+template<class PatchType>
+void Foam::AMIMethod<SourcePatch, TargetPatch>::triangulatePatch
+(
+    const PatchType& patch,
+    List<DynamicList<face>>& tris,
+    List<scalar>& magSf
+) const
+{
+    const pointField& points = patch.points();
+    tris.setSize(patch.size());
+    magSf.setSize(patch.size());
+
+    // Using methods that index into existing points
+    forAll(patch, facei)
+    {
+        switch (triMode_)
+        {
+            case faceAreaIntersect::tmFan:
+            {
+                faceAreaIntersect::triangleFan(patch[facei], tris[facei]);
+                break;
+            }
+            case faceAreaIntersect::tmMesh:
+            {
+                patch[facei].triangles(points, tris[facei]);
+                break;
+            }
+        }
+
+        const DynamicList<face>& triFaces = tris[facei];
+        magSf[facei] = 0;
+        for (const face& f : triFaces)
+        {
+            magSf[facei] +=
+                triPointRef
+                (
+                    points[f[0]],
+                    points[f[1]],
+                    points[f[2]]
+                ).mag();
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class SourcePatch, class TargetPatch>
@@ -317,8 +364,6 @@ Foam::AMIMethod<SourcePatch, TargetPatch>::AMIMethod
 (
     const SourcePatch& srcPatch,
     const TargetPatch& tgtPatch,
-    const scalarField& srcMagSf,
-    const scalarField& tgtMagSf,
     const faceAreaIntersect::triangulationMode& triMode,
     const bool reverseTarget,
     const bool requireMatch
@@ -328,18 +373,14 @@ Foam::AMIMethod<SourcePatch, TargetPatch>::AMIMethod
     tgtPatch_(tgtPatch),
     reverseTarget_(reverseTarget),
     requireMatch_(requireMatch),
-    srcMagSf_(srcMagSf),
+    srcMagSf_(srcPatch.size(), 1.0),
-    tgtMagSf_(tgtMagSf),
+    tgtMagSf_(tgtPatch.size(), 1.0),
     srcNonOverlap_(),
     triMode_(triMode)
-{}
+{
-
+    // Note: setting srcMagSf and tgtMagSf to 1 by default for 1-to-1 methods
-
+    // - others will need to overwrite as necessary
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+}
-
-template<class SourcePatch, class TargetPatch>
-Foam::AMIMethod<SourcePatch, TargetPatch>::~AMIMethod()
-{}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpcnCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -92,10 +92,10 @@ protected:
         const bool requireMatch_;
 
         //- Source face areas
-        const scalarField& srcMagSf_;
+        List<scalar> srcMagSf_;
 
         //- Target face areas
-        const scalarField& tgtMagSf_;
+        List<scalar> tgtMagSf_;
 
         //- Labels of faces that are not overlapped by any target faces
         //  (should be empty for correct functioning)
@@ -154,6 +154,15 @@ protected:
                 DynamicList<label>& faceIDs
             ) const;
 
+            //- Helper function to decompose a patch
+            template<class PatchType>
+            void triangulatePatch
+            (
+                const PatchType& patch,
+                List<DynamicList<face>>& tris,
+                List<scalar>& magSf
+            ) const;
+
 
 public:
 
@@ -169,8 +178,6 @@ public:
         (
             const SourcePatch& srcPatch,
             const TargetPatch& tgtPatch,
-            const scalarField& srcMagSf,
-            const scalarField& tgtMagSf,
             const faceAreaIntersect::triangulationMode& triMode,
             const bool reverseTarget,
             const bool requireMatch
@@ -178,8 +185,6 @@ public:
         (
             srcPatch,
             tgtPatch,
-            srcMagSf,
-            tgtMagSf,
             triMode,
             reverseTarget,
             requireMatch
@@ -191,8 +196,6 @@ public:
     (
         const SourcePatch& srcPatch,
         const TargetPatch& tgtPatch,
-        const scalarField& srcMagSf,
-        const scalarField& tgtMagSf,
         const faceAreaIntersect::triangulationMode& triMode,
         const bool reverseTarget,
         const bool requireMatch
@@ -204,8 +207,6 @@ public:
         const word& methodName,
         const SourcePatch& srcPatch,
         const TargetPatch& tgtPatch,
-        const scalarField& srcMagSf,
-        const scalarField& tgtMagSf,
         const faceAreaIntersect::triangulationMode& triMode,
         const bool reverseTarget,
         const bool requireMatch
@@ -213,7 +214,7 @@ public:
 
 
     //- Destructor
-    virtual ~AMIMethod();
+    virtual ~AMIMethod() = default;
 
 
     // Member Functions
@@ -227,6 +228,18 @@ public:
             //- Flag to indicate that interpolation patches are conformal
             virtual bool conformal() const;
 
+            //- Return const access to source patch face areas
+            inline const List<scalar>& srcMagSf() const;
+
+            //- Return access to source patch face areas
+            inline List<scalar>& srcMagSf();
+
+            //- Return const access to target patch face areas
+            inline const List<scalar>& tgtMagSf() const;
+
+            //- Return access to target patch face areas
+            inline List<scalar>& tgtMagSf();
+
 
         // Manipulation
 

src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethodI.H

@@ -23,6 +23,38 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+template<class SourcePatch, class TargetPatch>
+inline const Foam::List<Foam::scalar>&
+Foam::AMIMethod<SourcePatch, TargetPatch>::srcMagSf() const
+{
+    return srcMagSf_;
+}
+
+
+template<class SourcePatch, class TargetPatch>
+inline Foam::List<Foam::scalar>&
+Foam::AMIMethod<SourcePatch, TargetPatch>::srcMagSf()
+{
+    return srcMagSf_;
+}
+
+
+template<class SourcePatch, class TargetPatch>
+inline const Foam::List<Foam::scalar>&
+Foam::AMIMethod<SourcePatch, TargetPatch>::tgtMagSf() const
+{
+    return tgtMagSf_;
+}
+
+
+template<class SourcePatch, class TargetPatch>
+inline Foam::List<Foam::scalar>&
+Foam::AMIMethod<SourcePatch, TargetPatch>::tgtMagSf()
+{
+    return tgtMagSf_;
+}
+
+
 template<class SourcePatch, class TargetPatch>
 inline const Foam::labelList&
 Foam::AMIMethod<SourcePatch, TargetPatch>::srcNonOverlap() const

src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethodNew.C

@@ -32,8 +32,6 @@ Foam::AMIMethod<SourcePatch, TargetPatch>::New
     const word& methodName,
     const SourcePatch& srcPatch,
     const TargetPatch& tgtPatch,
-    const scalarField& srcMagSf,
-    const scalarField& tgtMagSf,
     const faceAreaIntersect::triangulationMode& triMode,
     const bool reverseTarget,
     const bool requireMatch
@@ -62,8 +60,6 @@ Foam::AMIMethod<SourcePatch, TargetPatch>::New
         (
             srcPatch,
             tgtPatch,
-            srcMagSf,
-            tgtMagSf,
             triMode,
             reverseTarget,
             requireMatch

src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C

@@ -46,10 +46,8 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
 
     const vectorField& srcCf = this->srcPatch_.faceCentres();
 
-    forAll(srcNbr, i)
+    for (const label srcI : srcNbr)
     {
-        label srcI = srcNbr[i];
-
         if ((mapFlag[srcI] == 0) && (srcTgtSeed[srcI] == -1))
         {
             // first attempt: match by comparing face centres
@@ -69,9 +67,8 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
             tol = max(SMALL, 0.0001*sqrt(tol));
 
             bool found = false;
-            forAll(tgtNbr, j)
+            for (const label tgtI : tgtNbr)
             {
-                label tgtI = tgtNbr[j];
                 const face& tgtF = this->tgtPatch_[tgtI];
                 const point tgtC = tgtF.centre(tgtPoints);
 
@@ -92,9 +89,8 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
             {
                 const vector srcN = srcF.normal(srcPoints);
 
-                forAll(tgtNbr, j)
+                for (const label tgtI : tgtNbr)
                 {
-                    label tgtI = tgtNbr[j];
                     const face& tgtF = this->tgtPatch_[tgtI];
                     pointHit ray = tgtF.ray(srcCf[srcI], srcN, tgtPoints);
 
@@ -126,9 +122,8 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
                         << endl;
 
                     Pout<< "target neighbours:" << nl;
-                    forAll(tgtNbr, j)
+                    for (const label tgtI : tgtNbr)
                     {
-                        label tgtI = tgtNbr[j];
                         const face& tgtF = this->tgtPatch_[tgtI];
 
                         Pout<< "face id: " << tgtI
@@ -192,8 +187,6 @@ Foam::directAMI<SourcePatch, TargetPatch>::directAMI
 (
     const SourcePatch& srcPatch,
     const TargetPatch& tgtPatch,
-    const scalarField& srcMagSf,
-    const scalarField& tgtMagSf,
     const faceAreaIntersect::triangulationMode& triMode,
     const bool reverseTarget,
     const bool requireMatch
@@ -203,8 +196,6 @@ Foam::directAMI<SourcePatch, TargetPatch>::directAMI
     (
         srcPatch,
         tgtPatch,
-        srcMagSf,
-        tgtMagSf,
         triMode,
         reverseTarget,
         requireMatch

src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H

@@ -109,8 +109,6 @@ public:
         (
             const SourcePatch& srcPatch,
             const TargetPatch& tgtPatch,
-            const scalarField& srcMagSf,
-            const scalarField& tgtMagSf,
             const faceAreaIntersect::triangulationMode& triMode,
             const bool reverseTarget = false,
             const bool requireMatch = true

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

@@ -143,6 +143,8 @@ bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
 
     bool faceProcessed = false;
 
+    label maxNeighbourFaces = nbrFaces.size();
+
     do
     {
         // process new target face
@@ -168,10 +170,17 @@ bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
             );
 
             faceProcessed = true;
+
+            maxNeighbourFaces = max(maxNeighbourFaces, nbrFaces.size());
         }
 
     } while (nbrFaces.size() > 0);
 
+    if (debug > 1)
+    {
+        DebugVar(maxNeighbourFaces);
+    }
+
     return faceProcessed;
 }
 
@@ -201,14 +210,13 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
         {
             for (label faceT : visitedFaces)
             {
-                scalar area = interArea(faceS, faceT);
+                const scalar threshold =
-                scalar areaTotal = this->srcMagSf_[srcFacei];
+                    this->srcMagSf_[faceS]*faceAreaIntersect::tolerance();
 
+                // store when intersection fractional area > tolerance
                 // Check that faces have enough overlap for robust walking
-                if (area/areaTotal > faceAreaIntersect::tolerance())
+                if (overlaps(faceS, faceT, threshold))
                 {
-                    // TODO - throwing area away - re-use in next iteration?
-
                     seedFaces[faceS] = faceT;
 
                     if (!valuesSet)
@@ -231,7 +239,7 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
     {
         // try to use existing seed
         bool foundNextSeed = false;
-        for (label facei = startSeedi; facei < mapFlag.size(); facei++)
+        for (label facei = startSeedi; facei < mapFlag.size(); ++facei)
         {
             if (mapFlag[facei])
             {
@@ -259,7 +267,7 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
         }
 
         foundNextSeed = false;
-        for (label facei = startSeedi; facei < mapFlag.size(); facei++)
+        for (label facei = startSeedi; facei < mapFlag.size(); ++facei)
         {
             if (mapFlag[facei])
             {
@@ -305,16 +313,25 @@ Foam::scalar Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
     const face& tgt = this->tgtPatch_[tgtFacei];
 
     // quick reject if either face has zero area
-    // Note: do not use stored face areas for target patch
+    if
-    const scalar tgtMag = tgt.mag(tgtPoints);
+    (
-    if ((this->srcMagSf_[srcFacei] < ROOTVSMALL) || (tgtMag < ROOTVSMALL))
+        (this->srcMagSf_[srcFacei] < ROOTVSMALL)
+     || (this->tgtMagSf_[tgtFacei] < ROOTVSMALL)
+    )
     {
         return area;
     }
 
     // create intersection object
-    bool cache = true;
+    faceAreaIntersect inter
-    faceAreaIntersect inter(srcPoints, tgtPoints, this->reverseTarget_, cache);
+    (
+        srcPoints,
+        tgtPoints,
+        this->srcTris_[srcFacei],
+        this->tgtTris_[tgtFacei],
+        this->reverseTarget_,
+        AMIInterpolation<SourcePatch, TargetPatch>::cacheIntersections_
+    );
 
     // crude resultant norm
     vector n(-this->srcPatch_.faceNormals()[srcFacei]);
@@ -330,9 +347,7 @@ Foam::scalar Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
 
     if (magN > ROOTVSMALL)
     {
-        area = inter.calc(src, tgt, n/magN, this->triMode_);
+        area = inter.calc(src, tgt, n/magN);
-        DebugVar(inter.triangles());
-        DebugVar(inter.triangles().size());
     }
     else
     {
@@ -354,6 +369,71 @@ Foam::scalar Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
 }
 
 
+template<class SourcePatch, class TargetPatch>
+bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::overlaps
+(
+    const label srcFacei,
+    const label tgtFacei,
+    const scalar threshold
+) const
+{
+    const pointField& srcPoints = this->srcPatch_.points();
+    const pointField& tgtPoints = this->tgtPatch_.points();
+
+    // references to candidate faces
+    const face& src = this->srcPatch_[srcFacei];
+    const face& tgt = this->tgtPatch_[tgtFacei];
+
+    // quick reject if either face has zero area
+    if
+    (
+        (this->srcMagSf_[srcFacei] < ROOTVSMALL)
+     || (this->tgtMagSf_[tgtFacei] < ROOTVSMALL)
+    )
+    {
+        return false;
+    }
+
+    faceAreaIntersect inter
+    (
+        srcPoints,
+        tgtPoints,
+        this->srcTris_[srcFacei],
+        this->tgtTris_[tgtFacei],
+        this->reverseTarget_,
+        AMIInterpolation<SourcePatch, TargetPatch>::cacheIntersections_
+    );
+
+    // crude resultant norm
+    vector n(-this->srcPatch_.faceNormals()[srcFacei]);
+    if (this->reverseTarget_)
+    {
+        n -= this->tgtPatch_.faceNormals()[tgtFacei];
+    }
+    else
+    {
+        n += this->tgtPatch_.faceNormals()[tgtFacei];
+    }
+    scalar magN = mag(n);
+
+    if (magN > ROOTVSMALL)
+    {
+        return inter.overlaps(src, tgt, n/magN, threshold);
+    }
+    else
+    {
+        WarningInFunction
+            << "Invalid normal for source face " << srcFacei
+            << " points " << UIndirectList<point>(srcPoints, src)
+            << " target face " << tgtFacei
+            << " points " << UIndirectList<point>(tgtPoints, tgt)
+            << endl;
+    }
+
+    return false;
+}
+
+
 template<class SourcePatch, class TargetPatch>
 void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::
 restartUncoveredSourceFace
@@ -458,8 +538,6 @@ Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::faceAreaWeightAMI
 (
     const SourcePatch& srcPatch,
     const TargetPatch& tgtPatch,
-    const scalarField& srcMagSf,
-    const scalarField& tgtMagSf,
     const faceAreaIntersect::triangulationMode& triMode,
     const bool reverseTarget,
     const bool requireMatch,
@@ -470,14 +548,17 @@ Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::faceAreaWeightAMI
     (
         srcPatch,
         tgtPatch,
-        srcMagSf,
-        tgtMagSf,
         triMode,
         reverseTarget,
         requireMatch
     ),
-    restartUncoveredSourceFace_(restartUncoveredSourceFace)
+    restartUncoveredSourceFace_(restartUncoveredSourceFace),
-{}
+    srcTris_(),
+    tgtTris_()
+{
+    this->triangulatePatch(srcPatch, srcTris_, this->srcMagSf_);
+    this->triangulatePatch(tgtPatch, tgtTris_, this->tgtMagSf_);
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //

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

@@ -58,6 +58,12 @@ private:
         //- Flag to restart uncovered source faces
         const bool restartUncoveredSourceFace_;
 
+        //- Storage for src-side triangle decomposition
+        List<DynamicList<face>> srcTris_;
+
+        //- Storage for tgt-side triangle decomposition
+        List<DynamicList<face>> tgtTris_;
+
 
 protected:
 
@@ -126,6 +132,14 @@ protected:
                 const label tgtFacei
             ) const;
 
+            //- Return true if faces overlap
+            virtual bool overlaps
+            (
+                const label srcFacei,
+                const label tgtFacei,
+                const scalar threshold
+            ) const;
+
 
 public:
 
@@ -140,8 +154,6 @@ public:
         (
             const SourcePatch& srcPatch,
             const TargetPatch& tgtPatch,
-            const scalarField& srcMagSf,
-            const scalarField& tgtMagSf,
             const faceAreaIntersect::triangulationMode& triMode,
             const bool reverseTarget = false,
             const bool requireMatch = true,

src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C

@@ -85,9 +85,8 @@ void Foam::mapNearestAMI<SourcePatch, TargetPatch>::setNextNearestFaces
 
     srcFacei = -1;
 
-    forAll(srcNbr, i)
+    for (const label facei : srcNbr)
     {
-        label facei = srcNbr[i];
         if (mapFlag[facei])
         {
             srcFacei = facei;
@@ -149,11 +148,11 @@ Foam::label Foam::mapNearestAMI<SourcePatch, TargetPatch>::findMappedSrcFace
             {
                 const labelList& nbrFaces = this->tgtPatch_.faceFaces()[tgtI];
 
-                forAll(nbrFaces, i)
+                for (const label nbrFacei : nbrFaces)
                 {
-                    if (!visitedFaces.found(nbrFaces[i]))
+                    if (findIndex(visitedFaces, nbrFacei) == -1)
                     {
-                        testFaces.append(nbrFaces[i]);
+                        testFaces.append(nbrFacei);
                     }
                 }
             }
@@ -172,8 +171,6 @@ Foam::mapNearestAMI<SourcePatch, TargetPatch>::mapNearestAMI
 (
     const SourcePatch& srcPatch,
     const TargetPatch& tgtPatch,
-    const scalarField& srcMagSf,
-    const scalarField& tgtMagSf,
     const faceAreaIntersect::triangulationMode& triMode,
     const bool reverseTarget,
     const bool requireMatch
@@ -183,8 +180,6 @@ Foam::mapNearestAMI<SourcePatch, TargetPatch>::mapNearestAMI
     (
         srcPatch,
         tgtPatch,
-        srcMagSf,
-        tgtMagSf,
         triMode,
         reverseTarget,
         requireMatch
@@ -280,7 +275,7 @@ void Foam::mapNearestAMI<SourcePatch, TargetPatch>::calculate
             label srcFacei = srcFaces[0];
             scalar d = magSqr(tgtC - srcCf[srcFacei]);
 
-            for (label i = 1; i < srcFaces.size(); i++)
+            for (label i = 1; i < srcFaces.size(); ++i)
             {
                 label srcI = srcFaces[i];
                 scalar dNew = magSqr(tgtC - srcCf[srcI]);
@@ -375,7 +370,7 @@ void Foam::mapNearestAMI<SourcePatch, TargetPatch>::normaliseWeights
                 label minFaceI = addr[0];
                 scalar minWeight = wghts[0];
 
-                for (label i = 0; i < addr.size(); i++)
+                for (label i = 0; i < addr.size(); ++i)
                 {
                     if (wghts[i] < minWeight)
                     {
@@ -408,7 +403,7 @@ void Foam::mapNearestAMI<SourcePatch, TargetPatch>::normaliseWeights
                 label minFaceI = addr[0];
                 scalar minWeight = wghts[0];
 
-                for (label i = 0; i < addr.size(); i++)
+                forAll(addr, i)
                 {
                     if (wghts[i] < minWeight)
                     {

src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.H

@@ -113,8 +113,6 @@ public:
         (
             const SourcePatch& srcPatch,
             const TargetPatch& tgtPatch,
-            const scalarField& srcMagSf,
-            const scalarField& tgtMagSf,
             const faceAreaIntersect::triangulationMode& triMode,
             const bool reverseTarget = false,
             const bool requireMatch = true

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

@@ -60,8 +60,6 @@ partialFaceAreaWeightAMI
 (
     const SourcePatch& srcPatch,
     const TargetPatch& tgtPatch,
-    const scalarField& srcMagSf,
-    const scalarField& tgtMagSf,
     const faceAreaIntersect::triangulationMode& triMode,
     const bool reverseTarget,
     const bool requireMatch
@@ -71,8 +69,6 @@ partialFaceAreaWeightAMI
     (
         srcPatch,
         tgtPatch,
-        srcMagSf,
-        tgtMagSf,
         triMode,
         reverseTarget,
         requireMatch

src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.H

@@ -90,8 +90,6 @@ public:
         (
             const SourcePatch& srcPatch,
             const TargetPatch& tgtPatch,
-            const scalarField& srcMagSf,
-            const scalarField& tgtMagSf,
             const faceAreaIntersect::triangulationMode& triMode,
             const bool reverseTarget = false,
             const bool requireMatch = true

src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,12 +34,11 @@ const Foam::Enum
 Foam::faceAreaIntersect::triangulationModeNames_
 {
     { triangulationMode::tmFan, "fan" },
-    { triangulationMode::tmMesh, "mesh" },
+    { triangulationMode::tmMesh, "mesh" }
 };
 
 Foam::scalar Foam::faceAreaIntersect::tol = 1e-6;
 
-
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 void Foam::faceAreaIntersect::triSliceWithPlane
@@ -49,7 +48,7 @@ void Foam::faceAreaIntersect::triSliceWithPlane
     FixedList<triPoints, 10>& tris,
     label& nTris,
     const scalar len
-)
+) const
 {
     // distance to cutting plane
     FixedList<scalar, 3> d;
@@ -220,9 +219,11 @@ void Foam::faceAreaIntersect::triSliceWithPlane
 Foam::scalar Foam::faceAreaIntersect::triangleIntersect
 (
     const triPoints& src,
-    const triPoints& tgt,
+    const point& tgt0,
+    const point& tgt1,
+    const point& tgt2,
     const vector& n
-)
+) const
 {
     // Work storage
     FixedList<triPoints, 10> workTris1;
@@ -241,8 +242,9 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
         // Cut triangle src with plane and put resulting sub-triangles in
         // workTris1 list
 
-        scalar s = mag(tgt[1] - tgt[0]);
+        scalar s = mag(tgt1 - tgt0);
-        plane pl0(tgt[0], tgt[1], tgt[1] + s*n);
+        //plane pl0(tgt0, tgt1, tgt1 + s*n);
+        plane pl0(tgt0, (tgt0 - tgt1)^(-s*n), true);
         triSliceWithPlane(src, pl0, workTris1, nWorkTris1, t);
     }
 
@@ -256,12 +258,13 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
         // Cut workTris1 with plane and put resulting sub-triangles in
         // workTris2 list (re-use tris storage)
 
-        scalar s = mag(tgt[2] - tgt[1]);
+        scalar s = mag(tgt2 - tgt1);
-        plane pl1(tgt[1], tgt[2], tgt[2] + s*n);
+        //plane pl1(tgt1, tgt2, tgt2 + s*n);
+        plane pl1(tgt1, (tgt1 - tgt2)^(-s*n), true);
 
         nWorkTris2 = 0;
 
-        for (label i = 0; i < nWorkTris1; i++)
+        for (label i = 0; i < nWorkTris1; ++i)
         {
             triSliceWithPlane(workTris1[i], pl1, workTris2, nWorkTris2, t);
         }
@@ -277,12 +280,13 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
         // Cut workTris2 with plane and put resulting sub-triangles in
         // workTris1 list (re-use workTris1 storage)
 
-        scalar s = mag(tgt[2] - tgt[0]);
+        scalar s = mag(tgt2 - tgt0);
-        plane pl2(tgt[2], tgt[0], tgt[0] + s*n);
+        //plane pl2(tgt2, tgt0, tgt0 + s*n);
+        plane pl2(tgt2, (tgt2 - tgt0)^(-s*n), true);
 
         nWorkTris1 = 0;
 
-        for (label i = 0; i < nWorkTris2; i++)
+        for (label i = 0; i < nWorkTris2; ++i)
         {
             triSliceWithPlane(workTris2[i], pl2, workTris1, nWorkTris1, t);
         }
@@ -295,7 +299,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
         {
             // Calculate area of sub-triangles
             scalar area = 0.0;
-            for (label i = 0; i < nWorkTris1; i++)
+            for (label i = 0; i < nWorkTris1; ++i)
             {
                 area += triArea(workTris1[i]);
 
@@ -317,12 +321,16 @@ Foam::faceAreaIntersect::faceAreaIntersect
 (
     const pointField& pointsA,
     const pointField& pointsB,
+    const DynamicList<face>& trisA,
+    const DynamicList<face>& trisB,
     const bool reverseB,
     const bool cacheTriangulation
 )
 :
     pointsA_(pointsA),
     pointsB_(pointsB),
+    trisA_(trisA),
+    trisB_(trisB),
     reverseB_(reverseB),
     cacheTriangulation_(cacheTriangulation),
     triangles_(cacheTriangulation ? 10 : 0)
@@ -343,7 +351,7 @@ void Foam::faceAreaIntersect::triangulate
 
     switch (triMode)
     {
-        case tmFan:
+        case triangulationMode::tmFan:
         {
             for (label i = 0; i < f.nTriangles(); ++i)
             {
@@ -355,7 +363,7 @@ void Foam::faceAreaIntersect::triangulate
 
             break;
         }
-        case tmMesh:
+        case triangulationMode::tmMesh:
         {
             const label nFaceTris = f.nTriangles();
 
@@ -374,37 +382,50 @@ void Foam::faceAreaIntersect::triangulate
 }
 
 
+
 Foam::scalar Foam::faceAreaIntersect::calc
 (
     const face& faceA,
     const face& faceB,
-    const vector& n,
+    const vector& n
-    const triangulationMode& triMode
+) const
-)
 {
-    // split faces into triangles
+    if (cacheTriangulation_)
-    faceList trisA, trisB;
+    {
-    triangulate(faceA, pointsA_, triMode, trisA);
+        triangles_.clear();
-    triangulate(faceB, pointsB_, triMode, trisB);
+    }
 
     // Intersect triangles
     scalar totalArea = 0.0;
-    triangles_.clear();
+    for (const face& triA : trisA_)
+    {
+        triPoints tpA = getTriPoints(pointsA_, triA, false);
 
-    forAll(trisA, tA)
+        for (const face& triB : trisB_)
         {
-        triPoints tpA = getTriPoints(pointsA_, trisA[tA], false);
+            if (reverseB_)
-
-//        if (triArea(tpA) > ROOTVSMALL)
             {
-            forAll(trisB, tB)
+                totalArea +=
-            {
+                    triangleIntersect
-                triPoints tpB = getTriPoints(pointsB_, trisB[tB], !reverseB_);
+                    (
-
+                        tpA,
-//                if (triArea(tpB) > ROOTVSMALL)
+                        pointsB_[triB[0]],
-                {
+                        pointsB_[triB[1]],
-                    totalArea += triangleIntersect(tpA, tpB, n);
+                        pointsB_[triB[2]],
+                        n
+                    );
             }
+            else
+            {
+                totalArea +=
+                    triangleIntersect
+                    (
+                        tpA,
+                        pointsB_[triB[2]],
+                        pointsB_[triB[1]],
+                        pointsB_[triB[0]],
+                        n
+                    );
             }
         }
     }
@@ -413,4 +434,56 @@ Foam::scalar Foam::faceAreaIntersect::calc
 }
 
 
+bool Foam::faceAreaIntersect::overlaps
+(
+    const face& faceA,
+    const face& faceB,
+    const vector& n,
+    const scalar threshold
+) const
+{
+    // Intersect triangles
+    scalar totalArea = 0.0;
+    for (const face& triA : trisA_)
+    {
+        const triPoints tpA = getTriPoints(pointsA_, triA, false);
+
+        for (const face& triB : trisB_)
+        {
+            if (reverseB_)
+            {
+                totalArea +=
+                    triangleIntersect
+                    (
+                        tpA,
+                        pointsB_[triB[0]],
+                        pointsB_[triB[1]],
+                        pointsB_[triB[2]],
+                        n
+                    );
+            }
+            else
+            {
+                totalArea +=
+                    triangleIntersect
+                    (
+                        tpA,
+                        pointsB_[triB[2]],
+                        pointsB_[triB[1]],
+                        pointsB_[triB[0]],
+                        n
+                    );
+            }
+
+            if (totalArea > threshold)
+            {
+                return true;
+            }
+        }
+    }
+
+    return false;
+}
+
+
 // ************************************************************************* //

src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -70,12 +70,18 @@ private:
 
     // Private data
 
-        //- Reference to the points of sideA
+        //- Reference to the points of face A
         const pointField& pointsA_;
 
-        //- Reference to the points of sideB
+        //- Reference to the points of face B
         const pointField& pointsB_;
 
+        //- Triangle decomposition of face A
+        const DynamicList<face>& trisA_;
+
+        //- Triangle decomposition of face B
+        const DynamicList<face>& trisB_;
+
         //- Flag to reverse B faces
         const bool reverseB_;
 
@@ -83,7 +89,7 @@ private:
         bool cacheTriangulation_;
 
         //- Final triangulation
-        DynamicList<triPoints> triangles_;
+        mutable DynamicList<triPoints> triangles_;
 
 
     // Static data members
@@ -123,7 +129,6 @@ private:
         //- Return triangle area
         inline scalar triArea(const triPoints& t) const;
 
-
         //- Slice triangle with plane and generate new cut sub-triangles
         void triSliceWithPlane
         (
@@ -132,15 +137,17 @@ private:
             FixedList<triPoints, 10>& tris,
             label& nTris,
             const scalar len
-        );
+        ) const;
 
         //- Return area of intersection of triangles src and tgt
         scalar triangleIntersect
         (
             const triPoints& src,
-            const triPoints& tgt,
+            const point& tgt0,
+            const point& tgt1,
+            const point& tgt2,
             const vector& n
-        );
+        ) const;
 
 
 public:
@@ -152,6 +159,8 @@ public:
         (
             const pointField& pointsA,
             const pointField& pointsB,
+            const DynamicList<face>& trisA,
+            const DynamicList<face>& trisB,
             const bool reverseB = false,
             const bool cacheTriangulation = false
         );
@@ -180,14 +189,29 @@ public:
         //- Non-const access to the triangulation
         inline DynamicList<triPoints>& triangles();
 
+        //- Decompose face into triangle fan
+        static inline void triangleFan
+        (
+            const face& f,
+            DynamicList<face>& faces
+        );
+
         //- Return area of intersection of faceA with faceB
         scalar calc
+        (
+            const face& faceA,
+            const face& faceB,
+            const vector& n
+        ) const;
+
+        //- Return area of intersection of faceA with faceB
+        bool overlaps
         (
             const face& faceA,
             const face& faceB,
             const vector& n,
-            const triangulationMode& triMode
+            const scalar threshold
-        );
+        ) const;
 };
 
 

src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,6 +25,29 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+inline void Foam::faceAreaIntersect::triangleFan
+(
+    const face& f,
+    DynamicList<face>& faces
+)
+{
+    if (f.size() > 2)
+    {
+        const label v0 = 0;
+
+        face indices(3);
+
+        for (label i = 1; i < f.size() - 1; ++i)
+        {
+            indices[0] = f[v0];
+            indices[1] = f[i];
+            indices[2] = f[i + 1];
+            faces.append(indices);
+        }
+    }
+}
+
+
 inline void Foam::faceAreaIntersect::setTriPoints
 (
     const point& a,
@@ -113,4 +136,5 @@ Foam::DynamicList<Foam::triPoints>& Foam::faceAreaIntersect::triangles()
     return triangles_;
 }
 
+
 // ************************************************************************* //

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C

@@ -526,17 +526,17 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     AMIPtr_(nullptr),
     AMIMethod_
     (
-        AMIPatchToPatchInterpolation::wordTointerpolationMethod
+        AMIPatchToPatchInterpolation::interpolationMethodNames_
-        (
+        [
             dict.lookupOrDefault
             (
                 "method",
-                AMIPatchToPatchInterpolation::interpolationMethodToWord
+                AMIPatchToPatchInterpolation::interpolationMethodNames_
-                (
+                [
                     AMIPatchToPatchInterpolation::imFaceAreaWeight
+                ]
             )
-            )
+        ]
-        )
     ),
     AMIReverse_(dict.lookupOrDefault("flipNormals", false)),
     AMIRequireMatch_(true),
@@ -1096,10 +1096,10 @@ void Foam::cyclicAMIPolyPatch::write(Ostream& os) const
         os.writeEntry
         (
             "method",
-            AMIPatchToPatchInterpolation::interpolationMethodToWord
+            AMIPatchToPatchInterpolation::interpolationMethodNames_
-            (
+            [
                 AMIMethod_
-            )
+            ]
         );
     }
 

src/sampling/meshToMesh/meshToMesh.C

@@ -850,10 +850,10 @@ Foam::meshToMesh::meshToMesh
     constructNoCuttingPatches
     (
         interpolationMethodNames_[method],
-        AMIPatchToPatchInterpolation::interpolationMethodToWord
+        AMIPatchToPatchInterpolation::interpolationMethodNames_
-        (
+        [
             interpolationMethodAMI(method)
-        ),
+        ],
         interpAllPatches
     );
 }
@@ -917,10 +917,10 @@ Foam::meshToMesh::meshToMesh
     constructFromCuttingPatches
     (
         interpolationMethodNames_[method],
-        AMIPatchToPatchInterpolation::interpolationMethodToWord
+        AMIPatchToPatchInterpolation::interpolationMethodNames_
-        (
+        [
             interpolationMethodAMI(method)
-        ),
+        ],
         patchMap,
         cuttingPatches,
         normalise