Fixing conflict in triSurfaceMesh.C, committing failed automatic merge of master.

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C

@@ -52,24 +52,13 @@ void Foam::mapDistribute::distribute
 
             if (domain != Pstream::myProcNo() && map.size())
             {
-                List<T> subField(map.size());
-                forAll(map, i)
-                {
-                    subField[i] = field[map[i]];
-                }
                 OPstream toNbr(Pstream::blocking, domain);
-                toNbr << subField;
+                toNbr << UIndirectList<T>(field, map);
             }
         }
 
         // Subset myself
-        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+        UIndirectList<T> subField(field, subMap[Pstream::myProcNo()]);
-
-        List<T> subField(mySubMap.size());
-        forAll(mySubMap, i)
-        {
-            subField[i] = field[mySubMap[i]];
-        }
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
@@ -126,13 +115,7 @@ void Foam::mapDistribute::distribute
         List<T> newField(constructSize);
 
         // Subset myself
-        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+        UIndirectList<T> subField(field, subMap[Pstream::myProcNo()]);
-
-        List<T> subField(mySubMap.size());
-        forAll(mySubMap, i)
-        {
-            subField[i] = field[mySubMap[i]];
-        }
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
@@ -152,16 +135,8 @@ void Foam::mapDistribute::distribute
             if (Pstream::myProcNo() == sendProc)
             {
                 // I am sender. Send to recvProc.
-                const labelList& map = subMap[recvProc];
-
-                List<T> subField(map.size());
-                forAll(map, i)
-                {
-                    subField[i] = field[map[i]];
-                }
-
                 OPstream toNbr(Pstream::scheduled, recvProc);
-                toNbr << subField;
+                toNbr << UIndirectList<T>(field, subMap[recvProc]);
             }
             else
             {
@@ -374,24 +349,13 @@ void Foam::mapDistribute::distribute
 
             if (domain != Pstream::myProcNo() && map.size())
             {
-                List<T> subField(map.size());
-                forAll(map, i)
-                {
-                    subField[i] = field[map[i]];
-                }
                 OPstream toNbr(Pstream::blocking, domain);
-                toNbr << subField;
+                toNbr << UIndirectList<T>(field, map);
             }
         }
 
         // Subset myself
-        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+        UIndirectList<T> subField(field, subMap[Pstream::myProcNo()]);
-
-        List<T> subField(mySubMap.size());
-        forAll(mySubMap, i)
-        {
-            subField[i] = field[mySubMap[i]];
-        }
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
@@ -449,13 +413,7 @@ void Foam::mapDistribute::distribute
         List<T> newField(constructSize, nullValue);
 
         // Subset myself
-        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+        UIndirectList<T> subField(field, subMap[Pstream::myProcNo()]);
-
-        List<T> subField(mySubMap.size());
-        forAll(mySubMap, i)
-        {
-            subField[i] = field[mySubMap[i]];
-        }
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
@@ -475,16 +433,8 @@ void Foam::mapDistribute::distribute
             if (Pstream::myProcNo() == sendProc)
             {
                 // I am sender. Send to recvProc.
-                const labelList& map = subMap[recvProc];
-
-                List<T> subField(map.size());
-                forAll(map, i)
-                {
-                    subField[i] = field[map[i]];
-                }
-
                 OPstream toNbr(Pstream::scheduled, recvProc);
-                toNbr << subField;
+                toNbr << UIndirectList<T>(field, subMap[recvProc]);
             }
             else
             {

src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H

@@ -166,12 +166,9 @@ public:
                 const intersection::direction dir = intersection::VECTOR
             ) const;
 
-            //- Fast intersection with a ray. Distance is normalised distance
+            //- Fast intersection with a ray.
-            //  to the intersection (normalised with the vector magnitude).
+            //  For a hit, the pointHit.distance() is the line parameter t :
-            //  For a hit, the distance is signed. Positive number
+            //  intersection=p+t*q. Only defined for VISIBLE, FULL_RAY or
-            //  represents the point in front of triangle. In case of miss
-            //  pointHit position is undefined.
-            //  Only defined for VISIBLE, FULL_RAY or
             //  HALF_RAY. tol increases the virtual size of the triangle
             // by a relative factor.
             inline pointHit intersection

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -45,6 +45,7 @@ Description
 #include "OFstream.H"
 #include "layerParameters.H"
 #include "combineFaces.H"
