Adding boolean to cellSizeFunctions to bail out quickly if told that operating

on a surface point.

Added tolerance static varible to cellSizeFunction and checks for points on or
very near to the surface and "snap" them to the surface to cell size.

Moved

    scalar maxProtrusionDistance = maxSurfaceProtrusion(vert);

out of the for loop in dualCellLargestSurfaceProtrusion, shouldn't have been in
- erroneous and expensive to calculate.

Stablised triSurfaceTools normalisation of c and added more info to FatalError
message.

src/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C

@@ -181,7 +181,8 @@ bool Foam::cellSizeControlSurfaces::samePriorityNext(label i) const
 
 Foam::scalar Foam::cellSizeControlSurfaces::cellSize
 (
-    const point& pt
+    const point& pt,
+    bool isSurfacePoint
 ) const
 {
     scalar sizeAccumulator = 0;
@@ -193,7 +194,7 @@ Foam::scalar Foam::cellSizeControlSurfaces::cellSize
 
         scalar sizeI;
 
-        if (cSF.cellSize(pt, sizeI))
+        if (cSF.cellSize(pt, sizeI, isSurfacePoint))
         {
             sizeAccumulator += sizeI;
             numberOfFunctions++;

src/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.H

@@ -126,7 +126,11 @@ public:
         // Query
 
             //- Return the cell size at the given location
-            scalar cellSize(const point& pt) const;
+            scalar cellSize
+            (
+                const point& pt,
+                bool isSurfacePoint = false
+            ) const;
 };
 
 

src/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/cellSizeFunction/cellSizeFunction.C

@@ -37,6 +37,8 @@ namespace Foam
 defineTypeNameAndDebug(cellSizeFunction, 0);
 defineRunTimeSelectionTable(cellSizeFunction, dictionary);
 
+scalar cellSizeFunction::snapToSurfaceTol_ = 1e-10;
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 

src/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/cellSizeFunction/cellSizeFunction.H

@@ -73,6 +73,13 @@ public:
 
 protected:
 
+    // Static data
+
+        //- Point closeness tolerance to a surface where the function "snaps" to
+        //  including the surface
+        static scalar snapToSurfaceTol_;
+
+
     // Protected data
 
         //- Reference to the conformalVoronoiMesh holding this cvs object
@@ -165,7 +172,12 @@ public:
         //  Return a boolean specifying if the function was used, i.e. false if
         //  the point was not in range of the surface for a spatially varying
         //  size.
-        virtual bool cellSize(const point& pt, scalar& size) const = 0;
+        virtual bool cellSize
+        (
+            const point& pt,
+            scalar& size,
+            bool isSurfacePoint = false
+        ) const = 0;
 };
 
 

src/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearDistance/linearDistance.C

@@ -37,6 +37,7 @@ namespace Foam
 defineTypeNameAndDebug(linearDistance, 0);
 addToRunTimeSelectionTable(cellSizeFunction, linearDistance, dictionary);
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 linearDistance::linearDistance
@@ -65,8 +66,20 @@ scalar linearDistance::sizeFunction(scalar d) const
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool linearDistance::cellSize(const point& pt, scalar& size) const
+bool linearDistance::cellSize
+(
+    const point& pt,
+    scalar& size,
+    bool isSurfacePoint
+) const
 {
+    if (isSurfacePoint)
+    {
+        size = surfaceCellSize_;
+
+        return true;
+    }
+
     size = 0;
 
     List<pointIndexHit> hits;
@@ -89,6 +102,15 @@ bool linearDistance::cellSize(const point& pt, scalar& size) const
             return true;
         }
 
+        // If the nearest point is essentially on the surface, do not do a
+        // getVolumeType calculation, as it will be prone to error.
+        if (mag(pt  - hitInfo.hitPoint()) < snapToSurfaceTol_)
+        {
+            size = sizeFunction(0);
+
+            return true;
+        }
+
         pointField ptF(1, pt);
         List<searchableSurface::volumeType> vTL;
 

