From a29da9dc885b0e59d08720ee928642567d5c234a Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Thu, 1 Sep 2011 15:45:20 +0100
Subject: [PATCH] ENH: Updated AMI functionality - re: projection surface

---
 src/AMIInterpolation/AMIInterpolation.C       | 156 +++++++++---------
 src/AMIInterpolation/AMIInterpolation.H       |  14 +-
 .../cyclicAMIPolyPatch/cyclicAMIPolyPatch.C   | 105 +++++++-----
 .../cyclicAMIPolyPatch/cyclicAMIPolyPatch.H   |   5 +-
 4 files changed, 143 insertions(+), 137 deletions(-)

diff --git a/src/AMIInterpolation/AMIInterpolation.C b/src/AMIInterpolation/AMIInterpolation.C
index c116ee2030f..7ffbb3e90a5 100644
--- a/src/AMIInterpolation/AMIInterpolation.C
+++ b/src/AMIInterpolation/AMIInterpolation.C
@@ -499,11 +499,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::projectPointsToSurface
     pointField& pts
 ) const
 {
-    if (!projectPoints_)
-    {
-        return;
-    }
-
     if (debug)
     {
         Info<< "AMI: projecting points to surface" << endl;
@@ -803,7 +798,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
 (
     const primitivePatch& srcPatch,
     const primitivePatch& tgtPatch,
-    const searchableSurface& surf,
     label srcFaceI,
     label tgtFaceI
 )
@@ -824,6 +818,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
     if (debug)
     {
         Info<< "AMI: calcAddressing" << endl;
+        writePatch(srcPatch, "VTK", "source");
+        writePatch(tgtPatch, "VTK", "target");
     }
 
     // temporary storage for addressing and weights
@@ -833,57 +829,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
     List<DynamicList<scalar> > tgtWght(tgtPatch.size());
 
 
-    // create new patches for source and target
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    pointField srcPoints = srcPatch.points();
-    primitivePatch srcPatch0
-    (
-        SubList<face>
-        (
-            srcPatch,
-            srcPatch.size(),
-            0
-        ),
-        srcPoints
-    );
-
-    if (debug)
-    {
-        OFstream os("amiSrcPoints.obj");
-        forAll(srcPoints, i)
-        {
-            meshTools::writeOBJ(os, srcPoints[i]);
-        }
-    }
-
-    pointField tgtPoints = tgtPatch.points();
-    primitivePatch tgtPatch0
-    (
-        SubList<face>
-        (
-            tgtPatch,
-            tgtPatch.size(),
-            0
-        ),
-        tgtPoints
-    );
-
-    if (debug)
-    {
-        OFstream os("amiTgtPoints.obj");
-        forAll(tgtPoints, i)
-        {
-            meshTools::writeOBJ(os, tgtPoints[i]);
-        }
-    }
-
-
-    // map source and target patches onto projection surface
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    projectPointsToSurface(surf, srcPoints);
-    projectPointsToSurface(surf, tgtPoints);
-
-
     // find initial face match using brute force/octree search
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     if ((srcFaceI == -1) || (tgtFaceI == -1))
@@ -891,9 +836,9 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
         srcFaceI = 0;
         tgtFaceI = 0;
         bool foundFace = false;
-        forAll(srcPatch0, faceI)
+        forAll(srcPatch, faceI)
         {
-            tgtFaceI = findTargetFace(srcFaceI, srcPatch0, tgtPatch0);
+            tgtFaceI = findTargetFace(faceI, srcPatch, tgtPatch);
             if (tgtFaceI >= 0)
             {
                 srcFaceI = faceI;
@@ -911,7 +856,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
                 "("
                     "const primitivePatch&, "
                     "const primitivePatch&, "
-                    "const searchableSurface&, "
                     "label, "
                     "label"
                 ")"
@@ -927,7 +871,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
 
     // construct weights and addressing
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    label nFacesRemaining = srcPatch0.size();
+    label nFacesRemaining = srcPatch.size();
 
     // list of tgt face neighbour faces
     DynamicList<label> nbrFaces(10);
@@ -953,7 +897,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
 
         // append initial target face and neighbours
         nbrFaces.append(tgtFaceI);
-        appendNbrFaces(tgtFaceI, tgtPatch0, visitedFaces, nbrFaces);
+        appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
 
         bool faceProcessed = false;
 
@@ -962,7 +906,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
             // process new target face
             tgtFaceI = nbrFaces.remove();
             visitedFaces.append(tgtFaceI);
-            scalar area = interArea(srcFaceI, tgtFaceI, srcPatch0, tgtPatch0);
+            scalar area = interArea(srcFaceI, tgtFaceI, srcPatch, tgtPatch);
 
             // store when intersection area > 0
             if (area > 0)
@@ -973,7 +917,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
                 tgtAddr[tgtFaceI].append(srcFaceI);
                 tgtWght[tgtFaceI].append(area);
 
-                appendNbrFaces(tgtFaceI, tgtPatch0, visitedFaces, nbrFaces);
+                appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
 
                 faceProcessed = true;
             }
@@ -996,8 +940,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
             (
                 srcFaceI,
                 tgtFaceI,
-                srcPatch0,
-                tgtPatch0,
+                srcPatch,
+                tgtPatch,
                 mapFlag,
                 seedFaces,
                 visitedFaces
@@ -1084,9 +1028,8 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
 (
     const SourcePatch& srcPatch,
     const TargetPatch& tgtPatch,
-    const searchableSurface& surf,
-    const faceAreaIntersect::triangulationMode& triMode,
-    const bool projectPoints
+    const autoPtr<searchableSurface>& surfPtr,
+    const faceAreaIntersect::triangulationMode& triMode
 )
 :
     srcAddress_(),
@@ -1095,24 +1038,74 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
     tgtWeights_(),
     startSeedI_(0),
     triMode_(triMode),
-    projectPoints_(projectPoints),
     srcMapPtr_(NULL),
     tgtMapPtr_(NULL)
 {
     label srcSize = returnReduce(srcPatch.size(), sumOp<label>());
     label tgtSize = returnReduce(tgtPatch.size(), sumOp<label>());
 
-    Info<< "AMI: Creating interpolation addressing and weights" << nl
-        << "    Source faces = " << srcSize << nl
-        << "    Target faces = " << tgtSize << endl;
+    Info<< "AMI: Creating addressing and weights between "
+        << srcSize << " source faces and " << tgtSize << " target faces"
+        << endl;
 
-    if (debug)
+    if (surfPtr.valid())
     {
-        writePatch(srcPatch, "VTK", "source");
-        writePatch(tgtPatch, "VTK", "target");
-    }
+        // create new patches for source and target
+        pointField srcPoints = srcPatch.points();
+        primitivePatch srcPatch0
+        (
+            SubList<face>
+            (
+                srcPatch,
+                srcPatch.size(),
+                0
+            ),
+            srcPoints
+        );
+
+        if (debug)
+        {
+            OFstream os("amiSrcPoints.obj");
+            forAll(srcPoints, i)
+            {
+                meshTools::writeOBJ(os, srcPoints[i]);
+            }
+        }
+
+        pointField tgtPoints = tgtPatch.points();
+        primitivePatch tgtPatch0
+        (
+            SubList<face>
+            (
+                tgtPatch,
+                tgtPatch.size(),
+                0
+            ),
+            tgtPoints
+        );
 
-    update(srcPatch, tgtPatch, surf);
+        if (debug)
+        {
+            OFstream os("amiTgtPoints.obj");
+            forAll(tgtPoints, i)
+            {
+                meshTools::writeOBJ(os, tgtPoints[i]);
+            }
+        }
+
+
+        // map source and target patches onto projection surface
+        projectPointsToSurface(surfPtr(), srcPoints);
+        projectPointsToSurface(surfPtr(), tgtPoints);
+
+
+        // calculate AMI interpolation
+        update(srcPatch0, tgtPatch0);
+    }
+    else
+    {
+        update(srcPatch, tgtPatch);
+    }
 }
 
 
@@ -1129,8 +1122,7 @@ template<class SourcePatch, class TargetPatch>
 void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
 (
     const primitivePatch& srcPatch,
-    const primitivePatch& tgtPatch,
-    const searchableSurface& surf
+    const primitivePatch& tgtPatch
 )
 {
     static label patchI = 0;
@@ -1180,7 +1172,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
 
 
         // calculate AMI interpolation
-        calcAddressing(srcPatch, newTgtPatch, surf);
+        calcAddressing(srcPatch, newTgtPatch);
 
         // Now
         // ~~~
@@ -1260,7 +1252,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
     {
         checkPatches(srcPatch, tgtPatch);
 
-        calcAddressing(srcPatch, tgtPatch, surf);
+        calcAddressing(srcPatch, tgtPatch);
 
         if (debug)
         {
diff --git a/src/AMIInterpolation/AMIInterpolation.H b/src/AMIInterpolation/AMIInterpolation.H
index 4173362f542..2142114661a 100644
--- a/src/AMIInterpolation/AMIInterpolation.H
+++ b/src/AMIInterpolation/AMIInterpolation.H
@@ -129,11 +129,10 @@ class AMIInterpolation
         //- Face triangulation mode
         const faceAreaIntersect::triangulationMode triMode_;
 
-        //- Flag to indicate whether or not to project points
-        const bool projectPoints_;
-
+        //- Source map pointer - parallel running only
         autoPtr<mapDistribute> srcMapPtr_;
 
+        //- Target map pointer - parallel running only
         autoPtr<mapDistribute> tgtMapPtr_;
 
 
@@ -260,7 +259,6 @@ class AMIInterpolation
             (
                 const primitivePatch& srcPatch,
                 const primitivePatch& tgtPatch,
-                const searchableSurface& surf,
                 label srcFaceI = -1,
                 label tgtFaceI = -1
             );
@@ -289,9 +287,8 @@ public:
         (
             const SourcePatch& srcPatch,
             const TargetPatch& tgtPatch,
-            const searchableSurface& surf,
-            const faceAreaIntersect::triangulationMode& triMode,
-            const bool projectPoints
+            const autoPtr<searchableSurface>& surf,
+            const faceAreaIntersect::triangulationMode& triMode
         );
 
 
@@ -327,8 +324,7 @@ public:
             void update
             (
                 const primitivePatch& srcPatch,
-                const primitivePatch& tgtPatch,
-                const searchableSurface& surf
+                const primitivePatch& tgtPatch
             );
 
 
diff --git a/src/AMIInterpolation/patches/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C b/src/AMIInterpolation/patches/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
index 8c36c919079..58cd9d8e2a6 100644
--- a/src/AMIInterpolation/patches/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
+++ b/src/AMIInterpolation/patches/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
@@ -221,36 +221,11 @@ void Foam::cyclicAMIPolyPatch::calcTransforms
 }
 
 
-const Foam::searchableSurface& Foam::cyclicAMIPolyPatch::surf()
+void Foam::cyclicAMIPolyPatch::resetAMI()
 {
-    word surfType(surfDict_.lookup("type"));
-    word surfName(surfDict_.lookupOrDefault("name", name()));
-
-    const polyMesh& mesh = boundaryMesh().mesh();
-
-    surfPtr_ =
-        searchableSurface::New
-        (
-            surfType,
-            IOobject
-            (
-                surfName,
-                mesh.time().constant(),
-                "triSurface",
-                mesh,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE
-            ),
-            surfDict_
-        );
+    AMIPtr_.clear();
 
-    return surfPtr_();
-}
-
-
-const Foam::AMIPatchToPatchInterpolation& Foam::cyclicAMIPolyPatch::AMI()
-{
-    if (!AMIPtr_.valid() && owner())
+    if (owner())
     {
         const polyPatch& nbr = nbrPatch();
         pointField nbrPoints = nbrPatch().localPoints();
@@ -282,8 +257,6 @@ const Foam::AMIPatchToPatchInterpolation& Foam::cyclicAMIPolyPatch::AMI()
             meshTools::writeOBJ(osO, this->localFaces(), localPoints());
         }
 
-        bool projectPoints(readBool(surfDict_.lookup("projectPoints")));
-
         // Construct/apply AMI interpolation to determine addressing and weights
         AMIPtr_.reset
         (
@@ -291,14 +264,11 @@ const Foam::AMIPatchToPatchInterpolation& Foam::cyclicAMIPolyPatch::AMI()
             (
                 *this,
                 nbrPatch0,
-                surf(),
-                faceAreaIntersect::tmMesh,
-                projectPoints
+                surfPtr_,
+                faceAreaIntersect::tmMesh
             )
         );
     }
-
-    return AMIPtr_();
 }
 
 
@@ -345,11 +315,7 @@ void Foam::cyclicAMIPolyPatch::movePoints
 
     calcTransforms();
 
-    if (owner())
-    {
-        AMIPtr_.clear();
-        AMI();
-    }
+    resetAMI();
 }
 
 
@@ -611,6 +577,59 @@ bool Foam::cyclicAMIPolyPatch::owner() const
 }
 
 
+const Foam::autoPtr<Foam::searchableSurface>&
+Foam::cyclicAMIPolyPatch::surfPtr()
+{
+    word surfType(surfDict_.lookup("type"));
+
+    if (!surfPtr_.valid() && owner() && surfType != "none")
+    {
+        word surfName(surfDict_.lookupOrDefault("name", name()));
+
+        const polyMesh& mesh = boundaryMesh().mesh();
+
+        surfPtr_ =
+            searchableSurface::New
+            (
+                surfType,
+                IOobject
+                (
+                    surfName,
+                    mesh.time().constant(),
+                    "triSurface",
+                    mesh,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                ),
+                surfDict_
+            );
+    }
+
+    return surfPtr_;
+}
+
+
+const Foam::AMIPatchToPatchInterpolation& Foam::cyclicAMIPolyPatch::AMI()
+{
+    if (!owner())
+    {
+        FatalErrorIn
+        (
+            "const AMIPatchToPatchInterpolation& cyclicAMIPolyPatch::AMI()"
+        )
+            << "AMI interpolator only available to owner patch"
+            << abort(FatalError);
+    }
+
+    if (!AMIPtr_.valid())
+    {
+        resetAMI();
+    }
+
+    return AMIPtr_();
+}
+
+
 void Foam::cyclicAMIPolyPatch::transformPosition(pointField& l) const
 {
     if (!parallel())
@@ -704,11 +723,7 @@ void Foam::cyclicAMIPolyPatch::calcGeometry
         nbrAreas
     );
 
-    if (owner())
-    {
-        AMIPtr_.clear();
-        AMI();
-    }
+    resetAMI();
 }
 
 
diff --git a/src/AMIInterpolation/patches/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H b/src/AMIInterpolation/patches/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H
index 1ad095036b0..352bbddb54a 100644
--- a/src/AMIInterpolation/patches/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H
+++ b/src/AMIInterpolation/patches/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H
@@ -109,6 +109,9 @@ private:
             const vectorField& half1Areas
         );
 
+        //- Reset the AMI interpolator
+        void resetAMI();
+
 
 protected:
 
@@ -264,7 +267,7 @@ public:
             inline const cyclicAMIPolyPatch& nbrPatch() const;
 
             //- Return a reference to the projection surface
-            const searchableSurface& surf();
+            const autoPtr<searchableSurface>& surfPtr();
 
             //- Return a reference to the AMI interpolator
             const AMIPatchToPatchInterpolation& AMI();
-- 
GitLab