+#include "IOmanip.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1414,8 +1415,11 @@ void Foam::autoLayerDriver::calculateLayerThickness
     const indirectPrimitivePatch& pp,
     const labelList& patchIDs,
     const scalarField& patchExpansionRatio,
-    const scalarField& patchFinalLayerRatio,
+
-    const scalarField& patchRelMinThickness,
+    const bool relativeSizes,
+    const scalarField& patchFinalLayerThickness,
+    const scalarField& patchMinThickness,
+
     const labelList& cellLevel,
     const labelList& patchNLayers,
     const scalar edge0Len,
@@ -1428,22 +1432,100 @@ void Foam::autoLayerDriver::calculateLayerThickness
     const fvMesh& mesh = meshRefiner_.mesh();
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-    if (min(patchRelMinThickness) < 0 || max(patchRelMinThickness) > 2)
+
+    // Rework patch-wise layer parameters into minimum per point
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Reuse input fields
+    expansionRatio.setSize(pp.nPoints());
+    expansionRatio = GREAT;
+    thickness.setSize(pp.nPoints());
+    thickness = GREAT;
+    minThickness.setSize(pp.nPoints());
+    minThickness = GREAT;
+
+    forAll(patchIDs, i)
+    {
+        label patchI = patchIDs[i];
+
+        const labelList& meshPoints = patches[patchI].meshPoints();
+
+        forAll(meshPoints, patchPointI)
+        {
+            label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
+
+            expansionRatio[ppPointI] = min
+            (
+                expansionRatio[ppPointI],
+                patchExpansionRatio[patchI]
+            );
+            thickness[ppPointI] = min
+            (
+                thickness[ppPointI],
+                patchFinalLayerThickness[patchI]
+            );
+            minThickness[ppPointI] = min
+            (
+                minThickness[ppPointI],
+                patchMinThickness[patchI]
+            );
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        expansionRatio,
+        minEqOp<scalar>(),
+        GREAT,              // null value
+        false               // no separation
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        thickness,
+        minEqOp<scalar>(),
+        GREAT,              // null value
+        false               // no separation
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        minThickness,
+        minEqOp<scalar>(),
+        GREAT,              // null value
+        false               // no separation
+    );
+
+
+    // Now the thicknesses are set according to the minimum of connected
+    // patches.
+
+
+    // Rework relative thickness into absolute
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // by multiplying with the internal cell size.
+
+    if (relativeSizes)
+    {
+        if (min(patchMinThickness) < 0 || max(patchMinThickness) > 2)
         {
             FatalErrorIn("calculateLayerThickness(..)")
                 << "Thickness should be factor of local undistorted cell size."
                 << " Valid values are [0..2]." << nl
-            << " minThickness:" << patchRelMinThickness
+                << " minThickness:" << patchMinThickness
                 << exit(FatalError);
         }
 
 
-    // Per point the max cell level of connected cells
+        // Determine per point the max cell level of connected cells
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
         labelList maxPointLevel(pp.nPoints(), labelMin);
 
-    {
         forAll(pp, i)
         {
             label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]];
@@ -1465,113 +1547,44 @@ void Foam::autoLayerDriver::calculateLayerThickness
             labelMin,           // null value
             false               // no separation
         );
-    }
 
 
-    // Rework patch-wise layer parameters into minimum per point
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    expansionRatio.setSize(pp.nPoints());
-    expansionRatio = GREAT;
-    scalarField finalLayerRatio(pp.nPoints(), GREAT);
-    scalarField relMinThickness(pp.nPoints(), GREAT);
-
-    {
-        forAll(patchIDs, i)
-        {
-            label patchI = patchIDs[i];
-
-            const labelList& meshPoints = patches[patchI].meshPoints();
-
-            forAll(meshPoints, patchPointI)
-            {
-                label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
-
-                expansionRatio[ppPointI] = min
-                (
-                    expansionRatio[ppPointI],
-                    patchExpansionRatio[patchI]
-                );
-                finalLayerRatio[ppPointI] = min
-                (
-                    finalLayerRatio[ppPointI],
-                    patchFinalLayerRatio[patchI]
-                );
-                relMinThickness[ppPointI] = min
-                (
-                    relMinThickness[ppPointI],
-                    patchRelMinThickness[patchI]
-                );
-            }
-        }
-
-        syncTools::syncPointList
-        (
-            mesh,
-            pp.meshPoints(),
-            expansionRatio,
-            minEqOp<scalar>(),
-            GREAT,              // null value
-            false               // no separation
-        );
-        syncTools::syncPointList
-        (
-            mesh,
-            pp.meshPoints(),
-            finalLayerRatio,
-            minEqOp<scalar>(),
-            GREAT,              // null value
-            false               // no separation
-        );
-        syncTools::syncPointList
-        (
-            mesh,
-            pp.meshPoints(),
-            relMinThickness,
-            minEqOp<scalar>(),
-            GREAT,              // null value
-            false               // no separation
-        );
-    }
-
-
-
-    // Per mesh point the expansion parameters
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    thickness.setSize(pp.nPoints());
-    minThickness.setSize(pp.nPoints());
-
         forAll(maxPointLevel, pointI)
         {
             // Find undistorted edge size for this level.
             scalar edgeLen = edge0Len/(1<<maxPointLevel[pointI]);
+            thickness[pointI] *= edgeLen;
+            minThickness[pointI] *= edgeLen;
+        }
+    }
 