src/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearDistance/linearDistance.H

@@ -104,7 +104,12 @@ public:
         //  Return a boolean specifying if the function was used, i.e. false if
         //  the point was not in range of the surface for a spatially varying
         //  size.
-        virtual bool cellSize(const point& pt, scalar& size) const;
+        virtual bool cellSize
+        (
+            const point& pt,
+            scalar& size,
+            bool isSurfacePoint = false
+        ) const;
 };
 
 

src/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/surfaceOffsetLinearDistance/surfaceOffsetLinearDistance.C

@@ -83,8 +83,20 @@ scalar surfaceOffsetLinearDistance::sizeFunction(scalar d) const
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool surfaceOffsetLinearDistance::cellSize(const point& pt, scalar& size) const
+bool surfaceOffsetLinearDistance::cellSize
+(
+    const point& pt,
+    scalar& size,
+    bool isSurfacePoint
+) const
 {
+    if (isSurfacePoint)
+    {
+        size = surfaceCellSize_;
+
+        return true;
+    }
+
     size = 0;
 
     List<pointIndexHit> hits;
@@ -107,6 +119,15 @@ bool surfaceOffsetLinearDistance::cellSize(const point& pt, scalar& size) const
             return true;
         }
 
+        // If the nearest point is essentially on the surface, do not do a
+        // getVolumeType calculation, as it will be prone to error.
+        if (mag(pt  - hitInfo.hitPoint()) < snapToSurfaceTol_)
+        {
+            size = sizeFunction(0);
+
+            return true;
+        }
+
         pointField ptF(1, pt);
         List<searchableSurface::volumeType> vTL;
 

src/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/surfaceOffsetLinearDistance/surfaceOffsetLinearDistance.H

@@ -111,7 +111,12 @@ public:
         //  Return a boolean specifying if the function was used, i.e. false if
         //  the point was not in range of the surface for a spatially varying
         //  size.
-        virtual bool cellSize(const point& pt, scalar& size) const;
+        virtual bool cellSize
+        (
+            const point& pt,
+            scalar& size,
+            bool isSurfacePoint = false
+        ) const;
 };
 
 

src/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniform/uniform.C

@@ -53,9 +53,19 @@ uniform::uniform
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool uniform::cellSize(const point& pt, scalar& size) const
+bool uniform::cellSize
+(
+    const point& pt,
+    scalar& size,
+    bool isSurfacePoint
+) const
 {
-    size = 0;
+    if (isSurfacePoint)
+    {
+        size = cellSize_;
+
+        return true;
+    }
 
     if (sideMode_ == BOTHSIDES)
     {
@@ -64,6 +74,28 @@ bool uniform::cellSize(const point& pt, scalar& size) const
         return true;
     }
 
+    size = 0;
+
+    List<pointIndexHit> hits;
+
+    surface_.findNearest
+    (
+        pointField(1, pt),
+        scalarField(1, sqr(snapToSurfaceTol_)),
+        hits
+    );
+
+    const pointIndexHit& hitInfo = hits[0];
+
+    // If the nearest point is essentially on the surface, do not do a
+    // getVolumeType calculation, as it will be prone to error.
+    if (hitInfo.hit())
+    {
+        size = cellSize_;
+
+        return true;
+    }
+
     pointField ptF(1, pt);
     List<searchableSurface::volumeType> vTL(1);
 

src/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniform/uniform.H

@@ -86,7 +86,12 @@ public:
         //  Return a boolean specifying if the function was used, i.e. false if
         //  the point was not in range of the surface for a spatially varying
         //  size.
-        virtual bool cellSize(const point& pt, scalar& size) const;
+        virtual bool cellSize
+        (
+            const point& pt,
+            scalar& size,
+            bool isSurfacePoint = false
+        ) const;
 };
 
 

