From 736b551ae2546c42cf31d36bafcf3d9cb1fbd2dd Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 30 Sep 2014 17:20:27 +0100
Subject: [PATCH] ENH: mappedPatch: survive neg-tet volume cells

In case of point not being found re-search for nearest
---
 .../mappedPolyPatch/mappedPatchBase.C         | 87 +++++++++++++------
 .../mappedPolyPatch/mappedPatchBase.H         | 12 ++-
 2 files changed, 70 insertions(+), 29 deletions(-)

diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
index 6ab8930a42..ea24a1c80e 100644
--- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
+++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
@@ -41,6 +41,7 @@ License
 #include "triPointRef.H"
 #include "syncTools.H"
 #include "treeDataCell.H"
+#include "DynamicField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,14 +53,15 @@ namespace Foam
     const char* Foam::NamedEnum
     <
         Foam::mappedPatchBase::sampleMode,
-        5
+        6
     >::names[] =
     {
         "nearestCell",
         "nearestPatchFace",
         "nearestPatchFaceAMI",
         "nearestPatchPoint",
-        "nearestFace"
+        "nearestFace",
+        "nearestOnlyCell"
     };
 
     template<>
@@ -76,7 +78,7 @@ namespace Foam
 }
 
 
-const Foam::NamedEnum<Foam::mappedPatchBase::sampleMode, 5>
+const Foam::NamedEnum<Foam::mappedPatchBase::sampleMode, 6>
     Foam::mappedPatchBase::sampleModeNames_;
 
 const Foam::NamedEnum<Foam::mappedPatchBase::offsetMode, 3>
@@ -187,6 +189,7 @@ void Foam::mappedPatchBase::collectSamples
 // for samples being found in two processors.
 void Foam::mappedPatchBase::findSamples
 (
+    const sampleMode mode,
     const pointField& samples,
     labelList& sampleProcs,
     labelList& sampleIndices,
@@ -199,7 +202,7 @@ void Foam::mappedPatchBase::findSamples
     // All the info for nearest. Construct to miss
     List<nearInfo> nearest(samples.size());
 
-    switch (mode_)
+    switch (mode)
     {
         case NEARESTCELL:
         {
@@ -210,7 +213,7 @@ void Foam::mappedPatchBase::findSamples
                     "mappedPatchBase::findSamples(const pointField&,"
                     " labelList&, labelList&, pointField&) const"
                 )   << "No need to supply a patch name when in "
-                    << sampleModeNames_[mode_] << " mode." << exit(FatalError);
+                    << sampleModeNames_[mode] << " mode." << exit(FatalError);
             }
 
             //- Note: face-diagonal decomposition
@@ -244,6 +247,36 @@ void Foam::mappedPatchBase::findSamples
             break;
         }
 
+        case NEARESTONLYCELL:
+        {
+            if (samplePatch_.size() && samplePatch_ != "none")
+            {
+                FatalErrorIn
+                (
+                    "mappedPatchBase::findSamples(const pointField&,"
+                    " labelList&, labelList&, pointField&) const"
+                )   << "No need to supply a patch name when in "
+                    << sampleModeNames_[mode] << " mode." << exit(FatalError);
+            }
+
+            //- Note: face-diagonal decomposition
+            const indexedOctree<Foam::treeDataCell>& tree = mesh.cellTree();
+
+            forAll(samples, sampleI)
+            {
+                const point& sample = samples[sampleI];
+
+                nearest[sampleI].first() = tree.findNearest(sample, sqr(GREAT));
+                nearest[sampleI].second().first() = magSqr
+                (
+                    nearest[sampleI].first().hitPoint()
+                   -sample
+                );
+                nearest[sampleI].second().second() = Pstream::myProcNo();
+            }
+            break;
+        }
+
         case NEARESTPATCHFACE:
         {
             Random rndGen(123456);
@@ -398,7 +431,7 @@ void Foam::mappedPatchBase::findSamples
                     "mappedPatchBase::findSamples(const pointField&,"
                     " labelList&, labelList&, pointField&) const"
                 )   << "No need to supply a patch name when in "