+
+
+    // Rework thickness (of final layer) into overall thickness of all layers
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    forAll(thickness, pointI)
+    {
         // Calculate layer thickness based on expansion ratio
         // and final layer height
         if (expansionRatio[pointI] == 1)
         {
-            thickness[pointI] =
+            thickness[pointI] *= patchNLayers[pointI];
-                finalLayerRatio[pointI]
-              * patchNLayers[pointI]
-              * edgeLen;
-            minThickness[pointI] = relMinThickness[pointI]*edgeLen;
         }
         else
         {
+
             scalar invExpansion = 1.0 / expansionRatio[pointI];
             label nLay = patchNLayers[pointI];
-            thickness[pointI] =
+            thickness[pointI] *=
-                finalLayerRatio[pointI]
+                (1.0 - pow(invExpansion, nLay))
-              * edgeLen
-              * (1.0 - pow(invExpansion, nLay))
               / (1.0 - invExpansion);
-            minThickness[pointI] = relMinThickness[pointI]*edgeLen;
         }
     }
 
-    Info<< "calculateLayerThickness : min:" << gMin(thickness)
+
-        << " max:" << gMax(thickness) << endl;
+    //Info<< "calculateLayerThickness : min:" << gMin(thickness)
+    //    << " max:" << gMax(thickness) << endl;
 }
 
 
@@ -2621,8 +2634,11 @@ void Foam::autoLayerDriver::addLayers
         pp,
         meshMover.adaptPatchIDs(),
         layerParams.expansionRatio(),
-        layerParams.finalLayerRatio(),
+
-        layerParams.minThickness(),
+        layerParams.relativeSizes(),        // thickness relative to cellsize?
+        layerParams.finalLayerThickness(),  // wanted thicknes
+        layerParams.minThickness(),         // minimum thickness
+
         cellLevel,
         patchNLayers,
         edge0Len,
@@ -2632,6 +2648,79 @@ void Foam::autoLayerDriver::addLayers
         expansionRatio
     );
 