src/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniformDistance/uniformDistance.C

@@ -55,8 +55,20 @@ uniformDistance::uniformDistance
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool uniformDistance::cellSize(const point& pt, scalar& size) const
+bool uniformDistance::cellSize
+(
+    const point& pt,
+    scalar& size,
+    bool isSurfacePoint
+) const
 {
+    if (isSurfacePoint)
+    {
+        size = cellSize_;
+
+        return true;
+    }
+
     size = 0;
 
     List<pointIndexHit> hits;
@@ -68,7 +80,9 @@ bool uniformDistance::cellSize(const point& pt, scalar& size) const
         hits
     );
 
-    if (hits[0].hit())
+    const pointIndexHit& hitInfo = hits[0];
+
+    if (hitInfo.hit())
     {
         if (sideMode_ == BOTHSIDES)
         {
@@ -77,6 +91,15 @@ bool uniformDistance::cellSize(const point& pt, scalar& size) const
             return true;
         }
 
+        // If the nearest point is essentially on the surface, do not do a
+        // getVolumeType calculation, as it will be prone to error.
+        if (mag(pt  - hitInfo.hitPoint()) < snapToSurfaceTol_)
+        {
+            size = cellSize_;
+
+            return true;
+        }
+
         pointField ptF(1, pt);
         List<searchableSurface::volumeType> vTL;
 

src/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniformDistance/uniformDistance.H

@@ -92,7 +92,12 @@ public:
         //  Return a boolean specifying if the function was used, i.e. false if
         //  the point was not in range of the surface for a spatially varying
         //  size.
-        virtual bool cellSize(const point& pt, scalar& size) const;
+        virtual bool cellSize
+        (
+            const point& pt,
+            scalar& size,
+            bool isSurfacePoint = false
+        ) const;
 };
 
 

src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C

@@ -32,8 +32,6 @@ License
 #include "ulong.H"
 #include "surfaceFeatures.H"
 
-#include "zeroGradientPointPatchField.H"
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::conformalVoronoiMesh::conformalVoronoiMesh
@@ -117,51 +115,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
     writeMesh();
     timeCheck();
 
-    // Test field creation and cell size functions
-
-    Info<< "Create fvMesh" << endl;
-
-    fvMesh fMesh
-    (
-        IOobject
-        (
-            Foam::polyMesh::defaultRegion,
-            runTime_.constant(),
-            runTime_,
-            IOobject::MUST_READ
-        )
-    );
-
-    timeCheck();
-
-    Info<< "Create volScalarField" << endl;
-
-    volScalarField cellSizeTest
-    (
-        IOobject
-        (
-            "cellSizeTest",
-            runTime_.timeName(),
-            runTime_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        fMesh,
-        dimensionedScalar("cellSize", dimLength, 0),
-        zeroGradientPointPatchField<scalar>::typeName
-    );
-
-    scalarField& cellSize = cellSizeTest.internalField();
-
-    const vectorField& C = fMesh.cellCentres();
-
-    forAll(cellSize, i)
-    {
-        cellSize[i] = cellSizeControl().cellSize(C[i]);
-    }
-
-    cellSizeTest.write();
-
+    writeTargetCellSize();
     timeCheck();
 }
 
