ENH: sample on nearest cell, not cell containing triangle centre.

Much more robust. Also now works on triSurfaces with non-compact point numbering.

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C

@@ -25,11 +25,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sampledTriSurfaceMesh.H"
-#include "dictionary.H"
+#include "treeDataPoint.H"
-#include "polyMesh.H"
-#include "polyPatch.H"
-#include "volFields.H"
 #include "meshSearch.H"
+#include "Tuple2.H"
+#include "globalIndex.H"
 
 #include "addToRunTimeSelectionTable.H"
 
@@ -44,6 +43,25 @@ namespace Foam
         sampledTriSurfaceMesh,
         word
     );
+
+    //- Private class for finding nearest
+    //  - global index
+    //  - sqr(distance)
+    typedef Tuple2<scalar, label> nearInfo;
+
+    class nearestEqOp
+    {
+
+    public:
+
+        void operator()(nearInfo& x, const nearInfo& y) const
+        {
+            if (y.first() < x.first())
+            {
+                x = y;
+            }
+        }
+    };
 }
 
 
@@ -63,16 +81,16 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
         (
             name,
             mesh.time().constant(), // instance
-            "triSurface",       // local
+            "triSurface",           // local
-            mesh,               // registry
+            mesh,                   // registry
             IOobject::MUST_READ,
-            IOobject::NO_WRITE
+            IOobject::NO_WRITE,
+            false
         )
     ),
     needsUpdate_(true),
     cellLabels_(0),
-    pointMap_(0),
+    pointToFace_(0)
-    faceMap_(0)
 {}
 
 
@@ -90,16 +108,16 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
         (
             dict.lookup("surface"),
             mesh.time().constant(), // instance
-            "triSurface",       // local
+            "triSurface",           // local
-            mesh,               // registry
+            mesh,                   // registry
             IOobject::MUST_READ,
-            IOobject::NO_WRITE
+            IOobject::NO_WRITE,
+            false
         )
     ),
     needsUpdate_(true),
     cellLabels_(0),
-    pointMap_(0),
+    pointToFace_(0)
-    faceMap_(0)
 {}
 
 
@@ -128,8 +146,7 @@ bool Foam::sampledTriSurfaceMesh::expire()
     sampledSurface::clearGeom();
     MeshStorage::clear();
     cellLabels_.clear();
-    pointMap_.clear();
+    pointToFace_.clear();
-    faceMap_.clear();
 
     needsUpdate_ = true;
     return true;
@@ -148,93 +165,137 @@ bool Foam::sampledTriSurfaceMesh::update()
     // Does approximation by looking at the face centres only
     const pointField& fc = surface_.faceCentres();
 
-    cellLabels_.setSize(fc.size());
-
     meshSearch meshSearcher(mesh(), false);
 
+    const indexedOctree<treeDataPoint>& cellCentreTree =
+        meshSearcher.cellCentreTree();
+
+
+    // Global numbering for cells - only used to uniquely identify local cells.
+    globalIndex globalCells(mesh().nCells());
+    List<nearInfo> nearest(fc.size());
+    forAll(nearest, i)
+    {
+        nearest[i].first() = GREAT;
+        nearest[i].second() = labelMax;
+    }
+
+    // Search triangles using nearest on local mesh
     forAll(fc, triI)
     {
-        cellLabels_[triI] = meshSearcher.findCell
+        pointIndexHit nearInfo = cellCentreTree.findNearest
         (
             fc[triI],
-            -1,                 // seed
+            sqr(GREAT)
-            true                // use tree
         );
-    }
+        if (nearInfo.hit())
-
-
-    boolList include(surface_.size(), true);
-
-    label nNotFound = 0;
-
-    forAll(cellLabels_, triI)
-    {
-        if (cellLabels_[triI] == -1)
         {
-            include[triI] = false;
+            nearest[triI].first() = magSqr(nearInfo.hitPoint()-fc[triI]);
-            nNotFound++;
+            nearest[triI].second() = globalCells.toGlobal(nearInfo.index());
         }
     }
-    label nTotalNotFound = returnReduce(nNotFound, sumOp<label>());
+
+    // See which processor has the nearest.
+    Pstream::listCombineGather(nearest, nearestEqOp());
+    Pstream::listCombineScatter(nearest);
+
+    boolList include(surface_.size(), false);
+
+    cellLabels_.setSize(fc.size());
+    cellLabels_ = -1;
+
+    label nFound = 0;
+    forAll(nearest, triI)
+    {
+        if (nearest[triI].second() == labelMax)
+        {
+            // Not found on any processor. How to map?
+        }
+        else if (globalCells.isLocal(nearest[triI].second()))
+        {
+            cellLabels_[triI] = globalCells.toLocal(nearest[triI].second());
+
+            include[triI] = true;
+            nFound++;
+        }
+    }
+
 
     if (debug)
     {
-        Pout<< "surface:" << surface_.size()
+        Pout<< "Local out of faces:" << cellLabels_.size()
-            << "  included:" << surface_.size()-nNotFound
+            << " keeping:" << nFound << endl;
-            << "  total not found" << nTotalNotFound << endl;
     }
 
+    // Now subset the surface. Do not use triSurface::subsetMesh since requires
+    // original surface to be in compact numbering.
+
+    const triSurface& s = surface_;
+
+    // Compact to original triangle
+    labelList faceMap(s.size());
+    // Compact to original points
+    labelList pointMap(s.points().size());
+    // From original point to compact points
+    labelList reversePointMap(s.points().size(), -1);
 
