From 145029c67f9ec5c6fa3d703bf7a5ca579958503d Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Wed, 2 Jul 2014 16:11:35 +0100
Subject: [PATCH] ENH: directAMI - migrating new functionlity from internal
 development line

---
 .../AMIMethod/directAMI/directAMI.C           | 123 ++++++++++++++++--
 .../AMIMethod/directAMI/directAMI.H           |  12 +-
 2 files changed, 123 insertions(+), 12 deletions(-)

diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C
index 70673d31d9..b63e0f7ae0 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C
@@ -30,7 +30,7 @@ License
 template<class SourcePatch, class TargetPatch>
 void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
 (
-    boolList& mapFlag,
+    labelList& mapFlag,
     labelList& srcTgtSeed,
     DynamicList<label>& srcSeeds,
     DynamicList<label>& nonOverlapFaces,
@@ -50,19 +50,32 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
     {
         label srcI = srcNbr[i];
 
-        if (mapFlag[srcI] && srcTgtSeed[srcI] == -1)
+        if ((mapFlag[srcI] == 0) && (srcTgtSeed[srcI] == -1))
         {
+            // first attempt: match by comparing face centres
             const face& srcF = this->srcPatch_[srcI];
-            const vector srcN = srcF.normal(srcPoints);
+            const point& srcC = srcCf[srcI];
+
+            scalar tol = GREAT;
+            forAll(srcF, fpI)
+            {
+                const point& p = srcPoints[srcF[fpI]];
+                scalar d2 = magSqr(p - srcC);
+                if (d2 < tol)
+                {
+                    tol = d2;
+                }
+            }
+            tol = max(SMALL, 0.0001*sqrt(tol));
 
             bool found = false;
             forAll(tgtNbr, j)
             {
                 label tgtI = tgtNbr[j];
                 const face& tgtF = this->tgtPatch_[tgtI];
-                pointHit ray = tgtF.ray(srcCf[srcI], srcN, tgtPoints);
+                const point tgtC = tgtF.centre(tgtPoints);
 
-                if (ray.hit())
+                if (mag(srcC - tgtC) < tol)
                 {
                     // new match - append to lists
                     found = true;
@@ -74,11 +87,57 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
                 }
             }
 
+            // second attempt: match by shooting a ray into the tgt face
             if (!found)
             {
-                // no match available for source face srcI
-                mapFlag[srcI] = false;
+                const vector srcN = srcF.normal(srcPoints);
+
+                forAll(tgtNbr, j)
+                {
+                    label tgtI = tgtNbr[j];
+                    const face& tgtF = this->tgtPatch_[tgtI];
+                    pointHit ray = tgtF.ray(srcCf[srcI], srcN, tgtPoints);
+
+                    if (ray.hit())
+                    {
+                        // new match - append to lists
+                        found = true;
+
+                        srcTgtSeed[srcI] = tgtI;
+                        srcSeeds.append(srcI);
+
+                        break;
+                    }
+                }
+            }
+
+            // no match available for source face srcI
+            if (!found)
+            {
+                mapFlag[srcI] = -1;
                 nonOverlapFaces.append(srcI);
+
+                if (debug)
+                {
+                    Pout<< "source face not found: id=" << srcI
+                        << " centre=" << srcCf[srcI]
+                        << " normal=" << srcF.normal(srcPoints)
+                        << " points=" << srcF.points(srcPoints)
+                        << endl;
+
+                    Pout<< "target neighbours:" << nl;
+                    forAll(tgtNbr, j)
+                    {
+                        label tgtI = tgtNbr[j];
+                        const face& tgtF = this->tgtPatch_[tgtI];
+
+                        Pout<< "face id: " << tgtI
+                            << " centre=" << tgtF.centre(tgtPoints)
+                            << " normal=" << tgtF.normal(tgtPoints)
+                            << " points=" << tgtF.points(tgtPoints)
+                            << endl;
+                    }
+                }
             }
         }
     }
@@ -96,6 +155,36 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
 }
 
 