@@ -174,6 +128,15 @@ Foam::conformalVoronoiMesh::~conformalVoronoiMesh()
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
+Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
+(
+    const point& pt
+) const
+{
+    return I;
+}
+
+
 void Foam::conformalVoronoiMesh::insertSurfacePointPairs
 (
     const List<pointIndexHit>& surfaceHits,
@@ -434,18 +397,9 @@ void Foam::conformalVoronoiMesh::insertMultipleEdgePointGroup
 
 void Foam::conformalVoronoiMesh::createFeaturePoints()
 {
-    const boundBox& bb = geometryToConformTo().bounds();
-
-    scalar span
-    (
-        max(mag(bb.max().x()), mag(bb.min().x()))
-      + max(mag(bb.max().y()), mag(bb.min().y()))
-      + max(mag(bb.max().z()), mag(bb.min().z()))
-    );
-
     Info<< nl << "Creating bounding points" << endl;
 
-    scalar bigSpan = 10*span;
+    scalar bigSpan = 10*cvMeshControls().span();
 
     insertPoint(point(-bigSpan, -bigSpan, -bigSpan), Vb::FAR_POINT);
     insertPoint(point(-bigSpan, -bigSpan,  bigSpan), Vb::FAR_POINT);
@@ -827,6 +781,10 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
     std::list<Facet> facets;
     incident_facets(vit, std::back_inserter(facets));
 
+    point vert(topoint(vit->point()));
+
+    scalar maxProtrusionDistance = maxSurfaceProtrusion(vert);
+
     for
     (
         std::list<Facet>::iterator fit=facets.begin();
@@ -840,10 +798,6 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
          && !is_infinite(fit->first->neighbor(fit->second))
         )
         {
-            point vert(topoint(vit->point()));
-
-            scalar maxProtrusionDistance = maxSurfaceProtrusion(vert);
-
             point edgeMid =
                 0.5
                *(
@@ -860,7 +814,7 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
                 edgeMid,
                 surfHit,
                 hitSurface
-           );
+            );
 
             if (surfHit.hit())
             {
@@ -901,7 +855,8 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
 
     scalar exclusionRangeCoeff = 0.2;
 
-    scalar exclusionRangeSqr = sqr(exclusionRangeCoeff*targetCellSize(pt));
+    scalar exclusionRangeSqr =
+        sqr(exclusionRangeCoeff*targetCellSize(pt, true));
 
     // 0.01 and 1000 determined from speed tests, varying the indexedOctree
     // rebuild frequency and balance of additions to queries.
@@ -1789,7 +1744,7 @@ void Foam::conformalVoronoiMesh::conformToSurface()
 
     label iterationNo = 0;
 
-    label maxIterations = 2;
+    label maxIterations = 10;
     Info << "    MAX ITERATIONS HARD CODED TO "<< maxIterations << endl;
 
     // Set totalHits to a positive value to enter the while loop on the first
@@ -1873,7 +1828,7 @@ void Foam::conformalVoronoiMesh::conformToSurface()
             (
                 featureEdgeHits,
                 featureEdgeFeaturesHit,
-                "edgeGenerationLocations_" + name(iterationNo) + ".obj"
+                "edgeConformationLocations_" + name(iterationNo) + ".obj"
             );
         }
 

src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H

@@ -131,8 +131,14 @@ class conformalVoronoiMesh
         //- Write the elapsedCpuTime
         inline void timeCheck() const;
 
-        //- Return the local target cell size at the given location
-        inline scalar targetCellSize(const point& pt) const;
+        //- Return the local target cell size at the given location.  Takes
+        //  boolean argument to allow speed-up of queries if the point is going
+        //  to be on a surface.
+        inline scalar targetCellSize
+        (
+            const point& pt,
+            bool isSurfacePoint = false
+        ) const;
 
         //- Return the local point pair separation at the given location
         inline scalar pointPairDistance(const point& pt) const;
@@ -153,7 +159,7 @@ class conformalVoronoiMesh
         inline scalar minimumEdgeLength(const point& pt) const;
 
         //- Return the required alignment directions at the given location
-        inline tensor requiredAlignment(const point& pt) const;
+        tensor requiredAlignment(const point& pt) const;
 
         //- Insert point and return it's index
         inline label insertPoint
@@ -412,6 +418,9 @@ public:
                 const fileName& fName
             ) const;
 
+            //- Read mesh from file as an fvMesh then calculate and write a
+            // field of the target cell size
+            void writeTargetCellSize() const;
 };
 
 

src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H

