From 83b357ab6ef8ccba520549d798dba1cd60ff0971 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 25 Jul 2012 12:09:09 +0100
Subject: [PATCH] ENH: AMIInterpolation: use stable guess for staritng seed.
 Add restart (under debug flag)

---
 .../AMIInterpolation/AMIInterpolation.C       | 241 +++++++++++++++---
 .../AMIInterpolation/AMIInterpolation.H       |  25 ++
 .../faceAreaIntersect/faceAreaIntersect.C     |  10 +-
 .../faceAreaIntersect/faceAreaIntersect.H     |  10 +-
 .../faceAreaIntersect/faceAreaIntersectI.H    |  12 +-
 5 files changed, 260 insertions(+), 38 deletions(-)

diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
index 9db5dc0e12..0b0e3800f7 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
@@ -782,6 +782,62 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::appendNbrFaces
 }
 
 
+template<class SourcePatch, class TargetPatch>
+bool Foam::AMIInterpolation<SourcePatch, TargetPatch>::processSourceFace
+(
+    const SourcePatch& srcPatch,
+    const TargetPatch& tgtPatch,
+    const label srcFaceI,
+    const label tgtStartFaceI,
+
+    // list of tgt face neighbour faces
+    DynamicList<label>& nbrFaces,
+    // list of faces currently visited for srcFaceI to avoid multiple hits
+    DynamicList<label>& visitedFaces,
+
+    // temporary storage for addressing and weights
+    List<DynamicList<label> >& srcAddr,
+    List<DynamicList<scalar> >& srcWght,
+    List<DynamicList<label> >& tgtAddr,
+    List<DynamicList<scalar> >& tgtWght
+)
+{
+    nbrFaces.clear();
+    visitedFaces.clear();
+
+    // append initial target face and neighbours
+    nbrFaces.append(tgtStartFaceI);
+    appendNbrFaces(tgtStartFaceI, tgtPatch, visitedFaces, nbrFaces);
+
+    bool faceProcessed = false;
+
+    do
+    {
+        // process new target face
+        label tgtFaceI = nbrFaces.remove();
+        visitedFaces.append(tgtFaceI);
+        scalar area = interArea(srcFaceI, tgtFaceI, srcPatch, tgtPatch);
+
+        // store when intersection area > 0
+        if (area > 0)
+        {
+            srcAddr[srcFaceI].append(tgtFaceI);
+            srcWght[srcFaceI].append(area);
+
+            tgtAddr[tgtFaceI].append(srcFaceI);
+            tgtWght[tgtFaceI].append(area);
+
+            appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
+
+            faceProcessed = true;
+        }
+
+    } while (nbrFaces.size() > 0);
+
+    return faceProcessed;
+}
+
+
 template<class SourcePatch, class TargetPatch>
 void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
 (
@@ -811,7 +867,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
                 label faceT = visitedFaces[j];
                 scalar area = interArea(faceS, faceT, srcPatch0, tgtPatch0);
 
-                if (area > 0)
+                // Check that faces have enough overlap for robust walking
+                if (area/srcMagSf_[srcFaceI] > faceAreaIntersect::tolerance())
                 {
                     // TODO - throwing area away - re-use in next iteration?
 
@@ -864,7 +921,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
                 << "target face" << endl;
         }
 
-//        foundNextSeed = false;
+        foundNextSeed = false;
         for (label faceI = startSeedI; faceI < mapFlag.size(); faceI++)
         {
             if (mapFlag[faceI])
@@ -887,7 +944,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
 
         FatalErrorIn
         (
-            "void Foam::cyclicAMIPolyPatch::setNextFaces"
+            "void Foam::AMIInterpolation<SourcePatch, TargetPatch>::"
+            "setNextFaces"
             "("
                 "label&, "
                 "label&, "
@@ -946,6 +1004,25 @@ Foam::scalar Foam::AMIInterpolation<SourcePatch, TargetPatch>::interArea
     {
         area = inter.calc(src, tgt, n, triMode_);
     }
+    else
+    {
+        WarningIn
+        (
+            "void Foam::AMIInterpolation<SourcePatch, TargetPatch>::"
+            "interArea"
+            "("
+                "const label, "
+                "const label, "
+                "const SourcePatch&, "
+                "const TargetPatch&"
+            ") const"
+        )   << "Invalid normal for source face " << srcFaceI
+            << " points " << UIndirectList<point>(srcPoints, src)
+            << " target face " << tgtFaceI
+            << " points " << UIndirectList<point>(tgtPoints, tgt)
+            << endl;
+    }
+
 
     if ((debug > 1) && (area > 0))
     {
@@ -956,6 +1033,105 @@ Foam::scalar Foam::AMIInterpolation<SourcePatch, TargetPatch>::interArea
 }
 
 
+template<class SourcePatch, class TargetPatch>
+void Foam::AMIInterpolation<SourcePatch, TargetPatch>::
+restartUncoveredSourceFace
+(
+    const SourcePatch& srcPatch,
+    const TargetPatch& tgtPatch,
+    List<DynamicList<label> >& srcAddr,
+    List<DynamicList<scalar> >& srcWght,
+    List<DynamicList<label> >& tgtAddr,
+    List<DynamicList<scalar> >& tgtWght
+)
+{
+    // Collect all src faces with a low weight
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelHashSet lowWeightFaces(100);
+    forAll(srcWght, srcFaceI)
+    {
+        scalar s = sum(srcWght[srcFaceI]);
+        scalar t = s/srcMagSf_[srcFaceI];
+
+        if (t < 0.5)
+        {
+            lowWeightFaces.insert(srcFaceI);
+        }
+    }
+
+    Info<< "AMIInterpolation : restarting search on "
+        << returnReduce(lowWeightFaces.size(), sumOp<label>())
+        << " faces since sum of weights < 0.5" << endl;
+
+    if (lowWeightFaces.size() > 0)
+    {
+        // Erase all the lowWeight source faces from the target
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        DynamicList<label> okSrcFaces(10);
+        DynamicList<scalar> okSrcWeights(10);
+        forAll(tgtAddr, tgtFaceI)
+        {
+            okSrcFaces.clear();
+            okSrcWeights.clear();
+            DynamicList<label>& srcFaces = tgtAddr[tgtFaceI];
+            DynamicList<scalar>& srcWeights = tgtWght[tgtFaceI];
+            forAll(srcFaces, i)
+            {
+                if (!lowWeightFaces.found(srcFaces[i]))
+                {
+                    okSrcFaces.append(srcFaces[i]);
+                    okSrcWeights.append(srcWeights[i]);
+                }
+            }
+            if (okSrcFaces.size() < srcFaces.size())
+            {
+                srcFaces.transfer(okSrcFaces);
+                srcWeights.transfer(okSrcWeights);
+            }
+        }
+
+
+
+        // Restart search from best hit
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        // list of tgt face neighbour faces
+        DynamicList<label> nbrFaces(10);
+
+        // list of faces currently visited for srcFaceI to avoid multiple hits
+        DynamicList<label> visitedFaces(10);
+
+        forAllConstIter(labelHashSet, lowWeightFaces, iter)
+        {
+            label srcFaceI = iter.key();
+            label tgtFaceI = findTargetFace(srcFaceI, srcPatch);
+            if (tgtFaceI != -1)
+            {
+                //bool faceProcessed =
+                processSourceFace
+                (
+                    srcPatch,
+                    tgtPatch,
+                    srcFaceI,
+                    tgtFaceI,
+
+                    nbrFaces,
+                    visitedFaces,
+
+                    srcAddr,
+                    srcWght,
+                    tgtAddr,
+                    tgtWght
+                );
+                // ? Check faceProcessed to see if restarting has worked.
+            }
+        }
+    }
+}
+
+
 template<class SourcePatch, class TargetPatch>
 void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
 (
@@ -1073,37 +1249,22 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
     DynamicList<label> nonOverlapFaces;
     do
     {
-        nbrFaces.clear();
-        visitedFaces.clear();
+        // Do advancing front starting from srcFaceI,tgtFaceI
+        bool faceProcessed = processSourceFace
+        (
+            srcPatch,
+            tgtPatch,
+            srcFaceI,
+            tgtFaceI,
 
-        // append initial target face and neighbours
-        nbrFaces.append(tgtFaceI);
-        appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
+            nbrFaces,
+            visitedFaces,
 
-        bool faceProcessed = false;
-
-        do
-        {
-            // process new target face
-            tgtFaceI = nbrFaces.remove();
-            visitedFaces.append(tgtFaceI);
-            scalar area = interArea(srcFaceI, tgtFaceI, srcPatch, tgtPatch);
-
-            // store when intersection area > 0
-            if (area > 0)
-            {
-                srcAddr[srcFaceI].append(tgtFaceI);
-                srcWght[srcFaceI].append(area);
-
-                tgtAddr[tgtFaceI].append(srcFaceI);
-                tgtWght[tgtFaceI].append(area);
-
-                appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
-
-                faceProcessed = true;
-            }
-
-        } while (nbrFaces.size() > 0);
+            srcAddr,
+            srcWght,
+            tgtAddr,
+            tgtWght
+        );
 
         mapFlag[srcFaceI] = false;
 
@@ -1140,6 +1301,22 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
         srcNonOverlap_.transfer(nonOverlapFaces);
     }
 
+
+    // Check for badly covered faces
+    if (debug)
+    {
+        restartUncoveredSourceFace
+        (
+            srcPatch,
+            tgtPatch,
+            srcAddr,
+            srcWght,
+            tgtAddr,
+            tgtWght
+        );
+    }
+
+
     // transfer data to persistent storage
     forAll(srcAddr, i)
     {
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
index b393030823..5b91afc38a 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
@@ -246,6 +246,31 @@ class AMIInterpolation
                 DynamicList<label>& faceIDs
             ) const;
 
+            bool processSourceFace
+            (
+                const SourcePatch& srcPatch,
+                const TargetPatch& tgtPatch,
+                const label srcFaceI,
+                const label tgtStartFaceI,
+
+                DynamicList<label>& nbrFaces,
+                DynamicList<label>& visitedFaces,
+                List<DynamicList<label> >& srcAddr,
+                List<DynamicList<scalar> >& srcWght,
+                List<DynamicList<label> >& tgtAddr,
+                List<DynamicList<scalar> >& tgtWght
+            );
+
+            void restartUncoveredSourceFace
+            (
+                const SourcePatch& srcPatch,
+                const TargetPatch& tgtPatch,
+                List<DynamicList<label> >& srcAddr,
+                List<DynamicList<scalar> >& srcWght,
+                List<DynamicList<label> >& tgtAddr,
+                List<DynamicList<scalar> >& tgtWght
+            );
+
             //- Set the source and target seed faces
             void setNextFaces
             (
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
index db8a0ba294..4135cfe041 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,10 @@ License
 
 #include "faceAreaIntersect.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+Foam::scalar Foam::faceAreaIntersect::tol = 1e-6;
+
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 void Foam::faceAreaIntersect::triSliceWithPlane
@@ -36,8 +40,6 @@ void Foam::faceAreaIntersect::triSliceWithPlane
     const scalar len
 )
 {
-    const scalar matchTol = 1e-6;
-
     // distance to cutting plane
     FixedList<scalar, 3> d;
 
@@ -51,7 +53,7 @@ void Foam::faceAreaIntersect::triSliceWithPlane
     {
         d[i] = ((tri[i] - p.refPoint()) & p.normal());
 
-        if (mag(d[i]) < matchTol*len)
+        if (mag(d[i]) < tol*len)
         {
             nCoPlanar++;
             copI = i;
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
index 2a27e30960..bfbe7c08a2 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,6 +78,11 @@ private:
         const bool reverseB_;
 
 
+    // Static data members
+
+        static scalar tol;
+
+
     // Private Member Functions
 
         //- Get triPoints from face
@@ -152,6 +157,9 @@ public:
 
     // Public Member Functions
 
+        //- Fraction of local length scale to use as intersection tolerance
+        inline static scalar& tolerance();
+
         //- Return area of intersection of faceA with faceB
         scalar calc
         (
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
index 0b8265aad6..87fab781e0 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,6 +23,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
 inline void Foam::faceAreaIntersect::setTriPoints
 (
     const point& a,
@@ -106,4 +108,12 @@ inline Foam::scalar Foam::faceAreaIntersect::triArea(const triPoints& t) const
 }
 
 
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::scalar& Foam::faceAreaIntersect::tolerance()
+{
+    return tol;
+}
+
+
 // ************************************************************************* //