-    // Add to master all triangles are outside all meshes.
     {
-        boolList onAnyProc(surface_.size(), false);
+        label newPointI = 0;
-        Pstream::listCombineGather(onAnyProc, orEqOp<bool>());
+        label newTriI = 0;
 
-        if (Pstream::master())
+        forAll(s, triI)
         {
-            label nAdded = 0;
+            if (include[triI])
-            forAll(onAnyProc, triI)
             {
-                if (!onAnyProc[triI])
+                faceMap[newTriI++] = triI;
+
+                const labelledTri& f = s[triI];
+                forAll(f, fp)
                 {
-                    include[triI] = true;
+                    if (reversePointMap[f[fp]] == -1)
-                    nAdded++;
+                    {
+                        pointMap[newPointI] = f[fp];
+                        reversePointMap[f[fp]] = newPointI++;
+                    }
                 }
             }
+        }
+        faceMap.setSize(newTriI);
+        pointMap.setSize(newPointI);
+    }
 
-            if (debug)
+    // Subset cellLabels
-            {
+    cellLabels_ = UIndirectList<label>(cellLabels_, faceMap)();
-                Pout<< "nAdded to master:" << nAdded << endl;
+
-            }
+    // Store any face per point
+    pointToFace_.setSize(pointMap.size());
+
+    // Create faces and points for subsetted surface
+    faceList& faces = this->storedFaces();
+    faces.setSize(faceMap.size());
+    forAll(faceMap, i)
+    {
+        const triFace& f = s[faceMap[i]];
+        triFace newF
+        (
+            reversePointMap[f[0]],
+            reversePointMap[f[1]],
+            reversePointMap[f[2]]
+        );
+        faces[i] = newF.triFaceFace();
+
+        forAll(newF, fp)
+        {
+            pointToFace_[newF[fp]] = i;
         }
     }
 
-
+    this->storedPoints() = pointField(s.points(), pointMap);
-    // Now subset the surface
-    triSurface localSurface
-    (
-        surface_.subsetMesh
-        (
-            include,
-            pointMap_,
-            faceMap_
-        )
-    );
-
-    // And convert into faces.
-    faceList& faces = this->storedFaces();
-    faces.setSize(localSurface.size());
-    forAll(localSurface, i)
-    {
-        faces[i] = localSurface[i].triFaceFace();
-    }
-
-    this->storedPoints() = localSurface.points();
 
     if (debug)
     {
         print(Pout);
-        Pout << endl;
+        Pout<< endl;
     }
 
     needsUpdate_ = false;

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H

@@ -29,15 +29,13 @@ Description
     A sampledSurface from a triSurfaceMesh. It samples on the points/triangles
     of the triSurface.
 
+    It samples using the cell nearest to the triangle centre so does
+    not check the cell the centre is actually in ...
+
     In parallel every processor just operates on the part of the surface
     where the face centres are inside the mesh. It is then up to the
     caller to stitch the partial surfaces together.
 
-    No check is done for triangle centres being on multiple processors
-    (should never happen but truncation errors ...)
-
-    Any value on a triangle outside the mesh will be set to 0.
-
 SourceFiles
     sampledTriSurfaceMesh.C
 
@@ -75,14 +73,11 @@ class sampledTriSurfaceMesh
         //- Track if the surface needs an update
         mutable bool needsUpdate_;
 
-        //- Local cell labels
+        //- From local surface triangle to mesh cell.
         labelList cellLabels_;
 
         //- From local surface back to surface_
-        labelList pointMap_;
+        labelList pointToFace_;
-
-        //- From local surface back to surface_
-        labelList faceMap_;
 
 
     // Private Member Functions

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C

@@ -36,20 +36,12 @@ Foam::sampledTriSurfaceMesh::sampleField
 ) const
 {
     // One value per face
-    tmp<Field<Type> > tvalues(new Field<Type>(faceMap_.size()));
+    tmp<Field<Type> > tvalues(new Field<Type>(cellLabels_.size()));
     Field<Type>& values = tvalues();
 
-    forAll(faceMap_, i)
+    forAll(cellLabels_, triI)
     {
-        label cellI = cellLabels_[faceMap_[i]];
+        values[triI] = vField[cellLabels_[triI]];
-        if (cellI != -1)
-        {
-            values[i] = vField[cellI];
-        }
-        else
-        {
-            values[i] = pTraits<Type>::zero;
-        }
     }
 
     return tvalues;
@@ -64,24 +56,15 @@ Foam::sampledTriSurfaceMesh::interpolateField
 ) const
 {
     // One value per vertex
-    tmp<Field<Type> > tvalues(new Field<Type>(pointMap_.size()));
+    tmp<Field<Type> > tvalues(new Field<Type>(pointToFace_.size()));
     Field<Type>& values = tvalues();
 
-    forAll(pointMap_, i)
+    forAll(pointToFace_, pointI)
     {
-        label origPointI = pointMap_[i];
+        label triI = pointToFace_[pointI];
-        label origTriI = surface_.pointFaces()[origPointI][0];
+        label cellI = cellLabels_[triI];
 
-        label cellI = cellLabels_[origTriI];
+        values[pointI] = interpolator.interpolate(points()[pointI], cellI);
-
-        if (cellI != -1)
-        {
-            values[i] = interpolator.interpolate(points()[i], cellI);
-        }
-        else
-        {
-            values[i] = pTraits<Type>::zero;
-        }
     }
 
     return tvalues;