-                    << sampleModeNames_[mode_] << " mode." << exit(FatalError);
+                    << sampleModeNames_[mode] << " mode." << exit(FatalError);
             }
 
             //- Note: face-diagonal decomposition
@@ -565,7 +598,7 @@ void Foam::mappedPatchBase::calcMapping() const
     labelList sampleProcs;
     labelList sampleIndices;
     pointField sampleLocations;
-    findSamples(samples, sampleProcs, sampleIndices, sampleLocations);
+    findSamples(mode_, samples, sampleProcs, sampleIndices, sampleLocations);
 
     // Check for samples that were not found. This will only happen for
     // NEARESTCELL since finds cell containing a location
@@ -607,32 +640,36 @@ void Foam::mappedPatchBase::calcMapping() const
                 hasWarned = true;
             }
 
-            // Reset the samples that cannot be found to the cell centres.
-            pointField patchCc;
-            {
-                List<pointField> globalCc(Pstream::nProcs());
-                globalCc[Pstream::myProcNo()] = patch_.faceCellCentres();
-                Pstream::gatherList(globalCc);
-                Pstream::scatterList(globalCc);
-                patchCc = ListListOps::combine<pointField>
-                (
-                    globalCc,
-                    accessOp<pointField>()
-                );
-            }
+            // Collect the samples that cannot be found
+            DynamicList<label> subMap;
+            DynamicField<point> subSamples;
 
             forAll(sampleProcs, sampleI)
             {
                 if (sampleProcs[sampleI] == -1)
                 {
-                    // Reset to cell centres
-                    samples[sampleI] = patchCc[sampleI];
+                    subMap.append(sampleI);
+                    subSamples.append(samples[sampleI]);
                 }
             }
 
-            // And re-search. Note: could be optimised to only search missing
-            // points.
-            findSamples(samples, sampleProcs, sampleIndices, sampleLocations);
+            // And re-search for pure nearest (should not fail)
+            labelList subSampleProcs;
+            labelList subSampleIndices;
+            pointField subSampleLocations;
+            findSamples
+            (
+                NEARESTONLYCELL,
+                subSamples,
+                subSampleProcs,
+                subSampleIndices,
+                subSampleLocations
+            );
+
+            // Insert
+            UIndirectList<label>(sampleProcs, subMap) = subSampleProcs;
+            UIndirectList<label>(sampleIndices, subMap) = subSampleIndices;
+            UIndirectList<point>(sampleLocations, subMap) = subSampleLocations;
         }
     }
 
diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H
index 7dd98c1806..e6fcbb540e 100644
--- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H
+++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,8 @@ Description
 
         // What to sample:
         // - nearestCell         : sample cell containing point
+        // - nearestOnlyCell     : nearest sample cell (even if not containing
+        //                         point)
         // - nearestPatchFace    : nearest face on selected patch
         // - nearestPatchFaceAMI : nearest face on selected patch
                                    - patches need not conform
@@ -107,11 +109,12 @@ public:
         //- Mesh items to sample
         enum sampleMode
         {
-            NEARESTCELL,         // nearest cell
+            NEARESTCELL,         // nearest cell containing sample
             NEARESTPATCHFACE,    // nearest face on selected patch
             NEARESTPATCHFACEAMI, // nearest patch face + AMI interpolation
             NEARESTPATCHPOINT,   // nearest point on selected patch
-            NEARESTFACE          // nearest face
+            NEARESTFACE,         // nearest face
+            NEARESTONLYCELL      // nearest cell (even if not containing cell)
         };
 
         //- How to project face centres
@@ -122,7 +125,7 @@ public:
             NORMAL              // use face normal + distance
         };
 
-        static const NamedEnum<sampleMode, 5> sampleModeNames_;
+        static const NamedEnum<sampleMode, 6> sampleModeNames_;
 
         static const NamedEnum<offsetMode, 3> offsetModeNames_;
 
@@ -255,6 +258,7 @@ protected:
         //- Find cells/faces containing samples
         void findSamples
         (
+            const sampleMode mode,      // search mode
             const pointField&,
             labelList& sampleProcs,     // processor containing sample
             labelList& sampleIndices,   // local index of cell/face