ENH: additional handling for out-of-range sampling (#1891)

- when sampling onto a meshed surface, the sampling surface may be
  outside of the mesh region, or simply too far away to be considered
  reasonable.

  Can now specify a max search distance and default values for samples
  that are too distant.
  If a default value is not specified, uses Type(Zero).

  Eg,

      maxDistance     0.005;
      defaultValue
      {
          "p.*"   1e5;
          T       273.15;
          U       (-100 -100 -100);
      }

src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C

@@ -170,7 +170,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
     // Does approximation by looking at the face centres only
     const pointField& fc = surface_.faceCentres();
 
-    List<nearInfo> nearest(fc.size(), nearInfo(GREAT, labelMax));
+    List<nearInfo> nearest(fc.size(), nearInfo(Foam::sqr(GREAT), labelMax));
 
     if (sampleSource_ == samplingSource::cells)
     {
@@ -181,13 +181,14 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
         forAll(fc, facei)
         {
             const point& pt = fc[facei];
+            auto& near = nearest[facei];
 
-            pointIndexHit info = cellTree.findNearest(pt, sqr(GREAT));
+            pointIndexHit info = cellTree.findNearest(pt, Foam::sqr(GREAT));
 
             if (info.hit())
             {
-                nearest[facei].first()  = magSqr(info.hitPoint()-pt);
+                near.first()  = magSqr(info.hitPoint()-pt);
-                nearest[facei].second() = globalCells.toGlobal(info.index());
+                near.second() = globalCells.toGlobal(info.index());
             }
         }
     }
@@ -200,14 +201,15 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
         forAll(fc, facei)
         {
             const point& pt = fc[facei];
+            auto& near = nearest[facei];
 
             if (cellTree.bb().contains(pt))
             {
                 const label index = cellTree.findInside(pt);
                 if (index != -1)
                 {
-                    nearest[facei].first()  = 0;
+                    near.first()  = 0;
-                    nearest[facei].second() = globalCells.toGlobal(index);
+                    near.second() = globalCells.toGlobal(index);
                 }
             }
         }
@@ -222,13 +224,14 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
         forAll(fc, facei)
         {
             const point& pt = fc[facei];
+            auto& near = nearest[facei];
 
             pointIndexHit info = bndTree.findNearest(pt, sqr(GREAT));
 
             if (info.hit())
             {
-                nearest[facei].first()  = magSqr(info.hitPoint()-pt);
+                near.first()  = magSqr(info.hitPoint()-pt);
-                nearest[facei].second() =
+                near.second() =
                     globalCells.toGlobal
                     (
                         bndTree.shapes().faceLabels()[info.index()]
@@ -250,16 +253,26 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
 
     forAll(nearest, facei)
     {
-        const label index = nearest[facei].second();
+        const auto& near = nearest[facei];
+
+        const label index = near.second();
 
         if (index == labelMax)
         {
-            // Not found on any processor. How to map?
+            // Not found on any processor, or simply too far.
+            // How to map?
         }
         else if (globalCells.isLocal(index))
         {
-            cellOrFaceLabels[facei] = globalCells.toLocal(index);
             facesToSubset.set(facei);
+
+            // Mark as special if too far away
+            cellOrFaceLabels[facei] =
+            (
+                (near.first() < maxDistanceSqr_)
+              ? globalCells.toLocal(index)
+              : -1
+            );
         }
     }
 
@@ -295,7 +308,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
     {
         // With point interpolation
 
-        samplePoints_.resize(pointMap.size());
+        samplePoints_ = points();
         sampleElements_.resize(pointMap.size(), -1);
 
         // Store any face per point (without using pointFaces())
@@ -314,32 +327,25 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
 
         if (sampleSource_ == samplingSource::cells)
         {
-            // samplePoints_   : per surface point a location inside the cell
+            // sampleElements_ : per surface point, the cell
-            // sampleElements_ : per surface point the cell
+            // samplePoints_   : per surface point, a location inside the cell
 
-            forAll(points(), pointi)
+            forAll(samplePoints_, pointi)
             {
-                const point& pt = points()[pointi];
+                // Use _copy_ of point
+                const point pt = samplePoints_[pointi];
 
                 const label celli = cellOrFaceLabels[pointToFace[pointi]];
 
                 sampleElements_[pointi] = celli;
 
-                // Check if point inside cell
                 if
                 (
-                    mesh().pointInCell
+                    celli >= 0
-                    (
+                 && !mesh().pointInCell(pt, celli, meshSearcher.decompMode())
-                        pt,
-                        sampleElements_[pointi],
-                        meshSearcher.decompMode()
-                    )
                 )
                 {
-                    samplePoints_[pointi] = pt;
+                    // Point not inside cell
-                }
-                else
-                {
                     // Find nearest point on faces of cell
 
                     scalar minDistSqr = VGREAT;
@@ -348,7 +354,12 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
                     {
                         const face& f = mesh().faces()[facei];
 
-                        pointHit info = f.nearestPoint(pt, mesh().points());
+                        pointHit info =
+                            f.nearestPoint
+                            (
+                                pt,
+                                mesh().points()
+                            );
 
                         if (info.distance() < minDistSqr)
                         {
@@ -361,50 +372,52 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
         }
         else if (sampleSource_ == samplingSource::insideCells)
         {
-            // samplePoints_   : per surface point a location inside the cell
             // sampleElements_ : per surface point the cell
+            // samplePoints_   : per surface point a location inside the cell
 
-            forAll(points(), pointi)
+            forAll(samplePoints_, pointi)
             {
-                const point& pt = points()[pointi];
-
                 const label celli = cellOrFaceLabels[pointToFace[pointi]];
 
                 sampleElements_[pointi] = celli;
-                samplePoints_[pointi] = pt;
             }
         }
         else  // samplingSource::boundaryFaces
         {
-            // samplePoints_   : per surface point a location on the boundary
+            // sampleElements_ : per surface point, the boundary face containing
-            // sampleElements_ : per surface point the boundary face containing
             //                   the location
+            // samplePoints_   : per surface point, a location on the boundary
 
-            forAll(points(), pointi)
+            forAll(samplePoints_, pointi)
             {
-                const point& pt = points()[pointi];
+                const point& pt = samplePoints_[pointi];
 
                 const label facei = cellOrFaceLabels[pointToFace[pointi]];
 
                 sampleElements_[pointi] = facei;
-                samplePoints_[pointi] = mesh().faces()[facei].nearestPoint
+
-                (
+                if (facei >= 0)
-                    pt,
+                {
-                    mesh().points()
+                    samplePoints_[pointi] =
-                ).rawPoint();
+                        mesh().faces()[facei].nearestPoint
+                        (
+                            pt,
+                            mesh().points()
+                        ).rawPoint();
+                }
             }
         }
     }
     else
     {
         // if sampleSource_ == cells:
-        //      sampleElements_ : per surface triangle the cell
+        //      sampleElements_ : per surface face, the cell
         //      samplePoints_   : n/a
         // if sampleSource_ == insideCells:
-        //      sampleElements_ : -1 or per surface triangle the cell
+        //      sampleElements_ : -1 or per surface face, the cell
         //      samplePoints_   : n/a
-        // else:
+        // if sampleSource_ == boundaryFaces:
-        //      sampleElements_ : per surface triangle the boundary face
+        //      sampleElements_ : per surface face, the boundary face
         //      samplePoints_   : n/a
 
         sampleElements_.transfer(cellOrFaceLabels);
@@ -432,14 +445,14 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
             forAll(samplePoints_, pointi)
             {
                 meshTools::writeOBJ(str, points()[pointi]);
-                vertI++;
+                ++vertI;
 
                 meshTools::writeOBJ(str, samplePoints_[pointi]);
-                vertI++;
+                ++vertI;
 
                 label elemi = sampleElements_[pointi];
                 meshTools::writeOBJ(str, centres[elemi]);
-                vertI++;
+                ++vertI;
                 str << "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI << nl;
             }
         }
@@ -449,11 +462,11 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
             forAll(sampleElements_, triI)
             {
                 meshTools::writeOBJ(str, faceCentres()[triI]);
-                vertI++;
+                ++vertI;
 
                 label elemi = sampleElements_[triI];
                 meshTools::writeOBJ(str, centres[elemi]);
-                vertI++;
+                ++vertI;
                 str << "l " << vertI-1 << ' ' << vertI << nl;
             }
         }
@@ -475,7 +488,7 @@ Foam::sampledMeshedSurface::sampledMeshedSurface
 )
 :
     sampledSurface(name, mesh),
-    MeshStorage(),
+    meshedSurface(),
     surfaceName_(surfaceName),
     surface_
     (
@@ -487,7 +500,9 @@ Foam::sampledMeshedSurface::sampledMeshedSurface
     keepIds_(true),
     zoneIds_(),
     sampleElements_(),
-    samplePoints_()
+    samplePoints_(),
+    maxDistanceSqr_(Foam::sqr(GREAT)),
+    defaultValues_()
 {}
 
 
@@ -499,7 +514,7 @@ Foam::sampledMeshedSurface::sampledMeshedSurface
 )
 :
     sampledSurface(name, mesh, dict),
-    MeshStorage(),
+    meshedSurface(),
     surfaceName_
     (
         meshedSurface::findFile
@@ -518,8 +533,25 @@ Foam::sampledMeshedSurface::sampledMeshedSurface
     keepIds_(dict.getOrDefault("keepIds", true)),
     zoneIds_(),
     sampleElements_(),
-    samplePoints_()
+    samplePoints_(),
+    maxDistanceSqr_(Foam::sqr(GREAT)),
+    defaultValues_(dict.subOrEmptyDict("defaultValue"))
 {
+    if (dict.readIfPresent("maxDistance", maxDistanceSqr_))
+    {
+        // Info<< "Limit to maxDistance " << maxDistanceSqr_ << nl;
+        // if (defaultValues_.empty())
+        // {
+        //     Info<< "defaultValues = zero" << nl;
+        // }
+        // else
+        // {
+        //     defaultValues_.writeEntry(Info);
+        // }
+
+        maxDistanceSqr_ = Foam::sqr(maxDistanceSqr_);
+    }
+
     wordRes includePatches;
     dict.readIfPresent("patches", includePatches);
     includePatches.uniq();
@@ -597,7 +629,7 @@ bool Foam::sampledMeshedSurface::expire()
     }
 
     sampledSurface::clearGeom();
-    MeshStorage::clear();
+    meshedSurface::clear();
     zoneIds_.clear();
 
     sampleElements_.clear();

src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H

@@ -35,33 +35,35 @@ Description
 
     - 6 different modes:
         - source=cells, interpolate=false:
-            finds per triangle centre the \em nearest cell centre and uses its value
+            finds per surface face centre the \em nearest cell centre
+            and uses its value
         - source=cells, interpolate=true:
-            finds per triangle centre the \em nearest cell centre.
+            finds per surface face centre the \em nearest cell centre.
             Per surface point checks if this nearest cell is the one containing
             point; otherwise projects the point onto the nearest point on
             the boundary of the cell (to make sure interpolateCellPoint
             gets a valid location)
 
         - source=insideCells, interpolate=false:
-            finds per triangle centre the cell containing it and uses its value.
+            finds per surface face centre the cell containing it
-            Trims triangles outside mesh.
+            and uses its value.
+            Trims surface faces outside of the mesh.
         - source=insideCells, interpolate=true:
             Per surface point interpolate cell containing it.
 
         - source=boundaryFaces, interpolate=false:
-            finds per triangle centre the nearest point on the boundary
+            finds per surface face centre the \em nearest point on the boundary
             (uncoupled faces only) and uses the value (or 0 if the nearest
             is on an empty boundary)
         - source=boundaryFaces, interpolate=true:
-            finds per triangle centre the nearest point on the boundary
+            finds per surface face centre the \em nearest point on the boundary
             (uncoupled faces only).
             Per surface point projects the point onto this boundary face
             (to make sure interpolateCellPoint gets a valid location)
 
-    - since it finds a nearest per triangle each triangle is guaranteed
+    - since it finds the nearest per surface face, each surface face
-    to be on one processor only. So after stitching (by sampledSurfaces)
+      is guaranteed to be on one processor only.
-    the original surface should be complete.
+      So after stitching the original surface should be complete.
 
     This is often embedded as part of a sampled surfaces function object.
 
@@ -75,6 +77,16 @@ Usage
             type    meshedSurface;
             surface something.obj;
             source  cells;
+
+            //- Max sampling distance
+            maxDistance 0.005;
+
+            //- Fallback for missed sampling in 'cells' mode
+            defaultValue
+            {
+                "p.*"   1e5;
+                T       273.15;
+            }
         }
     }
     \endverbatim
@@ -90,8 +102,11 @@ Usage
         file     | Alternative file name                    | no  |
         fileType | The surface format                       | no  | (extension)
         scale    | Surface scaling factor                   | no  | 0
+        maxDistance | Max search distance                   | no  | GREAT
+        defaultValue | Value beyond max distance (dictionary) | no  | empty
     \endtable
 
+
 SourceFiles
     sampledMeshedSurface.C
     sampledMeshedSurfaceTemplates.C
@@ -102,8 +117,7 @@ SourceFiles
 #define sampledMeshedSurface_H
 
 #include "sampledSurface.H"
-#include "MeshedSurface.H"
+#include "MeshedSurfaces.H"
-#include "MeshedSurfacesFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -123,18 +137,19 @@ class sampledMeshedSurface
     public meshedSurface
 {
 public:
+
         //- Types of sampling regions
-        enum class samplingSource
+        enum class samplingSource : unsigned char
         {
             cells,          //!< Use nearest cell value
-            insideCells,    //!< Face within a cell
+            insideCells,    //!< Surface face within a cell, or trim
-            boundaryFaces   //!< Use boundary face values
+            boundaryFaces   //!< Use nearest boundary face values
         };
 
 private:
 
-        //- Private typedefs for convenience
+    //- The concrete storage type for faces and points
-        typedef meshedSurface MeshStorage;
+    typedef meshedSurface Mesh;
 
 
     // Private Data
@@ -165,6 +180,12 @@ private:
         //- Local points to sample per point
         pointField samplePoints_;
 
+        //- Max search distance squared
+        scalar maxDistanceSqr_;
+
+        //- Out-of-range fallback values (when beyond the max search distance)
+        dictionary defaultValues_;
+
 
     // Private Member Functions
 
@@ -239,13 +260,13 @@ public:
         //- Points of surface
         virtual const pointField& points() const
         {
-            return MeshStorage::points();
+            return Mesh::points();
         }
 
         //- Faces of surface
         virtual const faceList& faces() const
         {
-            return MeshStorage::surfFaces();
+            return Mesh::surfFaces();
         }
 
         //- Per-face zone/region information
@@ -257,31 +278,31 @@ public:
         //- Per-face identifier (eg, element Id)
         virtual const labelList& faceIds() const
         {
-            return MeshStorage::faceIds();
+            return Mesh::faceIds();
         }
 
         //- Face area vectors
         virtual const vectorField& Sf() const
         {
-            return MeshStorage::Sf();
+            return Mesh::Sf();
         }
 
         //- Face area magnitudes
         virtual const scalarField& magSf() const
         {
-            return MeshStorage::magSf();
+            return Mesh::magSf();
         }
 
         //- Face centres
         virtual const vectorField& Cf() const
         {
-            return MeshStorage::Cf();
+            return Mesh::Cf();
         }
 
         //- If element ids/order of the original surface are kept
         virtual bool hasFaceIds() const
         {
-            return keepIds_ && !MeshStorage::faceIds().empty();
+            return keepIds_ && !Mesh::faceIds().empty();
         }
 
         //- Sampling boundary values instead of cell values

src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C

@@ -37,8 +37,14 @@ Foam::sampledMeshedSurface::sampleOnFaces
     const interpolation<Type>& sampler
 ) const
 {
+    const Type deflt
+    (
+        defaultValues_.getOrDefault<Type>(sampler.psi().name(), Zero)
+    );
+
     const labelList& elements = sampleElements_;
 
+
     if (!onBoundary())
     {
         // Sample cells
@@ -48,7 +54,8 @@ Foam::sampledMeshedSurface::sampleOnFaces
             sampler,
             elements,
             faces(),
-            points()
+            points(),
+            deflt
         );
     }
 
@@ -60,33 +67,33 @@ Foam::sampledMeshedSurface::sampleOnFaces
     auto tvalues = tmp<Field<Type>>::New(elements.size());
     auto& values = tvalues.ref();
 
-    const polyBoundaryMesh& pbm = mesh().boundaryMesh();
-    const label nBnd = mesh().nBoundaryFaces();
-
     // Create flat boundary field
 
-    Field<Type> bVals(nBnd, Zero);
+    const polyBoundaryMesh& pbm = mesh().boundaryMesh();
+
+    Field<Type> bVals(mesh().nBoundaryFaces(), Zero);
 
     const auto& bField = sampler.psi().boundaryField();
 
     forAll(bField, patchi)
     {
-        const label bFacei = (pbm[patchi].start() - mesh().nInternalFaces());
+        SubList<Type>(bVals, pbm[patchi].range()) = bField[patchi];
-
-        SubList<Type>
-        (
-            bVals,
-            bField[patchi].size(),
-            bFacei
-        ) = bField[patchi];
     }
 
-    // Sample in flat boundary field
+    // Sample within the flat boundary field
 
     forAll(elements, i)
     {
         const label bFacei = (elements[i] - mesh().nInternalFaces());
-        values[i] = bVals[bFacei];
+
+        if (bFacei < 0)
+        {
+            values[i] = deflt;
+        }
+        else
+        {
+            values[i] = bVals[bFacei];
+        }
     }
 
     return tvalues;
@@ -100,34 +107,60 @@ Foam::sampledMeshedSurface::sampleOnPoints
     const interpolation<Type>& interpolator
 ) const
 {
-    // One value per vertex
+    const Type deflt
-    auto tvalues = tmp<Field<Type>>::New(sampleElements_.size());
+    (
+        defaultValues_.getOrDefault<Type>(interpolator.psi().name(), Zero)
+    );
+
+    const labelList& elements = sampleElements_;
+
+
+    // One value per vertex.
+    // - sampleElements_ and samplePoints_ have the identical size
+    auto tvalues = tmp<Field<Type>>::New(elements.size());
     auto& values = tvalues.ref();
 
     if (!onBoundary())
     {
         // Sample cells
 
-        forAll(sampleElements_, pointi)
+        forAll(elements, i)
         {
-            values[pointi] = interpolator.interpolate
+            const label celli = elements[i];
-            (
+
-                samplePoints_[pointi],
+            if (celli < 0)
-                sampleElements_[pointi]
+            {
-            );
+                values[i] = deflt;
+            }
+            else
+            {
+                values[i] = interpolator.interpolate
+                (
+                    samplePoints_[i],
+                    celli
+                );
+            }
         }
     }
-    else
+
+
+    //
+    // Sample boundary faces
+    //
+
+    forAll(elements, i)
     {
-        // Sample boundary faces
+        const label facei = elements[i];
 
-        forAll(samplePoints_, pointi)
+        if (facei < 0)
         {
-            const label facei = sampleElements_[pointi];
+            values[i] = deflt;
-
+        }
-            values[pointi] = interpolator.interpolate
+        else
+        {
+            values[i] = interpolator.interpolate
             (
-                samplePoints_[pointi],
+                samplePoints_[i],
                 mesh().faceOwner()[facei],
                 facei
             );

src/sampling/sampledSurface/sampledSurface/sampledSurface.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -132,14 +132,16 @@ protected:
 
     // Protected Member Functions
 
-        //- General loop for sampling elements to faces
+        //- General loop for sampling elements to faces.
+        //  The defaultValue is used for invalid (negative) elements
         template<class Type>
         static tmp<Field<Type>> sampleOnFaces
         (
             const interpolation<Type>& sampler,
             const labelUList& elements,
             const faceList& fcs,
-            const pointField& pts
+            const pointField& pts,
+            const Type& defaultValue = Type(Zero)
         );
 
 

src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,7 +37,8 @@ Foam::sampledSurface::sampleOnFaces
     const interpolation<Type>& sampler,
     const labelUList& elements,
     const faceList& fcs,
-    const pointField& pts
+    const pointField& pts,
+    const Type& defaultValue
 )
 {
     const label len = elements.size();
@@ -57,9 +58,16 @@ Foam::sampledSurface::sampleOnFaces
     for (label i=0; i < len; ++i)
     {
         const label celli = elements[i];
-        const point pt = fcs[i].centre(pts);
+        if (celli < 0)
+        {
+            values[i] = defaultValue;
+        }
+        else
+        {
+            const point pt = fcs[i].centre(pts);
 
-        values[i] = sampler.interpolate(pt, celli);
+            values[i] = sampler.interpolate(pt, celli);
+        }
     }
 
     return tvalues;

tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling

@@ -121,6 +121,44 @@ sampled
 }
 
 
+// Sample and write for meshed surfaces
+sampleMesh
+{
+    type    surfaces;
+    libs    (sampling);
+    log     true;
+
+    writeControl    writeTime;
+    writeInterval   1;
+
+    surfaceFormat   vtk;
+
+    fields      (p rho U T);
+
+    surfaces
+    {
+        // Oversized sampling - for general testing
+        meshed_sample
+        {
+            type        meshedSurface;
+            surface     oversized_sample.obj;
+            source      cells;
+            maxDistance 0.0025;
+            defaultValue
+            {
+                T     273;
+            }
+        }
+
+        meshed_interpolate
+        {
+            $meshed_sample;
+            interpolate true;
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Calculations  * * * * * * * * * * * * * * * //
 
 massflow