+template<class SourcePatch, class TargetPatch>
+void Foam::directAMI<SourcePatch, TargetPatch>::restartAdvancingFront
+(
+    labelList& mapFlag,
+    DynamicList<label>& nonOverlapFaces,
+    label& srcFaceI,
+    label& tgtFaceI
+) const
+{
+    forAll(mapFlag, faceI)
+    {
+        if (mapFlag[faceI] == 0)
+        {
+            tgtFaceI = this->findTargetFace(faceI);
+
+            if (tgtFaceI < 0)
+            {
+                mapFlag[faceI] = -1;
+                nonOverlapFaces.append(faceI);
+            }
+            else
+            {
+                srcFaceI = faceI;
+                break;
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class SourcePatch, class TargetPatch>
@@ -176,15 +265,19 @@ void Foam::directAMI<SourcePatch, TargetPatch>::calculate
     srcTgtSeed[srcFaceI] = tgtFaceI;
 
     // list to keep track of whether src face can be mapped
-    boolList mapFlag(srcAddr.size(), true);
+    // 1 = mapped, 0 = untested, -1 = cannot map
+    labelList mapFlag(srcAddr.size(), 0);
 
+    label nTested = 0;
     DynamicList<label> nonOverlapFaces;
     do
     {
         srcAddr[srcFaceI].append(tgtFaceI);
         tgtAddr[tgtFaceI].append(srcFaceI);
 
-        mapFlag[srcFaceI] = false;
+        mapFlag[srcFaceI] = 1;
+
+        nTested++;
 
         // Do advancing front starting from srcFaceI, tgtFaceI
         appendToDirectSeeds
@@ -196,6 +289,12 @@ void Foam::directAMI<SourcePatch, TargetPatch>::calculate
             srcFaceI,
             tgtFaceI
         );
+
+        if (srcFaceI < 0 && nTested < this->srcPatch_.size())
+        {
+            restartAdvancingFront(mapFlag, nonOverlapFaces, srcFaceI, tgtFaceI);
+        }
+
     } while (srcFaceI >= 0);
 
     if (nonOverlapFaces.size() != 0)
@@ -211,14 +310,16 @@ void Foam::directAMI<SourcePatch, TargetPatch>::calculate
     forAll(srcAddr, i)
     {
         scalar magSf = this->srcMagSf_[i];
-        srcAddress[i].transfer(srcAddr[i]);
+//        srcWeights[i] = scalarList(srcAddr[i].size(), magSf);
         srcWeights[i] = scalarList(1, magSf);
+        srcAddress[i].transfer(srcAddr[i]);
     }
     forAll(tgtAddr, i)
     {
         scalar magSf = this->tgtMagSf_[i];
-        tgtAddress[i].transfer(tgtAddr[i]);
+//        tgtWeights[i] = scalarList(tgtAddr[i].size(), magSf);
         tgtWeights[i] = scalarList(1, magSf);
+        tgtAddress[i].transfer(tgtAddr[i]);
     }
 }
 
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H
index dc273c792d..d74094d0b5 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H
@@ -67,7 +67,7 @@ private:
             //- Append to list of src face seed indices
             void appendToDirectSeeds
             (
-                boolList& mapFlag,
+                labelList& mapFlag,
                 labelList& srcTgtSeed,
                 DynamicList<label>& srcSeeds,
                 DynamicList<label>& nonOverlapFaces,
@@ -75,6 +75,16 @@ private:
                 label& tgtFaceI
             ) const;
 
+            //- Restart the advancing front - typically happens for
+            //  disconnected regions
+            void restartAdvancingFront
+            (
+                labelList& mapFlag,
+                DynamicList<label>& nonOverlapFaces,
+                label& srcFaceI,
+                label& tgtFaceI
+            ) const;
+
 
         // Evaluation