+
+    // Print a bit
+    {
+        const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+        Info<< nl
+            << "patch               faces    layers avg thickness[m]" << nl
+            << "                                    near-wall overall" << nl
+            << "-----               -----    ------ --------- -------" << endl;
+
+        forAll(meshMover.adaptPatchIDs(), i)
+        {
+            label patchI = meshMover.adaptPatchIDs()[i];
+
+            const labelList& meshPoints = patches[patchI].meshPoints();
+
+            //scalar maxThickness = -VGREAT;
+            //scalar minThickness = VGREAT;
+            scalar sumThickness = 0;
+            scalar sumNearWallThickness = 0;
+
+            forAll(meshPoints, patchPointI)
+            {
+                label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
+
+                //maxThickness = max(maxThickness, thickness[ppPointI]);
+                //minThickness = min(minThickness, thickness[ppPointI]);
+                sumThickness += thickness[ppPointI];
+
+                label nLay = patchNLayers[ppPointI];
+                if (nLay > 0)
+                {
+                    if (expansionRatio[ppPointI] == 1)
+                    {
+                        sumNearWallThickness += thickness[ppPointI]/nLay;
+                    }
+                    else
+                    {
+                        scalar s =
+                            (1.0-pow(expansionRatio[ppPointI], nLay))
+                          / (1.0-expansionRatio[ppPointI]);
+                        sumNearWallThickness += thickness[ppPointI]/s;
+                    }
+                }
+            }
+
+            label totNPoints = returnReduce(meshPoints.size(), sumOp<label>());
+
+            //reduce(maxThickness, maxOp<scalar>());
+            //reduce(minThickness, minOp<scalar>());
+            scalar avgThickness =
+                returnReduce(sumThickness, sumOp<scalar>())
+              / totNPoints;
+            scalar avgNearWallThickness =
+                returnReduce(sumNearWallThickness, sumOp<scalar>())
+              / totNPoints;
+
+            Info<< setf(ios_base::left) << setw(19) << patches[patchI].name();
+            //Sout.unsetf(ios_base::left);
+            Info<< setprecision(3)
+                << " " << setw(8)
+                << returnReduce(patches[patchI].size(), sumOp<scalar>())
+                << " " << setw(6) << layerParams.numLayers()[patchI]
+                << " " << setw(8) << avgNearWallThickness
+                << "  " << setw(8) << avgThickness
+                //<< " " << setw(8) << minThickness
+                //<< " " << setw(8) << maxThickness
+                << endl;
+        }
+        Info<< endl;
+    }
+
+
     // Calculate wall to medial axis distance for smoothing displacement
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -2835,7 +2924,8 @@ void Foam::autoLayerDriver::addLayers
             nPatchFaceLayers
         );
 
-        // Calculate displacement for first layer for addPatchLayer
+        // Calculate displacement for first layer for addPatchLayer.
+        // (first layer = layer of cells next to the original mesh)
         vectorField firstDisp(patchNLayers.size(), vector::zero);
 
         forAll(patchNLayers, i)
@@ -2851,9 +2941,9 @@ void Foam::autoLayerDriver::addLayers
                     label nLay = nPatchPointLayers[i];
                     scalar h =
                         pow(expansionRatio[i], nLay - 1)
-                      * (mag(patchDisp[i])*(1.0 - expansionRatio[i]))
+                      * (1.0 - expansionRatio[i])
                       / (1.0 - pow(expansionRatio[i], nLay));
-                    firstDisp[i] = h/mag(patchDisp[i])*patchDisp[i];
+                    firstDisp[i] = h*patchDisp[i];
                 }
             }
         }
@@ -2868,7 +2958,7 @@ void Foam::autoLayerDriver::addLayers
             pp,
             nPatchFaceLayers,   // layers per face
             nPatchPointLayers,  // layers per point
-            firstDisp,          // thickness of first layer
+            firstDisp,          // thickness of layer nearest internal mesh
             meshMod
         );
 

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H

@@ -255,9 +255,12 @@ class autoLayerDriver
                 (
                     const indirectPrimitivePatch& pp,
                     const labelList& patchIDs,
+
                     const scalarField& patchExpansionRatio,
-                    const scalarField& patchFinalLayerRatio,
+                    const bool relativeSizes,
-                    const scalarField& patchRelMinThickness,
+                    const scalarField& patchFinalLayerThickness,
+                    const scalarField& patchMinThickness,
+
                     const labelList& cellLevel,
                     const labelList& patchNLayers,
                     const scalar edge0Len,

src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C

@@ -163,7 +163,8 @@ Foam::layerParameters::layerParameters
         numLayers_.size(),
         readScalar(dict.lookup("expansionRatio"))
     ),
-    finalLayerRatio_
+    relativeSizes_(false),
+    finalLayerThickness_
     (
         numLayers_.size(),
         readScalar(dict.lookup("finalLayerRatio"))
@@ -244,10 +245,11 @@ Foam::layerParameters::layerParameters
         boundaryMesh.size(),
         readScalar(dict.lookup("expansionRatio"))
     ),