@@ -35,10 +35,11 @@ inline void Foam::conformalVoronoiMesh::timeCheck() const
 
 inline Foam::scalar Foam::conformalVoronoiMesh::targetCellSize
 (
-    const point& pt
+    const point& pt,
+    bool isSurfacePoint
 ) const
 {
-    return cvMeshControls_.defaultCellSize();
+    return cellSizeControl().cellSize(pt, isSurfacePoint);
 }
 
 
@@ -47,7 +48,10 @@ inline Foam::scalar Foam::conformalVoronoiMesh::pointPairDistance
     const point& pt
 ) const
 {
-    return targetCellSize(pt)*cvMeshControls_.pointPairDistanceCoeff();
+    // Point pair distances are always going to be at the surface, so the
+    // targetCellSize can be told to do a quick, surface only check.
+
+    return targetCellSize(pt, true)*cvMeshControls_.pointPairDistanceCoeff();
 }
 
 
@@ -67,10 +71,13 @@ inline Foam::scalar Foam::conformalVoronoiMesh::featurePointExclusionDistanceSqr
     const point& pt
 ) const
 {
+    // Exclusion distance tests are always going to be at the surface, so the
+    // targetCellSize can be told to do a quick, surface only check.
+
     return
         sqr
         (
-            targetCellSize(pt)
+            targetCellSize(pt, true)
            *cvMeshControls_.featurePointExclusionDistanceCoeff()
         );
 }
@@ -103,16 +110,6 @@ inline Foam::scalar Foam::conformalVoronoiMesh::minimumEdgeLength
 }
 
 
-inline Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
-(
-    const point& pt
-) const
-{
-    Info<< "STILL TO IMPLEMENT requiredAlignment." << endl;
-    return I;
-}
-
-
 inline Foam::label Foam::conformalVoronoiMesh::insertPoint
 (
     const point& p,

src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C

@@ -27,6 +27,7 @@ License
 #include "conformalVoronoiMesh.H"
 #include "IOstreams.H"
 #include "OFstream.H"
+#include "zeroGradientPointPatchField.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -176,4 +177,51 @@ void Foam::conformalVoronoiMesh::writeDual
 }
 
 
