diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict
index 92b1c02e35..88c86bf482 100644
--- a/applications/utilities/postProcessing/sampling/sample/sampleDict
+++ b/applications/utilities/postProcessing/sampling/sample/sampleDict
@@ -31,18 +31,66 @@ setFormat raw;
 //      dx          : DX scalar or vector format
 //      vtk         : VTK ascii format
 //      raw         : x y z value format for use with e.g. gnuplot 'splot'.
+//      boundaryData: written in a form that can be used with timeVaryingMapped
+//                    boundary conditions (move to constant/boundaryData)
+//      starcd      : Star-CD format (geometry in .cel/.vrt, data in .usr)
+//      nastran     : Nastran (.nas) format
 //
 // Note:
-// other formats such as obj, stl, etc can also be written (by proxy)
-// but without any values!
+// - other formats such as obj, stl, etc can also be written (by proxy)
+//   but without any values!
+// - boundaryData format: typical use:
+//      - use a surfaces functionObject to sample a patch:
+//            type            surfaces;
+//            surfaceFormat   boundaryData;
+//            fields          ( p );
+//            surfaces
+//            (
+//                inlet
+//                {
+//                    type            patch;
+//                    patches         (inpletPatch);
+//                    interpolate     false;
+//                }
+//            );
+//      - the boundaryData writer will write postProcessing/surfaces/inlet/
+//      - move this to constant/boundaryData/inlet of the destination case
+//      - use a timeVaryingMappedFixedValue bc to read&interpolate and use
+//        as fixedValue.
+// - boundaryData: 'interpolate false' writes face centres, 'interpolate true'
+//   writes points. For 2D case the face centres might only be a single column
+//   so cannot be used in timeVaryingMapped with standard mapping (since
+//   uses triangulation to do interpolation).
+
 surfaceFormat vtk;
 
-// optionally define extra controls for the output formats
+// Optionally define extra controls for the output formats
 formatOptions
 {
     ensight
     {
+        // ascii/binary format
         format  ascii;
+        //collateTimes true;  // write single file containing multiple timesteps
+                              // (only for static surfaces)
+    }
+    nastran
+    {
+        // From OpenFOAM field name to Nastran field name
+        fields
+        (
+            (U      PLOAD4)
+            (p      PLOAD2)
+        );
+        // Optional scale
+        scale 1.0;
+        // Optional format
+        // short/long/free format
+        // short: scalar field width 8. default.
+        // long : scalar field width 16
+        // free : comma separated, field width according to writePrecision
+        //        setting
+        format  short;      // short/long/free
     }
 }
 
@@ -143,9 +191,19 @@ sets
         type        patchSeed;
         axis        xyz;
         patches     (".*Wall.*");
-        // Number of points to seed. Divided amongst all processors according
-        // to fraction of patches they hold.
+
+        // Subset patch faces by:
+
+        // 1. Number of points to seed. Divided amongst all processors
+        //   according to fraction of patches they hold.
         maxPoints   100;
+
+        // 2. Specified set of locations. This selects for every the specified
+        //    point the nearest patch face. (in addition the number of points
+        //    is also truncated by the maxPoints setting)
+        //    The difference with patchCloud is that this selects patch
+        //    face centres, not an arbitrary location on the face.
+        points   ((0.049 0.099 0.005)(0.051 0.054 0.005));
     }
 
 );
@@ -249,8 +307,10 @@ surfaces
         isoValue        0.5;
         interpolate     true;
 
-        //zone            ABC;          // Optional: zone only
-        //exposedPatchName fixedWalls;  // Optional: zone only
+        // zone            ABC;         // Optional: zone only
+        // exposedPatchName fixedWalls; // Optional: zone only
+
+        // bounds       (1 1 1)(2 2 2); // Optional: limit extent
 
         // regularise      false;    // Optional: do not simplify
         // mergeTol        1e-10;    // Optional: fraction of mesh bounding box
@@ -265,6 +325,9 @@ surfaces
         isoValue        0.5;
         interpolate     false;
         regularise      false;              // do not simplify
+
+        // bounds       (1 1 1)(2 2 2); // Optional: limit extent
+
         // mergeTol        1e-10;    // Optional: fraction of mesh bounding box
                                      // to merge points (default=1e-6)
     }
@@ -281,8 +344,10 @@ surfaces
         }
         interpolate     true;
 
-        //zone            ABC;          // Optional: zone only
-        //exposedPatchName fixedWalls;  // Optional: zone only
+        // zone             ABC;          // Optional: zone only
+        // exposedPatchName fixedWalls;  // Optional: zone only
+
+        // bounds       (1 1 1)(2 2 2); // Optional: limit extent
 
         // regularise      false;    // Optional: do not simplify
         // mergeTol        1e-10;    // Optional: fraction of mesh bounding box
@@ -305,6 +370,8 @@ surfaces
                                 // of isoSurfaceCell
         interpolate     false;
         regularise      false;  // Optional: do not simplify