-    finalLayerRatio_
+    relativeSizes_(dict.lookup("relativeSizes")),
+    finalLayerThickness_
     (
         boundaryMesh.size(),
-        readScalar(dict.lookup("finalLayerRatio"))
+        readScalar(dict.lookup("finalLayerThickness"))
     ),
     minThickness_
     (
@@ -323,14 +325,21 @@ Foam::layerParameters::layerParameters
             numLayers_[patchI] =
                 readLabel(layerDict.lookup("nSurfaceLayers"));
 
-            //- Patch-wise layer parameters disabled for now. Just remove
+            layerDict.readIfPresent
-            //  settings in initialiser list and uncomment below.
+            (
-            //expansionRatio_[patchI] =
+                "expansionRatio",
-            //    readScalar(layerDict.lookup("expansionRatio"));
+                expansionRatio_[patchI]
-            //finalLayerRatio_[patchI] =
+            );
-            //    readScalar(layerDict.lookup("finalLayerRatio"));
+            layerDict.readIfPresent
-            //minThickness_[patchI] =
+            (
-            //    readScalar(layerDict.lookup("minThickness"));
+                "finalLayerThickness",
+                finalLayerThickness_[patchI]
+            );
+            layerDict.readIfPresent
+            (
+                "minThickness",
+                minThickness_[patchI]
+            );
         }
     }
 

src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H

@@ -39,6 +39,7 @@ SourceFiles
 #include "dictionary.H"
 #include "scalarField.H"
 #include "labelList.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -70,14 +71,10 @@ class layerParameters
 
             scalarField expansionRatio_;
 
-            //- Wanted thickness of final added cell layer. If multiple layers
+            Switch relativeSizes_;
-            //  is the thickness of the layer furthest away from the wall.
+
-            //  Relative to undistorted size of cell outside layer.
+            scalarField finalLayerThickness_;
-            scalarField finalLayerRatio_;
 
-            //- Minimum thickness of cell layer. If for any reason layer
-            //  cannot be above minThickness do not add layer.
-            //  Relative to undistorted size of cell outside layer.
             scalarField minThickness_;
 
 
@@ -164,18 +161,27 @@ public:
                     return expansionRatio_;
                 }
 
-                //- Wanted thickness of final added cell layer. If multiple
+                //- Are size parameters relative to inner cell size or
-                //  layers
+                //  absolute distances.
-                //  is the thickness of the layer furthest away from the wall.
+                bool relativeSizes() const
-                //  Relative to undistorted size of cell outside layer.
-                const scalarField& finalLayerRatio() const
                 {
-                    return finalLayerRatio_;
+                    return relativeSizes_;
+                }
+
+                //- Wanted thickness of final added cell layer. If multiple
+                //  layers is the thickness of the layer furthest away
+                //  from the wall (i.e. nearest the original mesh)
+                //  If relativeSize() this number is relative to undistorted
+                //  size of the cell outside layer.
+                const scalarField& finalLayerThickness() const
+                {
+                    return finalLayerThickness_;
                 }
 
                 //- Minimum thickness of cell layer. If for any reason layer
                 //  cannot be above minThickness do not add layer.
-                //  Relative to undistorted size of cell outside layer.
+                //  If relativeSize() this number is relative to undistorted
+                //  size of the cell outside layer.
                 const scalarField& minThickness() const
                 {
                     return minThickness_;

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H

@@ -294,9 +294,10 @@ public:
             //  - nPointLayers : number of layers per (patch)point
             //  - nFaceLayers : number of layers per (patch) face
             //  - firstDisplacement : displacement per point for first
-            //    layer of points. If zero do not add point.
+            //    layer of points (i.e. nearest to original mesh). If zero
+            //    do not add point.
             //  Layer thicknesses are calculated to constant geometric
-            //  expansion. Use expansionRatio for constant size.
+            //  expansion. Use expansionRatio 1 for constant size.
             //  Sets addedPoints_ which is per pp point a list of points
             //  added.
             //  Note: firstDisplacement has to be parallel synchronised before

src/meshTools/searchableSurface/triSurfaceMesh.C

@@ -237,6 +237,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
         )
     ),
     triSurface(s),
+    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
     surfaceClosed_(-1)
 {
     bounds() = boundBox(points());
@@ -281,6 +282,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
             searchableSurface::objectPath()
         )
     ),
+    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
     surfaceClosed_(-1)
 {
     bounds() = boundBox(points());
@@ -328,6 +330,7 @@ Foam::triSurfaceMesh::triSurfaceMesh
             searchableSurface::objectPath()
         )
     ),
+    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
     surfaceClosed_(-1)
 {
     scalar scaleFactor = 0;
@@ -336,10 +339,19 @@ Foam::triSurfaceMesh::triSurfaceMesh
     // eg, CAD geometries are often done in millimeters
     if (dict.readIfPresent("scale", scaleFactor) && scaleFactor > 0)
     {
+        Info<< searchableSurface::name() << " : using scale " << scaleFactor
+            << endl;
         triSurface::scalePoints(scaleFactor);
     }
 
     bounds() = boundBox(points());
+
+    // Have optional non-standard search tolerance for gappy surfaces.
+    if (dict.readIfPresent("tolerance", tolerance_) && tolerance_ > 0)
+    {
+        Info<< searchableSurface::name() << " : using intersection tolerance "
+            << tolerance_ << endl;
+    }
 }
 
 