+void Foam::conformalVoronoiMesh::writeTargetCellSize() const
+{
+    Info<< nl << "Create fvMesh" << endl;
+
+    fvMesh fMesh
+    (
+        IOobject
+        (
+            Foam::polyMesh::defaultRegion,
+            runTime_.constant(),
+            runTime_,
+            IOobject::MUST_READ
+        )
+    );
+
+    timeCheck();
+
+    Info<< nl << "Create targetCellSize volScalarField" << endl;
+
+    volScalarField targetCellSize
+    (
+        IOobject
+        (
+            "targetCellSize",
+            runTime_.timeName(),
+            runTime_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        fMesh,
+        dimensionedScalar("cellSize", dimLength, 0),
+        zeroGradientPointPatchField<scalar>::typeName
+    );
+
+    scalarField& cellSize = targetCellSize.internalField();
+
+    const vectorField& C = fMesh.cellCentres();
+
+    forAll(cellSize, i)
+    {
+        cellSize[i] = cellSizeControl().cellSize(C[i]);
+    }
+
+    targetCellSize.write();
+}
+
+
 // ************************************************************************* //

src/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H

@@ -53,8 +53,8 @@ class indexedCell
 {
     // Private data
 
-        //- The index for this triangle Cell
-        //  -1: triangle and changed and associated data needs updating
+        //- The index for this tetrahedral cell
+        //  -1: tetrahedral and changed and associated data needs updating
         // >=0: index of cells, set by external update algorithm
         int index_;
 

src/conformalVoronoiMesh/cvControls/cvControls.C

@@ -38,7 +38,17 @@ Foam::cvControls::cvControls
     cvMesh_(cvMesh),
     cvMeshDict_(cvMeshDict)
 {
+    // General parameters
+
+    const boundBox& bb = cvMesh_.geometryToConformTo().bounds();
+
+    span_ =
+        max(mag(bb.max().x()), mag(bb.min().x()))
+      + max(mag(bb.max().y()), mag(bb.min().y()))
+      + max(mag(bb.max().z()), mag(bb.min().z()));
+
     // Surface conformation controls
+
     const dictionary& surfDict(cvMeshDict_.subDict("surfaceConformation"));
 
     pointPairDistanceCoeff_ = readScalar
@@ -69,9 +79,8 @@ Foam::cvControls::cvControls
     maxQuadAngle_= readScalar(surfDict.lookup("maxQuadAngle"));
 
     // Motion control controls
-    const dictionary& motionDict(cvMeshDict_.subDict("motionControl"));
 
-    defaultCellSize_ = readScalar(motionDict.lookup("defaultCellSize"));
+    const dictionary& motionDict(cvMeshDict_.subDict("motionControl"));
 
     const dictionary& insertionDict
     (
@@ -104,6 +113,7 @@ Foam::cvControls::cvControls
     );
 
     // polyMesh filtering controls
+
     const dictionary& filteringDict(cvMeshDict_.subDict("polyMeshFiltering"));
 
     minimumEdgeLengthCoeff_ = readScalar

src/conformalVoronoiMesh/cvControls/cvControls.H

@@ -60,6 +60,11 @@ class cvControls
         //- Reference to the cvMeshDict
         const IOdictionary& cvMeshDict_;
 
+        // General parameters
+
+            //- Span of the domain
+            scalar span_;
+
         // Surface conformation controls
 
             //- Point pair spacing coefficient - fraction of the local target
@@ -91,10 +96,6 @@ class cvControls
 
     // Motion control controls
 
-        // Storing the default cell size value for the moment until a method for
-        // storing the spatially resolved size requirement is implemented
-        scalar defaultCellSize_;
-
         // Point insertion criteria
 
             scalar cellCentreInsertionDistCoeff_;
@@ -138,8 +139,8 @@ public:
             //- Return the cvMeshDict
             inline const IOdictionary& cvMeshDict() const;
 
-            //- Return the defaultCellSize
-            inline scalar defaultCellSize() const;
+            //- Return the span
+            inline scalar span() const;
 
             //- Return the pointPairDistanceCoeff
             inline scalar pointPairDistanceCoeff() const;

src/conformalVoronoiMesh/cvControls/cvControlsI.H

@@ -32,9 +32,9 @@ inline const Foam::IOdictionary& Foam::cvControls::cvMeshDict() const
 }
 
 
-inline Foam::scalar Foam::cvControls::defaultCellSize() const
+inline Foam::scalar Foam::cvControls::span() const
 {
-    return defaultCellSize_;
+    return span_;
 }
 
 

src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C

@@ -2225,7 +2225,11 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
         // Nearest to face interior. Use faceNormal to determine side
         scalar c = (sample - nearestPoint) & surf.faceNormals()[nearestFaceI];
 
-        c /= (mag(sample - nearestPoint)*mag(surf.faceNormals()[nearestFaceI]));
+        c /= max
+        (
+            VSMALL,
+            mag(sample - nearestPoint)*mag(surf.faceNormals()[nearestFaceI])
+        );
 
         if (mag(c) < 0.9)
         {
@@ -2235,6 +2239,10 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
                 << "perpendicular to the normal." << nl
                 << "sample: " << sample << nl
                 << "nearestPoint: " << nearestPoint << nl
+                << "sample - nearestPoint: " << sample - nearestPoint << nl
+                << "normal: " << surf.faceNormals()[nearestFaceI] << nl
+                << "mag(sample - nearestPoint): "
+                << mag(sample - nearestPoint) << nl
                 << "normalised dot product: " << c << nl
                 << "triangle vertices: " << nl
                 << "    " << points[f[0]] << nl