+        // bounds       (1 1 1)(2 2 2); // Optional: limit extent
+
         // mergeTol 1e-10;      // Optional: fraction of mesh bounding box
                                 // to merge points (default=1e-6)
     }
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C
new file mode 100644
index 0000000000..8085ccd8f0
--- /dev/null
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C
@@ -0,0 +1,242 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "triangle.H"
+#include "triPoints.H"
+#include "plane.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Point, class PointRef>
+template<class AboveOp, class BelowOp>
+inline void Foam::triangle<Point, PointRef>::triSliceWithPlane
+(
+    const plane& pl,
+    const triPoints& tri,
+    AboveOp& aboveOp,
+    BelowOp& belowOp
+)
+{
+    // distance to cutting plane
+    FixedList<scalar, 3> d;
+
+    // determine how many of the points are above the cutting plane
+    label nPos = 0;
+    label posI = -1;
+    label negI = -1;
+    forAll(tri, i)
+    {
+        d[i] = ((tri[i] - pl.refPoint()) & pl.normal());
+
+        if (d[i] > 0)
+        {
+            nPos++;
+            posI = i;
+        }
+        else
+        {
+            negI = i;
+        }
+    }
+
+    if (nPos == 3)
+    {
+        aboveOp(tri);
+    }
+    else if (nPos == 2)
+    {
+        // point under the plane
+        label i0 = negI;
+
+        // indices of remaining points
+        label i1 = d.fcIndex(i0);
+        label i2 = d.fcIndex(i1);
+
+        // determine the two intersection points
+        point p01 = planeIntersection(d, tri, i0, i1);
+        point p02 = planeIntersection(d, tri, i0, i2);
+
+        aboveOp(triPoints(tri[i1], tri[i2], p02));
+        aboveOp(triPoints(tri[i1], p02, p01));
+        belowOp(triPoints(tri[i0], p01, p02));
+    }
+    else if (nPos == 1)
+    {
+        // point above the plane
+        label i0 = posI;
+
+        // indices of remaining points
+        label i1 = d.fcIndex(i0);
+        label i2 = d.fcIndex(i1);
+
+        // determine the two intersection points
+        point p01 = planeIntersection(d, tri, i0, i1);
+        point p02 = planeIntersection(d, tri, i0, i2);
+
+        belowOp(triPoints(tri[i1], tri[i2], p02));
+        belowOp(triPoints(tri[i1], p02, p01));
+        aboveOp(triPoints(tri[i0], p01, p02));
+    }
+    else
+    {
+        // All below
+        belowOp(tri);
+    }
+}
+
+
+template<class Point, class PointRef>
+template<class AboveOp, class BelowOp>
+inline void Foam::triangle<Point, PointRef>::sliceWithPlane
+(
+    const plane& pl,
+    AboveOp& aboveOp,
+    BelowOp& belowOp
+) const
+{
+    triSliceWithPlane(pl, triPoints(a_, b_, c_), aboveOp, belowOp);
+}
+
+
+template<class Point, class PointRef>
+template<class InsideOp, class OutsideOp>
+inline void Foam::triangle<Point, PointRef>::triangleOverlap
+(
+    const vector& n,
+    const triangle<Point, PointRef>& tgt,
+    InsideOp& insideOp,
+    OutsideOp& outsideOp
+) const
+{
+    // There are two possibilities with this algorithm - we either cut
+    // the outside triangles with all the edges or not (and keep them
+    // as disconnected triangles). In the first case
+    // we cannot do any evaluation short cut so we've chosen not to re-cut
+    // the outside triangles.
+
+
+    triIntersectionList insideTrisA;
+    label nInsideA = 0;
+    storeOp insideOpA(insideTrisA, nInsideA);
+
+    triIntersectionList outsideTrisA;
+    label nOutsideA = 0;
+    storeOp outsideOpA(outsideTrisA, nOutsideA);
+
+
+    const triPoints thisTri(a_, b_, c_);
+
+
+    // Cut original triangle with tgt edge 0.
+    // From *this to insideTrisA, outsideTrisA.
+    {
+        scalar s = Foam::mag(tgt.b() - tgt.a());
+        const plane pl0(tgt.a(), tgt.b(), tgt.b() + s*n);
+        triSliceWithPlane(pl0, thisTri, insideOpA, outsideOpA);
+    }
+
+    // Shortcut if nothing cut
+
+    if (insideOpA.nTris_ == 0)
+    {
+        outsideOp(thisTri);
+        return;
+    }
+
+    if (outsideOpA.nTris_ == 0)
+    {
+        insideOp(thisTri);
+        return;
+    }
+
+
+    // Cut all triangles with edge 1.
+    // From insideTrisA to insideTrisB, outsideTrisA
+
+    triIntersectionList insideTrisB;
+    label nInsideB = 0;
+    storeOp insideOpB(insideTrisB, nInsideB);
+
+    //triIntersectionList outsideTrisB;
+    //label nOutsideB = 0;
+    //storeOp outsideOpB(outsideTrisB, nOutsideB);
+
+    {
+        scalar s = Foam::mag(tgt.c() - tgt.b());
+        const plane pl0(tgt.b(), tgt.c(), tgt.c() + s*n);
+
+        for (label i = 0; i < insideOpA.nTris_; i++)
+        {
+            const triPoints& tri = insideOpA.tris_[i];
+            triSliceWithPlane(pl0, tri, insideOpB, outsideOpA);
+        }
+
+        //// Recut outside triangles (not necessary if only interested in
+        //// intersection properties)
+        //for (label i = 0; i < outsideOpA.nTris_; i++)
+        //{
+        //    const triPoints& tri = outsideOpA.tris_[i];
+        //    triSliceWithPlane(pl0, tri, outsideOpB, outsideOpB);
+        //}
+    }
+
+
+    // Cut all triangles with edge 2.
+    // From insideTrisB to insideTrisA, outsideTrisA
+    {
+        scalar s = Foam::mag(tgt.a() - tgt.c());
+        const plane pl0(tgt.c(), tgt.a(), tgt.a() + s*n);
+
+        insideOpA.nTris_ = 0;
+        //outsideOpA.nTris_ = 0;
+        for (label i = 0; i < insideOpB.nTris_; i++)
+        {
+            const triPoints& tri = insideOpB.tris_[i];
+            triSliceWithPlane(pl0, tri, insideOpA, outsideOpA);
+        }
+
+        //// Recut outside triangles (not necessary if only interested in
+        //// intersection properties)
+        //for (label i = 0; i < outsideOpB.nTris_; i++)
+        //{
+        //    const triPoints& tri = outsideOpB.tris_[i];
+        //    triSliceWithPlane(pl0, tri, outsideOpA, outsideOpA);
+        //}
+    }
+
+    // Transfer from A to argument
+    for (label i = 0; i < insideOpA.nTris_; i++)
+    {
+        insideOp(insideOpA.tris_[i]);
+    }
+
+    for (label i = 0; i < outsideOpA.nTris_; i++)
+    {
+        outsideOp(outsideOpA.tris_[i]);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
index 6e8030c797..de477d20e0 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.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-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,6 +52,8 @@ namespace Foam
 
 class Istream;
 class Ostream;
+class triPoints;
+class plane;
 
 // Forward declaration of friend functions and operators
 
@@ -71,6 +73,8 @@ inline Ostream& operator<<
     const triangle<Point, PointRef>&
 );
 
+typedef triangle<point, const point&> triPointRef;
+
 
 /*---------------------------------------------------------------------------*\
                           Class triangle Declaration
@@ -79,27 +83,98 @@ inline Ostream& operator<<
 template<class Point, class PointRef>
 class triangle
 {
+public:
+
+    // Public typedefs
+
+        //- Storage type for triangles originating from intersecting triangle
+        //  with another triangle
+        typedef FixedList<triPoints, 27> triIntersectionList;
+
+        //- Return types for classify
+        enum proxType
+        {
+            NONE,
+            POINT,  // Close to point
+            EDGE    // Close to edge
+        };
+
+
+    // Public classes
+
+        // Classes for use in sliceWithPlane. What to do with decomposition
+        // of triangle.
+
+            //- Dummy
+            class dummyOp
+            {
+            public:
+                inline void operator()(const triPoints&);
+            };
+
+            //- Sum resulting areas
+            class sumAreaOp
+            {
+            public:
+                scalar area_;
+
+                inline sumAreaOp();
+
+                inline void operator()(const triPoints&);
+            };
+
+            //- Store resulting tris
+            class storeOp
+            {
+            public:
+                triIntersectionList& tris_;
+                label& nTris_;
+
+                inline storeOp(triIntersectionList&, label&);
+
+                inline void operator()(const triPoints&);
+            };
+
+
+private:
+
     // Private data
 
         PointRef a_, b_, c_;
 
 
+   // Private Member Functions
+
+        //- Helper: calculate intersection point
+        inline static point planeIntersection
+        (
+            const FixedList<scalar, 3>& d,
+            const triPoints& t,
+            const label negI,
+            const label posI
+        );
+
+        //- Helper: slice triangle with plane
+        template<class AboveOp, class BelowOp>
+        inline static void triSliceWithPlane
+        (
+            const plane& pl,
+            const triPoints& tri,
+            AboveOp& aboveOp,
+            BelowOp& belowOp
+        );
+
+
 public:
 
-    //- Return types for classify
-    enum proxType
-    {
-        NONE,
-        POINT,  // Close to point
-        EDGE    // Close to edge
-    };
-
-
     // Constructors
 
         //- Construct from three points
         inline triangle(const Point& a, const Point& b, const Point& c);
 
+        //- Construct from three points
+        inline triangle(const FixedList<Point, 3>&);
+
         //- Construct from three points in the list of points
         //  The indices could be from triFace etc.
         inline triangle
@@ -237,6 +312,27 @@ public:
                 pointHit& edgePoint
             ) const;
 
+            //- Decompose triangle into triangles above and below plane
+            template<class AboveOp, class BelowOp>
+            inline void sliceWithPlane
+            (
+                const plane& pl,
+                AboveOp& aboveOp,
+                BelowOp& belowOp
+            ) const;
+
+            //- Decompose triangle into triangles inside and outside
+            //  (with respect to user provided normal) other
+            //  triangle.
+            template<class InsideOp, class OutsideOp>
+            inline void triangleOverlap
+            (
+                const vector& n,
+                const triangle<Point, PointRef>& tri,
+                InsideOp& insideOp,
+                OutsideOp& outsideOp
+            ) const;
+
 
     // IOstream operators
 
@@ -264,6 +360,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "triangle.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
index fbc5a14286..e3f8993e2e 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.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-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,7 @@ License
 
 #include "IOstreams.H"
 #include "pointHit.H"
+#include "triPoints.H"
 #include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@@ -43,6 +44,18 @@ inline Foam::triangle<Point, PointRef>::triangle
 {}
 
 
+template<class Point, class PointRef>
+inline Foam::triangle<Point, PointRef>::triangle
+(
+    const FixedList<Point, 3>& tri
+)
+:
+    a_(tri[0]),
+    b_(tri[1]),
+    c_(tri[2])
+{}
+
+
 template<class Point, class PointRef>
 inline Foam::triangle<Point, PointRef>::triangle
 (
@@ -804,6 +817,66 @@ inline Foam::pointHit Foam::triangle<Point, PointRef>::nearestPoint
 }
 
 
+template<class Point, class PointRef>
+inline void Foam::triangle<Point, PointRef>::dummyOp::operator()
+(
+    const triPoints&
+)
+{}
+
+
+template<class Point, class PointRef>
+inline Foam::triangle<Point, PointRef>::sumAreaOp::sumAreaOp()
+:
+    area_(0.0)
+{}
+
+
+template<class Point, class PointRef>
+inline void Foam::triangle<Point, PointRef>::sumAreaOp::operator()
+(
+    const triPoints& tri
+)
+{
+    area_ += triangle<Point, const Point&>(tri).mag();
+}
+
+
+template<class Point, class PointRef>
+inline Foam::triangle<Point, PointRef>::storeOp::storeOp
+(
+    triIntersectionList& tris,
+    label& nTris
+)
+:
+    tris_(tris),
+    nTris_(nTris)
+{}
+
+
+template<class Point, class PointRef>
+inline void Foam::triangle<Point, PointRef>::storeOp::operator()
+(
+    const triPoints& tri
+)
+{
+    tris_[nTris_++] = tri;
+}
+
+
+template<class Point, class PointRef>
+inline Foam::point Foam::triangle<Point, PointRef>::planeIntersection
+(
+    const FixedList<scalar, 3>& d,
+    const triPoints& t,
+    const label negI,
+    const label posI
+)
+{
+    return (d[posI]*t[negI] - d[negI]*t[posI])/(-d[negI] + d[posI]);
+}
+
+
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 template<class Point, class PointRef>
diff --git a/src/mesh/extrudeModel/planeExtrusion/planeExtrusion.H b/src/mesh/extrudeModel/planeExtrusion/planeExtrusion.H
index 9fea06cf2c..b832f7e2be 100644
--- a/src/mesh/extrudeModel/planeExtrusion/planeExtrusion.H
+++ b/src/mesh/extrudeModel/planeExtrusion/planeExtrusion.H
@@ -33,8 +33,8 @@ SeeAlso
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef plane_H
-#define plane_H
+#ifndef planeExtrusion_H
+#define planeExtrusion_H
 
 #include "linearNormal.H"
 
diff --git a/src/sampling/sampledSet/patchSeed/patchSeedSet.C b/src/sampling/sampledSet/patchSeed/patchSeedSet.C
index 3560fe4204..9cdcc7d230 100644
--- a/src/sampling/sampledSet/patchSeed/patchSeedSet.C
+++ b/src/sampling/sampledSet/patchSeed/patchSeedSet.C
@@ -30,9 +30,8 @@ License
 #include "treeDataFace.H"
 #include "Time.H"
 #include "meshTools.H"
-//#include "Random.H"
-// For 'facePoint' helper function only
 #include "mappedPatchBase.H"
+#include "indirectPrimitivePatch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -85,19 +84,145 @@ void Foam::patchSeedSet::calcSamples
     }
 
 
-    label totalSize = returnReduce(sz, sumOp<label>());
+    if (!rndGenPtr_.valid())
+    {
+        rndGenPtr_.reset(new Random(0));
+    }
+    Random& rndGen = rndGenPtr_();
+
+
+    if (selectedLocations_.size())
+    {
+        DynamicList<label> newPatchFaces(patchFaces.size());
+
+        // Find the nearest patch face
+        {
+            // 1. All processors find nearest local patch face for all
+            //    selectedLocations
+
+            // All the info for nearest. Construct to miss
+            List<mappedPatchBase::nearInfo> nearest(selectedLocations_.size());
+
+            const indirectPrimitivePatch pp
+            (
+                IndirectList<face>(mesh().faces(), patchFaces),
+                mesh().points()
+            );
+
+            treeBoundBox patchBb
+            (
+                treeBoundBox(pp.points(), pp.meshPoints()).extend
+                (
+                    rndGen,
+                    1e-4
+                )
+            );
+            patchBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+            patchBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+            indexedOctree<treeDataFace> boundaryTree
+            (
+                treeDataFace    // all information needed to search faces
+                (
+                    false,      // do not cache bb
+                    mesh(),
+                    patchFaces  // boundary faces only
+                ),
+                patchBb,        // overall search domain
+                8,              // maxLevel
+                10,             // leafsize
+                3.0             // duplicity
+            );
+
+            // Get some global dimension so all points are equally likely
+            // to be found
+            const scalar globalDistSqr
+            (
+                //magSqr
+                //(
+                //    boundBox
+                //    (
+                //        pp.points(),
+                //        pp.meshPoints(),
+                //        true
+                //    ).span()
+                //)
+                GREAT
+            );
+
+            forAll(selectedLocations_, sampleI)
+            {
+                const point& sample = selectedLocations_[sampleI];
+
+                pointIndexHit& nearInfo = nearest[sampleI].first();
+                nearInfo = boundaryTree.findNearest
+                (
+                    sample,
+                    globalDistSqr
+                );
+
+                if (!nearInfo.hit())
+                {
+                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                    nearest[sampleI].second().second() =
+                        Pstream::myProcNo();
+                }
+                else
+                {
+                    point fc(pp[nearInfo.index()].centre(pp.points()));
+                    nearInfo.setPoint(fc);
+                    nearest[sampleI].second().first() = magSqr(fc-sample);
+                    nearest[sampleI].second().second() =
+                        Pstream::myProcNo();
+                }
+            }
+
+
+            // 2. Reduce on master. Select nearest processor.
+
+            // Find nearest. Combine on master.
+            Pstream::listCombineGather(nearest, mappedPatchBase::nearestEqOp());
+            Pstream::listCombineScatter(nearest);
+
+
+            // 3. Pick up my local faces that have won
+
+            forAll(nearest, sampleI)
+            {
+                if (nearest[sampleI].first().hit())
+                {
+                    label procI = nearest[sampleI].second().second();
+                    label index = nearest[sampleI].first().index();
+
+                    if (procI == Pstream::myProcNo())
+                    {
+                        newPatchFaces.append(pp.addressing()[index]);
+                    }
+                }
+            }
+        }
+
+        if (debug)
+        {
+            Pout<< "Found " << newPatchFaces.size()
+                << " out of " << selectedLocations_.size()
+                << " on local processor" << endl;
+        }
+
+        patchFaces.transfer(newPatchFaces);
+    }
 
 
     // Shuffle and truncate if in random mode
+    label totalSize = returnReduce(patchFaces.size(), sumOp<label>());
+
     if (maxPoints_ < totalSize)
     {
         // Check what fraction of maxPoints_ I need to generate locally.
-        label myMaxPoints = label(scalar(sz)/totalSize*maxPoints_);
+        label myMaxPoints =
+            label(scalar(patchFaces.size())/totalSize*maxPoints_);
 
-        rndGenPtr_.reset(new Random(123456));
-        Random& rndGen = rndGenPtr_();
-
-        labelList subset = identity(sz);
+        labelList subset = identity(patchFaces.size());
         for (label iter = 0; iter < 4; iter++)
         {
             forAll(subset, i)
@@ -115,7 +240,7 @@ void Foam::patchSeedSet::calcSamples
         if (debug)
         {
             Pout<< "In random mode : selected " << patchFaces.size()
-                << " faces out of " << sz << endl;
+                << " faces out of " << patchFaces.size() << endl;
         }
     }
 
@@ -135,6 +260,9 @@ void Foam::patchSeedSet::calcSamples
     forAll(patchFaces, i)
     {
         label faceI = patchFaces[i];
+
+        // Slightly shift point in since on warped face face-diagonal
+        // decomposition might be outside cell for face-centre decomposition!
         pointIndexHit info = mappedPatchBase::facePoint
         (
             mesh(),
@@ -217,9 +345,15 @@ Foam::patchSeedSet::patchSeedSet
             wordReList(dict.lookup("patches"))
         )
     ),
-    //searchDist_(readScalar(dict.lookup("maxDistance"))),
-    //offsetDist_(readScalar(dict.lookup("offsetDist"))),
-    maxPoints_(readLabel(dict.lookup("maxPoints")))
+    maxPoints_(readLabel(dict.lookup("maxPoints"))),
+    selectedLocations_
+    (
+        dict.lookupOrDefault<pointField>
+        (
+            "points",
+            pointField(0)
+        )
+    )
 {
     genSamples();
 
diff --git a/src/sampling/sampledSet/patchSeed/patchSeedSet.H b/src/sampling/sampledSet/patchSeed/patchSeedSet.H
index a6e795ad2b..c2143a2098 100644
--- a/src/sampling/sampledSet/patchSeed/patchSeedSet.H
+++ b/src/sampling/sampledSet/patchSeed/patchSeedSet.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,18 +58,15 @@ class patchSeedSet
         //- Patches to sample
         const labelHashSet patchSet_;
 
-//        //- Maximum distance to look for nearest
-//        const scalar searchDist_;
-//
-//        //- Offset distance
-//        const scalar offsetDist_;
-
-        //- Maximum number of patch faces to seed
+        //- Maximum number of patch faces to seed (if in random subset mode)
         const label maxPoints_;
 
         //- Random number generator (if maxPoints < num patch faces)
         autoPtr<Random> rndGenPtr_;
 
+        //- Patch faces to seed selected based on nearness to supplied points
+        const pointField selectedLocations_;
+
 
     // Private Member Functions
 
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
index 2dafe3e98a..fd2927ca71 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
@@ -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-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -278,7 +278,8 @@ void Foam::distanceSurface::createGeometry()
                 cellDistance,
                 pointDistance_,
                 distance_,
-                regularise_
+                regularise_,
+                bounds_
             )
         );
     }
@@ -291,7 +292,8 @@ void Foam::distanceSurface::createGeometry()
                 cellDistance,
                 pointDistance_,
                 distance_,
-                regularise_
+                regularise_,
+                bounds_
             )
         );
     }
@@ -336,6 +338,7 @@ Foam::distanceSurface::distanceSurface
     cell_(dict.lookupOrDefault("cell", true)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", false)),
+    bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
     zoneKey_(keyType::null),
     needsUpdate_(true),
     isoSurfCellPtr_(NULL),
@@ -364,7 +367,8 @@ Foam::distanceSurface::distanceSurface
     const bool signedDistance,
     const bool cell,
     const Switch regularise,
-    const Switch average
+    const Switch average,
+    const boundBox& bounds
 )
 :
     sampledSurface(name, mesh, interpolate),
@@ -390,6 +394,7 @@ Foam::distanceSurface::distanceSurface
     cell_(cell),
     regularise_(regularise),
     average_(average),
+    bounds_(bounds),
     zoneKey_(keyType::null),
     needsUpdate_(true),
     isoSurfCellPtr_(NULL),
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
index 73cd555e05..aa3f632511 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
@@ -75,6 +75,9 @@ class distanceSurface
         //- Whether to recalculate cell values as average of point values
         const Switch average_;
 
+        //- Optional bounding box to trim triangles against
+        const boundBox bounds_;
+
         //- If restricted to zones, name of this zone or a regular expression
         keyType zoneKey_;
 
@@ -144,7 +147,8 @@ public:
             const bool signedDistance,
             const bool cell,
             const Switch regularise,
-            const Switch average
+            const Switch average,
+            const boundBox& bounds = boundBox::greatBox
         );
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index e3105b35a5..09d0a819d4 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -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-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,12 +32,36 @@ License
 #include "surfaceFields.H"
 #include "OFstream.H"
 #include "meshTools.H"
+#include "triSurfaceSearch.H"
+#include "surfaceIntersection.H"
+#include "intersectedSurface.H"
+#include "searchableBox.H"
+#include "triSurfaceMesh.H"
+#include "triPoints.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
 defineTypeNameAndDebug(isoSurface, 0);
+
+// Helper class for slicing triangles
+class storeOp
+{
+public:
+    DynamicList<triPoints>& tris_;
+
+    inline storeOp(DynamicList<triPoints>& tris)
+    :
+        tris_(tris)
+    {}
+
+    inline void operator()(const triPoints& tri)
+    {
+        tris_.append(tri);
+    }
+};
+
 }
 
 
@@ -473,41 +497,41 @@ void Foam::isoSurface::calcCutTypes
 }
 
 
-// Return the two common points between two triangles
-Foam::labelPair Foam::isoSurface::findCommonPoints
-(
-    const labelledTri& tri0,
-    const labelledTri& tri1
-)
-{
-    labelPair common(-1, -1);
-
-    label fp0 = 0;
-    label fp1 = findIndex(tri1, tri0[fp0]);
-
-    if (fp1 == -1)
-    {
-        fp0 = 1;
-        fp1 = findIndex(tri1, tri0[fp0]);
-    }
-
-    if (fp1 != -1)
-    {
-        // So tri0[fp0] is tri1[fp1]
-
-        // Find next common point
-        label fp0p1 = tri0.fcIndex(fp0);
-        label fp1p1 = tri1.fcIndex(fp1);
-        label fp1m1 = tri1.rcIndex(fp1);
-
-        if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
-        {
-            common[0] = tri0[fp0];
-            common[1] = tri0[fp0p1];
-        }
-    }
-    return common;
-}
+//// Return the two common points between two triangles
+//Foam::labelPair Foam::isoSurface::findCommonPoints
+//(
+//    const labelledTri& tri0,
+//    const labelledTri& tri1
+//)
+//{
+//    labelPair common(-1, -1);
+//
+//    label fp0 = 0;
+//    label fp1 = findIndex(tri1, tri0[fp0]);
+//
+//    if (fp1 == -1)
+//    {
+//        fp0 = 1;
+//        fp1 = findIndex(tri1, tri0[fp0]);
+//    }
+//
+//    if (fp1 != -1)
+//    {
+//        // So tri0[fp0] is tri1[fp1]
+//
+//        // Find next common point
+//        label fp0p1 = tri0.fcIndex(fp0);
+//        label fp1p1 = tri1.fcIndex(fp1);
+//        label fp1m1 = tri1.rcIndex(fp1);
+//
+//        if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
+//        {
+//            common[0] = tri0[fp0];
+//            common[1] = tri0[fp0p1];
+//        }
+//    }
+//    return common;
+//}
 
 
 // Caculate centre of surface.
@@ -523,71 +547,71 @@ Foam::point Foam::isoSurface::calcCentre(const triSurface& s)
 }
 
 
-// Replace surface (localPoints, localTris) with single point. Returns
-// point. Destructs arguments.
-Foam::pointIndexHit Foam::isoSurface::collapseSurface
-(
-    pointField& localPoints,
-    DynamicList<labelledTri, 64>& localTris
-)
-{
-    pointIndexHit info(false, vector::zero, localTris.size());
-
-    if (localTris.size() == 1)
-    {
-        const labelledTri& tri = localTris[0];
-        info.setPoint(tri.centre(localPoints));
-        info.setHit();
-    }
-    else if (localTris.size() == 2)
-    {
-        // Check if the two triangles share an edge.
-        const labelledTri& tri0 = localTris[0];
-        const labelledTri& tri1 = localTris[0];
-
-        labelPair shared = findCommonPoints(tri0, tri1);
-
-        if (shared[0] != -1)
-        {
-            info.setPoint
-            (
-                0.5
-              * (
-                    tri0.centre(localPoints)
-                  + tri1.centre(localPoints)
-                )
-            );
-            info.setHit();
-        }
-    }
-    else if (localTris.size())
-    {
-        // Check if single region. Rare situation.
-        triSurface surf
-        (
-            localTris,
-            geometricSurfacePatchList(0),
-            localPoints,
-            true
-        );
-        localTris.clearStorage();
-
-        labelList faceZone;
-        label nZones = surf.markZones
-        (
-            boolList(surf.nEdges(), false),
-            faceZone
-        );
-
-        if (nZones == 1)
-        {
-            info.setPoint(calcCentre(surf));
-            info.setHit();
-        }
-    }
-
-    return info;
-}
+//// Replace surface (localPoints, localTris) with single point. Returns
+//// point. Destructs arguments.
+//Foam::pointIndexHit Foam::isoSurface::collapseSurface
+//(
+//    pointField& localPoints,
+//    DynamicList<labelledTri, 64>& localTris
+//)
+//{
+//    pointIndexHit info(false, vector::zero, localTris.size());
+//
+//    if (localTris.size() == 1)
+//    {
+//        const labelledTri& tri = localTris[0];
+//        info.setPoint(tri.centre(localPoints));
+//        info.setHit();
+//    }
+//    else if (localTris.size() == 2)
+//    {
+//        // Check if the two triangles share an edge.
+//        const labelledTri& tri0 = localTris[0];
+//        const labelledTri& tri1 = localTris[0];
+//
+//        labelPair shared = findCommonPoints(tri0, tri1);
+//
+//        if (shared[0] != -1)
+//        {
+//            info.setPoint
+//            (
+//                0.5
+//              * (
+//                    tri0.centre(localPoints)
+//                  + tri1.centre(localPoints)
+//                )
+//            );
+//            info.setHit();
+//        }
+//    }
+//    else if (localTris.size())
+//    {
+//        // Check if single region. Rare situation.
+//        triSurface surf
+//        (
+//            localTris,
+//            geometricSurfacePatchList(0),
+//            localPoints,
+//            true
+//        );
+//        localTris.clearStorage();
+//
+//        labelList faceZone;
+//        label nZones = surf.markZones
+//        (
+//            boolList(surf.nEdges(), false),
+//            faceZone
+//        );
+//
+//        if (nZones == 1)
+//        {
+//            info.setPoint(calcCentre(surf));
+//            info.setHit();
+//        }
+//    }
+//
+//    return info;
+//}
 
 
 // Determine per cell centre whether all the intersections get collapsed
@@ -1133,6 +1157,220 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
 }
 
 
+void Foam::isoSurface::trimToPlanes
+(
+    const PtrList<plane>& planes,
+    const triPointRef& tri,
+    DynamicList<point>& newTriPoints
+)
+{
+    // Buffer for generated triangles
+    DynamicList<triPoints> insideTrisA;
+    storeOp insideOpA(insideTrisA);
+
+    // Buffer for generated triangles
+    DynamicList<triPoints> insideTrisB;
+    storeOp insideOpB(insideTrisB);
+
+    triPointRef::dummyOp dop;
+
+    // Store starting triangle in insideTrisA
+    insideOpA(triPoints(tri.a(), tri.b(), tri.c()));
+
+
+    bool useA = true;
+
+    forAll(planes, faceI)
+    {
+        const plane& pl = planes[faceI];
+
+        if (useA)
+        {
+            insideOpB.tris_.clear();
+            forAll(insideOpA.tris_, i)
+            {
+                const triPoints& tri = insideOpA.tris_[i];
+                triPointRef(tri).sliceWithPlane(pl, insideOpB, dop);
+            }
+        }
+        else
+        {
+            insideOpA.tris_.clear();
+            forAll(insideOpB.tris_, i)
+            {
+                const triPoints& tri = insideOpB.tris_[i];
+                triPointRef(tri).sliceWithPlane(pl, insideOpA, dop);
+            }
+        }
+        useA = !useA;
+    }
+
+
+    // Transfer
+    if (useA)
+    {
+        forAll(insideOpA.tris_, i)
+        {
+            const triPoints& tri = insideOpA.tris_[i];
+            newTriPoints.append(tri[0]);
+            newTriPoints.append(tri[1]);
+            newTriPoints.append(tri[2]);
+        }
+    }
+    else
+    {
+        forAll(insideOpB.tris_, i)
+        {
+            const triPoints& tri = insideOpB.tris_[i];
+            newTriPoints.append(tri[0]);
+            newTriPoints.append(tri[1]);
+            newTriPoints.append(tri[2]);
+        }
+    }
+}
+
+
+void Foam::isoSurface::trimToBox
+(
+    const treeBoundBox& bb,
+    DynamicList<point>& triPoints,
+    DynamicList<label>& triMap      // trimmed to original tris
+)
+{
+    if (debug)
+    {
+        Pout<< "isoSurface : trimming to " << bb << endl;
+    }
+
+    // Generate inwards pointing planes
+    PtrList<plane> planes(6);
+    const pointField pts(bb.treeBoundBox::points());
+    forAll(treeBoundBox::faces, faceI)
+    {
+        const face& f = treeBoundBox::faces[faceI];
+        const vector& n = treeBoundBox::faceNormals[faceI];
+        planes.set(faceI, new plane(pts[f[0]], -n));
+    }
+
+    label nTris = triPoints.size()/3;
+
+    DynamicList<point> newTriPoints(triPoints.size()/16);
+    triMap.setCapacity(nTris/16);
+
+    label vertI = 0;
+    for (label triI = 0; triI < nTris; triI++)
+    {
+        const point& p0 = triPoints[vertI++];
+        const point& p1 = triPoints[vertI++];
+        const point& p2 = triPoints[vertI++];
+
+        label oldNPoints = newTriPoints.size();
+        trimToPlanes
+        (
+            planes,
+            triPointRef(p0, p1, p2),
+            newTriPoints
+        );
+
+        label nCells = (newTriPoints.size()-oldNPoints)/3;
+        for (label i = 0; i < nCells; i++)
+        {
+            triMap.append(triI);
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "isoSurface : trimmed from " << nTris
+            << " down to " << triMap.size()
+            << " triangles." << endl;
+    }
+
+    triPoints.transfer(newTriPoints);
+}
+
+
+void Foam::isoSurface::trimToBox
+(
+    const treeBoundBox& bb,
+    DynamicList<point>& triPoints,  // new points
+    DynamicList<label>& triMap,     // map from (new) triangle to original
+    labelList& triPointMap,         // map from (new) point to original
+    labelList& interpolatedPoints,  // labels of newly introduced points
+    List<FixedList<label, 3> >& interpolatedOldPoints,// and their interpolation
+    List<FixedList<scalar, 3> >& interpolationWeights
+)
+{
+    const List<point> oldTriPoints(triPoints);
+
+    // Trim triPoints, return map
+    trimToBox(bb, triPoints, triMap);
+
+
+    // Find point correspondence:
+    // - one-to-one for preserved original points (triPointMap)
+    // - interpolation for newly introduced points
+    //   (interpolatedOldPoints)
+    label sz = oldTriPoints.size()/100;
+    DynamicList<label> dynInterpolatedPoints(sz);
+    DynamicList<FixedList<label, 3> > dynInterpolatedOldPoints(sz);
+    DynamicList<FixedList<scalar, 3> > dynInterpolationWeights(sz);
+
+
+    triPointMap.setSize(triPoints.size());
+    forAll(triMap, triI)
+    {
+        label oldTriI = triMap[triI];
+
+        // Find point correspondence. Assumes coordinate is bit-copy.
+        for (label i = 0; i < 3; i++)
+        {
+            label pointI = 3*triI+i;
+            const point& pt = triPoints[pointI];
+
+            // Compare to old-triangle's points
+            label matchPointI = -1;
+            for (label j = 0; j < 3; j++)
+            {
+                label oldPointI = 3*oldTriI+j;
+                if (pt == oldTriPoints[oldPointI])
+                {
+                    matchPointI = oldPointI;
+                    break;
+                }
+            }
+
+            triPointMap[pointI] = matchPointI;
+
+            // If new point: calculate and store interpolation
+            if (matchPointI == -1)
+            {
+                dynInterpolatedPoints.append(pointI);
+
+                FixedList<label, 3> oldPoints;
+                oldPoints[0] = 3*oldTriI;
+                oldPoints[1] = 3*oldTriI+1;
+                oldPoints[2] = 3*oldTriI+2;
+                dynInterpolatedOldPoints.append(oldPoints);
+
+                triPointRef tri(oldTriPoints, oldPoints);
+                scalarList bary;
+                tri.barycentric(pt, bary);
+                FixedList<scalar, 3> weights;
+                weights[0] = bary[0];
+                weights[1] = bary[1];
+                weights[2] = bary[2];
+                dynInterpolationWeights.append(weights);
+            }
+        }
+    }
+
+    interpolatedPoints.transfer(dynInterpolatedPoints);
+    interpolatedOldPoints.transfer(dynInterpolatedOldPoints);
+    interpolationWeights.transfer(dynInterpolationWeights);
+}
+
+
 // Does face use valid vertices?
 bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
 {
@@ -1204,435 +1442,435 @@ bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
 }
 
 
-void Foam::isoSurface::calcAddressing
-(
-    const triSurface& surf,
-    List<FixedList<label, 3> >& faceEdges,
-    labelList& edgeFace0,
-    labelList& edgeFace1,
-    Map<labelList>& edgeFacesRest
-) const
-{
-    const pointField& points = surf.points();
-
-    pointField edgeCentres(3*surf.size());
-    label edgeI = 0;
-    forAll(surf, triI)
-    {
-        const labelledTri& tri = surf[triI];
-        edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
-        edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
-        edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
-    }
-
-    pointField mergedCentres;
-    labelList oldToMerged;
-    bool hasMerged = mergePoints
-    (
-        edgeCentres,
-        mergeDistance_,
-        false,
-        oldToMerged,
-        mergedCentres
-    );
-
-    if (debug)
-    {
-        Pout<< "isoSurface : detected "
-            << mergedCentres.size()
-            << " geometric edges on " << surf.size() << " triangles." << endl;
-    }
-
-    if (!hasMerged)
-    {
-        return;
-    }
-
-
-    // Determine faceEdges
-    faceEdges.setSize(surf.size());
-    edgeI = 0;
-    forAll(surf, triI)
-    {
-        faceEdges[triI][0] = oldToMerged[edgeI++];
-        faceEdges[triI][1] = oldToMerged[edgeI++];
-        faceEdges[triI][2] = oldToMerged[edgeI++];
-    }
-
-
-    // Determine edgeFaces
-    edgeFace0.setSize(mergedCentres.size());
-    edgeFace0 = -1;
-    edgeFace1.setSize(mergedCentres.size());
-    edgeFace1 = -1;
-    edgeFacesRest.clear();
-
-    // Overflow edge faces for geometric shared edges that turned
-    // out to be different anyway.
-    EdgeMap<labelList> extraEdgeFaces(mergedCentres.size()/100);
-
-    forAll(oldToMerged, oldEdgeI)
-    {
-        label triI = oldEdgeI / 3;
-        label edgeI = oldToMerged[oldEdgeI];
-
-        if (edgeFace0[edgeI] == -1)
-        {
-            // First triangle for edge
-            edgeFace0[edgeI] = triI;
-        }
-        else
-        {
-            //- Check that the two triangles actually topologically
-            //  share an edge
-            const labelledTri& prevTri = surf[edgeFace0[edgeI]];
-            const labelledTri& tri = surf[triI];
-
-            label fp = oldEdgeI % 3;
-
-            edge e(tri[fp], tri[tri.fcIndex(fp)]);
-
-            label prevTriIndex = -1;
-
-            forAll(prevTri, i)
-            {
-                if (edge(prevTri[i], prevTri[prevTri.fcIndex(i)]) == e)
-                {
-                    prevTriIndex = i;
-                    break;
-                }
-            }
-
-            if (prevTriIndex == -1)
-            {
-                // Different edge. Store for later.
-                EdgeMap<labelList>::iterator iter = extraEdgeFaces.find(e);
-
-                if (iter != extraEdgeFaces.end())
-                {
-                    labelList& eFaces = iter();
-                    label sz = eFaces.size();
-                    eFaces.setSize(sz+1);
-                    eFaces[sz] = triI;
-                }
-                else
-                {
-                    extraEdgeFaces.insert(e, labelList(1, triI));
-                }
-            }
-            else if (edgeFace1[edgeI] == -1)
-            {
-                edgeFace1[edgeI] = triI;
-            }
-            else
-            {
-                //WarningIn("orientSurface(triSurface&)")
-                //    << "Edge " << edgeI << " with centre "
-                //    << mergedCentres[edgeI]
-                //    << " used by more than two triangles: "
-                //    << edgeFace0[edgeI] << ", "
-                //    << edgeFace1[edgeI] << " and " << triI << endl;
-                Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
-
-                if (iter != edgeFacesRest.end())
-                {
-                    labelList& eFaces = iter();
-                    label sz = eFaces.size();
-                    eFaces.setSize(sz+1);
-                    eFaces[sz] = triI;
-                }
-                else
-                {
-                    edgeFacesRest.insert(edgeI, labelList(1, triI));
-                }
-            }
-        }
-    }
-
-    // Add extraEdgeFaces
-    edgeI = edgeFace0.size();
-
-    edgeFace0.setSize(edgeI + extraEdgeFaces.size());
-    edgeFace1.setSize(edgeI + extraEdgeFaces.size(), -1);
-
-    forAllConstIter(EdgeMap<labelList>, extraEdgeFaces, iter)
-    {
-        const labelList& eFaces = iter();
-
-        // The current edge will become edgeI. Replace all occurrences in
-        // faceEdges
-        forAll(eFaces, i)
-        {
-            label triI = eFaces[i];
-            const labelledTri& tri = surf[triI];
-
-            FixedList<label, 3>& fEdges = faceEdges[triI];
-            forAll(tri, fp)
-            {
-                edge e(tri[fp], tri[tri.fcIndex(fp)]);
-
-                if (e == iter.key())
-                {
-                    fEdges[fp] = edgeI;
-                    break;
-                }
-            }
-        }
-
-
-        // Add face to edgeFaces
-
-        edgeFace0[edgeI] = eFaces[0];
-
-        if (eFaces.size() >= 2)
-        {
-            edgeFace1[edgeI] = eFaces[1];
-
-            if (eFaces.size() > 2)
-            {
-                edgeFacesRest.insert
-                (
-                    edgeI,
-                    SubList<label>(eFaces, eFaces.size()-2, 2)
-                );
-            }
-        }
-
-        edgeI++;
-    }
-}
-
-
-void Foam::isoSurface::walkOrientation
-(
-    const triSurface& surf,
-    const List<FixedList<label, 3> >& faceEdges,
-    const labelList& edgeFace0,
-    const labelList& edgeFace1,
-    const label seedTriI,
-    labelList& flipState
-)
-{
-    // Do walk for consistent orientation.
-    DynamicList<label> changedFaces(surf.size());
-
-    changedFaces.append(seedTriI);
-
-    while (changedFaces.size())
-    {
-        DynamicList<label> newChangedFaces(changedFaces.size());
-
-        forAll(changedFaces, i)
-        {
-            label triI = changedFaces[i];
-            const labelledTri& tri = surf[triI];
-            const FixedList<label, 3>& fEdges = faceEdges[triI];
-
-            forAll(fEdges, fp)
-            {
-                label edgeI = fEdges[fp];
-
-                // my points:
-                label p0 = tri[fp];
-                label p1 = tri[tri.fcIndex(fp)];
-
-                label nbrI =
-                (
-                    edgeFace0[edgeI] != triI
-                  ? edgeFace0[edgeI]
-                  : edgeFace1[edgeI]
-                );
-
-                if (nbrI != -1 && flipState[nbrI] == -1)
-                {
-                    const labelledTri& nbrTri = surf[nbrI];
-
-                    // nbr points
-                    label nbrFp = findIndex(nbrTri, p0);
-
-                    if (nbrFp == -1)
-                    {
-                        FatalErrorIn("isoSurface::walkOrientation(..)")
-                            << "triI:" << triI
-                            << " tri:" << tri
-                            << " p0:" << p0
-                            << " p1:" << p1
-                            << " fEdges:" << fEdges
-                            << " edgeI:" << edgeI
-                            << " edgeFace0:" << edgeFace0[edgeI]
-                            << " edgeFace1:" << edgeFace1[edgeI]
-                            << " nbrI:" << nbrI
-                            << " nbrTri:" << nbrTri
-                            << abort(FatalError);
-                    }
-
-
-                    label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
-
-                    bool sameOrientation = (p1 == nbrP1);
-
-                    if (flipState[triI] == 0)
-                    {
-                        flipState[nbrI] = (sameOrientation ? 0 : 1);
-                    }
-                    else
-                    {
-                        flipState[nbrI] = (sameOrientation ? 1 : 0);
-                    }
-                    newChangedFaces.append(nbrI);
-                }
-            }
-        }
-
-        changedFaces.transfer(newChangedFaces);
-    }
-}
-
-
-void Foam::isoSurface::orientSurface
-(
-    triSurface& surf,
-    const List<FixedList<label, 3> >& faceEdges,
-    const labelList& edgeFace0,
-    const labelList& edgeFace1,
-    const Map<labelList>& edgeFacesRest
-)
-{
-    // -1 : unvisited
-    //  0 : leave as is
-    //  1 : flip
-    labelList flipState(surf.size(), -1);
-
-    label seedTriI = 0;
-
-    while (true)
-    {
-        // Find first unvisited triangle
-        for
-        (
-            ;
-            seedTriI < surf.size() && flipState[seedTriI] != -1;
-            seedTriI++
-        )
-        {}
-
-        if (seedTriI == surf.size())
-        {
-            break;
-        }
-
-        // Note: Determine orientation of seedTriI?
-        // for now assume it is ok
-        flipState[seedTriI] = 0;
-
-        walkOrientation
-        (
-            surf,
-            faceEdges,
-            edgeFace0,
-            edgeFace1,
-            seedTriI,
-            flipState
-        );
-    }
-
-    // Do actual flipping
-    surf.clearOut();
-    forAll(surf, triI)
-    {
-        if (flipState[triI] == 1)
-        {
-            labelledTri tri(surf[triI]);
-
-            surf[triI][0] = tri[0];
-            surf[triI][1] = tri[2];
-            surf[triI][2] = tri[1];
-        }
-        else if (flipState[triI] == -1)
-        {
-            FatalErrorIn
-            (
-                "isoSurface::orientSurface(triSurface&, const label)"
-            )   << "problem" << abort(FatalError);
-        }
-    }
-}
-
-
-// Checks if triangle is connected through edgeI only.
-bool Foam::isoSurface::danglingTriangle
-(
-    const FixedList<label, 3>& fEdges,
-    const labelList& edgeFace1
-)
-{
-    label nOpen = 0;
-    forAll(fEdges, i)
-    {
-        if (edgeFace1[fEdges[i]] == -1)
-        {
-            nOpen++;
-        }
-    }
-
-    if (nOpen == 1 || nOpen == 2 || nOpen == 3)
-    {
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-// Mark triangles to keep. Returns number of dangling triangles.
-Foam::label Foam::isoSurface::markDanglingTriangles
-(
-    const List<FixedList<label, 3> >& faceEdges,
-    const labelList& edgeFace0,
-    const labelList& edgeFace1,
-    const Map<labelList>& edgeFacesRest,
-    boolList& keepTriangles
-)
-{
-    keepTriangles.setSize(faceEdges.size());
-    keepTriangles = true;
-
-    label nDangling = 0;
-
-    // Remove any dangling triangles
-    forAllConstIter(Map<labelList>,  edgeFacesRest, iter)
-    {
-        // These are all the non-manifold edges. Filter out all triangles
-        // with only one connected edge (= this edge)
-
-        label edgeI = iter.key();
-        const labelList& otherEdgeFaces = iter();
-
-        // Remove all dangling triangles
-        if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
-        {
-            keepTriangles[edgeFace0[edgeI]] = false;
-            nDangling++;
-        }
-        if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
-        {
-            keepTriangles[edgeFace1[edgeI]] = false;
-            nDangling++;
-        }
-        forAll(otherEdgeFaces, i)
-        {
-            label triI = otherEdgeFaces[i];
-            if (danglingTriangle(faceEdges[triI], edgeFace1))
-            {
-                keepTriangles[triI] = false;
-                nDangling++;
-            }
-        }
-    }
-    return nDangling;
-}
+//void Foam::isoSurface::calcAddressing
+//(
+//    const triSurface& surf,
+//    List<FixedList<label, 3> >& faceEdges,
+//    labelList& edgeFace0,
+//    labelList& edgeFace1,
+//    Map<labelList>& edgeFacesRest
+//) const
+//{
+//    const pointField& points = surf.points();
+//
+//    pointField edgeCentres(3*surf.size());
+//    label edgeI = 0;
+//    forAll(surf, triI)
+//    {
+//        const labelledTri& tri = surf[triI];
+//        edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
+//        edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
+//        edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
+//    }
+//
+//    pointField mergedCentres;
+//    labelList oldToMerged;
+//    bool hasMerged = mergePoints
+//    (
+//        edgeCentres,
+//        mergeDistance_,
+//        false,
+//        oldToMerged,
+//        mergedCentres
+//    );
+//
+//    if (debug)
+//    {
+//        Pout<< "isoSurface : detected "
+//            << mergedCentres.size()
+//            << " geometric edges on " << surf.size() << " triangles." << endl;
+//    }
+//
+//    if (!hasMerged)
+//    {
+//        return;
+//    }
+//
+//
+//    // Determine faceEdges
+//    faceEdges.setSize(surf.size());
+//    edgeI = 0;
+//    forAll(surf, triI)
+//    {
+//        faceEdges[triI][0] = oldToMerged[edgeI++];
+//        faceEdges[triI][1] = oldToMerged[edgeI++];
+//        faceEdges[triI][2] = oldToMerged[edgeI++];
+//    }
+//
+//
+//    // Determine edgeFaces
+//    edgeFace0.setSize(mergedCentres.size());
+//    edgeFace0 = -1;
+//    edgeFace1.setSize(mergedCentres.size());
+//    edgeFace1 = -1;
+//    edgeFacesRest.clear();
+//
+//    // Overflow edge faces for geometric shared edges that turned
+//    // out to be different anyway.
+//    EdgeMap<labelList> extraEdgeFaces(mergedCentres.size()/100);
+//
+//    forAll(oldToMerged, oldEdgeI)
+//    {
+//        label triI = oldEdgeI / 3;
+//        label edgeI = oldToMerged[oldEdgeI];
+//
+//        if (edgeFace0[edgeI] == -1)
+//        {
+//            // First triangle for edge
+//            edgeFace0[edgeI] = triI;
+//        }
+//        else
+//        {
+//            //- Check that the two triangles actually topologically
+//            //  share an edge
+//            const labelledTri& prevTri = surf[edgeFace0[edgeI]];
+//            const labelledTri& tri = surf[triI];
+//
+//            label fp = oldEdgeI % 3;
+//
+//            edge e(tri[fp], tri[tri.fcIndex(fp)]);
+//
+//            label prevTriIndex = -1;
+//
+//            forAll(prevTri, i)
+//            {
+//                if (edge(prevTri[i], prevTri[prevTri.fcIndex(i)]) == e)
+//                {
+//                    prevTriIndex = i;
+//                    break;
+//                }
+//            }
+//
+//            if (prevTriIndex == -1)
+//            {
+//                // Different edge. Store for later.
+//                EdgeMap<labelList>::iterator iter = extraEdgeFaces.find(e);
+//
+//                if (iter != extraEdgeFaces.end())
+//                {
+//                    labelList& eFaces = iter();
+//                    label sz = eFaces.size();
+//                    eFaces.setSize(sz+1);
+//                    eFaces[sz] = triI;
+//                }
+//                else
+//                {
+//                    extraEdgeFaces.insert(e, labelList(1, triI));
+//                }
+//            }
+//            else if (edgeFace1[edgeI] == -1)
+//            {
+//                edgeFace1[edgeI] = triI;
+//            }
+//            else
+//            {
+//                //WarningIn("orientSurface(triSurface&)")
+//                //    << "Edge " << edgeI << " with centre "
+//                //    << mergedCentres[edgeI]
+//                //    << " used by more than two triangles: "
+//                //    << edgeFace0[edgeI] << ", "
+//                //    << edgeFace1[edgeI] << " and " << triI << endl;
+//                Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
+//
+//                if (iter != edgeFacesRest.end())
+//                {
+//                    labelList& eFaces = iter();
+//                    label sz = eFaces.size();
+//                    eFaces.setSize(sz+1);
+//                    eFaces[sz] = triI;
+//                }
+//                else
+//                {
+//                    edgeFacesRest.insert(edgeI, labelList(1, triI));
+//                }
+//            }
+//        }
+//    }
+//
+//    // Add extraEdgeFaces
+//    edgeI = edgeFace0.size();
+//
+//    edgeFace0.setSize(edgeI + extraEdgeFaces.size());
+//    edgeFace1.setSize(edgeI + extraEdgeFaces.size(), -1);
+//
+//    forAllConstIter(EdgeMap<labelList>, extraEdgeFaces, iter)
+//    {
+//        const labelList& eFaces = iter();
+//
+//        // The current edge will become edgeI. Replace all occurrences in
+//        // faceEdges
+//        forAll(eFaces, i)
+//        {
+//            label triI = eFaces[i];
+//            const labelledTri& tri = surf[triI];
+//
+//            FixedList<label, 3>& fEdges = faceEdges[triI];
+//            forAll(tri, fp)
+//            {
+//                edge e(tri[fp], tri[tri.fcIndex(fp)]);
+//
+//                if (e == iter.key())
+//                {
+//                    fEdges[fp] = edgeI;
+//                    break;
+//                }
+//            }
+//        }
+//
+//
+//        // Add face to edgeFaces
+//
+//        edgeFace0[edgeI] = eFaces[0];
+//
+//        if (eFaces.size() >= 2)
+//        {
+//            edgeFace1[edgeI] = eFaces[1];
+//
+//            if (eFaces.size() > 2)
+//            {
+//                edgeFacesRest.insert
+//                (
+//                    edgeI,
+//                    SubList<label>(eFaces, eFaces.size()-2, 2)
+//                );
+//            }
+//        }
+//
+//        edgeI++;
+//    }
+//}
+//
+//
+//void Foam::isoSurface::walkOrientation
+//(
+//    const triSurface& surf,
+//    const List<FixedList<label, 3> >& faceEdges,
+//    const labelList& edgeFace0,
+//    const labelList& edgeFace1,
+//    const label seedTriI,
+//    labelList& flipState
+//)
+//{
+//    // Do walk for consistent orientation.
+//    DynamicList<label> changedFaces(surf.size());
+//
+//    changedFaces.append(seedTriI);
+//
+//    while (changedFaces.size())
+//    {
+//        DynamicList<label> newChangedFaces(changedFaces.size());
+//
+//        forAll(changedFaces, i)
+//        {
+//            label triI = changedFaces[i];
+//            const labelledTri& tri = surf[triI];
+//            const FixedList<label, 3>& fEdges = faceEdges[triI];
+//
+//            forAll(fEdges, fp)
+//            {
+//                label edgeI = fEdges[fp];
+//
+//                // my points:
+//                label p0 = tri[fp];
+//                label p1 = tri[tri.fcIndex(fp)];
+//
+//                label nbrI =
+//                (
+//                    edgeFace0[edgeI] != triI
+//                  ? edgeFace0[edgeI]
+//                  : edgeFace1[edgeI]
+//                );
+//
+//                if (nbrI != -1 && flipState[nbrI] == -1)
+//                {
+//                    const labelledTri& nbrTri = surf[nbrI];
+//
+//                    // nbr points
+//                    label nbrFp = findIndex(nbrTri, p0);
+//
+//                    if (nbrFp == -1)
+//                    {
+//                        FatalErrorIn("isoSurface::walkOrientation(..)")
+//                            << "triI:" << triI
+//                            << " tri:" << tri
+//                            << " p0:" << p0
+//                            << " p1:" << p1
+//                            << " fEdges:" << fEdges
+//                            << " edgeI:" << edgeI
+//                            << " edgeFace0:" << edgeFace0[edgeI]
+//                            << " edgeFace1:" << edgeFace1[edgeI]
+//                            << " nbrI:" << nbrI
+//                            << " nbrTri:" << nbrTri
+//                            << abort(FatalError);
+//                    }
+//
+//
+//                    label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
+//
+//                    bool sameOrientation = (p1 == nbrP1);
+//
+//                    if (flipState[triI] == 0)
+//                    {
+//                        flipState[nbrI] = (sameOrientation ? 0 : 1);
+//                    }
+//                    else
+//                    {
+//                        flipState[nbrI] = (sameOrientation ? 1 : 0);
+//                    }
+//                    newChangedFaces.append(nbrI);
+//                }
+//            }
+//        }
+//
+//        changedFaces.transfer(newChangedFaces);
+//    }
+//}
+//
+//
+//void Foam::isoSurface::orientSurface
+//(
+//    triSurface& surf,
+//    const List<FixedList<label, 3> >& faceEdges,
+//    const labelList& edgeFace0,
+//    const labelList& edgeFace1,
+//    const Map<labelList>& edgeFacesRest
+//)
+//{
+//    // -1 : unvisited
+//    //  0 : leave as is
+//    //  1 : flip
+//    labelList flipState(surf.size(), -1);
+//
+//    label seedTriI = 0;
+//
+//    while (true)
+//    {
+//        // Find first unvisited triangle
+//        for
+//        (
+//            ;
+//            seedTriI < surf.size() && flipState[seedTriI] != -1;
+//            seedTriI++
+//        )
+//        {}
+//
+//        if (seedTriI == surf.size())
+//        {
+//            break;
+//        }
+//
+//        // Note: Determine orientation of seedTriI?
+//        // for now assume it is ok
+//        flipState[seedTriI] = 0;
+//
+//        walkOrientation
+//        (
+//            surf,
+//            faceEdges,
+//            edgeFace0,
+//            edgeFace1,
+//            seedTriI,
+//            flipState
+//        );
+//    }
+//
+//    // Do actual flipping
+//    surf.clearOut();
+//    forAll(surf, triI)
+//    {
+//        if (flipState[triI] == 1)
+//        {
+//            labelledTri tri(surf[triI]);
+//
+//            surf[triI][0] = tri[0];
+//            surf[triI][1] = tri[2];
+//            surf[triI][2] = tri[1];
+//        }
+//        else if (flipState[triI] == -1)
+//        {
+//            FatalErrorIn
+//            (
+//                "isoSurface::orientSurface(triSurface&, const label)"
+//            )   << "problem" << abort(FatalError);
+//        }
+//    }
+//}
+//
+//
+//// Checks if triangle is connected through edgeI only.
+//bool Foam::isoSurface::danglingTriangle
+//(
+//    const FixedList<label, 3>& fEdges,
+//    const labelList& edgeFace1
+//)
+//{
+//    label nOpen = 0;
+//    forAll(fEdges, i)
+//    {
+//        if (edgeFace1[fEdges[i]] == -1)
+//        {
+//            nOpen++;
+//        }
+//    }
+//
+//    if (nOpen == 1 || nOpen == 2 || nOpen == 3)
+//    {
+//        return true;
+//    }
+//    else
+//    {
+//        return false;
+//    }
+//}
+//
+//
+//// Mark triangles to keep. Returns number of dangling triangles.
+//Foam::label Foam::isoSurface::markDanglingTriangles
+//(
+//    const List<FixedList<label, 3> >& faceEdges,
+//    const labelList& edgeFace0,
+//    const labelList& edgeFace1,
+//    const Map<labelList>& edgeFacesRest,
+//    boolList& keepTriangles
+//)
+//{
+//    keepTriangles.setSize(faceEdges.size());
+//    keepTriangles = true;
+//
+//    label nDangling = 0;
+//
+//    // Remove any dangling triangles
+//    forAllConstIter(Map<labelList>,  edgeFacesRest, iter)
+//    {
+//        // These are all the non-manifold edges. Filter out all triangles
+//        // with only one connected edge (= this edge)
+//
+//        label edgeI = iter.key();
+//        const labelList& otherEdgeFaces = iter();
+//
+//        // Remove all dangling triangles
+//        if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
+//        {
+//            keepTriangles[edgeFace0[edgeI]] = false;
+//            nDangling++;
+//        }
+//        if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
+//        {
+//            keepTriangles[edgeFace1[edgeI]] = false;
+//            nDangling++;
+//        }
+//        forAll(otherEdgeFaces, i)
+//        {
+//            label triI = otherEdgeFaces[i];
+//            if (danglingTriangle(faceEdges[triI], edgeFace1))
+//            {
+//                keepTriangles[triI] = false;
+//                nDangling++;
+//            }
+//        }
+//    }
+//    return nDangling;
+//}
 
 
 Foam::triSurface Foam::isoSurface::subsetMesh
@@ -1715,6 +1953,7 @@ Foam::isoSurface::isoSurface
     const scalarField& pVals,
     const scalar iso,
     const bool regularise,
+    const boundBox& bounds,
     const scalar mergeTol
 )
 :
@@ -1722,6 +1961,7 @@ Foam::isoSurface::isoSurface
     pVals_(pVals),
     iso_(iso),
     regularise_(regularise),
+    bounds_(bounds),
     mergeDistance_(mergeTol*mesh_.bounds().mag())
 {
     if (debug)
@@ -1957,57 +2197,102 @@ Foam::isoSurface::isoSurface
     }
 
 
-
-    DynamicList<point> triPoints(nCutCells_);
-    DynamicList<label> triMeshCells(nCutCells_);
-
-    generateTriPoints
-    (
-        cValsPtr_(),
-        pVals_,
-
-        meshC,
-        mesh_.points(),
-
-        snappedPoints,
-        snappedCc,
-        snappedPoint,
-
-        triPoints,
-        triMeshCells
-    );
-
-    if (debug)
     {
-        Pout<< "isoSurface : generated " << triMeshCells.size()
-            << " unmerged triangles from " << triPoints.size()
-            << " unmerged points." << endl;
-    }
+        DynamicList<point> triPoints(3*nCutCells_);
+        DynamicList<label> triMeshCells(nCutCells_);
 
-
-    // Merge points and compact out non-valid triangles
-    labelList triMap;           // merged to unmerged triangle
-    triSurface::operator=
-    (
-        stitchTriPoints
+        generateTriPoints
         (
-            true,               // check for duplicate tris
-            triPoints,
-            triPointMergeMap_,  // unmerged to merged point
-            triMap
-        )
-    );
+            cValsPtr_(),
+            pVals_,
 
-    if (debug)
-    {
-        Pout<< "isoSurface : generated " << triMap.size()
-            << " merged triangles." << endl;
-    }
+            meshC,
+            mesh_.points(),
 
-    meshCells_.setSize(triMap.size());
-    forAll(triMap, i)
-    {
-        meshCells_[i] = triMeshCells[triMap[i]];
+            snappedPoints,
+            snappedCc,
+            snappedPoint,
+
+            triPoints,      // 3 points of the triangle
+            triMeshCells    // per triangle the originating cell
+        );
+
+        if (debug)
+        {
+            Pout<< "isoSurface : generated " << triMeshCells.size()
+                << " unmerged triangles from " << triPoints.size()
+                << " unmerged points." << endl;
+        }
+
+        label nOldPoints = triPoints.size();
+
+        // Trimmed to original triangle
+        DynamicList<label> trimTriMap;
+        // Trimmed to original point
+        labelList trimTriPointMap;
+        if (bounds_ != boundBox::greatBox)
+        {
+            trimToBox
+            (
+                treeBoundBox(bounds_),
+                triPoints,              // new points
+                trimTriMap,             // map from (new) triangle to original
+                trimTriPointMap,        // map from (new) point to original
+                interpolatedPoints_,    // labels of newly introduced points
+                interpolatedOldPoints_, // and their interpolation
+                interpolationWeights_
+            );
+            triMeshCells = labelField(triMeshCells, trimTriMap);
+        }
+
+
+        // Merge points and compact out non-valid triangles
+        labelList triMap;           // merged to unmerged triangle
+        triSurface::operator=
+        (
+            stitchTriPoints
+            (
+                true,               // check for duplicate tris
+                triPoints,
+                triPointMergeMap_,  // unmerged to merged point
+                triMap
+            )
+        );
+
+        if (debug)
+        {
+            Pout<< "isoSurface : generated " << triMap.size()
+                << " merged triangles." << endl;
+        }
+
+
+        if (bounds_ != boundBox::greatBox)
+        {
+            // Adjust interpolatedPoints_
+            inplaceRenumber(triPointMergeMap_, interpolatedPoints_);
+
+            // Adjust triPointMergeMap_
+            labelList newTriPointMergeMap(nOldPoints, -1);
+            forAll(trimTriPointMap, trimPointI)
+            {
+                label oldPointI = trimTriPointMap[trimPointI];
+                if (oldPointI >= 0)
+                {
+                    label pointI = triPointMergeMap_[trimPointI];
+                    if (pointI >= 0)
+                    {
+                        newTriPointMergeMap[oldPointI] = pointI;
+                    }
+                }
+            }
+            triPointMergeMap_.transfer(newTriPointMergeMap);
+        }
+
+        meshCells_.setSize(triMap.size());
+        forAll(triMap, i)
+        {
+            meshCells_[i] = triMeshCells[triMap[i]];
+        }
     }
 
     if (debug)
@@ -2023,68 +2308,68 @@ Foam::isoSurface::isoSurface
     }
 
 
-    if (false)
-    {
-        List<FixedList<label, 3> > faceEdges;
-        labelList edgeFace0, edgeFace1;
-        Map<labelList> edgeFacesRest;
-
-
-        while (true)
-        {
-            // Calculate addressing
-            calcAddressing
-            (
-                *this,
-                faceEdges,
-                edgeFace0,
-                edgeFace1,
-                edgeFacesRest
-            );
-
-            // See if any dangling triangles
-            boolList keepTriangles;
-            label nDangling = markDanglingTriangles
-            (
-                faceEdges,
-                edgeFace0,
-                edgeFace1,
-                edgeFacesRest,
-                keepTriangles
-            );
-
-            if (debug)
-            {
-                Pout<< "isoSurface : detected " << nDangling
-                    << " dangling triangles." << endl;
-            }
-
-            if (nDangling == 0)
-            {
-                break;
-            }
-
-            // Create face map (new to old)
-            labelList subsetTriMap(findIndices(keepTriangles, true));
-
-            labelList subsetPointMap;
-            labelList reversePointMap;
-            triSurface::operator=
-            (
-                subsetMesh
-                (
-                    *this,
-                    subsetTriMap,
-                    reversePointMap,
-                    subsetPointMap
-                )
-            );
-            meshCells_ = labelField(meshCells_, subsetTriMap);
-            inplaceRenumber(reversePointMap, triPointMergeMap_);
-        }
-
-        orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
-    }
+    //if (false)
+    //{
+    //    List<FixedList<label, 3> > faceEdges;
+    //    labelList edgeFace0, edgeFace1;
+    //    Map<labelList> edgeFacesRest;
+    //
+    //
+    //    while (true)
+    //    {
+    //        // Calculate addressing
+    //        calcAddressing
+    //        (
+    //            *this,
+    //            faceEdges,
+    //            edgeFace0,
+    //            edgeFace1,
+    //            edgeFacesRest
+    //        );
+    //
+    //        // See if any dangling triangles
+    //        boolList keepTriangles;
+    //        label nDangling = markDanglingTriangles
+    //        (
+    //            faceEdges,
+    //            edgeFace0,
+    //            edgeFace1,
+    //            edgeFacesRest,
+    //            keepTriangles
+    //        );
+    //
+    //        if (debug)
+    //        {
+    //            Pout<< "isoSurface : detected " << nDangling
+    //                << " dangling triangles." << endl;
+    //        }
+    //
+    //        if (nDangling == 0)
+    //        {
+    //            break;
+    //        }
+    //
+    //        // Create face map (new to old)
+    //        labelList subsetTriMap(findIndices(keepTriangles, true));
+    //
+    //        labelList subsetPointMap;
+    //        labelList reversePointMap;
+    //        triSurface::operator=
+    //        (
+    //            subsetMesh
+    //            (
+    //                *this,
+    //                subsetTriMap,
+    //                reversePointMap,
+    //                subsetPointMap
+    //            )
+    //        );
+    //        meshCells_ = labelField(meshCells_, subsetTriMap);
+    //        inplaceRenumber(reversePointMap, triPointMergeMap_);
+    //    }
+    //
+    //    orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
+    //}
 
 
     if (debug)
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H
index 9f6dc7d76e..c74586a2ce 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.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-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -76,6 +76,8 @@ namespace Foam
 {
 
 class fvMesh;
+class plane;
+class treeBoundBox;
 
 /*---------------------------------------------------------------------------*\
                        Class isoSurface Declaration
@@ -116,10 +118,12 @@ class isoSurface
         //- Regularise?
         const Switch regularise_;
 
+        //- Optional bounds
+        const boundBox bounds_;
+
         //- When to merge points
         const scalar mergeDistance_;
 
-
         //- Whether face might be cut
         List<cellCutType> faceCutType_;
 
@@ -135,6 +139,15 @@ class isoSurface
         //- For every unmerged triangle point the point in the triSurface
         labelList triPointMergeMap_;
 
+        //- triSurface points that have weighted interpolation
+        DynamicList<label> interpolatedPoints_;
+
+        //- corresponding original, unmerged points
+        DynamicList<FixedList<label, 3> > interpolatedOldPoints_;
+
+        //- corresponding weights
+        DynamicList<FixedList<scalar, 3> > interpolationWeights_;
+
 
     // Private Member Functions
 
@@ -193,20 +206,8 @@ class isoSurface
             const scalarField& pVals
         );
 
-        static labelPair findCommonPoints
-        (
-            const labelledTri&,
-            const labelledTri&
-        );
-
         static point calcCentre(const triSurface&);
 
-        static pointIndexHit collapseSurface
-        (
-            pointField& localPoints,
-            DynamicList<labelledTri, 64>& localTris
-        );
-
         //- Determine per cc whether all near cuts can be snapped to single
         //  point.
         void calcSnappedCc
@@ -323,6 +324,17 @@ class isoSurface
             DynamicList<label>& triMeshCells
         ) const;
 
+        template<class Type>
+        static tmp<Field<Type> > interpolate
+        (
+            const label nPoints,
+            const labelList& triPointMergeMap,
+            const labelList& interpolatedPoints,
+            const List<FixedList<label, 3> >& interpolatedOldPoints,
+            const List<FixedList<scalar, 3> >& interpolationWeights,
+            const DynamicList<Type>& unmergedValues
+        );
+
         triSurface stitchTriPoints
         (
             const bool checkDuplicates,
@@ -331,57 +343,38 @@ class isoSurface
             labelList& triMap               // merged to unmerged triangle
         ) const;
 
+        //- Trim triangle to planes
+        static void trimToPlanes
+        (
+            const PtrList<plane>& planes,
+            const triPointRef& tri,
+            DynamicList<point>& newTriPoints
+        );
+
+        //- Trim all triangles to box
+        static void trimToBox
+        (
+            const treeBoundBox& bb,
+            DynamicList<point>& triPoints,
+            DynamicList<label>& triMeshCells
+        );
+
+        //- Trim all triangles to box. Determine interpolation
+        //  for existing and new points
+        static void trimToBox
+        (
+            const treeBoundBox& bb,
+            DynamicList<point>& triPoints,
+            DynamicList<label>& triMap,
+            labelList& triPointMap,
+            labelList& interpolatedPoints,
+            List<FixedList<label, 3> >& interpolatedOldPoints,
+            List<FixedList<scalar, 3> >& interpolationWeights
+        );
+
         //- Check single triangle for (topological) validity
         static bool validTri(const triSurface&, const label);
 
-        //- Determine edge-face addressing
-        void calcAddressing
-        (
-            const triSurface& surf,
-            List<FixedList<label, 3> >& faceEdges,
-            labelList& edgeFace0,
-            labelList& edgeFace1,
-            Map<labelList>& edgeFacesRest
-        ) const;
-
-        //- Determine orientation
-        static void walkOrientation
-        (
-            const triSurface& surf,
-            const List<FixedList<label, 3> >& faceEdges,
-            const labelList& edgeFace0,
-            const labelList& edgeFace1,
-            const label seedTriI,
-            labelList& flipState
-        );
-
-        //- Orient surface
-        static void orientSurface
-        (
-            triSurface&,
-            const List<FixedList<label, 3> >& faceEdges,
-            const labelList& edgeFace0,
-            const labelList& edgeFace1,
-            const Map<labelList>& edgeFacesRest
-        );
-
-        //- Is triangle (given by 3 edges) not fully connected?
-        static bool danglingTriangle
-        (
-            const FixedList<label, 3>& fEdges,
-            const labelList& edgeFace1
-        );
-
-        //- Mark all non-fully connected triangles
-        static label markDanglingTriangles
-        (
-            const List<FixedList<label, 3> >& faceEdges,
-            const labelList& edgeFace0,
-            const labelList& edgeFace1,
-            const Map<labelList>& edgeFacesRest,
-            boolList& keepTriangles
-        );
-
         static triSurface subsetMesh
         (
             const triSurface& s,
@@ -392,6 +385,10 @@ class isoSurface
 
 public:
 
+    //- Declare friendship with isoSurfaceCell to share some functionality
+    friend class isoSurfaceCell;
+
+
     //- Runtime type information
     TypeName("isoSurface");
 
@@ -407,24 +404,19 @@ public:
             const scalarField& pointIsoVals,
             const scalar iso,
             const bool regularise,
+            const boundBox& bounds = boundBox::greatBox,
             const scalar mergeTol = 1e-6    // fraction of bounding box
         );
 
 
     // Member Functions
 
-        //- For every face original cell in mesh
+        //- For every triangle the original cell in mesh
         const labelList& meshCells() const
         {
             return meshCells_;
         }
 
-        //- For every unmerged triangle point the point in the triSurface
-        const labelList& triPointMergeMap() const
-        {
-            return triPointMergeMap_;
-        }
-
         //- Interpolates cCoords,pCoords. Uses the references to the original
         //  fields used to create the iso surface.
         template<class Type>
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
index 8123541ebe..b32edd1595 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,6 +30,9 @@ License
 #include "tetMatcher.H"
 #include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
+#include "Time.H"
+#include "triPoints.H"
+#include "isoSurface.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1222,144 +1225,6 @@ void Foam::isoSurfaceCell::calcAddressing
 }
 
 
-//void Foam::isoSurfaceCell::walkOrientation
-//(
-//    const triSurface& surf,
-//    const List<FixedList<label, 3> >& faceEdges,
-//    const labelList& edgeFace0,
-//    const labelList& edgeFace1,
-//    const label seedTriI,
-//    labelList& flipState
-//)
-//{
-//    // Do walk for consistent orientation.
-//    DynamicList<label> changedFaces(surf.size());
-//
-//    changedFaces.append(seedTriI);
-//
-//    while (changedFaces.size())
-//    {
-//        DynamicList<label> newChangedFaces(changedFaces.size());
-//
-//        forAll(changedFaces, i)
-//        {
-//            label triI = changedFaces[i];
-//            const labelledTri& tri = surf[triI];
-//            const FixedList<label, 3>& fEdges = faceEdges[triI];
-//
-//            forAll(fEdges, fp)
-//            {
-//                label edgeI = fEdges[fp];
-//
-//                // my points:
-//                label p0 = tri[fp];
-//                label p1 = tri[tri.fcIndex(fp)];
-//
-//                label nbrI =
-//                (
-//                    edgeFace0[edgeI] != triI
-//                  ? edgeFace0[edgeI]
-//                  : edgeFace1[edgeI]
-//                );
-//
-//                if (nbrI != -1 && flipState[nbrI] == -1)
-//                {
-//                    const labelledTri& nbrTri = surf[nbrI];
-//
-//                    // nbr points
-//                    label nbrFp = findIndex(nbrTri, p0);
-//                    label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
-//
-//                    bool sameOrientation = (p1 == nbrP1);
-//
-//                    if (flipState[triI] == 0)
-//                    {
-//                        flipState[nbrI] = (sameOrientation ? 0 : 1);
-//                    }
-//                    else
-//                    {
-//                        flipState[nbrI] = (sameOrientation ? 1 : 0);
-//                    }
-//                    newChangedFaces.append(nbrI);
-//                }
-//            }
-//        }
-//
-//        changedFaces.transfer(newChangedFaces);
-//    }
-//}
-//
-//
-//void Foam::isoSurfaceCell::orientSurface
-//(
-//    triSurface& surf,
-//    const List<FixedList<label, 3> >& faceEdges,
-//    const labelList& edgeFace0,
-//    const labelList& edgeFace1,
-//    const Map<labelList>& edgeFacesRest
-//)
-//{
-//    // -1 : unvisited
-//    //  0 : leave as is
-//    //  1 : flip
-//    labelList flipState(surf.size(), -1);
-//
-//    label seedTriI = 0;
-//
-//    while (true)
-//    {
-//        // Find first unvisited triangle
-//        for
-//        (
-//            ;
-//            seedTriI < surf.size() && flipState[seedTriI] != -1;
-//            seedTriI++
-//        )
-//        {}
-//
-//        if (seedTriI == surf.size())
-//        {
-//            break;
-//        }
-//
-//        // Note: Determine orientation of seedTriI?
-//        // for now assume it is ok
-//        flipState[seedTriI] = 0;
-//
-//        walkOrientation
-//        (
-//            surf,
-//            faceEdges,
-//            edgeFace0,
-//            edgeFace1,
-//            seedTriI,
-//            flipState
-//        );
-//    }
-//
-//    // Do actual flipping
-//    surf.clearOut();
-//    forAll(surf, triI)
-//    {
-//        if (flipState[triI] == 1)
-//        {
-//            labelledTri tri(surf[triI]);
-//
-//            surf[triI][0] = tri[0];
-//            surf[triI][1] = tri[2];
-//            surf[triI][2] = tri[1];
-//        }
-//        else if (flipState[triI] == -1)
-//        {
-//            FatalErrorIn
-//            (
-//                "isoSurfaceCell::orientSurface(triSurface&, const label)"
-//            )   << "problem" << abort(FatalError);
-//        }
-//    }
-//}
-
-
 // Checks if triangle is connected through edgeI only.
 bool Foam::isoSurfaceCell::danglingTriangle
 (
@@ -1517,6 +1382,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
     const scalarField& pVals,
     const scalar iso,
     const bool regularise,
+    const boundBox& bounds,
     const scalar mergeTol
 )
 :
@@ -1524,6 +1390,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
     cVals_(cVals),
     pVals_(pVals),
     iso_(iso),
+    bounds_(bounds),
     mergeDistance_(mergeTol*mesh.bounds().mag())
 {
     if (debug)
@@ -1607,56 +1474,104 @@ Foam::isoSurfaceCell::isoSurfaceCell
     }
 
 
-
-    DynamicList<point> triPoints(nCutCells_);
-    DynamicList<label> triMeshCells(nCutCells_);
-
-    generateTriPoints
-    (
-        cVals,
-        pVals,
-
-        mesh_.cellCentres(),
-        mesh_.points(),
-
-        snappedPoints,
-        snappedCc,
-        snappedPoint,
-
-        triPoints,
-        triMeshCells
-    );
-
-    if (debug)
     {
-        Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
-            << " unmerged triangles." << endl;
-    }
+        DynamicList<point> triPoints(nCutCells_);
+        DynamicList<label> triMeshCells(nCutCells_);
 
-    // Merge points and compact out non-valid triangles
-    labelList triMap;           // merged to unmerged triangle
-    triSurface::operator=
-    (
-        stitchTriPoints
+        generateTriPoints
         (
-            regularise,         // check for duplicate tris
+            cVals,
+            pVals,
+
+            mesh_.cellCentres(),
+            mesh_.points(),
+
+            snappedPoints,
+            snappedCc,
+            snappedPoint,
+
             triPoints,
-            triPointMergeMap_,  // unmerged to merged point
-            triMap
-        )
-    );
+            triMeshCells
+        );
 
-    if (debug)
-    {
-        Pout<< "isoSurfaceCell : generated " << triMap.size()
-            << " merged triangles." << endl;
+        if (debug)
+        {
+            Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
+                << " unmerged triangles." << endl;
+        }
+
+
+        label nOldPoints = triPoints.size();
+
+        // Trimmed to original triangle
+        DynamicList<label> trimTriMap;
+        // Trimmed to original point
+        labelList trimTriPointMap;
+        if (bounds_ != boundBox::greatBox)
+        {
+            isoSurface::trimToBox
+            (
+                treeBoundBox(bounds_),
+                triPoints,              // new points
+                trimTriMap,             // map from (new) triangle to original
+                trimTriPointMap,        // map from (new) point to original
+                interpolatedPoints_,    // labels of newly introduced points
+                interpolatedOldPoints_, // and their interpolation
+                interpolationWeights_
+            );
+            triMeshCells = labelField(triMeshCells, trimTriMap);
+        }
+
+
+
+        // Merge points and compact out non-valid triangles
+        labelList triMap;
+        triSurface::operator=
+        (
+            stitchTriPoints
+            (
+                regularise,         // check for duplicate tris
+                triPoints,
+                triPointMergeMap_,  // unmerged to merged point
+                triMap              // merged to unmerged triangle
+            )
+        );
+
+        if (debug)
+        {
+            Pout<< "isoSurfaceCell : generated " << triMap.size()
+                << " merged triangles." << endl;
+        }
+
+        if (bounds_ != boundBox::greatBox)
+        {
+            // Adjust interpolatedPoints_
+            inplaceRenumber(triPointMergeMap_, interpolatedPoints_);
+
+            // Adjust triPointMergeMap_
+            labelList newTriPointMergeMap(nOldPoints, -1);
+            forAll(trimTriPointMap, trimPointI)
+            {
+                label oldPointI = trimTriPointMap[trimPointI];
+                if (oldPointI >= 0)
+                {
+                    label pointI = triPointMergeMap_[trimPointI];
+                    if (pointI >= 0)
+                    {
+                        newTriPointMergeMap[oldPointI] = pointI;
+                    }
+                }
+            }
+            triPointMergeMap_.transfer(newTriPointMergeMap);
+        }
+
+        meshCells_.setSize(triMap.size());
+        forAll(triMap, i)
+        {
+            meshCells_[i] = triMeshCells[triMap[i]];
+        }
     }
 
-    meshCells_.setSize(triMap.size());
-    forAll(triMap, i)
-    {
-        meshCells_[i] = triMeshCells[triMap[i]];
-    }
 
     if (debug)
     {
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
index 7e12efc7bb..910f3fe46e 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.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-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,6 +48,7 @@ SourceFiles
 #include "labelPair.H"
 #include "pointIndexHit.H"
 #include "PackedBoolList.H"
+#include "boundBox.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -55,6 +56,7 @@ namespace Foam
 {
 
 class polyMesh;
+class plane;
 
 /*---------------------------------------------------------------------------*\
                        Class isoSurfaceCell Declaration
@@ -91,6 +93,9 @@ class isoSurfaceCell
         //- isoSurfaceCell value
         const scalar iso_;
 
+        //- Optional bounds
+        const boundBox bounds_;
+
         //- When to merge points
         const scalar mergeDistance_;
 
@@ -106,6 +111,15 @@ class isoSurfaceCell
         //- For every unmerged triangle point the point in the triSurface
         labelList triPointMergeMap_;
 
+        //- triSurface points that have weighted interpolation
+        DynamicList<label> interpolatedPoints_;
+
+        //- corresponding original, unmerged points
+        DynamicList<FixedList<label, 3> > interpolatedOldPoints_;
+
+        //- corresponding weights
+        DynamicList<FixedList<scalar, 3> > interpolationWeights_;
+
 
     // Private Member Functions
 
@@ -214,19 +228,15 @@ class isoSurfaceCell
         void generateTriPoints
         (
             const DynamicList<Type>& snapped,
-            const scalar isoVal0,
             const scalar s0,
             const Type& p0,
             const label p0Index,
-            const scalar isoVal1,
             const scalar s1,
             const Type& p1,
             const label p1Index,
-            const scalar isoVal2,
             const scalar s2,
             const Type& p2,
             const label p2Index,
-            const scalar isoVal3,
             const scalar s3,
             const Type& p3,
             const label p3Index,
@@ -271,27 +281,6 @@ class isoSurfaceCell
             Map<labelList>& edgeFacesRest
         ) const;
 
-        ////- Determine orientation
-        //static void walkOrientation
-        //(
-        //    const triSurface& surf,
-        //    const List<FixedList<label, 3> >& faceEdges,
-        //    const labelList& edgeFace0,
-        //    const labelList& edgeFace1,
-        //    const label seedTriI,
-        //    labelList& flipState
-        //);
-
-        ////- Orient surface
-        //static void orientSurface
-        //(
-        //    triSurface&,
-        //    const List<FixedList<label, 3> >& faceEdges,
-        //    const labelList& edgeFace0,
-        //    const labelList& edgeFace1,
-        //    const Map<labelList>& edgeFacesRest
-        //);
-
         //- Is triangle (given by 3 edges) not fully connected?
         static bool danglingTriangle
         (
@@ -317,8 +306,6 @@ class isoSurfaceCell
             labelList& newToOldPoints
         );
 
-        //- Combine all triangles inside a cell into a minimal triangulation
-        void combineCellTriangles();
 
 public:
 
@@ -336,6 +323,7 @@ public:
             const scalarField& pointValues,
             const scalar iso,
             const bool regularise,
+            const boundBox& bounds = boundBox::greatBox,
             const scalar mergeTol = 1e-6    // fraction of bounding box
         );
 
@@ -348,24 +336,6 @@ public:
             return meshCells_;
         }
 
-        //- For every unmerged triangle point the point in the triSurface
-        const labelList triPointMergeMap() const
-        {
-            return triPointMergeMap_;
-        }
-
-
-        //- Interpolates cCoords,pCoords. Takes the original fields
-        //  used to create the iso surface.
-        template<class Type>
-        tmp<Field<Type> > interpolate
-        (
-            const scalarField& cVals,
-            const scalarField& pVals,
-            const Field<Type>& cCoords,
-            const Field<Type>& pCoords
-        ) const;
-
         //- Interpolates cCoords,pCoords.
         template<class Type>
         tmp<Field<Type> > interpolate
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
index c78aafeee5..8fdee4831b 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
@@ -26,6 +26,7 @@ License
 #include "isoSurfaceCell.H"
 #include "polyMesh.H"
 #include "tetMatcher.H"
+#include "isoSurface.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -76,22 +77,18 @@ void Foam::isoSurfaceCell::generateTriPoints
 (
     const DynamicList<Type>& snapped,
 
-    const scalar isoVal0,
     const scalar s0,
     const Type& p0,
     const label p0Index,
 
-    const scalar isoVal1,
     const scalar s1,
     const Type& p1,
     const label p1Index,
 
-    const scalar isoVal2,
     const scalar s2,
     const Type& p2,
     const label p2Index,
 
-    const scalar isoVal3,
     const scalar s3,
     const Type& p3,
     const label p3Index,
@@ -124,160 +121,196 @@ void Foam::isoSurfaceCell::generateTriPoints
         case 0x0F:
         break;
 
-        case 0x0E:
         case 0x01:
+        case 0x0E:
         {
-            // 0 is common point. Orient such that normal points in positive
-            // gradient direction
-            if (isoVal0 >= isoVal1)
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+            );
+            if (triIndex == 0x0E)
             {
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-            }
-            else
-            {
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+                // Flip normals
+                label sz = pts.size();
+                Swap(pts[sz-2], pts[sz-1]);
             }
         }
         break;
 
-        case 0x0D:
         case 0x02:
+        case 0x0D:
         {
-            // 1 is common point
-            if (isoVal1 >= isoVal0)
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+            );
+
+            if (triIndex == 0x0D)
             {
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-            }
-            else
-            {
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+                // Flip normals
+                label sz = pts.size();
+                Swap(pts[sz-2], pts[sz-1]);
             }
         }
         break;
 
-        case 0x0C:
         case 0x03:
+        case 0x0C:
         {
-            Type s02 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
-            Type s13 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
+            Type p0p2 =
+                generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
+            Type p1p3 =
+                generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
 
-            if (isoVal0 >= isoVal3)
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+            );
+            pts.append(p1p3);
+            pts.append(p0p2);
+
+            pts.append(p1p3);
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+            );
+            pts.append(p0p2);
+
+            if (triIndex == 0x0C)
             {
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-                pts.append(s02);
-                pts.append(s13);
-                pts.append(s13);
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-                pts.append(s02);
-            }
-            else
-            {
-                pts.append(s02);
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-                pts.append(s13);
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-                pts.append(s13);
-                pts.append(s02);
+                // Flip normals
+                label sz = pts.size();
+                Swap(pts[sz-5], pts[sz-4]);
+                Swap(pts[sz-2], pts[sz-1]);
             }
         }
         break;
 
-        case 0x0B:
         case 0x04:
+        case 0x0B:
         {
-            // 2 is common point
-            if (isoVal2 >= isoVal0)
+            pts.append
+            (
+                generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index)
+            );
+
+            if (triIndex == 0x0B)
             {
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
-            }
-            else
-            {
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+                // Flip normals
+                label sz = pts.size();
+                Swap(pts[sz-2], pts[sz-1]);
             }
         }
         break;
 
-        case 0x0A:
         case 0x05:
+        case 0x0A:
         {
-            Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+            Type p0p1 =
+                generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type p2p3 =
+                generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
 
-            if (isoVal3 >= isoVal0)
+            pts.append(p0p1);
+            pts.append(p2p3);
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+            );
+
+            pts.append(p0p1);
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+            );
+            pts.append(p2p3);
+
+            if (triIndex == 0x0A)
             {
-                pts.append(s01);
-                pts.append(s23);
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-                pts.append(s01);
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-                pts.append(s23);
-            }
-            else
-            {
-                pts.append(s23);
-                pts.append(s01);
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-                pts.append(s01);
-                pts.append(s23);
+                // Flip normals
+                label sz = pts.size();
+                Swap(pts[sz-5], pts[sz-4]);
+                Swap(pts[sz-2], pts[sz-1]);
             }
         }
         break;
 
-        case 0x09:
         case 0x06:
+        case 0x09:
         {
-            Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+            Type p0p1 =
+                generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type p2p3 =
+                generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
 
-            if (isoVal3 >= isoVal1)
+            pts.append(p0p1);
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)
+            );
+            pts.append(p2p3);
+
+            pts.append(p0p1);
+            pts.append(p2p3);
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)
+            );
+
+            if (triIndex == 0x09)
             {
-                pts.append(s01);
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-                pts.append(s23);
-                pts.append(s01);
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-                pts.append(s23);
-            }
-            else
-            {
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-                pts.append(s01);
-                pts.append(s23);
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-                pts.append(s01);
-                pts.append(s23);
+                // Flip normals
+                label sz = pts.size();
+                Swap(pts[sz-5], pts[sz-4]);
+                Swap(pts[sz-2], pts[sz-1]);
             }
         }
         break;
 
-        case 0x07:
         case 0x08:
+        case 0x07:
         {
-            // 3 is common point
-            if (isoVal3 >= isoVal0)
+            pts.append
+            (
+                generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index)
+            );
+            if (triIndex == 0x07)
             {
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
-            }
-            else
-            {
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+                // Flip normals
+                label sz = pts.size();
+                Swap(pts[sz-2], pts[sz-1]);
             }
         }
         break;
@@ -341,22 +374,18 @@ void Foam::isoSurfaceCell::generateTriPoints
                     (
                         snappedPoints,
 
-                        pVals_[f0[1]],
                         pVals[f0[1]],
                         pCoords[f0[1]],
                         snappedPoint[f0[1]],
 
-                        pVals_[f0[0]],
                         pVals[f0[0]],
                         pCoords[f0[0]],
                         snappedPoint[f0[0]],
 
-                        pVals_[f0[2]],
                         pVals[f0[2]],
                         pCoords[f0[2]],
                         snappedPoint[f0[2]],
 
-                        pVals_[oppositeI],
                         pVals[oppositeI],
                         pCoords[oppositeI],
                         snappedPoint[oppositeI],
@@ -370,22 +399,18 @@ void Foam::isoSurfaceCell::generateTriPoints
                     (
                         snappedPoints,
 
-                        pVals_[f0[0]],
                         pVals[f0[0]],
                         pCoords[f0[0]],
                         snappedPoint[f0[0]],
 
-                        pVals_[f0[1]],
                         pVals[f0[1]],
                         pCoords[f0[1]],
                         snappedPoint[f0[1]],
 
-                        pVals_[f0[2]],
                         pVals[f0[2]],
                         pCoords[f0[2]],
                         snappedPoint[f0[2]],
 
-                        pVals_[oppositeI],
                         pVals[oppositeI],
                         pCoords[oppositeI],
                         snappedPoint[oppositeI],
@@ -411,7 +436,6 @@ void Foam::isoSurfaceCell::generateTriPoints
                     }
 
                     label fp = f.fcIndex(fp0);
-
                     for (label i = 2; i < f.size(); i++)
                     {
                         label nextFp = f.fcIndex(fp);
@@ -425,22 +449,18 @@ void Foam::isoSurfaceCell::generateTriPoints
                             (
                                 snappedPoints,
 
-                                pVals_[tri[1]],
                                 pVals[tri[1]],
                                 pCoords[tri[1]],
                                 snappedPoint[tri[1]],
 
-                                pVals_[tri[0]],
                                 pVals[tri[0]],
                                 pCoords[tri[0]],
                                 snappedPoint[tri[0]],
 
-                                pVals_[tri[2]],
                                 pVals[tri[2]],
                                 pCoords[tri[2]],
                                 snappedPoint[tri[2]],
 
-                                cVals_[cellI],
                                 cVals[cellI],
                                 cCoords[cellI],
                                 snappedCc[cellI],
@@ -454,22 +474,18 @@ void Foam::isoSurfaceCell::generateTriPoints
                             (
                                 snappedPoints,
 
-                                pVals_[tri[0]],
                                 pVals[tri[0]],
                                 pCoords[tri[0]],
                                 snappedPoint[tri[0]],
 
-                                pVals_[tri[1]],
                                 pVals[tri[1]],
                                 pCoords[tri[1]],
                                 snappedPoint[tri[1]],
 
-                                pVals_[tri[2]],
                                 pVals[tri[2]],
                                 pCoords[tri[2]],
                                 snappedPoint[tri[2]],
 
-                                cVals_[cellI],
                                 cVals[cellI],
                                 cCoords[cellI],
                                 snappedCc[cellI],
@@ -495,7 +511,7 @@ void Foam::isoSurfaceCell::generateTriPoints
 
     if (countNotFoundTets > 0)
     {
-        WarningIn("Foam::isoSurfaceCell::generateTriPoints")
+        WarningIn("Foam::isoSurfaceCell::generateTriPoints(..)")
             << "Could not find " << countNotFoundTets
             << " tet base points, which may lead to inverted triangles."
             << endl;
@@ -510,13 +526,11 @@ template<class Type>
 Foam::tmp<Foam::Field<Type> >
 Foam::isoSurfaceCell::interpolate
 (
-    const scalarField& cVals,
-    const scalarField& pVals,
     const Field<Type>& cCoords,
     const Field<Type>& pCoords
 ) const
 {
-    DynamicList<Type> triPoints(nCutCells_);
+    DynamicList<Type> triPoints(3*nCutCells_);
     DynamicList<label> triMeshCells(nCutCells_);
 
     // Dummy snap data
@@ -524,59 +538,6 @@ Foam::isoSurfaceCell::interpolate
     labelList snappedCc(mesh_.nCells(), -1);
     labelList snappedPoint(mesh_.nPoints(), -1);
 
-
-    generateTriPoints
-    (
-        cVals,
-        pVals,
-
-        cCoords,
-        pCoords,
-
-        snappedPoints,
-        snappedCc,
-        snappedPoint,
-
-        triPoints,
-        triMeshCells
-    );
-
-
-    // One value per point
-    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
-    Field<Type>& values = tvalues();
-
-    forAll(triPoints, i)
-    {
-        label mergedPointI = triPointMergeMap_[i];
-
-        if (mergedPointI >= 0)
-        {
-            values[mergedPointI] = triPoints[i];
-        }
-    }
-
-    return tvalues;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type> >
-Foam::isoSurfaceCell::interpolate
-(
-    const Field<Type>& cCoords,
-    const Field<Type>& pCoords
-) const
-{
-    DynamicList<Type> triPoints(nCutCells_);
-    DynamicList<label> triMeshCells(nCutCells_);
-
-    // Dummy snap data
-    DynamicList<Type> snappedPoints;
-    labelList snappedCc(mesh_.nCells(), -1);
-    labelList snappedPoint(mesh_.nPoints(), -1);
-
-
     generateTriPoints
     (
         cVals_,
@@ -593,22 +554,15 @@ Foam::isoSurfaceCell::interpolate
         triMeshCells
     );
 
-
-    // One value per point
-    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
-    Field<Type>& values = tvalues();
-
-    forAll(triPoints, i)
-    {
-        label mergedPointI = triPointMergeMap_[i];
-
-        if (mergedPointI >= 0)
-        {
-            values[mergedPointI] = triPoints[i];
-        }
-    }
-
-    return tvalues;
+    return isoSurface::interpolate
+    (
+        points().size(),
+        triPointMergeMap_,
+        interpolatedPoints_,
+        interpolatedOldPoints_,
+        interpolationWeights_,
+        triPoints
+    );
 }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
index 6f64c9fe54..0c4f0f7811 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
@@ -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-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -120,7 +120,7 @@ Foam::isoSurface::adaptPatchFields
         {
             fvPatchField<Type>& pfld = const_cast<fvPatchField<Type>&>
             (
-                fld.boundaryField()[patchI]
+                sliceFld.boundaryField()[patchI]
             );
 
             const scalarField& w = mesh.weights().boundaryField()[patchI];
@@ -189,6 +189,8 @@ Type Foam::isoSurface::generatePoint
 }
 
 
+// Note: cannot use simpler isoSurfaceCell::generateTriPoints since
+//       the need here to sometimes pass in remote 'snappedPoints'
 template<class Type>
 void Foam::isoSurface::generateTriPoints
 (
@@ -240,8 +242,8 @@ void Foam::isoSurface::generateTriPoints
         case 0x0F:
         break;
 
-        case 0x0E:
         case 0x01:
+        case 0x0E:
             points.append
             (
                 generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1)
@@ -254,10 +256,16 @@ void Foam::isoSurface::generateTriPoints
             (
                 generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
             );
+            if (triIndex == 0x0E)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x0D:
         case 0x02:
+        case 0x0D:
             points.append
             (
                 generatePoint(s1,p1,hasSnap1,snapP1,s0,p0,hasSnap0,snapP0)
@@ -270,97 +278,133 @@ void Foam::isoSurface::generateTriPoints
             (
                 generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
             );
+            if (triIndex == 0x0D)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x0C:
         case 0x03:
-        {
-            Type tp1 =
-                generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
-            Type tp2 =
-                generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
+        case 0x0C:
+            {
+                Type p0p2 =
+                    generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
+                Type p1p3 =
+                    generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
 
-            points.append
-            (
-                generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
-            );
-            points.append(tp1);
-            points.append(tp2);
-            points.append(tp2);
-            points.append
-            (
-                generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
-            );
-            points.append(tp1);
-        }
+                points.append
+                (
+                    generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
+                );
+                points.append(p1p3);
+                points.append(p0p2);
+
+                points.append(p1p3);
+                points.append
+                (
+                    generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
+                );
+                points.append(p0p2);
+            }
+            if (triIndex == 0x0C)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-5], points[sz-4]);
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x0B:
         case 0x04:
-        {
-            points.append
-            (
-                generatePoint(s2,p2,hasSnap2,snapP2,s0,p0,hasSnap0,snapP0)
-            );
-            points.append
-            (
-                generatePoint(s2,p2,hasSnap2,snapP2,s1,p1,hasSnap1,snapP1)
-            );
-            points.append
-            (
-                generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3)
-            );
-        }
+        case 0x0B:
+            {
+                points.append
+                (
+                    generatePoint(s2,p2,hasSnap2,snapP2,s0,p0,hasSnap0,snapP0)
+                );
+                points.append
+                (
+                    generatePoint(s2,p2,hasSnap2,snapP2,s1,p1,hasSnap1,snapP1)
+                );
+                points.append
+                (
+                    generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3)
+                );
+            }
+            if (triIndex == 0x0B)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x0A:
         case 0x05:
-        {
-            Type tp0 =
-                generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
-            Type tp1 =
-                generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
+        case 0x0A:
+            {
+                Type p0p1 =
+                    generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
+                Type p2p3 =
+                    generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
 
-            points.append(tp0);
-            points.append(tp1);
-            points.append
-            (
-                generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
-            );
-            points.append(tp0);
-            points.append
-            (
-                generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
-            );
-            points.append(tp1);
-        }
+                points.append(p0p1);
+                points.append(p2p3);
+                points.append
+                (
+                    generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
+                );
+
+                points.append(p0p1);
+                points.append
+                (
+                    generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
+                );
+                points.append(p2p3);
+            }
+            if (triIndex == 0x0A)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-5], points[sz-4]);
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x09:
         case 0x06:
-        {
-            Type tp0 =
-                generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
-            Type tp1 =
-                generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
+        case 0x09:
+            {
+                Type p0p1 =
+                    generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
+                Type p2p3 =
+                    generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
 
-            points.append(tp0);
-            points.append
-            (
-                generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
-            );
-            points.append(tp1);
-            points.append(tp0);
-            points.append
-            (
-                generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
-            );
-            points.append(tp1);
-        }
+                points.append(p0p1);
+                points.append
+                (
+                    generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
+                );
+                points.append(p2p3);
+
+                points.append(p0p1);
+                points.append(p2p3);
+                points.append
+                (
+                    generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
+                );
+            }
+            if (triIndex == 0x09)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-5], points[sz-4]);
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x07:
         case 0x08:
+        case 0x07:
             points.append
             (
                 generatePoint(s3,p3,hasSnap3,snapP3,s0,p0,hasSnap0,snapP0)
@@ -373,6 +417,12 @@ void Foam::isoSurface::generateTriPoints
             (
                 generatePoint(s3,p3,hasSnap3,snapP3,s1,p1,hasSnap1,snapP1)
             );
+            if (triIndex == 0x07)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
     }
 }
@@ -689,6 +739,69 @@ void Foam::isoSurface::generateTriPoints
 //}
 
 
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::isoSurface::interpolate
+(
+    const label nPoints,
+    const labelList& triPointMergeMap,
+    const labelList& interpolatedPoints,
+    const List<FixedList<label, 3> >& interpolatedOldPoints,
+    const List<FixedList<scalar, 3> >& interpolationWeights,
+    const DynamicList<Type>& unmergedValues
+)
+{
+    // One value per point
+    tmp<Field<Type> > tvalues(new Field<Type>(nPoints, pTraits<Type>::zero));
+    Field<Type>& values = tvalues();
+
+
+    // Pass1: unweighted average of merged point values
+    {
+        labelList nValues(values.size(), 0);
+
+        forAll(unmergedValues, i)
+        {
+            label mergedPointI = triPointMergeMap[i];
+
+            if (mergedPointI >= 0)
+            {
+                values[mergedPointI] += unmergedValues[i];
+                nValues[mergedPointI]++;
+            }
+        }
+
+        forAll(values, i)
+        {
+            if (nValues[i] > 0)
+            {
+                values[i] /= scalar(nValues[i]);
+            }
+        }
+    }
+
+
+    // Pass2: weighted average for remaining values (from clipped triangles)
+
+    forAll(interpolatedPoints, i)
+    {
+        label pointI = interpolatedPoints[i];
+        const FixedList<label, 3>& oldPoints = interpolatedOldPoints[i];
+        const FixedList<scalar, 3>& w = interpolationWeights[i];
+
+        // Note: zeroing should not be necessary if interpolation only done
+        //       for newly introduced points (i.e. not in triPointMergeMap)
+        values[pointI] = pTraits<Type>::zero;
+        forAll(oldPoints, j)
+        {
+            values[pointI] = w[j]*unmergedValues[oldPoints[j]];
+        }
+    }
+
+    return tvalues;
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type> >
 Foam::isoSurface::interpolate
@@ -707,7 +820,7 @@ Foam::isoSurface::interpolate
     > > c2(adaptPatchFields(cCoords));
 
 
-    DynamicList<Type> triPoints(nCutCells_);
+    DynamicList<Type> triPoints(3*nCutCells_);
     DynamicList<label> triMeshCells(nCutCells_);
 
     // Dummy snap data
@@ -731,52 +844,15 @@ Foam::isoSurface::interpolate
         triMeshCells
     );
 
-
-    // One value per point
-    tmp<Field<Type> > tvalues
+    return interpolate
     (
-        new Field<Type>(points().size(), pTraits<Type>::zero)
+        points().size(),
+        triPointMergeMap_,
+        interpolatedPoints_,
+        interpolatedOldPoints_,
+        interpolationWeights_,
+        triPoints
     );
-    Field<Type>& values = tvalues();
-    labelList nValues(values.size(), 0);
-
-    forAll(triPoints, i)
-    {
-        label mergedPointI = triPointMergeMap_[i];
-
-        if (mergedPointI >= 0)
-        {
-            values[mergedPointI] += triPoints[i];
-            nValues[mergedPointI]++;
-        }
-    }
-
-    if (debug)
-    {
-        Pout<< "nValues:" << values.size() << endl;
-        label nMult = 0;
-        forAll(nValues, i)
-        {
-            if (nValues[i] == 0)
-            {
-                FatalErrorIn("isoSurface::interpolate(..)")
-                    << "point:" << i << " nValues:" << nValues[i]
-                    << abort(FatalError);
-            }
-            else if (nValues[i] > 1)
-            {
-                nMult++;
-            }
-        }
-        Pout<< "Of which mult:" << nMult << endl;
-    }
-
-    forAll(values, i)
-    {
-        values[i] /= scalar(nValues[i]);
-    }
-
-    return tvalues;
 }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index bb171c2c0e..0399eeae06 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -122,18 +122,52 @@ void Foam::sampledIsoSurface::getIsoFields() const
     // Get pointField
     // ~~~~~~~~~~~~~~
 
+    // In case of multiple iso values we don't want to calculate multiple e.g.
+    // "volPointInterpolate(p)" so register it and re-use it. This is the
+    // same as the 'cache' functionality from volPointInterpolate but
+    // unfortunately that one does not guarantee that the field pointer
+    // remain: e.g. some other
+    // functionObject might delete the cached version.
+    // (volPointInterpolation::interpolate with cache=false deletes any
+    //  registered one or if mesh.changing())
+
     if (!subMeshPtr_.valid())
     {
-        word pointFldName = "volPointInterpolate(" + isoField_ + ')';
+        const word pointFldName =
+            "volPointInterpolate_"
+          + type()
+          + "("
+          + isoField_
+          + ')';
 
         if (fvm.foundObject<pointScalarField>(pointFldName))
         {
             if (debug)
             {
-                Info<< "sampledIsoSurface::getIsoFields() : lookup pointField "
-                    << pointFldName << endl;
+                Info<< "sampledIsoSurface::getIsoFields() :"
+                    << " lookup pointField " << pointFldName << endl;
             }
-            pointFieldPtr_ = &fvm.lookupObject<pointScalarField>(pointFldName);
+            const pointScalarField& pfld = fvm.lookupObject<pointScalarField>
+            (
+                pointFldName
+            );
+
+            if (!pfld.upToDate(*volFieldPtr_))
+            {
+                if (debug)
+                {
+                    Info<< "sampledIsoSurface::getIsoFields() :"
+                        << " updating pointField " << pointFldName << endl;
+                }
+                // Update the interpolated value
+                volPointInterpolation::New(fvm).interpolate
+                (
+                    *volFieldPtr_,
+                    const_cast<pointScalarField&>(pfld)
+                );
+            }
+
+            pointFieldPtr_ = &pfld;
         }
         else
         {
@@ -141,33 +175,23 @@ void Foam::sampledIsoSurface::getIsoFields() const
 
             if (debug)
             {
-                Info<< "sampledIsoSurface::getIsoFields() : "
-                    << "checking pointField " << pointFldName
-                    << " for same time " << fvm.time().timeName()
-                    << endl;
+                Info<< "sampledIsoSurface::getIsoFields() :"
+                    << " creating and storing pointField "
+                    << pointFldName << " for time "
+                    << fvm.time().timeName() << endl;
             }
 
-            if
+            tmp<pointScalarField> tpfld
             (
-                storedPointFieldPtr_.empty()
-             || (fvm.time().timeName() != storedPointFieldPtr_().instance())
-            )
-            {
-                if (debug)
-                {
-                    Info<< "sampledIsoSurface::getIsoFields() :"
-                        << " interpolating volField " << volFieldPtr_->name()
-                        << " to get pointField " << pointFldName << endl;
-                }
-
-                storedPointFieldPtr_.reset
+                volPointInterpolation::New(fvm).interpolate
                 (
-                    volPointInterpolation::New(fvm)
-                    .interpolate(*volFieldPtr_).ptr()
-                );
-                storedPointFieldPtr_->checkOut();
-                pointFieldPtr_ = storedPointFieldPtr_.operator->();
-            }
+                    *volFieldPtr_,
+                    pointFldName,
+                    false
+                )
+            );
+            pointFieldPtr_ = tpfld.ptr();
+            const_cast<pointScalarField*>(pointFieldPtr_)->store();
         }
 
 
@@ -233,8 +257,10 @@ void Foam::sampledIsoSurface::getIsoFields() const
 
         // Pointfield on submesh
 
-        word pointFldName =
-            "volPointInterpolate("
+        const word pointFldName =
+            "volPointInterpolate_"
+          + type()
+          + "("
           + volSubFieldPtr_->name()
           + ')';
 
@@ -245,11 +271,28 @@ void Foam::sampledIsoSurface::getIsoFields() const
                 Info<< "sampledIsoSurface::getIsoFields() :"
                     << " submesh lookup pointField " << pointFldName << endl;
             }
-            storedPointSubFieldPtr_.clear();
-            pointSubFieldPtr_ = &subFvm.lookupObject<pointScalarField>
+            const pointScalarField& pfld = subFvm.lookupObject<pointScalarField>
             (
                 pointFldName
             );
+
+            if (!pfld.upToDate(*volSubFieldPtr_))
+            {
+                if (debug)
+                {
+                    Info<< "sampledIsoSurface::getIsoFields() :"
+                        << " updating submesh pointField "
+                        << pointFldName << endl;
+                }
+                // Update the interpolated value
+                volPointInterpolation::New(subFvm).interpolate
+                (
+                    *volSubFieldPtr_,
+                    const_cast<pointScalarField&>(pfld)
+                );
+            }
+
+            pointFieldPtr_ = &pfld;
         }
         else
         {
@@ -260,15 +303,15 @@ void Foam::sampledIsoSurface::getIsoFields() const
                     << volSubFieldPtr_->name()
                     << " to get submesh pointField " << pointFldName << endl;
             }
-            storedPointSubFieldPtr_.reset
+            tmp<pointScalarField> tpfld
             (
                 volPointInterpolation::New
                 (
                     subFvm
-                ).interpolate(*volSubFieldPtr_).ptr()
+                ).interpolate(*volSubFieldPtr_)
             );
-            storedPointSubFieldPtr_->checkOut();
-            pointSubFieldPtr_ = storedPointSubFieldPtr_.operator->();
+            pointSubFieldPtr_ = tpfld.ptr();
+            const_cast<pointScalarField*>(pointSubFieldPtr_)->store();
         }
 
 
@@ -349,28 +392,34 @@ bool Foam::sampledIsoSurface::updateGeometry() const
 
     if (subMeshPtr_.valid())
     {
+        const volScalarField& vfld = *volSubFieldPtr_;
+
         surfPtr_.reset
         (
             new isoSurface
             (
-                *volSubFieldPtr_,
+                vfld,
                 *pointSubFieldPtr_,
                 isoVal_,
                 regularise_,
+                bounds_,
                 mergeTol_
             )
         );
     }
     else
     {
+        const volScalarField& vfld = *volFieldPtr_;
+
         surfPtr_.reset
         (
             new isoSurface
             (
-                *volFieldPtr_,
+                vfld,
                 *pointFieldPtr_,
                 isoVal_,
                 regularise_,
+                bounds_,
                 mergeTol_
             )
         );
@@ -412,6 +461,7 @@ Foam::sampledIsoSurface::sampledIsoSurface
     sampledSurface(name, mesh, dict),
     isoField_(dict.lookup("isoField")),
     isoVal_(readScalar(dict.lookup("isoValue"))),
+    bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
     mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", false)),
@@ -422,7 +472,6 @@ Foam::sampledIsoSurface::sampledIsoSurface
     prevTimeIndex_(-1),
     storedVolFieldPtr_(NULL),
     volFieldPtr_(NULL),
-    storedPointFieldPtr_(NULL),
     pointFieldPtr_(NULL)
 {
     if (!sampledSurface::interpolate())
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
index 8cfd14c3c1..488ba56452 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
@@ -63,6 +63,9 @@ class sampledIsoSurface
         //- Iso value
         const scalar isoVal_;
 
+        //- Optional bounding box to trim triangles against
+        const boundBox bounds_;
+
         //- Merge tolerance
         const scalar mergeTol_;
 
@@ -94,7 +97,6 @@ class sampledIsoSurface
             mutable const volScalarField* volFieldPtr_;
 
             //- Cached pointfield
-            mutable autoPtr<pointScalarField> storedPointFieldPtr_;
             mutable const pointScalarField* pointFieldPtr_;
 
             // And on subsetted mesh
@@ -107,7 +109,6 @@ class sampledIsoSurface
                 mutable const volScalarField* volSubFieldPtr_;
 
                 //- Cached pointfield
-                mutable autoPtr<pointScalarField> storedPointSubFieldPtr_;
                 mutable const pointScalarField* pointSubFieldPtr_;
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
index 239217a352..4c21257ce7 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -148,7 +148,8 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
             cellAvg,
             pointFld().internalField(),
             isoVal_,
-            regularise_
+            regularise_,
+            bounds_
         );
 
         const_cast<sampledIsoSurfaceCell&>
@@ -166,7 +167,8 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
             cellFld.internalField(),
             pointFld().internalField(),
             isoVal_,
-            regularise_
+            regularise_,
+            bounds_
         );
 
         const_cast<sampledIsoSurfaceCell&>
@@ -185,6 +187,7 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
             << "    average        : " << average_ << nl
             << "    isoField       : " << isoField_ << nl
             << "    isoValue       : " << isoVal_ << nl
+            << "    bounds         : " << bounds_ << nl
             << "    points         : " << points().size() << nl
             << "    tris           : " << triSurface::size() << nl
             << "    cut cells      : " << meshCells_.size() << endl;
@@ -206,6 +209,7 @@ Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
     sampledSurface(name, mesh, dict),
     isoField_(dict.lookup("isoField")),
     isoVal_(readScalar(dict.lookup("isoValue"))),
+    bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", true)),
     zoneKey_(keyType::null),
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
index 7021bf0b97..0cfa6a7461 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
@@ -62,6 +62,9 @@ class sampledIsoSurfaceCell
         //- Iso value
         const scalar isoVal_;
 
+        //- Optional bounding box to trim triangles against
+        const boundBox bounds_;
+
         //- Whether to coarse
         const Switch regularise_;
 
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
index 7da10d7a57..7030b20d5e 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -230,6 +230,7 @@ void Foam::sampledCuttingPlane::createGeometry()
             pointDistance_,
             0.0,
             regularise_,
+            bounds_,
             mergeTol_
         )
         //new isoSurfaceCell
@@ -262,6 +263,7 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
 :
     sampledSurface(name, mesh, dict),
     plane_(dict),
+    bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
     mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", false)),
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
index a8aa477407..308bdcd46a 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
@@ -60,6 +60,9 @@ class sampledCuttingPlane
         //- Plane
         const plane plane_;
 
+        //- Optional bounding box to trim triangles against
+        const boundBox bounds_;
+
         //- Merge tolerance
         const scalar mergeTol_;
 
diff --git a/src/sampling/sampledSurface/sampledPatch/sampledPatch.C b/src/sampling/sampledSurface/sampledPatch/sampledPatch.C
index d07a6391f1..faeccd69a3 100644
--- a/src/sampling/sampledSurface/sampledPatch/sampledPatch.C
+++ b/src/sampling/sampledSurface/sampledPatch/sampledPatch.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,11 +84,7 @@ const Foam::labelList& Foam::sampledPatch::patchIDs() const
 {
     if (patchIDs_.empty())
     {
-        patchIDs_ = mesh().boundaryMesh().patchSet
-        (
-            patchNames_,
-            false
-        ).sortedToc();
+        patchIDs_ = mesh().boundaryMesh().patchSet(patchNames_).sortedToc();
     }
     return patchIDs_;
 }
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
index 47a9f9ebbc..4e12874c25 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -162,19 +162,7 @@ void Foam::sampledSurfaces::write()
 
         const label nFields = classifyFields();
 
-        if (Pstream::master())
-        {
-            if (debug)
-            {
-                Pout<< "Creating directory "
-                    << outputPath_/obr_.time().timeName() << nl << endl;
-
-            }
-
-            mkDir(outputPath_/obr_.time().timeName());
-        }
-
-        // Write geometry first if required,
+        // write geometry first if required,
         // or when no fields would otherwise be written
         if (nFields == 0 || formatter_->separateGeometry())
         {
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
index 094ae45cfc..654b3d438f 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
@@ -644,8 +644,25 @@ bool Foam::sampledTriSurfaceMesh::update()
         return false;
     }
 
+    // Calculate surface and mesh overlap bounding box
+    treeBoundBox bb
+    (
+        surface_.triSurface::points(),
+        surface_.triSurface::meshPoints()
+    );
+    bb.min() = max(bb.min(), mesh().bounds().min());
+    bb.max() = min(bb.max(), mesh().bounds().max());
+
+    // Extend a bit
+    const vector span(bb.span());
+
+    bb.min() -= 0.5*span;
+    bb.max() += 0.5*span;
+
+    bb.inflate(1e-6);
+
     // Mesh search engine, no triangulation of faces.
-    meshSearch meshSearcher(mesh(), polyMesh::FACE_PLANES);
+    meshSearch meshSearcher(mesh(), bb, polyMesh::FACE_PLANES);
 
     return update(meshSearcher);
 }