@@ -386,6 +398,9 @@ const Foam::indexedOctree<Foam::treeDataTriSurface>&
         bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
         bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 
+        scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+        indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
+
         tree_.reset
         (
             new indexedOctree<treeDataTriSurface>
@@ -397,6 +412,8 @@ const Foam::indexedOctree<Foam::treeDataTriSurface>&
                 3.0     // duplicity
             )
         );
+
+        indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
     }
 
     return tree_();
@@ -497,6 +514,9 @@ void Foam::triSurfaceMesh::findNearest
 
     info.setSize(samples.size());
 
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
+
     forAll(samples, i)
     {
         static_cast<pointIndexHit&>(info[i]) = octree.findNearest
@@ -505,6 +525,8 @@ void Foam::triSurfaceMesh::findNearest
             nearestDistSqr[i]
         );
     }
+
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 }
 
 
@@ -519,6 +541,9 @@ void Foam::triSurfaceMesh::findLine
 
     info.setSize(start.size());
 
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
+
     forAll(start, i)
     {
         static_cast<pointIndexHit&>(info[i]) = octree.findLine
@@ -527,6 +552,8 @@ void Foam::triSurfaceMesh::findLine
             end[i]
         );
     }
+
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 }
 
 
@@ -541,6 +568,9 @@ void Foam::triSurfaceMesh::findLineAny
 
     info.setSize(start.size());
 
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
+
     forAll(start, i)
     {
         static_cast<pointIndexHit&>(info[i]) = octree.findLineAny
@@ -549,6 +579,8 @@ void Foam::triSurfaceMesh::findLineAny
             end[i]
         );
     }
+
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 }
 
 
@@ -563,6 +595,9 @@ void Foam::triSurfaceMesh::findLineAll
 
     info.setSize(start.size());
 
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
+
     // Work array
     DynamicList<pointIndexHit, 1, 1> hits;
 
@@ -576,7 +611,7 @@ void Foam::triSurfaceMesh::findLineAll
     const scalarField magSqrDirVec(magSqr(dirVec));
     const vectorField smallVec
     (
-        Foam::sqrt(SMALL)*dirVec
+        indexedOctree<treeDataTriSurface>::perturbTol()*dirVec
       + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
     );
 
@@ -619,6 +654,8 @@ void Foam::triSurfaceMesh::findLineAll
             info[pointI].clear();
         }
     }
+
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 }
 
 
@@ -655,7 +692,13 @@ void Foam::triSurfaceMesh::getNormal
     {
         if (info[i].hit())
         {
-            normal[i] = faceNormals()[info[i].index()];
+            label triI = info[i].index();
+            //- Cached:
+            //normal[i] = faceNormals()[triI];
+
+            //- Uncached
+            normal[i] = triSurface::operator[](triI).normal(points());
+            normal[i] /= mag(normal[i]) + VSMALL;
         }
         else
         {
@@ -697,6 +740,9 @@ void Foam::triSurfaceMesh::getVolumeType
 {
     volType.setSize(points.size());
 
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
+
     forAll(points, pointI)
     {
         const point& pt = points[pointI];
@@ -705,6 +751,8 @@ void Foam::triSurfaceMesh::getVolumeType
         // - cheat conversion since same values
         volType[pointI] = static_cast<volumeType>(tree().getVolumeType(pt));
     }
+
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 }
 
 

src/meshTools/searchableSurface/triSurfaceMesh.H

@@ -28,6 +28,11 @@ Class
 Description
     IOoject and searching on triSurface
 
+    Note: when constructing from dictionary has optional parameters:
+        - scale     : scaling factor.
+        - tolerance : relative tolerance for doing intersections
+                      (see triangle::intersection)
+
 SourceFiles
     triSurfaceMesh.C
 
@@ -63,6 +68,9 @@ private:
 
     // Private member data
 
+        //- Optional tolerance to use in searches
+        scalar tolerance_;
+
         //- Search tree (triangles)
         mutable autoPtr<indexedOctree<treeDataTriSurface> > tree_;