diff --git a/applications/utilities/surface/surfaceAutoPatch/Make/files b/applications/utilities/surface/surfaceAutoPatch/Make/files
deleted file mode 100644
index 38a0ebbbf2..0000000000
--- a/applications/utilities/surface/surfaceAutoPatch/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-surfaceAutoPatch.C
-
-EXE = $(FOAM_APPBIN)/surfaceAutoPatch
diff --git a/applications/utilities/surface/surfaceAutoPatch/Make/options b/applications/utilities/surface/surfaceAutoPatch/Make/options
deleted file mode 100644
index 9f08e8d2a8..0000000000
--- a/applications/utilities/surface/surfaceAutoPatch/Make/options
+++ /dev/null
@@ -1,7 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude
-
-EXE_LIBS = \
-    -lmeshTools \
-    -ltriSurface
diff --git a/applications/utilities/surface/surfaceAutoPatch/surfaceAutoPatch.C b/applications/utilities/surface/surfaceAutoPatch/surfaceAutoPatch.C
deleted file mode 100644
index f53ec82963..0000000000
--- a/applications/utilities/surface/surfaceAutoPatch/surfaceAutoPatch.C
+++ /dev/null
@@ -1,125 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 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/>.
-
-Application
-    surfaceAutoPatch
-
-Description
-    Patches surface according to feature angle. Like autoPatch.
-
-\*---------------------------------------------------------------------------*/
-
-#include "triangle.H"
-#include "triSurface.H"
-#include "argList.H"
-#include "surfaceFeatures.H"
-#include "treeBoundBox.H"
-#include "meshTools.H"
-#include "OFstream.H"
-
-using namespace Foam;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
-int main(int argc, char *argv[])
-{
-    argList::noParallel();
-    argList::validArgs.append("input surfaceFile");
-    argList::validArgs.append("output surfaceFile");
-    argList::validArgs.append("includedAngle [0..180]");
-    argList args(argc, argv);
-
-    const fileName inFileName  = args[1];
-    const fileName outFileName = args[2];
-    const scalar includedAngle = args.argRead<scalar>(3);
-
-    Info<< "Surface        : " << inFileName << nl << endl;
-
-
-    // Read
-    // ~~~~
-
-    Info<< "Reading : " << inFileName << endl;
-    triSurface surf(inFileName);
-
-    Info<< "Read surface:" << endl;
-    surf.writeStats(Info);
-    Info<< endl;
-
-
-
-    // Construct features from surface&featureangle
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    Info<< "Constructing feature set from included angle " << includedAngle
-        << endl;
-
-    surfaceFeatures set(surf, includedAngle);
-
-    Info<< nl
-        << "Feature set:" << nl
-        << "    feature points : " << set.featurePoints().size() << nl
-        << "    feature edges  : " << set.featureEdges().size() << nl
-        << "    of which" << nl
-        << "        region edges   : " << set.nRegionEdges() << nl
-        << "        external edges : " << set.nExternalEdges() << nl
-        << "        internal edges : " << set.nInternalEdges() << nl
-        << endl;
-
-    // Get per-edge status.
-    boolList borderEdge(surf.nEdges(), false);
-    forAll(set.featureEdges(), i)
-    {
-        borderEdge[set.featureEdges()[i]] = true;
-    }
-
-    labelList faceRegion(surf.size());
-    label nRegions = surf.markZones(borderEdge, faceRegion);
-
-    // Reregion triangles.
-    forAll(surf, i)
-    {
-        surf[i].region() = faceRegion[i];
-    }
-
-    // Create some patches
-    surf.patches().setSize(nRegions);
-
-    forAll(surf.patches(), patchI)
-    {
-        surf.patches()[patchI].name() = "patch" + Foam::name(patchI);
-        surf.patches()[patchI].geometricType() = "empty";
-    }
-
-
-    Info<< "Writing : " << outFileName << endl;
-    surf.write(outFileName, true);
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/Allwmake b/applications/utilities/surface/surfaceBooleanFeatures/Allwmake
new file mode 100755
index 0000000000..50d640f2b2
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/Allwmake
@@ -0,0 +1,15 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+set -x
+
+if [ -z "$CGAL_ARCH_PATH" ]
+then
+    export COMPILE_FLAGS="-DNO_CGAL"
+else
+    wmake PolyhedronReader
+    export COMPILE_FLAGS='-IPolyhedronReader'
+    export LINK_FLAGS='${CGAL_LIBS} -lPolyhedronReader'
+fi
+wmake
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/CGAL3DKernel.H b/applications/utilities/surface/surfaceBooleanFeatures/CGAL3DKernel.H
new file mode 100644
index 0000000000..e753060ce0
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/CGAL3DKernel.H
@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Typedefs
+    CGAL3DKernel
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef CGAL3DKernel_H
+#define CGAL3DKernel_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef CGAL_INEXACT
+
+    // Fast kernel using a double as the storage type
+    #include "CGAL/Exact_predicates_inexact_constructions_kernel.h"
+    typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+#else
+
+    // Very robust but expensive kernel
+    #include "CGAL/Exact_predicates_exact_constructions_kernel.h"
+    typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+
+#endif
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/CGALIndexedPolyhedron.H b/applications/utilities/surface/surfaceBooleanFeatures/CGALIndexedPolyhedron.H
new file mode 100644
index 0000000000..0eba1bee2f
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/CGALIndexedPolyhedron.H
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Typedefs
+    IndexedPolyhedron
+
+Description
+    CGAL data structures used for triSurface handling
+
+    Define CGAL_INEXACT to use Exact_predicates_inexact_constructions kernel
+    otherwise the more robust but much less efficient
+    Exact_predicates_exact_constructions will be used.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IndexedPolyhedron_H
+#define IndexedPolyhedron_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Include uint.H before CGAL headers to define __STDC_LIMIT_MACROS
+#include "uint.H"
+
+#include "CGAL3DKernel.H"
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Nef_polyhedron_3.h>
+
+#include "label.H"
+
+
+typedef CGAL::Point_3<K>  Point;
+typedef CGAL::Segment_3<K> Segment;
+typedef CGAL::Direction_3<K>  Direction;
+typedef CGAL::Plane_3<K>  Plane;
+typedef CGAL::Triangle_3<K>  Triangle;
+
+// Define new class with color and define the polyhedron types
+template<class Refs>
+struct IndexedFace
+:
+    public CGAL::HalfedgeDS_face_base<Refs>
+{
+    Foam::label index;
+    Foam::label region;
+};
+struct My_items
+:
+    public CGAL::Polyhedron_items_3
+{
+    template<class Refs, class Traits>
+    struct Face_wrapper
+    {
+        typedef IndexedFace<Refs> Face;
+    };
+};
+
+
+//typedef CGAL::Polyhedron_3<K>  Polyhedron;
+typedef CGAL::Polyhedron_3<K, My_items>  Polyhedron;
+
+typedef Polyhedron::HalfedgeDS HalfedgeDS;
+typedef Polyhedron::Edge_iterator Edge_iterator;
+typedef Polyhedron::Vertex Vertex;
+typedef Polyhedron::Vertex_iterator Vertex_iterator;
+typedef Polyhedron::Halfedge_handle Halfedge_handle;
+typedef Polyhedron::Edge_iterator Edge_iterator;
+typedef Polyhedron::Facet_iterator Facet_iterator;
+typedef Polyhedron::Halfedge_around_facet_const_circulator HFCC;
+typedef Polyhedron::Vertex_const_iterator                  VCI;
+
+typedef CGAL::Nef_polyhedron_3<K>  Nef_polyhedron;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/Make/options b/applications/utilities/surface/surfaceBooleanFeatures/Make/options
index 727dbd1c2b..ed05fef039 100644
--- a/applications/utilities/surface/surfaceBooleanFeatures/Make/options
+++ b/applications/utilities/surface/surfaceBooleanFeatures/Make/options
@@ -1,9 +1,29 @@
+EXE_NDEBUG = -DNDEBUG
+/* EXE_NDEBUG = -g -O0 -DFULLDEBUG */
+
+
+c++CGALWARN = -Wno-old-style-cast
+
+/*-- Define NO_CGAL to avoid using CGAL altogether */
+/*-- Define CGAL_INEXACT to use inexact CGAL constructions */
+
+include $(GENERAL_RULES)/CGAL
+
 EXE_INC = \
-    -I$(LIB_SRC)/triSurface/lnInclude \
+    ${ROUNDING_MATH} \
+    ${EXE_NDEBUG} \
+    ${CGAL_INC} \
+    ${c++CGALWARN} \
+    $(COMPILE_FLAGS) \
+    -IPolyhedronReader \
+    -I$(FOAM_SRC)/surfMesh/lnInclude \
+    -I$(FOAM_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
+    -lsurfMesh \
     -ltriSurface \
     -ledgeMesh \
-    -lmeshTools
+    -lmeshTools \
+    $(LINK_FLAGS)
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/Make/files b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/Make/files
new file mode 100644
index 0000000000..d0e01d2747
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/Make/files
@@ -0,0 +1,3 @@
+PolyhedronReader.C
+
+LIB = $(FOAM_LIBBIN)/libPolyhedronReader
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/Make/options b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/Make/options
new file mode 100644
index 0000000000..6603adda0f
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/Make/options
@@ -0,0 +1,25 @@
+EXE_NDEBUG = -DNDEBUG
+/* EXE_NDEBUG = -g -O0 -DFULLDEBUG */
+
+
+c++CGALWARN = -Wno-old-style-cast
+
+/*-- Define CGAL_INEXACT to use inexact CGAL constructions */
+
+include $(GENERAL_RULES)/CGAL
+
+EXE_INC = \
+    ${ROUNDING_MATH} \
+    ${EXE_NDEBUG} \
+    ${CGAL_INC} \
+    ${c++CGALWARN} \
+    -I.. \
+    -I$(FOAM_SRC)/surfMesh/lnInclude \
+    -I$(FOAM_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/edgeMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I/usr/include/Qt
+
+LIB_LIBS = \
+    -L${CGAL_ARCH_PATH}/lib \
+    -ltriSurface
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReader.C b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReader.C
new file mode 100644
index 0000000000..462b9a6f05
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReader.C
@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "PolyhedronReader.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::PolyhedronReader::PolyhedronReader(const triSurface& s, Polyhedron& p)
+{
+    Build_triangle<HalfedgeDS> triangle(s);
+    p.delegate(triangle);
+    // Populate index and region
+    Foam::label nTris = 0;
+    for
+    (
+        Facet_iterator fi = p.facets_begin();
+        fi != p.facets_end();
+        ++fi
+    )
+    {
+        fi->index = nTris++;
+        fi->region = s[fi->index].region();
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReader.H b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReader.H
new file mode 100644
index 0000000000..ef3e88a0ac
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReader.H
@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::PolyhedronReader
+
+Description
+
+SourceFiles
+    PolyhedronReader.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PolyhedronReader_H
+#define PolyhedronReader_H
+
+#include "CGALIndexedPolyhedron.H"
+#include "triSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class PolyhedronReader Declaration
+\*---------------------------------------------------------------------------*/
+
+class PolyhedronReader
+{
+    // Private Classes
+
+        template<class HDS>
+        class Build_triangle
+        :
+            public CGAL::Modifier_base<HDS>
+        {
+            const triSurface& s_;
+
+        public:
+
+            Build_triangle(const triSurface& s);
+
+            void operator()(HDS& hds);
+        };
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        PolyhedronReader(const PolyhedronReader&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const PolyhedronReader&);
+
+
+public:
+
+    // Constructors
+
+        //- Populate polyhedron from surface
+        PolyhedronReader(const triSurface& s, Polyhedron& p);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "PolyhedronReaderTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReaderTemplates.C b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReaderTemplates.C
new file mode 100644
index 0000000000..53ea3b1b87
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReaderTemplates.C
@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "PolyhedronReader.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class HDS>
+Foam::PolyhedronReader::Build_triangle<HDS>::Build_triangle
+(
+    const triSurface& s
+)
+:
+    s_(s)
+{}
+
+
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+template<class HDS>
+void Foam::PolyhedronReader::Build_triangle<HDS>::operator()(HDS& hds)
+{
+    // Postcondition: hds is a valid polyhedral surface.
+    CGAL::Polyhedron_incremental_builder_3<HDS> B(hds, true);
+
+    B.begin_surface(s_.nPoints(), s_.size());
+
+    typedef typename HDS::Vertex Vertex;
+    typedef typename Vertex::Point Point;
+
+    const Foam::pointField& pts = s_.points();
+    forAll(pts, i)
+    {
+        const Foam::point& pt = pts[i];
+        B.add_vertex(Point(pt.x(), pt.y(), pt.z()));
+    }
+    forAll(s_, i)
+    {
+        const Foam::labelledTri& t = s_[i];
+        B.begin_facet();
+        B.add_vertex_to_facet(t[0]);
+        B.add_vertex_to_facet(t[1]);
+        B.add_vertex_to_facet(t[2]);
+        B.end_facet();
+    }
+    B.end_surface();
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C b/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
index 6c4e67dd7d..e948cd1fdb 100644
--- a/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
+++ b/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -74,11 +74,34 @@ Description
 #include "featureEdgeMesh.H"
 #include "extendedFeatureEdgeMesh.H"
 #include "triSurfaceSearch.H"
+#include "triSurfaceMesh.H"
 #include "OFstream.H"
+#include "OBJstream.H"
 #include "booleanSurface.H"
 #include "edgeIntersections.H"
 #include "meshTools.H"
-#include "labelPair.H"
+#include "DynamicField.H"
+
+
+#ifndef NO_CGAL
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+#include "CGALIndexedPolyhedron.H"
+#include "PolyhedronReader.H"
+typedef CGAL::AABB_face_graph_triangle_primitive
+<
+    Polyhedron, CGAL::Default, CGAL::Tag_false
+> Primitive;
+typedef CGAL::AABB_traits<K, Primitive> Traits;
+typedef CGAL::AABB_tree<Traits> Tree;
+
+typedef boost::optional<Tree::Intersection_and_primitive_id<Segment>::Type >
+Segment_intersection;
+
+#endif // NO_CGAL
+
 
 using namespace Foam;
 
@@ -105,7 +128,7 @@ bool intersectSurfaces
 
             scalarField surfPointTol
             (
-                1e-3*edgeIntersections::minEdgeLength(surf)
+                max(1e-3*edgeIntersections::minEdgeLength(surf), SMALL)
             );
 
             // Determine raw intersections
@@ -178,7 +201,7 @@ bool intersectSurfaces
 
             scalarField surf1PointTol
             (
-                1e-3*edgeIntersections::minEdgeLength(surf1)
+                max(1e-3*edgeIntersections::minEdgeLength(surf1), SMALL)
             );
 
             // Determine raw intersections
@@ -221,7 +244,7 @@ bool intersectSurfaces
 
             scalarField surf2PointTol
             (
-                1e-3*edgeIntersections::minEdgeLength(surf2)
+                max(1e-3*edgeIntersections::minEdgeLength(surf2), SMALL)
             );
 
             // Determine raw intersections
@@ -311,8 +334,8 @@ void calcEdgeCuts
     triSurface& surf1,
     triSurface& surf2,
     const bool perturb,
-    edgeIntersections& edge1Cuts,
-    edgeIntersections& edge2Cuts
+    edgeIntersections& edgeCuts1,
+    edgeIntersections& edgeCuts2
 )
 {
     if (perturb)
@@ -320,9 +343,9 @@ void calcEdgeCuts
         intersectSurfaces
         (
             surf1,
-            edge1Cuts,
+            edgeCuts1,
             surf2,
-            edge2Cuts
+            edgeCuts2
         );
     }
     else
@@ -332,11 +355,11 @@ void calcEdgeCuts
         Info<< "Determining intersections of surf1 edges with surf2 faces"
             << endl;
 
-        edge1Cuts = edgeIntersections
+        edgeCuts1 = edgeIntersections
         (
             surf1,
             querySurf2,
-            1e-3*edgeIntersections::minEdgeLength(surf1)
+            max(1e-3*edgeIntersections::minEdgeLength(surf1), SMALL)
         );
 
         triSurfaceSearch querySurf1(surf1);
@@ -344,161 +367,875 @@ void calcEdgeCuts
         Info<< "Determining intersections of surf2 edges with surf1 faces"
             << endl;
 
-        edge2Cuts = edgeIntersections
+        edgeCuts2 = edgeIntersections
         (
             surf2,
             querySurf1,
-            1e-3*edgeIntersections::minEdgeLength(surf2)
+            max(1e-3*edgeIntersections::minEdgeLength(surf2), SMALL)
         );
     }
 }
 
 
-void calcFeaturePoints(const pointField& points, const edgeList& edges)
+// CGAL variants
+
+#ifndef NO_CGAL
+
+void visitPointRegion
+(
+    const triSurface& s,
+    const label zoneI,
+    const label pointI,
+    const label startEdgeI,
+    const label startFaceI,
+    labelList& pFacesZone
+)
 {
-    edgeMesh eMesh(points, edges);
+    const labelList& eFaces = s.edgeFaces()[startEdgeI];
 
-    const labelListList& pointEdges = eMesh.pointEdges();
-
-    // Get total number of feature points
-    label nFeaturePoints = 0;
-    forAll(pointEdges, pI)
+    if (eFaces.size() == 2)
     {
-        const labelList& pEdges = pointEdges[pI];
-
-        if (pEdges.size() == 1)
+        label nextFaceI;
+        if (eFaces[0] == startFaceI)
         {
-            nFeaturePoints++;
+            nextFaceI = eFaces[1];
         }
-    }
-
-
-    // Calculate addressing from feature point to cut point and cut edge
-    labelList featurePointToCutPoint(nFeaturePoints);
-    labelList featurePointToCutEdge(nFeaturePoints);
-
-    label nFeatPts = 0;
-    forAll(pointEdges, pI)
-    {
-        const labelList& pEdges = pointEdges[pI];
-
-        if (pEdges.size() == 1)
+        else if (eFaces[1] == startFaceI)
         {
-            featurePointToCutPoint[nFeatPts] = pI;
-            featurePointToCutEdge[nFeatPts] = pEdges[0];
-            nFeatPts++;
+            nextFaceI = eFaces[0];
+        }
+        else
+        {
+            FatalErrorIn("visitPointRegion(..)")
+                << "problem" << exit(FatalError);
+            nextFaceI = -1;
+        }
+
+
+
+        label index = findIndex(s.pointFaces()[pointI], nextFaceI);
+
+        if (pFacesZone[index] == -1)
+        {
+            // Mark face as been visited.
+            pFacesZone[index] = zoneI;
+
+            // Step to next edge on face which is still using pointI
+            const labelList& fEdges = s.faceEdges()[nextFaceI];
+
+            label nextEdgeI = -1;
+
+            forAll(fEdges, i)
+            {
+                label edgeI = fEdges[i];
+                const edge& e = s.edges()[edgeI];
+
+                if (edgeI != startEdgeI && (e[0] == pointI || e[1] == pointI))
+                {
+                    nextEdgeI = edgeI;
+
+                    break;
+                }
+            }
+
+            if (nextEdgeI == -1)
+            {
+                FatalErrorIn("visitPointRegion()")
+                    << "Problem: cannot find edge out of " << fEdges
+                    << "on face " << nextFaceI << " that uses point " << pointI
+                    << " and is not edge " << startEdgeI << abort(FatalError);
+            }
+
+
+            visitPointRegion
+            (
+                s,
+                zoneI,
+                pointI,
+                nextEdgeI,
+                nextFaceI,
+                pFacesZone
+            );
         }
     }
 }
 
 
-int main(int argc, char *argv[])
+label dupNonManifoldPoints(triSurface& s, labelList& pointMap)
 {
-    argList::noParallel();
-    argList::validArgs.append("action");
-    argList::validArgs.append("surface file");
-    argList::validArgs.append("surface file");
+    const labelListList& pf = s.pointFaces();
+    const labelListList& fe = s.faceEdges();
+    const edgeList& edges = s.edges();
 
-    argList::addBoolOption
-    (
-        "surf1Baffle",
-        "Mark surface 1 as a baffle"
-    );
 
-    argList::addBoolOption
-    (
-        "surf2Baffle",
-        "Mark surface 2 as a baffle"
-    );
+    DynamicField<point> newPoints(s.points());
+    // From dupSurf back to s.pointa
+    DynamicList<label> newPointMap(identity(newPoints.size()));
+    List<labelledTri> newFaces(s);
+    label nNonManifold = 0;
 
-    argList::addBoolOption
-    (
-        "perturb",
-        "Perturb surface points to escape degenerate intersections"
-    );
 
-    argList::addBoolOption
-    (
-        "invertedSpace",
-        "do the surfaces have inverted space orientation, "
-        "i.e. a point at infinity is considered inside. "
-        "This is only sensible for union and intersection."
-    );
-
-    #include "setRootCase.H"
-    #include "createTime.H"
-
-    word action(args.args()[1]);
-
-    HashTable<booleanSurface::booleanOpType> validActions;
-    validActions.insert("intersection", booleanSurface::INTERSECTION);
-    validActions.insert("union", booleanSurface::UNION);
-    validActions.insert("difference", booleanSurface::DIFFERENCE);
-
-    if (!validActions.found(action))
+    forAll(pf, pointI)
     {
-        FatalErrorIn(args.executable())
-            << "Unsupported action " << action << endl
-            << "Supported actions:" << validActions.toc() << abort(FatalError);
+        const labelList& pFaces = pf[pointI];
+
+        // Visited faces (as indices into pFaces)
+        labelList pFacesZone(pFaces.size(), -1);
+
+        label nZones = 0;
+        label index = 0;
+        for (; index < pFacesZone.size(); index++)
+        {
+            if (pFacesZone[index] == -1)
+            {
+                label zoneI = nZones++;
+                pFacesZone[index] = zoneI;
+
+                label faceI = pFaces[index];
+                const labelList& fEdges = fe[faceI];
+
+                // Find starting edge
+                forAll(fEdges, fEdgeI)
+                {
+                    label edgeI = fEdges[fEdgeI];
+                    const edge& e = edges[edgeI];
+
+                    if (e[0] == pointI || e[1] == pointI)
+                    {
+                        visitPointRegion
+                        (
+                            s,
+                            zoneI,
+                            pointI,
+                            edgeI,
+                            faceI,
+                            pFacesZone
+                        );
+                    }
+                }
+            }
+        }
+
+
+        // Subset
+        if (nZones > 1)
+        {
+            for (label zoneI = 1; zoneI < nZones; zoneI++)
+            {
+                label newPointI = newPoints.size();
+                newPointMap.append(s.meshPoints()[pointI]);
+                newPoints.append(s.points()[s.meshPoints()[pointI]]);
+
+                forAll(pFacesZone, index)
+                {
+                    if (pFacesZone[index] == zoneI)
+                    {
+                        label faceI = pFaces[index];
+                        const labelledTri& localF = s.localFaces()[faceI];
+                        forAll(localF, fp)
+                        {
+                            if (localF[fp] == pointI)
+                            {
+                                newFaces[faceI][fp] = newPointI;
+                            }
+                        }
+                    }
+                }
+            }
+            nNonManifold++;
+        }
     }
 
-    fileName surf1Name(args.args()[2]);
-    Info<< "Reading surface " << surf1Name << endl;
-    triSurface surf1(surf1Name);
 
-    Info<< surf1Name << " statistics:" << endl;
-    surf1.writeStats(Info);
-    Info<< endl;
-
-    fileName surf2Name(args[3]);
-    Info<< "Reading surface " << surf2Name << endl;
-    triSurface surf2(surf2Name);
-
-    Info<< surf2Name << " statistics:" << endl;
-    surf2.writeStats(Info);
-    Info<< endl;
-
-    const bool surf1Baffle = args.optionFound("surf1Baffle");
-    const bool surf2Baffle = args.optionFound("surf2Baffle");
-
-    edgeIntersections edge1Cuts;
-    edgeIntersections edge2Cuts;
-
-    bool invertedSpace = args.optionFound("invertedSpace");
-
-    if (invertedSpace && validActions[action] == booleanSurface::DIFFERENCE)
+    Info<< "Duplicating " << nNonManifold << " points out of " << s.nPoints()
+        << endl;
+    if (nNonManifold > 0)
     {
-        FatalErrorIn(args.executable())
-            << "Inverted space only makes sense for union or intersection."
-            << exit(FatalError);
+        triSurface dupSurf(newFaces, s.patches(), newPoints, true);
+
+        // Create map from dupSurf localPoints to s.localPoints
+        const Map<label> mpm = s.meshPointMap();
+
+        const labelList& dupMp = dupSurf.meshPoints();
+
+        labelList dupPointMap(dupMp.size());
+        forAll(dupMp, pointI)
+        {
+            label dupMeshPointI = dupMp[pointI];
+            label meshPointI = newPointMap[dupMeshPointI];
+            dupPointMap[pointI] = mpm[meshPointI];
+        }
+
+
+        forAll(dupPointMap, pointI)
+        {
+            const point& dupPt = dupSurf.points()[dupMp[pointI]];
+            label sLocalPointI = dupPointMap[pointI];
+            label sMeshPointI = s.meshPoints()[sLocalPointI];
+            const point& sPt = s.points()[sMeshPointI];
+
+            if (mag(dupPt-sPt) > SMALL)
+            {
+                FatalErrorIn("dupNonManifoldPoints(..)")
+                    << "dupPt:" << dupPt
+                    << " sPt:" << sPt
+                    << exit(FatalError);
+            }
+        }
+
+
+        //s.transfer(dupSurf);
+        s = dupSurf;
+        pointMap = UIndirectList<label>(pointMap, dupPointMap)();
     }
 
-    // Calculate the points where the edges are cut by the other surface
-    calcEdgeCuts
-    (
-        surf1,
-        surf2,
-        args.optionFound("perturb"),
-        edge1Cuts,
-        edge2Cuts
-    );
+    return nNonManifold;
+}
 
+
+// Find intersections of surf1 by edges of surf2. Return number of degenerate
+// cuts (cuts through faces/edges/points)
+labelPair edgeIntersectionsCGAL
+(
+    const Tree& tree,
+    const triSurface& surf,
+    const pointField& points,
+    edgeIntersections& edgeCuts
+)
+{
+    const edgeList& edges = surf.edges();
+    const labelList& meshPoints = surf.meshPoints();
+
+    //Info<< "Intersecting CGAL surface ..." << endl;
+    List<List<pointIndexHit> > intersections(edges.size());
+    labelListList classifications(edges.size());
+
+    label nPoints = 0;
+    label nSegments = 0;
+
+    std::vector<Segment_intersection> segments;
+    forAll(edges, edgeI)
+    {
+        const edge& e = edges[edgeI];
+
+        const point& a = points[meshPoints[e[0]]];
+        const point& b = points[meshPoints[e[1]]];
+
+        K::Segment_3 segment_query
+        (
+            Point(a.x(), a.y(), a.z()),
+            Point(b.x(), b.y(), b.z())
+        );
+
+        segments.clear();
+        tree.all_intersections(segment_query, std::back_inserter(segments));
+
+        for
+        (
+            std::vector<Segment_intersection>::const_iterator iter =
+                segments.begin(),
+            end = segments.end();
+            iter != end;
+            ++iter
+        )
+        {
+            // Get intersection object
+            if (const Point* ptPtr = boost::get<Point>(&((*iter)->first)))
+            {
+                point pt
+                (
+                    CGAL::to_double(ptPtr->x()),
+                    CGAL::to_double(ptPtr->y()),
+                    CGAL::to_double(ptPtr->z())
+                );
+
+                Polyhedron::Face_handle f = (*iter)->second;
+
+                intersections[edgeI].append
+                (
+                    pointIndexHit
+                    (
+                        true,
+                        pt,
+                        f->index
+                    )
+                );
+                // Intersection on edge interior
+                classifications[edgeI].append(-1);
+                nPoints++;
+            }
+            else if
+            (
+                const Segment* sPtr = boost::get<Segment>(&((*iter)->first))
+            )
+            {
+                //std::cout
+                //    << "intersection object is a segment:" << sPtr->source()
+                //    << " " << sPtr->target() << std::endl;
+
+                Polyhedron::Face_handle f = (*iter)->second;
+                //std::cout<< "triangle:" << f->index
+                //    << " region:" << f->region << std::endl;
+
+                const point source
+                (
+                    CGAL::to_double(sPtr->source().x()),
+                    CGAL::to_double(sPtr->source().y()),
+                    CGAL::to_double(sPtr->source().z())
+                );
+
+                const point target
+                (
+                    CGAL::to_double(sPtr->target().x()),
+                    CGAL::to_double(sPtr->target().y()),
+                    CGAL::to_double(sPtr->target().z())
+                );
+
+                // Make up some intersection point
+                intersections[edgeI].append
+                (
+                    pointIndexHit
+                    (
+                        true,
+                        0.5*(source+target),
+                        f->index
+                    )
+                );
+                // Intersection aligned with face. Tbd: enums
+                classifications[edgeI].append(2);
+                nSegments++;
+            }
+        }
+    }
+
+    edgeCuts = edgeIntersections(intersections, classifications);
+
+    return labelPair(nPoints, nSegments);
+}
+
+
+// Intersect edges of surf1 until no more degenerate intersections. Return
+// number of degenerates
+labelPair edgeIntersectionsAndShuffleCGAL
+(
+    Random& rndGen,
+    const triSurface& surf2,
+    const scalarField& surf1PointTol,
+    triSurface& surf1,
+    edgeIntersections& edgeCuts1
+)
+{
+    //Info<< "Constructing CGAL surface ..." << endl;
+    Polyhedron p;
+    PolyhedronReader(surf2, p);
+
+
+    //Info<< "Constructing CGAL tree ..." << endl;
+    const Tree tree(p.facets_begin(), p.facets_end(), p);
+
+
+    labelPair cuts(0, 0);
+
+    // Update surface1 points until no longer intersecting
+    pointField surf1Points(surf1.points());
+
+    for (label iter = 0; iter < 5; iter++)
+    {
+        // See which edges of 1 intersect 2
+        Info<< "Determining intersections of " << surf1.nEdges()
+            << " edges with surface of " << label(tree.size()) << " triangles"
+            << endl;
+        cuts = edgeIntersectionsCGAL
+        (
+            tree,
+            surf1,
+            surf1Points,
+            edgeCuts1
+        );
+        Info<< "Determined intersections:" << nl
+            << "    edges               : " << surf1.nEdges() << nl
+            << "    non-degenerate cuts : " << cuts.first() << nl
+            << "    degenerate cuts     : " << cuts.second() << nl
+            << endl;
+
+        if (cuts.second() == 0)
+        {
+            break;
+        }
+
+        Info<< "Shuffling conflicting points" << endl;
+
+        const labelListList& edgeStat = edgeCuts1.classification();
+        const edgeList& edges = surf1.edges();
+        const labelList& mp = surf1.meshPoints();
+        const point p05(0.5, 0.5, 0.5);
+
+        forAll(edgeStat, edgeI)
+        {
+            const labelList& stat = edgeStat[edgeI];
+            forAll(stat, i)
+            {
+                if (stat[i] == 2)
+                {
+                    const edge& e = edges[edgeI];
+                    forAll(e, eI)
+                    {
+                        vector d = rndGen.vector01()-p05;
+                        surf1Points[mp[e[eI]]] += surf1PointTol[e[eI]]*d;
+                    }
+                }
+            }
+        }
+    }
+    surf1.movePoints(surf1Points);
+
+    return cuts;
+}
+
+
+// Return map from subSurf.edges() back to surf.edges()
+labelList matchEdges
+(
+    const triSurface& surf,
+    const triSurface& subSurf,
+    const labelList& pointMap
+)
+{
+    if (pointMap.size() != subSurf.nPoints())
+    {
+        FatalErrorIn("findEdges(..)") << "problem" << exit(FatalError);
+    }
+
+    labelList edgeMap(subSurf.nEdges(), -1);
+
+    const edgeList& edges = surf.edges();
+    const labelListList& pointEdges = surf.pointEdges();
+
+    const edgeList& subEdges = subSurf.edges();
+
+
+    forAll(subEdges, subEdgeI)
+    {
+        const edge& subE = subEdges[subEdgeI];
+
+        // Match points on edge to those on surf
+        label start = pointMap[subE[0]];
+        label end = pointMap[subE[1]];
+        const labelList& pEdges = pointEdges[start];
+        forAll(pEdges, pEdgeI)
+        {
+            label edgeI = pEdges[pEdgeI];
+            const edge& e = edges[edgeI];
+
+            if (e.otherVertex(start) == end)
+            {
+                if (edgeMap[subEdgeI] == -1)
+                {
+                    edgeMap[subEdgeI] = edgeI;
+                }
+                else if (edgeMap[subEdgeI] != edgeI)
+                {
+                    WarningIn("findEdges(..)") << "sub edge "
+                        << subE << " points:"
+                        << subE.line(subSurf.localPoints())
+                        << " matches to " << edgeI
+                        << " points:" << e.line(surf.localPoints())
+                        << " but also " << edgeMap[subEdgeI]
+                        << " points:"
+                        << edges[edgeMap[subEdgeI]].line(surf.localPoints())
+                        << exit(FatalError);
+                }
+                break;
+            }
+        }
+
+        if (edgeMap[subEdgeI] == -1)
+        {
+            FatalErrorIn("findEdges(..)") << "did not find edge matching "
+                << subE << " at:" << subSurf.localPoints()[subE[0]]
+                << subSurf.localPoints()[subE[1]]
+                << exit(FatalError);
+        }
+    }
+
+    return edgeMap;
+}
+
+
+void calcEdgeCutsCGAL
+(
+    triSurface& surf1,
+    triSurface& surf2,
+    const bool perturb,
+    edgeIntersections& edgeCuts1,
+    edgeIntersections& edgeCuts2
+)
+{
+    if (!perturb)
+    {
+        // See which edges of 1 intersect 2
+        {
+            Info<< "Constructing CGAL surface ..." << endl;
+            Polyhedron p;
+            PolyhedronReader(surf2, p);
+
+            Info<< "Constructing CGAL tree ..." << endl;
+            const Tree tree(p.facets_begin(), p.facets_end(), p);
+
+            edgeIntersectionsCGAL
+            (
+                tree,
+                surf1,
+                surf1.points(),
+                edgeCuts1
+            );
+        }
+        // See which edges of 2 intersect 1
+        {
+            Info<< "Constructing CGAL surface ..." << endl;
+            Polyhedron p;
+            PolyhedronReader(surf1, p);
+
+            Info<< "Constructing CGAL tree ..." << endl;
+            const Tree tree(p.facets_begin(), p.facets_end(), p);
+
+            edgeIntersectionsCGAL
+            (
+                tree,
+                surf2,
+                surf2.points(),
+                edgeCuts2
+            );
+        }
+    }
+    else
+    {
+        const scalarField surf1PointTol
+        (
+            max(1e-8*edgeIntersections::minEdgeLength(surf1), SMALL)
+        );
+        const scalarField surf2PointTol
+        (
+            max(1e-8*edgeIntersections::minEdgeLength(surf2), SMALL)
+        );
+
+
+        Random rndGen(0);
+
+        labelPair cuts1;
+        labelPair cuts2;
+
+        for (label iter = 0; iter < 10; iter++)
+        {
+            // Find intersections of surf1 edges with surf2 triangles
+            cuts1 = edgeIntersectionsAndShuffleCGAL
+            (
+                rndGen,
+                surf2,
+                surf1PointTol,
+                surf1,
+                edgeCuts1
+            );
+
+            // Find intersections of surf2 edges with surf1 triangles
+            cuts2 = edgeIntersectionsAndShuffleCGAL
+            (
+                rndGen,
+                surf1,
+                surf2PointTol,
+                surf2,
+                edgeCuts2
+            );
+
+            if (cuts1.second() + cuts2.second() == 0)
+            {
+                break;
+            }
+        }
+    }
+}
+
+
+void calcEdgeCutsBitsCGAL
+(
+    triSurface& surf1,
+    triSurface& surf2,
+    const bool perturb,
+    edgeIntersections& edgeCuts1,
+    edgeIntersections& edgeCuts2
+)
+{
+    label nZones1 = 1;
+    labelList zone1;
+    {
+        labelHashSet orientationEdge(surf1.size()/1000);
+        PatchTools::checkOrientation(surf1, false, &orientationEdge);
+        nZones1 = PatchTools::markZones(surf1, orientationEdge, zone1);
+
+        Info<< "Surface triangles " << surf1.size()
+            << " number of zones (orientation or non-manifold):"
+            << nZones1 << endl;
+    }
+
+    label nZones2 = 1;
+    labelList zone2;
+    {
+        labelHashSet orientationEdge(surf2.size()/1000);
+        PatchTools::checkOrientation(surf2, false, &orientationEdge);
+        nZones2 = PatchTools::markZones(surf2, orientationEdge, zone2);
+
+        Info<< "Surface triangles " << surf2.size()
+            << " number of zones (orientation or non-manifold):"
+            << nZones2 << endl;
+    }
+
+
+    if (nZones1 == 1 && nZones2 == 1)
+    {
+        calcEdgeCutsCGAL
+        (
+            surf1,
+            surf2,
+            perturb,
+            edgeCuts1,
+            edgeCuts2
+        );
+    }
+    else
+    {
+        edgeCuts1 = edgeIntersections
+        (
+            List<List<pointIndexHit> >(surf1.nEdges()),
+            labelListList(surf1.nEdges())
+        );
+        edgeCuts2 = edgeIntersections
+        (
+            List<List<pointIndexHit> >(surf2.nEdges()),
+            labelListList(surf2.nEdges())
+        );
+
+
+        for (label zone1I = 0; zone1I < nZones1; zone1I++)
+        {
+            // Generate sub surface for zone1I
+            boolList includeMap1(surf1.size(), false);
+
+            forAll(zone1, faceI)
+            {
+                if (zone1[faceI] == zone1I)
+                {
+                    includeMap1[faceI] = true;
+                }
+            }
+
+            // Subset. Map from local points on subset to local points on
+            // original
+            labelList pointMap1;
+            labelList faceMap1;
+            triSurface subSurf1
+            (
+                surf1.subsetMesh
+                (
+                    includeMap1,
+                    pointMap1,
+                    faceMap1
+                )
+            );
+
+            // Split non-manifold points; update pointMap
+            dupNonManifoldPoints(subSurf1, pointMap1);
+
+            const boundBox subBb1
+            (
+                subSurf1.points(),
+                subSurf1.meshPoints(),
+                false
+            );
+
+            const labelList edgeMap1
+            (
+                matchEdges
+                (
+                    surf1,
+                    subSurf1,
+                    pointMap1
+                )
+            );
+
+
+            for (label zone2I = 0; zone2I < nZones2; zone2I++)
+            {
+                // Generate sub surface for zone2I
+                boolList includeMap2(surf2.size(), false);
+
+                forAll(zone2, faceI)
+                {
+                    if (zone2[faceI] == zone2I)
+                    {
+                        includeMap2[faceI] = true;
+                    }
+                }
+
+                labelList pointMap2;
+                labelList faceMap2;
+                triSurface subSurf2
+                (
+                    surf2.subsetMesh
+                    (
+                        includeMap2,
+                        pointMap2,
+                        faceMap2
+                    )
+                );
+
+
+                const boundBox subBb2
+                (
+                    subSurf2.points(),
+                    subSurf2.meshPoints(),
+                    false
+                );
+
+                // Short-circuit expensive calculations
+                if (!subBb2.overlaps(subBb1))
+                {
+                    continue;
+                }
+
+
+                // Split non-manifold points; update pointMap
+                dupNonManifoldPoints(subSurf2, pointMap2);
+
+                const labelList edgeMap2
+                (
+                    matchEdges
+                    (
+                        surf2,
+                        subSurf2,
+                        pointMap2
+                    )
+                );
+
+
+                // Do cuts
+                edgeIntersections subEdgeCuts1;
+                edgeIntersections subEdgeCuts2;
+                calcEdgeCutsCGAL
+                (
+                    subSurf1,
+                    subSurf2,
+                    perturb,
+                    subEdgeCuts1,
+                    subEdgeCuts2
+                );
+
+                // Move original surface
+                {
+                    pointField points2(surf2.points());
+                    forAll(pointMap2, i)
+                    {
+                        label subMeshPointI = subSurf2.meshPoints()[i];
+                        label meshPointI = surf2.meshPoints()[pointMap2[i]];
+                        points2[meshPointI] = subSurf2.points()[subMeshPointI];
+                    }
+                    surf2.movePoints(points2);
+                }
+
+                // Insert into main structure
+                edgeCuts1.merge(subEdgeCuts1, edgeMap1, faceMap2);
+                edgeCuts2.merge(subEdgeCuts2, edgeMap2, faceMap1);
+            }
+
+
+            // Move original surface
+            {
+                pointField points1(surf1.points());
+                forAll(pointMap1, i)
+                {
+                    label subMeshPointI = subSurf1.meshPoints()[i];
+                    label meshPointI = surf1.meshPoints()[pointMap1[i]];
+                    points1[meshPointI] = subSurf1.points()[subMeshPointI];
+                }
+                surf1.movePoints(points1);
+            }
+        }
+    }
+}
+
+
+#endif // NO_CGAL
+
+
+//void calcFeaturePoints(const pointField& points, const edgeList& edges)
+//{
+//    edgeMesh eMesh(points, edges);
+//
+//    const labelListList& pointEdges = eMesh.pointEdges();
+//
+//
+//    // Get total number of feature points
+//    label nFeaturePoints = 0;
+//    forAll(pointEdges, pI)
+//    {
+//        const labelList& pEdges = pointEdges[pI];
+//
+//        if (pEdges.size() == 1)
+//        {
+//            nFeaturePoints++;
+//        }
+//    }
+//
+//
+//    // Calculate addressing from feature point to cut point and cut edge
+//    labelList featurePointToCutPoint(nFeaturePoints);
+//    labelList featurePointToCutEdge(nFeaturePoints);
+//
+//    label nFeatPts = 0;
+//    forAll(pointEdges, pI)
+//    {
+//        const labelList& pEdges = pointEdges[pI];
+//
+//        if (pEdges.size() == 1)
+//        {
+//            featurePointToCutPoint[nFeatPts] = pI;
+//            featurePointToCutEdge[nFeatPts] = pEdges[0];
+//            nFeatPts++;
+//        }
+//    }
+//
+//
+//
+//    label concaveStart = 0;
+//    label mixedStart = 0;
+//    label nonFeatureStart = nFeaturePoints;
+//
+//
+//    labelListList featurePointNormals(nFeaturePoints);
+//    labelListList featurePointEdges(nFeaturePoints);
+//    labelList regionEdges;
+//}
+
+
+autoPtr<extendedFeatureEdgeMesh> createEdgeMesh
+(
+    const IOobject& io,
+    const booleanSurface::booleanOpType action,
+    const bool surf1Baffle,
+    const bool surf2Baffle,
+    const bool invertedSpace,
+    const triSurface& surf1,
+    const edgeIntersections& edgeCuts1,
+    const triSurface& surf2,
+    const edgeIntersections& edgeCuts2
+)
+{
     // Determine intersection edges from the edge cuts
     surfaceIntersection inter
     (
         surf1,
-        edge1Cuts,
+        edgeCuts1,
         surf2,
-        edge2Cuts
-    );
-
-    fileName sFeatFileName
-    (
-        surf1Name.lessExt().name()
-      + "_"
-      + surf2Name.lessExt().name()
-      + "_"
-      + action
+        edgeCuts2
     );
 
     label nFeatEds = inter.cutEdges().size();
@@ -512,6 +1249,10 @@ int main(int argc, char *argv[])
     List<DynamicList<label> > edgeNormals(nFeatEds);
     List<DynamicList<label> > normalDirections(nFeatEds);
 
+
+    const triSurface& s1 = surf1;
+    const triSurface& s2 = surf2;
+
     forAllConstIter(labelPairLookup, inter.facePairToEdge(), iter)
     {
         const label& cutEdgeI = iter();
@@ -519,8 +1260,8 @@ int main(int argc, char *argv[])
 
         const edge& fE = inter.cutEdges()[cutEdgeI];
 
-        const vector& norm1 = surf1.faceNormals()[facePair.first()];
-        const vector& norm2 = surf2.faceNormals()[facePair.second()];
+        const vector& norm1 = s1.faceNormals()[facePair.first()];
+        const vector& norm2 = s2.faceNormals()[facePair.second()];
 
         DynamicList<label>& eNormals = edgeNormals[cutEdgeI];
         DynamicList<label>& nDirections = normalDirections[cutEdgeI];
@@ -546,7 +1287,7 @@ int main(int argc, char *argv[])
                     norm1,
                     norm2,
                     edgeDirections[cutEdgeI],
-                    surf1[facePair.first()].centre(surf1.points()),
+                    s1[facePair.first()].centre(s1.points()),
                     inter.cutPoints()[fE.start()]
                 )
             );
@@ -572,7 +1313,7 @@ int main(int argc, char *argv[])
                     norm2,
                     norm1,
                     edgeDirections[cutEdgeI],
-                    surf2[facePair.second()].centre(surf2.points()),
+                    s2[facePair.second()].centre(s2.points()),
                     inter.cutPoints()[fE.start()]
                 )
             );
@@ -600,7 +1341,7 @@ int main(int argc, char *argv[])
                         norm2,
                         norm1,
                         edgeDirections[cutEdgeI],
-                        surf2[facePair.second()].centre(surf2.points()),
+                        s2[facePair.second()].centre(s2.points()),
                         inter.cutPoints()[fE.start()]
                     )
                 );
@@ -630,7 +1371,7 @@ int main(int argc, char *argv[])
                         norm1,
                         norm2,
                         edgeDirections[cutEdgeI],
-                        surf1[facePair.first()].centre(surf1.points()),
+                        s1[facePair.first()].centre(s1.points()),
                         inter.cutPoints()[fE.start()]
                     )
                 );
@@ -675,7 +1416,7 @@ int main(int argc, char *argv[])
         }
     }
 
-    if (validActions[action] == booleanSurface::UNION)
+    if (action == booleanSurface::UNION)
     {
         if (!invertedSpace)
         {
@@ -688,7 +1429,7 @@ int main(int argc, char *argv[])
             internalStart = nIntOrExt;
         }
     }
-    else if (validActions[action] == booleanSurface::INTERSECTION)
+    else if (action == booleanSurface::INTERSECTION)
     {
         if (!invertedSpace)
         {
@@ -701,14 +1442,14 @@ int main(int argc, char *argv[])
             internalStart = 0;
         }
     }
-    else if (validActions[action] == booleanSurface::DIFFERENCE)
+    else if (action == booleanSurface::DIFFERENCE)
     {
         // All edges are external
         internalStart = nIntOrExt;
     }
     else
     {
-        FatalErrorIn(args.executable())
+        FatalErrorIn("createEdgeMesh(..)")
             << "Unsupported booleanSurface:booleanOpType and space "
             << action << " " << invertedSpace
             << abort(FatalError);
@@ -734,42 +1475,257 @@ int main(int argc, char *argv[])
         normalDirectionsTmp[i] = normalDirections[i];
     }
 
-    calcFeaturePoints(inter.cutPoints(), inter.cutEdges());
+    //calcFeaturePoints(inter.cutPoints(), inter.cutEdges());
 
-    extendedFeatureEdgeMesh feMesh
+    return autoPtr<extendedFeatureEdgeMesh>
+    (
+        new extendedFeatureEdgeMesh
+        (
+            io,
+            inter.cutPoints(),
+            inter.cutEdges(),
+
+            0,                  // concaveStart,
+            0,                  // mixedStart,
+            0,                  // nonFeatureStart,
+
+            internalStart,      // internalStart,
+            nIntOrExt,           // flatStart,
+            nIntOrExt + nFlat,   // openStart,
+            nIntOrExt + nFlat + nOpen,   // multipleStart,
+
+            normalsTmp,
+            normalVolumeTypesTmp,
+            edgeDirections,
+            normalDirectionsTmp,
+            edgeNormalsTmp,
+
+            labelListList(0),   // featurePointNormals,
+            labelListList(0),   // featurePointEdges,
+            labelList(0)        // regionEdges
+        )
+    );
+}
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+    argList::validArgs.append("action");
+    argList::validArgs.append("surface file");
+    argList::validArgs.append("surface file");
+
+    argList::addBoolOption
+    (
+        "surf1Baffle",
+        "Mark surface 1 as a baffle"
+    );
+
+    argList::addBoolOption
+    (
+        "surf2Baffle",
+        "Mark surface 2 as a baffle"
+    );
+
+    argList::addBoolOption
+    (
+        "perturb",
+        "Perturb surface points to escape degenerate intersections"
+    );
+
+    argList::addBoolOption
+    (
+        "invertedSpace",
+        "do the surfaces have inverted space orientation, "
+        "i.e. a point at infinity is considered inside. "
+        "This is only sensible for union and intersection."
+    );
+
+    argList::addOption
+    (
+        "trim",
+        "((surface1 volumeType) .. (surfaceN volumeType))",
+        "Trim resulting intersection with additional surfaces;"
+        " volumeType is 'inside' (keep (parts of) edges that are inside)"
+        ", 'outside' (keep (parts of) edges that are outside) or"
+        " 'mixed' (keep all)"
+    );
+
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+    word action(args.args()[1]);
+
+    HashTable<booleanSurface::booleanOpType> validActions;
+    validActions.insert("intersection", booleanSurface::INTERSECTION);
+    validActions.insert("union", booleanSurface::UNION);
+    validActions.insert("difference", booleanSurface::DIFFERENCE);
+
+    if (!validActions.found(action))
+    {
+        FatalErrorIn(args.executable())
+            << "Unsupported action " << action << endl
+            << "Supported actions:" << validActions.toc() << abort(FatalError);
+    }
+
+
+    List<Pair<word> > surfaceAndSide;
+    if (args.optionReadIfPresent("trim", surfaceAndSide))
+    {
+        Info<< "Trimming edges with " << surfaceAndSide << endl;
+    }
+
+
+    const word surf1Name(args.args()[2]);
+    Info<< "Reading surface " << surf1Name << endl;
+    triSurfaceMesh surf1
     (
         IOobject
         (
-            sFeatFileName + ".extendedFeatureEdgeMesh",
+            surf1Name,
             runTime.constant(),
-            "extendedFeatureEdgeMesh",
-            runTime,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        inter.cutPoints(),
-        inter.cutEdges(),
-
-        0,                  // concaveStart,
-        0,                  // mixedStart,
-        0,                  // nonFeatureStart,
-
-        internalStart,      // internalStart,
-        nIntOrExt,           // flatStart,
-        nIntOrExt + nFlat,   // openStart,
-        nIntOrExt + nFlat + nOpen,   // multipleStart,
-
-        normalsTmp,
-        normalVolumeTypesTmp,
-        edgeDirections,
-        normalDirectionsTmp,
-        edgeNormalsTmp,
-
-        labelListList(0),   // featurePointNormals,
-        labelListList(0),   // featurePointEdges,
-        labelList(0)        // regionEdges
+            triSurfaceMesh::meshSubDir,
+            runTime
+        )
     );
 
+    Info<< surf1Name << " statistics:" << endl;
+    surf1.writeStats(Info);
+    Info<< endl;
+
+    const word surf2Name(args[3]);
+    Info<< "Reading surface " << surf2Name << endl;
+    triSurfaceMesh surf2
+    (
+        IOobject
+        (
+            surf2Name,
+            runTime.constant(),
+            triSurfaceMesh::meshSubDir,
+            runTime
+        )
+    );
+
+    Info<< surf2Name << " statistics:" << endl;
+    surf2.writeStats(Info);
+    Info<< endl;
+
+    const bool surf1Baffle = args.optionFound("surf1Baffle");
+    const bool surf2Baffle = args.optionFound("surf2Baffle");
+
+    edgeIntersections edgeCuts1;
+    edgeIntersections edgeCuts2;
+
+    bool invertedSpace = args.optionFound("invertedSpace");
+
+    if (invertedSpace && validActions[action] == booleanSurface::DIFFERENCE)
+    {
+        FatalErrorIn(args.executable())
+            << "Inverted space only makes sense for union or intersection."
+            << exit(FatalError);
+    }
+
+
+#ifdef NO_CGAL
+    // Calculate the points where the edges are cut by the other surface
+    calcEdgeCuts
+    (
+        surf1,
+        surf2,
+        args.optionFound("perturb"),
+        edgeCuts1,
+        edgeCuts2
+    );
+#else
+    //calcEdgeCutsCGAL
+    calcEdgeCutsBitsCGAL
+    (
+        surf1,
+        surf2,
+        args.optionFound("perturb"),
+        edgeCuts1,
+        edgeCuts2
+    );
+#endif // NO_CGAL
+
+
+    const fileName sFeatFileName
+    (
+        fileName(surf1Name).lessExt().name()
+      + "_"
+      + fileName(surf2Name).lessExt().name()
+      + "_"
+      + action
+    );
+
+    autoPtr<extendedFeatureEdgeMesh> feMeshPtr
+    (
+        createEdgeMesh
+        (
+            IOobject
+            (
+                sFeatFileName + ".extendedFeatureEdgeMesh",
+                runTime.constant(),
+                "extendedFeatureEdgeMesh",
+                runTime,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            booleanSurface::booleanOpTypeNames[action],
+            surf1Baffle,
+            surf2Baffle,
+            invertedSpace,
+            surf1,
+            edgeCuts1,
+            surf2,
+            edgeCuts2
+        )
+    );
+
+
+    // Trim intersections
+    forAll(surfaceAndSide, trimI)
+    {
+        const word& trimName = surfaceAndSide[trimI].first();
+        const volumeType trimType
+        (
+            volumeType::names[surfaceAndSide[trimI].second()]
+        );
+
+        Info<< "Reading trim surface " << trimName << endl;
+        const triSurfaceMesh trimSurf
+        (
+            IOobject
+            (
+                trimName,
+                runTime.constant(),
+                triSurfaceMesh::meshSubDir,
+                runTime,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+        );
+
+        Info<< trimName << " statistics:" << endl;
+        trimSurf.writeStats(Info);
+        Info<< endl;
+
+        labelList pointMap;
+        labelList edgeMap;
+        feMeshPtr().trim
+        (
+            trimSurf,
+            trimType,
+            pointMap,
+            edgeMap
+        );
+    }
+
+
+    const extendedFeatureEdgeMesh& feMesh = feMeshPtr();
+
+    feMesh.writeStats(Info);
+
     feMesh.write();
 
     feMesh.writeObj(feMesh.path()/sFeatFileName);
@@ -782,7 +1738,7 @@ int main(int argc, char *argv[])
             (
                 sFeatFileName + ".eMesh",   // name
                 runTime.constant(),                         // instance
-                "triSurface",
+                triSurfaceMesh::meshSubDir,
                 runTime,                                    // registry
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
diff --git a/applications/utilities/surface/surfaceCheck/surfaceCheck.C b/applications/utilities/surface/surfaceCheck/surfaceCheck.C
index 6ec067aa84..c292ebe6b5 100644
--- a/applications/utilities/surface/surfaceCheck/surfaceCheck.C
+++ b/applications/utilities/surface/surfaceCheck/surfaceCheck.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -170,6 +170,101 @@ labelList countBins
 }
 
 
+
+void writeZoning
+(
+    const triSurface& surf,
+    const labelList& faceZone,
+    const word& fieldName,
+    const fileName& surfFilePath,
+    const fileName& surfFileNameBase
+)
+{
+    Info<< "Writing zoning to "
+        <<  fileName
+            (
+                surfFilePath
+              / fieldName
+              + '_'
+              + surfFileNameBase
+              + '.'
+              + vtkSurfaceWriter::typeName
+            )
+        << "..." << endl << endl;
+
+    // Convert data
+    scalarField scalarFaceZone(faceZone.size());
+    forAll(faceZone, i)
+    {
+        scalarFaceZone[i] = faceZone[i];
+    }
+    faceList faces(surf.size());
+    forAll(surf, i)
+    {
+        faces[i] = surf[i].triFaceFace();
+    }
+
+    vtkSurfaceWriter().write
+    (
+        surfFilePath,
+        surfFileNameBase,
+        surf.points(),
+        faces,
+        fieldName,
+        scalarFaceZone,
+        false               // face based data
+    );
+}
+
+
+void writeParts
+(
+    const triSurface& surf,
+    const label nFaceZones,
+    const labelList& faceZone,
+    const fileName& surfFilePath,
+    const fileName& surfFileNameBase
+)
+{
+    for (label zone = 0; zone < nFaceZones; zone++)
+    {
+        boolList includeMap(surf.size(), false);
+
+        forAll(faceZone, faceI)
+        {
+            if (faceZone[faceI] == zone)
+            {
+                includeMap[faceI] = true;
+            }
+        }
+
+        labelList pointMap;
+        labelList faceMap;
+
+        triSurface subSurf
+        (
+            surf.subsetMesh
+            (
+                includeMap,
+                pointMap,
+                faceMap
+            )
+        );
+
+        fileName subName
+        (
+            surfFilePath
+           /surfFileNameBase + "_" + name(zone) + ".obj"
+        );
+
+        Info<< "writing part " << zone << " size " << subSurf.size()
+            << " to " << subName << endl;
+
+        subSurf.write(subName);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -218,6 +313,20 @@ int main(int argc, char *argv[])
     surf.writeStats(Info);
     Info<< endl;
 
+
+    // Determine path and extension
+    fileName surfFileNameBase(surfFileName.name());
+    const word fileType = surfFileNameBase.ext();
+    // Strip extension
+    surfFileNameBase = surfFileNameBase.lessExt();
+    // If extension was .gz strip original extension
+    if (fileType == "gz")
+    {
+        surfFileNameBase = surfFileNameBase.lessExt();
+    }
+    const fileName surfFilePath(surfFileName.path());
+
+
     // write bounding box corners
     if (args.optionFound("blockMesh"))
     {
@@ -515,12 +624,12 @@ int main(int argc, char *argv[])
 
     DynamicList<label> problemFaces(surf.size()/100 + 1);
 
-    const labelListList& eFaces = surf.edgeFaces();
+    const labelListList& edgeFaces = surf.edgeFaces();
 
     label nSingleEdges = 0;
-    forAll(eFaces, edgeI)
+    forAll(edgeFaces, edgeI)
     {
-        const labelList& myFaces = eFaces[edgeI];
+        const labelList& myFaces = edgeFaces[edgeI];
 
         if (myFaces.size() == 1)
         {
@@ -531,9 +640,9 @@ int main(int argc, char *argv[])
     }
 
     label nMultEdges = 0;
-    forAll(eFaces, edgeI)
+    forAll(edgeFaces, edgeI)
     {
-        const labelList& myFaces = eFaces[edgeI];
+        const labelList& myFaces = edgeFaces[edgeI];
 
         if (myFaces.size() > 2)
         {
@@ -549,7 +658,8 @@ int main(int argc, char *argv[])
 
     if ((nSingleEdges != 0) || (nMultEdges != 0))
     {
-        Info<< "Surface is not closed since not all edges connected to "
+        Info<< "Surface is not closed since not all edges ("
+            << edgeFaces.size() << ") connected to "
             << "two faces:" << endl
             << "    connected to one face : " << nSingleEdges << endl
             << "    connected to >2 faces : " << nMultEdges << endl;
@@ -578,10 +688,9 @@ int main(int argc, char *argv[])
         boolList borderEdge(surf.nEdges(), false);
         if (splitNonManifold)
         {
-            const labelListList& eFaces = surf.edgeFaces();
-            forAll(eFaces, edgeI)
+            forAll(edgeFaces, edgeI)
             {
-                if (eFaces[edgeI].size() > 2)
+                if (edgeFaces[edgeI].size() > 2)
                 {
                     borderEdge[edgeI] = true;
                 }
@@ -597,85 +706,15 @@ int main(int argc, char *argv[])
         {
             Info<< "Splitting surface into parts ..." << endl << endl;
 
-            fileName surfFileNameBase(surfFileName.name());
-            const word fileType = surfFileNameBase.ext();
-            // Strip extension
-            surfFileNameBase = surfFileNameBase.lessExt();
-            // If extension was .gz strip original extension
-            if (fileType == "gz")
-            {
-                surfFileNameBase = surfFileNameBase.lessExt();
-            }
-
-
-            {
-                Info<< "Writing zoning to "
-                    <<  fileName
-                        (
-                            "zone_"
-                          + surfFileNameBase
-                          + '.'
-                          + vtkSurfaceWriter::typeName
-                        )
-                    << "..." << endl << endl;
-
-                // Convert data
-                scalarField scalarFaceZone(faceZone.size());
-                forAll(faceZone, i)
-                {
-                    scalarFaceZone[i] = faceZone[i];
-                }
-                faceList faces(surf.size());
-                forAll(surf, i)
-                {
-                    faces[i] = surf[i].triFaceFace();
-                }
-
-                vtkSurfaceWriter().write
-                (
-                    surfFileName.path(),
-                    surfFileNameBase,
-                    surf.points(),
-                    faces,
-                    "zone",
-                    scalarFaceZone,
-                    false               // face based data
-                );
-            }
-
-
-            for (label zone = 0; zone < numZones; zone++)
-            {
-                boolList includeMap(surf.size(), false);
-
-                forAll(faceZone, faceI)
-                {
-                    if (faceZone[faceI] == zone)
-                    {
-                        includeMap[faceI] = true;
-                    }
-                }
-
-                labelList pointMap;
-                labelList faceMap;
-
-                triSurface subSurf
-                (
-                    surf.subsetMesh
-                    (
-                        includeMap,
-                        pointMap,
-                        faceMap
-                    )
-                );
-
-                fileName subName(surfFileNameBase + "_" + name(zone) + ".obj");
-
-                Info<< "writing part " << zone << " size " << subSurf.size()
-                    << " to " << subName << endl;
-
-                subSurf.write(subName);
-            }
+            writeZoning(surf, faceZone, "zone", surfFilePath, surfFileNameBase);
+            writeParts
+            (
+                surf,
+                numZones,
+                faceZone,
+                surfFilePath,
+                surfFileNameBase
+            );
         }
     }
 
@@ -700,6 +739,15 @@ int main(int argc, char *argv[])
     if (numNormalZones > 1)
     {
         Info<< "More than one normal orientation." << endl;
+        writeZoning(surf, normalZone, "normal", surfFilePath, surfFileNameBase);
+        writeParts
+        (
+            surf,
+            numNormalZones,
+            normalZone,
+            surfFilePath,
+            surfFileNameBase + "_normal"
+        );
     }
     Info<< endl;
 
diff --git a/applications/utilities/surface/surfaceFeatureConvert/surfaceFeatureConvert.C b/applications/utilities/surface/surfaceFeatureConvert/surfaceFeatureConvert.C
index fadf9a7ec7..e62bee9a68 100644
--- a/applications/utilities/surface/surfaceFeatureConvert/surfaceFeatureConvert.C
+++ b/applications/utilities/surface/surfaceFeatureConvert/surfaceFeatureConvert.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -99,7 +99,6 @@ int main(int argc, char *argv[])
 
     mesh.write(exportName);
     mesh.writeStats(Info);
-    Info<< endl;
 
     Info<< "\nEnd\n" << endl;
 
diff --git a/applications/utilities/surface/surfaceInertia/surfaceInertia.C b/applications/utilities/surface/surfaceInertia/surfaceInertia.C
index 6b30f610b1..14ff3cc5c8 100644
--- a/applications/utilities/surface/surfaceInertia/surfaceInertia.C
+++ b/applications/utilities/surface/surfaceInertia/surfaceInertia.C
@@ -3,7 +3,7 @@
  \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
   \\    /    O peration      |
    \\  /     A nd            | Copyright (C) 2011-2013 OpenFOAM Foundation
-    \\/      M anipulation   |
+    \\/      M anipulation   | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -406,7 +406,7 @@ int main(int argc, char *argv[])
         str << "l " << 1 << ' ' << i + 1 << endl;
     }
 
-    Info<< nl << "End" << nl << endl;
+    Info<< "\nEnd\n" << endl;
 
     return 0;
 }
diff --git a/applications/utilities/surface/surfaceInflate/Make/files b/applications/utilities/surface/surfaceInflate/Make/files
new file mode 100644
index 0000000000..8677a01c39
--- /dev/null
+++ b/applications/utilities/surface/surfaceInflate/Make/files
@@ -0,0 +1,3 @@
+surfaceInflate.C
+
+EXE = $(FOAM_APPBIN)/surfaceInflate
diff --git a/applications/utilities/surface/surfaceInflate/Make/options b/applications/utilities/surface/surfaceInflate/Make/options
new file mode 100644
index 0000000000..6c6cd12287
--- /dev/null
+++ b/applications/utilities/surface/surfaceInflate/Make/options
@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(FOAM_SRC)/triSurface/lnInclude \
+    -I$(FOAM_SRC)/surfMesh/lnInclude \
+    -I$(FOAM_SRC)/meshTools/lnInclude
+
+
+EXE_LIBS = \
+    -ltriSurface \
+    -lsurfMesh \
+    -lmeshTools
diff --git a/applications/utilities/surface/surfaceInflate/surfaceInflate.C b/applications/utilities/surface/surfaceInflate/surfaceInflate.C
new file mode 100644
index 0000000000..dfc4154321
--- /dev/null
+++ b/applications/utilities/surface/surfaceInflate/surfaceInflate.C
@@ -0,0 +1,879 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Application
+    surfaceInflate
+
+Description
+    Inflates surface. WIP. Checks for overlaps and locally lowers
+    inflation distance
+
+Usage
+    - surfaceInflate [OPTION]
+
+    \param -checkSelfIntersection \n
+    Includes checks for self-intersection.
+
+    \param -nSmooth
+    Specify number of smoothing iterations
+
+    \param -featureAngle
+    Specify a feature angle
+
+
+    E.g. inflate surface by 2cm with 1.5 safety factor:
+        surfaceInflate DTC-scaled.obj 0.02 1.5 -featureAngle 45 -nSmooth 2
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "triSurfaceFields.H"
+#include "triSurfaceMesh.H"
+#include "triSurfaceGeoMesh.H"
+#include "PackedBoolList.H"
+#include "OBJstream.H"
+#include "surfaceFeatures.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scalar calcVertexNormalWeight
+(
+    const triFace& f,
+    const label pI,
+    const vector& fN,
+    const pointField& points
+)
+{
+    label index = findIndex(f, pI);
+
+    if (index == -1)
+    {
+        FatalErrorIn("calcVertexNormals()")
+            << "Point not in face" << abort(FatalError);
+    }
+
+    const vector e1 = points[f[index]] - points[f[f.fcIndex(index)]];
+    const vector e2 = points[f[index]] - points[f[f.rcIndex(index)]];
+
+    return mag(fN)/(magSqr(e1)*magSqr(e2) + VSMALL);
+}
+
+
+tmp<vectorField> calcVertexNormals(const triSurface& surf)
+{
+    // Weighted average of normals of faces attached to the vertex
+    // Weight = fA / (mag(e0)^2 * mag(e1)^2);
+    tmp<vectorField> tpointNormals
+    (
+        new pointField(surf.nPoints(), vector::zero)
+    );
+    vectorField& pointNormals = tpointNormals();
+
+    const pointField& points = surf.points();
+    const labelListList& pointFaces = surf.pointFaces();
+    const labelList& meshPoints = surf.meshPoints();
+
+    forAll(pointFaces, pI)
+    {
+        const labelList& pFaces = pointFaces[pI];
+
+        forAll(pFaces, fI)
+        {
+            const label faceI = pFaces[fI];
+            const triFace& f = surf[faceI];
+
+            vector fN = f.normal(points);
+
+            scalar weight = calcVertexNormalWeight
+            (
+                f,
+                meshPoints[pI],
+                fN,
+                points
+            );
+
+            pointNormals[pI] += weight*fN;
+        }
+
+        pointNormals[pI] /= mag(pointNormals[pI]) + VSMALL;
+    }
+
+    return tpointNormals;
+}
+
+
+tmp<vectorField> calcPointNormals
+(
+    const triSurface& s,
+    const PackedBoolList& isFeaturePoint,
+    const List<surfaceFeatures::edgeStatus>& edgeStat
+)
+{
+    //const pointField pointNormals(s.pointNormals());
+    tmp<vectorField> tpointNormals(calcVertexNormals(s));
+    vectorField& pointNormals = tpointNormals();
+
+
+    // feature edges: create edge normals from edgeFaces only.
+    {
+        const labelListList& edgeFaces = s.edgeFaces();
+
+        labelList nNormals(s.nPoints(), 0);
+        forAll(edgeStat, edgeI)
+        {
+            if (edgeStat[edgeI] != surfaceFeatures::NONE)
+            {
+                const edge& e = s.edges()[edgeI];
+                forAll(e, i)
+                {
+                    if (!isFeaturePoint[e[i]])
+                    {
+                        pointNormals[e[i]] = vector::zero;
+                    }
+                }
+            }
+        }
+
+        forAll(edgeStat, edgeI)
+        {
+            if (edgeStat[edgeI] != surfaceFeatures::NONE)
+            {
+                const labelList& eFaces = edgeFaces[edgeI];
+
+                // Get average edge normal
+                vector n = vector::zero;
+                forAll(eFaces, i)
+                {
+                    n += s.faceNormals()[eFaces[i]];
+                }
+                n /= eFaces.size();
+
+
+                // Sum to points
+                const edge& e = s.edges()[edgeI];
+                forAll(e, i)
+                {
+                    if (!isFeaturePoint[e[i]])
+                    {
+                        pointNormals[e[i]] += n;
+                        nNormals[e[i]]++;
+                    }
+                }
+            }
+        }
+
+        forAll(nNormals, pointI)
+        {
+            if (nNormals[pointI] > 0)
+            {
+                pointNormals[pointI] /= mag(pointNormals[pointI]);
+            }
+        }
+    }
+
+
+    forAll(pointNormals, pointI)
+    {
+        if (mag(mag(pointNormals[pointI])-1) > SMALL)
+        {
+            FatalErrorIn("calcPointNormals()") << "unitialised"
+                << exit(FatalError);
+        }
+    }
+
+    return tpointNormals;
+}
+
+
+void detectSelfIntersections
+(
+    const triSurfaceMesh& s,
+    PackedBoolList& isEdgeIntersecting
+)
+{
+    const edgeList& edges = s.edges();
+    const indexedOctree<treeDataTriSurface>& tree = s.tree();
+    const labelList& meshPoints = s.meshPoints();
+    const pointField& points = s.points();
+
+    isEdgeIntersecting.setSize(edges.size());
+    isEdgeIntersecting = false;
+
+    forAll(edges, edgeI)
+    {
+        const edge& e = edges[edgeI];
+
+        pointIndexHit hitInfo
+        (
+            tree.findLine
+            (
+                points[meshPoints[e[0]]],
+                points[meshPoints[e[1]]],
+                treeDataTriSurface::findSelfIntersectOp(tree, edgeI)
+            )
+        );
+
+        if (hitInfo.hit())
+        {
+            isEdgeIntersecting[edgeI] = true;
+        }
+    }
+}
+
+
+label detectIntersectionPoints
+(
+    const scalar tolerance,
+    const triSurfaceMesh& s,
+
+    const scalar extendFactor,
+    const pointField& initialPoints,
+    const vectorField& displacement,
+
+    const bool checkSelfIntersect,
+    const PackedBoolList& initialIsEdgeIntersecting,
+
+    PackedBoolList& isPointOnHitEdge,
+    scalarField& scale
+)
+{
+    const pointField initialLocalPoints(initialPoints, s.meshPoints());
+    const List<labelledTri>& localFaces = s.localFaces();
+
+
+    label nHits = 0;
+    isPointOnHitEdge.setSize(s.nPoints());
+    isPointOnHitEdge = false;
+
+
+    // 1. Extrusion offset vectors intersecting new surface location
+    {
+        scalar tol = max(tolerance, 10*s.tolerance());
+
+        // Collect all the edge vectors. Slightly shorten the edges to prevent
+        // finding lots of intersections. The fast triangle intersection routine
+        // has problems with rays passing through a point of the
+        // triangle.
+        pointField start(initialLocalPoints+tol*displacement);
+        pointField end(initialLocalPoints+extendFactor*displacement);
+
+        List<pointIndexHit> hits;
+        s.findLineAny(start, end, hits);
+
+        forAll(hits, pointI)
+        {
+            if
+            (
+                hits[pointI].hit()
+            &&  findIndex(localFaces[hits[pointI].index()], pointI) == -1
+            )
+            {
+                scale[pointI] = max(0.0, scale[pointI]-0.2);
+
+                isPointOnHitEdge[pointI] = true;
+                nHits++;
+            }
+        }
+    }
+
+
+    // 2. (new) surface self intersections
+    if (checkSelfIntersect)
+    {
+        PackedBoolList isEdgeIntersecting;
+        detectSelfIntersections(s, isEdgeIntersecting);
+
+        const edgeList& edges = s.edges();
+        const pointField& points = s.points();
+
+        forAll(edges, edgeI)
+        {
+            const edge& e = edges[edgeI];
+
+            if (isEdgeIntersecting[edgeI] && !initialIsEdgeIntersecting[edgeI])
+            {
+                if (isPointOnHitEdge.set(e[0]))
+                {
+                    label start = s.meshPoints()[e[0]];
+                    const point& pt = points[start];
+
+                    Pout<< "Additional self intersection at "
+                        << pt
+                        << endl;
+
+                    scale[e[0]] = max(0.0, scale[e[0]]-0.2);
+                    nHits++;
+                }
+                if (isPointOnHitEdge.set(e[1]))
+                {
+                    label end = s.meshPoints()[e[1]];
+                    const point& pt = points[end];
+
+                    Pout<< "Additional self intersection at "
+                        << pt
+                        << endl;
+
+                    scale[e[1]] = max(0.0, scale[e[1]]-0.2);
+                    nHits++;
+                }
+            }
+        }
+    }
+
+    return nHits;
+}
+
+
+tmp<scalarField> avg
+(
+    const triSurface& s,
+    const scalarField& fld,
+    const scalarField& edgeWeights
+)
+{
+    tmp<scalarField> tres(new scalarField(s.nPoints(), 0.0));
+    scalarField& res = tres();
+
+    scalarField sumWeight(s.nPoints(), 0.0);
+
+    const edgeList& edges = s.edges();
+
+    forAll(edges, edgeI)
+    {
+        const edge& e = edges[edgeI];
+        const scalar w = edgeWeights[edgeI];
+
+        res[e[0]] += w*fld[e[1]];
+        sumWeight[e[0]] += w;
+
+        res[e[1]] += w*fld[e[0]];
+        sumWeight[e[1]] += w;
+    }
+
+    // Average
+    // ~~~~~~~
+
+    forAll(res, pointI)
+    {
+        if (mag(sumWeight[pointI]) < VSMALL)
+        {
+            // Unconnected point. Take over original value
+            res[pointI] = fld[pointI];
+        }
+        else
+        {
+            res[pointI] /= sumWeight[pointI];
+        }
+    }
+
+    return tres;
+}
+
+
+void minSmooth
+(
+    const triSurface& s,
+    const PackedBoolList& isAffectedPoint,
+    const scalarField& fld,
+    scalarField& newFld
+)
+{
+
+    const edgeList& edges = s.edges();
+    const pointField& points = s.points();
+    const labelList& mp = s.meshPoints();
+
+    scalarField edgeWeights(edges.size());
+    forAll(edges, edgeI)
+    {
+        const edge& e = edges[edgeI];
+        scalar w = mag(points[mp[e[0]]]-points[mp[e[1]]]);
+
+        edgeWeights[edgeI] = 1.0/(max(w, SMALL));
+    }
+
+    tmp<scalarField> tavgFld = avg(s, fld, edgeWeights);
+
+    const scalarField& avgFld = tavgFld();
+
+    forAll(fld, pointI)
+    {
+        if (isAffectedPoint.get(pointI) == 1)
+        {
+            newFld[pointI] = min
+            (
+                fld[pointI],
+                0.5*fld[pointI] + 0.5*avgFld[pointI]
+                //avgFld[pointI]
+            );
+        }
+    }
+}
+
+
+void lloydsSmoothing
+(
+    const label nSmooth,
+    triSurface& s,
+    const PackedBoolList& isFeaturePoint,
+    const List<surfaceFeatures::edgeStatus>& edgeStat,
+    const PackedBoolList& isAffectedPoint
+)
+{
+    const labelList& meshPoints = s.meshPoints();
+    const edgeList& edges = s.edges();
+
+
+    PackedBoolList isSmoothPoint(isAffectedPoint);
+    // Extend isSmoothPoint
+    {
+        PackedBoolList newIsSmoothPoint(isSmoothPoint);
+        forAll(edges, edgeI)
+        {
+            const edge& e = edges[edgeI];
+            if (isSmoothPoint[e[0]])
+            {
+                newIsSmoothPoint[e[1]] = true;
+            }
+            if (isSmoothPoint[e[1]])
+            {
+                newIsSmoothPoint[e[0]] = true;
+            }
+        }
+        Info<< "From points-to-smooth " << isSmoothPoint.count()
+            << " to " << newIsSmoothPoint.count()
+            << endl;
+        isSmoothPoint.transfer(newIsSmoothPoint);
+    }
+
+    // Do some smoothing (Lloyds algorithm) around problematic points
+    for (label i = 0; i < nSmooth; i++)
+    {
+        const labelListList& pointFaces = s.pointFaces();
+        const pointField& faceCentres = s.faceCentres();
+
+        pointField newPoints(s.points());
+        forAll(isSmoothPoint, pointI)
+        {
+            if (isSmoothPoint[pointI] && !isFeaturePoint[pointI])
+            {
+                const labelList& pFaces = pointFaces[pointI];
+
+                point avg = point::zero;
+                forAll(pFaces, pFaceI)
+                {
+                    avg += faceCentres[pFaces[pFaceI]];
+                }
+                avg /= pFaces.size();
+                newPoints[meshPoints[pointI]] = avg;
+            }
+        }
+
+
+        // Move points on feature edges only according to feature edges.
+
+        const pointField& points = s.points();
+
+        vectorField pointSum(s.nPoints(), vector::zero);
+        labelList nPointSum(s.nPoints(), 0);
+
+        forAll(edges, edgeI)
+        {
+            if (edgeStat[edgeI] != surfaceFeatures::NONE)
+            {
+                const edge& e = edges[edgeI];
+                const point& start = points[meshPoints[e[0]]];
+                const point& end = points[meshPoints[e[1]]];
+                point eMid = 0.5*(start+end);
+                pointSum[e[0]] += eMid;
+                nPointSum[e[0]]++;
+                pointSum[e[1]] += eMid;
+                nPointSum[e[1]]++;
+            }
+        }
+
+        forAll(pointSum, pointI)
+        {
+            if
+            (
+                isSmoothPoint[pointI]
+             && isFeaturePoint[pointI]
+             && nPointSum[pointI] >= 2
+            )
+            {
+                newPoints[meshPoints[pointI]] =
+                    pointSum[pointI]/nPointSum[pointI];
+            }
+        }
+
+
+        s.movePoints(newPoints);
+
+
+        // Extend isSmoothPoint
+        {
+            PackedBoolList newIsSmoothPoint(isSmoothPoint);
+            forAll(edges, edgeI)
+            {
+                const edge& e = edges[edgeI];
+                if (isSmoothPoint[e[0]])
+                {
+                    newIsSmoothPoint[e[1]] = true;
+                }
+                if (isSmoothPoint[e[1]])
+                {
+                    newIsSmoothPoint[e[0]] = true;
+                }
+            }
+            Info<< "From points-to-smooth " << isSmoothPoint.count()
+                << " to " << newIsSmoothPoint.count()
+                << endl;
+            isSmoothPoint.transfer(newIsSmoothPoint);
+        }
+    }
+}
+
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::addNote("Inflates surface according to point normals.");
+
+    argList::noParallel();
+    argList::addNote
+    (
+        "Creates inflated version of surface using point normals."
+        " Takes surface, distance to inflate and additional safety factor"
+    );
+    argList::addBoolOption
+    (
+        "checkSelfIntersection",
+        "also check for self-intersection"
+    );
+    argList::addOption
+    (
+        "nSmooth",
+        "integer",
+        "number of smoothing iterations (default 20)"
+    );
+    argList::addOption
+    (
+        "featureAngle",
+        "scalar",
+        "feature angle"
+    );
+    argList::addBoolOption
+    (
+        "debug",
+        "switch on additional debug information"
+    );
+
+    argList::validArgs.append("inputFile");
+    argList::validArgs.append("distance");
+    argList::validArgs.append("safety factor [1..]");
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    runTime.functionObjects().off();
+
+    const word inputName(args.args()[1]);
+    const scalar distance(args.argRead<scalar>(2));
+    const scalar extendFactor(args.argRead<scalar>(3));
+    const bool checkSelfIntersect = args.optionFound("checkSelfIntersection");
+    const label nSmooth = args.optionLookupOrDefault("nSmooth", 10);
+    const scalar featureAngle = args.optionLookupOrDefault<scalar>
+    (
+        "featureAngle",
+        180
+    );
+    const bool debug = args.optionFound("debug");
+
+
+    Info<< "Inflating surface " << inputName
+        << " according to point normals" << nl
+        << "    distance               : " << distance << nl
+        << "    safety factor          : " << extendFactor << nl
+        << "    self intersection test : " << checkSelfIntersect << nl
+        << "    smoothing iterations   : " << nSmooth << nl
+        << "    feature angle          : " << featureAngle << nl
+        << endl;
+
+
+    if (extendFactor < 1 || extendFactor > 10)
+    {
+        FatalErrorIn(args.executable())
+            << "Illegal safety factor " << extendFactor
+            << ". It is usually 1..2"
+            << exit(FatalError);
+    }
+
+
+
+    // Load triSurfaceMesh
+    triSurfaceMesh s
+    (
+        IOobject
+        (
+            inputName,                          // name
+            runTime.constant(),                 // instance
+            "triSurface",                       // local
+            runTime,                            // registry
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        )
+    );
+
+    // Mark features
+    const surfaceFeatures features(s, 180.0-featureAngle);
+
+    Info<< "Detected features:" << nl
+        << "    feature points : " << features.featurePoints().size()
+        << " out of " << s.nPoints() << nl
+        << "    feature edges : " << features.featureEdges().size()
+        << " out of " << s.nEdges() << nl
+        << endl;
+
+    PackedBoolList isFeaturePoint(s.nPoints());
+    forAll(features.featurePoints(), i)
+    {
+        label pointI = features.featurePoints()[i];
+        isFeaturePoint[pointI] = true;
+    }
+    const List<surfaceFeatures::edgeStatus> edgeStat(features.toStatus());
+
+
+
+
+    const pointField initialPoints(s.points());
+
+
+    // Construct scale
+    Info<< "Constructing field scale\n" << endl;
+    triSurfacePointScalarField scale
+    (
+        IOobject
+        (
+            "scale",                            // name
+            runTime.timeName(),                 // instance
+            s,                                  // registry
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        s,
+        dimensionedScalar("scale", dimLength, 1.0)
+    );
+
+
+    // Construct unit normals
+
+    Info<< "Calculating vertex normals\n" << endl;
+    const pointField pointNormals
+    (
+        calcPointNormals
+        (
+            s,
+            isFeaturePoint,
+            edgeStat
+        )
+    );
+
+
+    // Construct pointDisplacement
+    Info<< "Constructing field pointDisplacement\n" << endl;
+    triSurfacePointVectorField pointDisplacement
+    (
+        IOobject
+        (
+            "pointDisplacement",                // name
+            runTime.timeName(),                 // instance
+            s,                                  // registry
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        s,
+        dimLength,
+        distance*scale*pointNormals
+    );
+
+
+    const labelList& meshPoints = s.meshPoints();
+
+
+    // Any point on any intersected edge in any of the iterations
+    PackedBoolList isScaledPoint(s.nPoints());
+
+
+    // Detect any self intersections on initial mesh
+    PackedBoolList initialIsEdgeIntersecting;
+    if (checkSelfIntersect)
+    {
+        detectSelfIntersections(s, initialIsEdgeIntersecting);
+    }
+    else
+    {
+        // Mark all edges as already self intersecting so avoid detecting any
+        // new ones
+        initialIsEdgeIntersecting.setSize(s.nEdges(), true);
+    }
+
+
+    // Inflate
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // Move to new location
+        pointField newPoints(initialPoints);
+        forAll(meshPoints, i)
+        {
+            newPoints[meshPoints[i]] += pointDisplacement[i];
+        }
+        s.movePoints(newPoints);
+        Info<< "Bounding box : " << s.bounds() << endl;
+
+
+        // Work out scale from pointDisplacement
+        forAll(scale, pointI)
+        {
+            if (s.pointFaces()[pointI].size())
+            {
+                scale[pointI] = mag(pointDisplacement[pointI])/distance;
+            }
+            else
+            {
+                scale[pointI] = 1.0;
+            }
+        }
+
+
+        Info<< "Scale        : " << gAverage(scale) << endl;
+        Info<< "Displacement : " << gAverage(pointDisplacement) << endl;
+
+
+
+        // Detect any intersections and scale back
+        PackedBoolList isAffectedPoint;
+        label nIntersections = detectIntersectionPoints
+        (
+            1e-9,       // intersection tolerance
+            s,
+            extendFactor,
+            initialPoints,
+            pointDisplacement,
+
+            checkSelfIntersect,
+            initialIsEdgeIntersecting,
+
+            isAffectedPoint,
+            scale
+        );
+        Info<< "Detected " << nIntersections << " intersections" << nl
+            << endl;
+
+        if (nIntersections == 0)
+        {
+            runTime.write();
+            break;
+        }
+
+
+        // Accumulate all affected points
+        forAll(isAffectedPoint, pointI)
+        {
+            if (isAffectedPoint[pointI])
+            {
+                isScaledPoint[pointI] = true;
+            }
+        }
+
+        // Smear out lowering of scale so any edges not found are
+        // still included
+        for (label i = 0; i < nSmooth; i++)
+        {
+            triSurfacePointScalarField oldScale(scale);
+            oldScale.rename("oldScale");
+            minSmooth
+            (
+                s,
+                PackedBoolList(s.nPoints(), true),
+                oldScale,
+                scale
+            );
+        }
+
+
+        // From scale update the pointDisplacement
+        pointDisplacement *= distance*scale/mag(pointDisplacement);
+
+
+        // Do some smoothing (Lloyds algorithm)
+        lloydsSmoothing(nSmooth, s, isFeaturePoint, edgeStat, isAffectedPoint);
+
+
+        // Update pointDisplacement
+        const pointField& pts = s.points();
+        forAll(meshPoints, i)
+        {
+            label meshPointI = meshPoints[i];
+            pointDisplacement[i] = pts[meshPointI]-initialPoints[meshPointI];
+        }
+
+
+        // Write
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+
+    Info<< "Detected overall intersecting points " << isScaledPoint.count()
+        << " out of " << s.nPoints() << nl << endl;
+
+
+    if (debug)
+    {
+        OBJstream str(runTime.path()/"isScaledPoint.obj");
+        forAll(isScaledPoint, pointI)
+        {
+            if (isScaledPoint[pointI])
+            {
+                str.write(initialPoints[meshPoints[pointI]]);
+            }
+        }
+    }
+
+
+    Info << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C b/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C
index d1e725faee..93f29a2283 100644
--- a/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C
+++ b/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -189,8 +189,6 @@ int main(int argc, char *argv[])
             << endl;
     }
 
-    pointField newPoints(surf1.localPoints());
-
     for (label iter = 0; iter < iters; iter++)
     {
         // Lambda
diff --git a/applications/utilities/surface/surfacePatch/Make/files b/applications/utilities/surface/surfacePatch/Make/files
new file mode 100644
index 0000000000..d983714ab3
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/Make/files
@@ -0,0 +1,7 @@
+searchableSurfaceModifier/searchableSurfaceModifier.C
+searchableSurfaceModifier/autoPatch.C
+searchableSurfaceModifier/cut.C
+
+surfacePatch.C
+
+EXE = $(FOAM_APPBIN)/surfacePatch
diff --git a/applications/utilities/surface/surfacePatch/Make/options b/applications/utilities/surface/surfacePatch/Make/options
new file mode 100644
index 0000000000..a8173dc045
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/Make/options
@@ -0,0 +1,10 @@
+EXE_INC = \
+    -IsearchableSurfaceModifier \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+
+EXE_LIBS = \
+    -ltriSurface \
+    -lmeshTools
diff --git a/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/autoPatch.C b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/autoPatch.C
new file mode 100644
index 0000000000..e8ca9a0125
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/autoPatch.C
@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "autoPatch.H"
+#include "addToRunTimeSelectionTable.H"
+#include "surfaceFeatures.H"
+#include "triSurfaceMesh.H"
+#include "PatchTools.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace searchableSurfaceModifiers
+{
+    defineTypeNameAndDebug(autoPatch, 0);
+    addToRunTimeSelectionTable
+    (
+        searchableSurfaceModifier,
+        autoPatch,
+        dictionary
+    );
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceModifiers::autoPatch::autoPatch
+(
+    const searchableSurfaces& geometry,
+    const dictionary& dict
+)
+:
+    searchableSurfaceModifier(geometry, dict),
+    featureAngle_(readScalar(dict.lookup("featureAngle")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceModifiers::autoPatch::~autoPatch()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::searchableSurfaceModifiers::autoPatch::modify
+(
+    const labelList& regions,
+    searchableSurface& geom
+) const
+{
+    triSurface& surf = refCast<triSurfaceMesh>(geom);
+
+    surfaceFeatures set(surf, 180.0-featureAngle_);
+
+    Info<< nl
+        << "Feature set:" << nl
+        << "    feature points : " << set.featurePoints().size() << nl
+        << "    feature edges  : " << set.featureEdges().size() << nl
+        << "    of which" << nl
+        << "        region edges   : " << set.nRegionEdges() << nl
+        << "        external edges : " << set.nExternalEdges() << nl
+        << "        internal edges : " << set.nInternalEdges() << nl
+        << endl;
+
+
+    // Get per-edge status.
+    boolList borderEdge(surf.nEdges(), false);
+    forAll(set.featureEdges(), i)
+    {
+        borderEdge[set.featureEdges()[i]] = true;
+    }
+
+
+    bool changed = false;
+
+    forAll(regions, i)
+    {
+        label regionI = regions[i];
+
+        labelList subRegion(surf.size(), -1);
+        label nSub = 0;
+
+        forAll(surf, faceI)
+        {
+            if (surf[faceI].region() == regionI && subRegion[faceI] == -1)
+            {
+                PatchTools::markZone(surf, borderEdge, faceI, nSub, subRegion);
+                nSub++;
+            }
+        }
+
+        Info<< "For region " << surf.patches()[regionI].name()
+            << " detected sub regions " << nSub << endl;
+
+
+        if (nSub > 1)
+        {
+            label nOldPatches = surf.patches().size();
+
+            for (label subI = 0; subI < nSub; subI++)
+            {
+                geometricSurfacePatch patch
+                (
+                    surf.patches()[regionI].geometricType(),
+                    surf.patches()[regionI].name() + Foam::name(subI),
+                    surf.patches().size()
+                );
+
+
+                Info<< "For region " << surf.patches()[regionI].name()
+                    << " creating region " << patch.name() << endl;
+
+                surf.patches().append(patch);
+            }
+
+            forAll(subRegion, faceI)
+            {
+                if (subRegion[faceI] != -1)
+                {
+                    surf[faceI].region() = nOldPatches + subRegion[faceI];
+                }
+            }
+            changed = true;
+        }
+    }
+
+    return changed;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/autoPatch.H b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/autoPatch.H
new file mode 100644
index 0000000000..4a7e0459b3
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/autoPatch.H
@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::searchableSurfaceModifier::autoPatch
+
+Description
+    Patchify triangles based on a feature angle. Replacement of surfaceAutoPatch
+
+SourceFiles
+    autoPatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef autoPatch_H
+#define autoPatch_H
+
+#include "searchableSurfaceModifier.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace searchableSurfaceModifiers
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class autoPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class autoPatch
+:
+    public searchableSurfaceModifier
+{
+    // Private data
+
+        //- Feature angle
+        const scalar featureAngle_;
+
+
+   // Private Member Functions
+
+public:
+
+    //- Runtime type information
+    TypeName("autoPatch");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        autoPatch(const searchableSurfaces&, const dictionary&);
+
+        //- Clone
+        autoPtr<searchableSurfaceModifier> clone() const
+        {
+            notImplemented("autoPtr<searchableSurfaceModifier> clone() const");
+            return autoPtr<searchableSurfaceModifier>(NULL);
+        }
+
+
+    //- Destructor
+    virtual ~autoPatch();
+
+
+    // Member Functions
+
+        //- Apply this selector
+        virtual bool modify(const labelList& regions, searchableSurface&) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace searchableSurfaceModifiers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/cut.C b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/cut.C
new file mode 100644
index 0000000000..c04c47da0b
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/cut.C
@@ -0,0 +1,390 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "cut.H"
+#include "addToRunTimeSelectionTable.H"
+#include "searchableSurfaces.H"
+#include "triSurfaceMesh.H"
+#include "searchableBox.H"
+#include "searchableRotatedBox.H"
+#include "surfaceIntersection.H"
+#include "intersectedSurface.H"
+#include "edgeIntersections.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace searchableSurfaceModifiers
+{
+    defineTypeNameAndDebug(cut, 0);
+    addToRunTimeSelectionTable(searchableSurfaceModifier, cut, dictionary);
+}
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::searchableSurfaceModifiers::cut::triangulate
+(
+    const faceList& fcs,
+    pointField& pts,
+    triSurface& cutSurf
+) const
+{
+    label nTris = 0;
+    forAll(fcs, i)
+    {
+        nTris += fcs[i].size()-2;
+    }
+
+    DynamicList<labelledTri> tris(nTris);
+
+    forAll(fcs, i)
+    {
+        const face& f = fcs[i];
+        // Triangulate around vertex 0
+        for (label fp = 1; fp < f.size()-1; fp++)
+        {
+            tris.append(labelledTri(f[0], f[fp], f[f.fcIndex(fp)], i));
+        }
+    }
+    geometricSurfacePatchList patches(fcs.size());
+    forAll(patches, patchI)
+    {
+        patches[patchI] = geometricSurfacePatch
+        (
+            "",
+            "patch" + Foam::name(patchI),
+            patchI
+        );
+    }
+    cutSurf = triSurface(tris.xfer(), patches, pts.xfer());
+}
+
+
+Foam::triSurface& Foam::searchableSurfaceModifiers::cut::triangulate
+(
+    const searchableSurface& cutter,
+    triSurface& cutSurf
+) const
+{
+    if (isA<searchableBox>(cutter))
+    {
+        const searchableBox& bb = refCast<const searchableBox>(cutter);
+
+        pointField pts(bb.points());
+        triangulate(treeBoundBox::faces, pts, cutSurf);
+
+        return cutSurf;
+    }
+    else if (isA<searchableRotatedBox>(cutter))
+    {
+        const searchableRotatedBox& bb =
+            refCast<const searchableRotatedBox>(cutter);
+
+        pointField pts(bb.points());
+        triangulate(treeBoundBox::faces, pts, cutSurf);
+
+        return cutSurf;
+    }
+    else if (isA<triSurfaceMesh>(cutter))
+    {
+        return const_cast<triSurfaceMesh&>
+        (
+            refCast<const triSurfaceMesh>(cutter)
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "searchableSurfaceModifiers::cut::triangulate\n"
+            "(\n"
+            "    const searchableSurface&,\n"
+            "    triSurface&\n"
+            ")"
+        )   << "Triangulation only supported for triSurfaceMesh, searchableBox"
+            << ", not for surface " << cutter.name()
+            << " of type " << cutter.type()
+            << exit(FatalError);
+        return const_cast<triSurfaceMesh&>
+        (
+            refCast<const triSurfaceMesh>(cutter)
+        );
+    }
+}
+
+
+// Keep on shuffling surface points until no more degenerate intersections.
+// Moves both surfaces and updates set of edge cuts.
+bool Foam::searchableSurfaceModifiers::cut::intersectSurfaces
+(
+    triSurface& surf1,
+    edgeIntersections& edgeCuts1,
+    triSurface& surf2,
+    edgeIntersections& edgeCuts2
+) const
+{
+    bool hasMoved1 = false;
+    bool hasMoved2 = false;
+
+    for (label iter = 0; iter < 10; iter++)
+    {
+        Info<< "Determining intersections of surf1 edges with surf2"
+            << " faces" << endl;
+
+        // Determine surface1 edge intersections. Allow surface to be moved.
+
+        // Number of iterations needed to resolve degenerates
+        label nIters1 = 0;
+        {
+            triSurfaceSearch querySurf2(surf2);
+
+            scalarField surf1PointTol
+            (
+                1E-6*edgeIntersections::minEdgeLength(surf1)
+            );
+
+            // Determine raw intersections
+            edgeCuts1 = edgeIntersections
+            (
+                surf1,
+                querySurf2,
+                surf1PointTol
+            );
+
+            // Shuffle a bit to resolve degenerate edge-face hits
+            {
+                pointField points1(surf1.points());
+
+                nIters1 =
+                    edgeCuts1.removeDegenerates
+                    (
+                        5,              // max iterations
+                        surf1,
+                        querySurf2,
+                        surf1PointTol,
+                        points1         // work array
+                    );
+
+                if (nIters1 != 0)
+                {
+                    // Update geometric quantities
+                    surf1.movePoints(points1);
+                    hasMoved1 = true;
+                }
+            }
+        }
+
+        Info<< "Determining intersections of surf2 edges with surf1"
+            << " faces" << endl;
+
+        label nIters2 = 0;
+        {
+            triSurfaceSearch querySurf1(surf1);
+
+            scalarField surf2PointTol
+            (
+                1E-6*edgeIntersections::minEdgeLength(surf2)
+            );
+
+            // Determine raw intersections
+            edgeCuts2 = edgeIntersections
+            (
+                surf2,
+                querySurf1,
+                surf2PointTol
+            );
+
+            // Shuffle a bit to resolve degenerate edge-face hits
+            {
+                pointField points2(surf2.points());
+
+                nIters2 =
+                    edgeCuts2.removeDegenerates
+                    (
+                        5,              // max iterations
+                        surf2,
+                        querySurf1,
+                        surf2PointTol,
+                        points2         // work array
+                    );
+
+                if (nIters2 != 0)
+                {
+                    // Update geometric quantities
+                    surf2.movePoints(points2);
+                    hasMoved2 = true;
+                }
+            }
+        }
+
+
+        if (nIters1 == 0 && nIters2 == 0)
+        {
+            //Info<< "** Resolved all intersections to be proper"
+            //    << "edge-face pierce" << endl;
+            break;
+        }
+    }
+
+    //if (hasMoved1)
+    //{
+    //    fileName newFile("surf1.ftr");
+    //    Info<< "Surface 1 has been moved. Writing to " << newFile
+    //        << endl;
+    //    surf1.write(newFile);
+    //}
+    //
+    //if (hasMoved2)
+    //{
+    //    fileName newFile("surf2.ftr");
+    //    Info<< "Surface 2 has been moved. Writing to " << newFile
+    //        << endl;
+    //    surf2.write(newFile);
+    //}
+
+    return hasMoved1 || hasMoved2;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceModifiers::cut::cut
+(
+    const searchableSurfaces& geometry,
+    const dictionary& dict
+)
+:
+    searchableSurfaceModifier(geometry, dict),
+    cutterNames_(dict_.lookup("cutters"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceModifiers::cut::~cut()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::searchableSurfaceModifiers::cut::modify
+(
+    const labelList& regions,
+    searchableSurface& geom
+) const
+{
+    triSurface& surf = refCast<triSurfaceMesh>(geom);
+
+    bool changed = false;
+
+    // Find the surfaces to cut with
+    forAll(cutterNames_, cutNameI)
+    {
+        labelList geomIDs =
+            findStrings(cutterNames_[cutNameI], geometry_.names());
+
+        forAll(geomIDs, j)
+        {
+            label geomI = geomIDs[j];
+            const searchableSurface& cutter = geometry_[geomI];
+
+            // Triangulate
+            triSurface work;
+            triSurface& cutSurf = triangulate(cutter, work);
+
+            // Determine intersection (with perturbation)
+            edgeIntersections edge1Cuts;
+            edgeIntersections edge2Cuts;
+            intersectSurfaces
+            (
+                surf,
+                edge1Cuts,
+                cutSurf,
+                edge2Cuts
+            );
+
+
+            // Determine intersection edges
+            surfaceIntersection inter(surf, edge1Cuts, cutSurf, edge2Cuts);
+
+
+            // Use intersection edges to cut up faces. (does all the hard work)
+            intersectedSurface surf3(surf, true, inter);
+
+
+            // Mark triangles based on whether they are inside or outside
+            List<volumeType> volTypes;
+            cutter.getVolumeType(surf3.faceCentres(), volTypes);
+
+            label nInside = 0;
+            forAll(volTypes, i)
+            {
+                if (volTypes[i] == volumeType::INSIDE)
+                {
+                    nInside++;
+                }
+            }
+
+            // Add a patch for inside the box
+            if (nInside > 0 && surf3.patches().size() > 0)
+            {
+                geometricSurfacePatchList newPatches(surf3.patches());
+                label sz = newPatches.size();
+                newPatches.setSize(sz+1);
+                newPatches[sz] = geometricSurfacePatch
+                (
+                    newPatches[sz-1].geometricType(),
+                    newPatches[sz-1].name() + "_inside",
+                    newPatches[sz-1].index()
+                );
+
+                Info<< "Moving " << nInside << " out of " << surf3.size()
+                    << " triangles to region "
+                    << newPatches[sz].name() << endl;
+
+
+                List<labelledTri> newTris(surf3);
+                forAll(volTypes, i)
+                {
+                    if (volTypes[i] == volumeType::INSIDE)
+                    {
+                        newTris[i].region() = sz;
+                    }
+                }
+                pointField newPoints(surf3.points());
+                surf = triSurface(newTris.xfer(), newPatches, newPoints.xfer());
+
+                changed = true;
+            }
+        }
+    }
+
+    return changed;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/cut.H b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/cut.H
new file mode 100644
index 0000000000..4081b753b3
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/cut.H
@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::searchableSurfaceModifier::cut
+
+Description
+    Patchify triangles based on orientation w.r.t other (triangulated or
+    triangulatable) surfaces
+
+SourceFiles
+    cut.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cut_H
+#define cut_H
+
+#include "searchableSurfaceModifier.H"
+#include "wordReList.H"
+#include "faceList.H"
+#include "pointField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class edgeIntersections;
+
+namespace searchableSurfaceModifiers
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class cut Declaration
+\*---------------------------------------------------------------------------*/
+
+class cut
+:
+    public searchableSurfaceModifier
+{
+    // Private data
+
+        //- Name of surfaces to cut with
+        const wordReList cutterNames_;
+
+
+   // Private Member Functions
+
+        //- Triangulate faces around 0th vertex
+        void triangulate(const faceList&, pointField&, triSurface&) const;
+
+        //- Triangulate searchableSurface (currently only supported for
+        //  searchableBox and triSurfaceMesh)
+        triSurface& triangulate(const searchableSurface&, triSurface&) const;
+
+        //- Intersect surfaces. Perturb to avoid degenerates.
+        bool intersectSurfaces
+        (
+            triSurface& surf1,
+            edgeIntersections& edgeCuts1,
+            triSurface& surf2,
+            edgeIntersections& edgeCuts2
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("cut");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        cut(const searchableSurfaces&, const dictionary&);
+
+        //- Clone
+        autoPtr<searchableSurfaceModifier> clone() const
+        {
+            notImplemented("autoPtr<searchableSurfaceModifier> clone() const");
+            return autoPtr<searchableSurfaceModifier>(NULL);
+        }
+
+
+    //- Destructor
+    virtual ~cut();
+
+
+    // Member Functions
+
+        //- Apply this selector
+        virtual bool modify(const labelList& regions, searchableSurface&) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace searchableSurfaceModifiers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/searchableSurfaceModifier.C b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/searchableSurfaceModifier.C
new file mode 100644
index 0000000000..e28036012c
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/searchableSurfaceModifier.C
@@ -0,0 +1,87 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "searchableSurfaceModifier.H"
+#include "triSurface.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(searchableSurfaceModifier, 0);
+    defineRunTimeSelectionTable(searchableSurfaceModifier, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceModifier::searchableSurfaceModifier
+(
+    const searchableSurfaces& geometry,
+    const dictionary& dict
+)
+:
+    geometry_(geometry),
+    dict_(dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceModifier::~searchableSurfaceModifier()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::searchableSurfaceModifier>
+Foam::searchableSurfaceModifier::New
+(
+    const word& type,
+    const searchableSurfaces& geometry,
+    const dictionary& dict
+)
+{
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(type);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "searchableSurfaceModifier::New"
+            "(const word&, const searchableSurfaces&, const dictionary&)"
+        )   << "Unknown searchableSurfaceModifier type "
+            << type << nl << nl
+            << "Valid searchableSurfaceModifier types : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<searchableSurfaceModifier>(cstrIter()(geometry, dict));
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/searchableSurfaceModifier.H b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/searchableSurfaceModifier.H
new file mode 100644
index 0000000000..8955dea741
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/searchableSurfaceModifier.H
@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::searchableSurfaceModifier
+
+Description
+    Changing a surface
+
+SourceFiles
+    searchableSurfaceModifier.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchableSurfaceModifier_H
+#define searchableSurfaceModifier_H
+
+#include "dictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "labelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class triSurface;
+class searchableSurface;
+class searchableSurfaces;
+
+/*---------------------------------------------------------------------------*\
+                     Class searchableSurfaceModifier Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchableSurfaceModifier
+{
+protected:
+
+    // Protected data
+
+        const searchableSurfaces& geometry_;
+
+        //- Input dictionary
+        const dictionary dict_;
+
+public:
+
+    //- Runtime type information
+    TypeName("searchableSurfaceModifier");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            searchableSurfaceModifier,
+            dictionary,
+            (
+                const searchableSurfaces& geometry,
+                const dictionary& dict
+            ),
+            (geometry, dict)
+        );
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        searchableSurfaceModifier(const searchableSurfaces&, const dictionary&);
+
+        //- Clone
+        autoPtr<searchableSurfaceModifier> clone() const
+        {
+            notImplemented("autoPtr<searchableSurfaceModifier> clone() const");
+            return autoPtr<searchableSurfaceModifier>(NULL);
+        }
+
+
+    // Selectors
+
+        //- Return a reference to the selected searchableSurfaceModifier
+        static autoPtr<searchableSurfaceModifier> New
+        (
+            const word& type,
+            const searchableSurfaces&,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~searchableSurfaceModifier();
+
+
+    // Member Functions
+
+        //- Do any operation on surface. Return true if anything changed
+        virtual bool modify(const labelList&, searchableSurface&) const = 0;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfacePatch/surfacePatch.C b/applications/utilities/surface/surfacePatch/surfacePatch.C
new file mode 100644
index 0000000000..3df36f6ec6
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/surfacePatch.C
@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Application
+    surfacePatch
+
+Description
+    Patches (regionises) a surface using a user-selectable method.
+
+    The current methods are either based on a geometric feature angle
+    (a replacement for the surfaceAutoPatch functionality) or on intersecting
+    a set of searchableSurfaces - this will re-triangulate the intersecting
+    triangles and regonise the resulting triangles according to being
+    inside/outside the surface.
+
+\*---------------------------------------------------------------------------*/
+
+#include "triSurface.H"
+#include "argList.H"
+#include "Time.H"
+#include "triSurfaceMesh.H"
+#include "searchableSurfaces.H"
+#include "searchableSurfaceModifier.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::validArgs.clear();
+    argList::noParallel();
+    #include "addDictOption.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+    const word dictName("surfacePatchDict");
+    #include "setSystemRunTimeDictionaryIO.H"
+    const IOdictionary meshDict(dictIO);
+
+
+    // Read geometry
+    // ~~~~~~~~~~~~~
+
+    searchableSurfaces allGeometry
+    (
+        IOobject
+        (
+            "abc",                      // dummy name
+            runTime.constant(),         // instance
+            triSurfaceMesh::meshSubDir, // local
+            runTime,                    // registry
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        ),
+        meshDict.subDict("geometry"),
+        meshDict.lookupOrDefault("singleRegionName", true)
+    );
+
+
+    const dictionary& surfacesDict = meshDict.subDict("surfaces");
+
+    forAllConstIter(dictionary, surfacesDict, surfacesIter)
+    {
+        if (surfacesIter().isDict())
+        {
+            const word& surfName = surfacesIter().keyword();
+            const dictionary& surfDict = surfacesIter().dict();
+
+            // Look up surface
+            searchableSurface& surf = allGeometry[surfName];
+
+
+            bool changed = false;
+
+            // Check for modifier on global level
+            if (surfDict.found("type"))
+            {
+                autoPtr<searchableSurfaceModifier> modifier
+                (
+                    searchableSurfaceModifier::New
+                    (
+                        surfDict.lookup("type"),
+                        allGeometry,
+                        surfDict
+                    )
+                );
+
+                if (modifier().modify(identity(surf.regions().size()), surf))
+                {
+                    changed = true;
+                }
+            }
+
+            // Check for modifier on a per-region level
+            if (surfDict.found("regions"))
+            {
+                const dictionary& regionsDict = surfDict.subDict("regions");
+                forAllConstIter(dictionary, regionsDict, regionsIter)
+                {
+                    const dictionary& regionDict = regionsIter().dict();
+                    const keyType& regionName = regionsIter().keyword();
+
+                    autoPtr<searchableSurfaceModifier> modifier
+                    (
+                        searchableSurfaceModifier::New
+                        (
+                            regionDict.lookup("type"),
+                            allGeometry,
+                            regionDict
+                        )
+                    );
+
+                    labelList regionIDs =
+                        findStrings(regionName, surf.regions());
+                    if (modifier().modify(regionIDs, surf))
+                    {
+                        changed = true;
+                    }
+                }
+            }
+
+
+            if (changed)
+            {
+                const fileName name(surf.name());
+                surf.rename(name.lessExt() + "_patched." + name.ext());
+
+                Info<< "Writing repatched surface " << name << " to "
+                    << surf.name() << nl << endl;
+
+                surf.write();
+            }
+            else
+            {
+                Info<< "Surface " << surf.name() << " unchanged" << nl << endl;
+            }
+        }
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfacePatch/surfacePatchDict b/applications/utilities/surface/surfacePatch/surfacePatchDict
new file mode 100644
index 0000000000..5e0f5f19a5
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/surfacePatchDict
@@ -0,0 +1,69 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      createPatchDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+geometry
+{
+    box
+    {
+        type    searchableRotatedBox;
+
+        // box rotated 45 degrees around z axis
+        span    (1 1 1);
+        origin  (0 0 0);
+        e1      (1 1 0);
+        e3      (0 0 1);
+    }
+
+    singleTri.obj
+    {
+        type triSurfaceMesh;
+    }
+
+    unitCube_me.stl
+    {
+        type triSurfaceMesh;
+    }
+};
+
+
+
+surfaces
+{
+    unitCube_me.stl
+    {
+        regions
+        {
+            // Divide region 'maxZ' into multiple regions according to
+            // a geometric feature angle
+            maxZ
+            {
+                type            autoPatch;
+                featureAngle    60;
+            }
+        }
+    }
+
+    singleTri.obj
+    {
+        // Regionise surface according to triangles (after cutting) being
+        // inside or outside the 'box' geometry
+        type        cut;
+        cutters     (box);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C b/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C
index 6e4db35a0b..001745fbca 100644
--- a/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C
+++ b/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,21 +31,17 @@ Description
 
 Note
     - best decomposition option is hierarchGeomDecomp since
-    guarantees square decompositions.
+      guarantees square decompositions.
     - triangles might be present on multiple processors.
     - merging uses geometric tolerance so take care with writing precision.
 
 \*---------------------------------------------------------------------------*/
 
-#include "treeBoundBox.H"
-#include "FixedList.H"
 #include "argList.H"
 #include "Time.H"
 #include "polyMesh.H"
 #include "distributedTriSurfaceMesh.H"
 #include "mapDistribute.H"
-#include "triSurfaceFields.H"
-#include "Pair.H"
 
 using namespace Foam;
 
@@ -172,7 +168,7 @@ int main(int argc, char *argv[])
         "triSurface",         // local
         runTime,              // registry
         IOobject::MUST_READ,
-        IOobject::NO_WRITE
+        IOobject::AUTO_WRITE
     );
 
     const fileName actualPath(io.filePath());
@@ -215,6 +211,7 @@ int main(int argc, char *argv[])
         Info<< "Writing dummy bounds dictionary to " << ioDict.name()
             << nl << endl;
 
+        // Force writing in ascii
         ioDict.regIOobject::writeObject
         (
             IOstream::ASCII,
@@ -239,23 +236,18 @@ int main(int argc, char *argv[])
             (
                 IOobject
                 (
-                    surfMesh.searchableSurface::name(),     // name
-                    surfMesh.searchableSurface::instance(), // instance
+                    "faceCentres",                                  // name
+                    surfMesh.searchableSurface::time().timeName(),  // instance
                     surfMesh.searchableSurface::local(),    // local
                     surfMesh,
                     IOobject::NO_READ,
                     IOobject::AUTO_WRITE
                 ),
                 surfMesh,
-                dimLength
+                dimLength,
+                s.faceCentres()
             )
         );
-        triSurfaceVectorField& fc = fcPtr();
-
-        forAll(fc, triI)
-        {
-            fc[triI] = s[triI].centre(s.points());
-        }
 
         // Steal pointer and store object on surfMesh
         fcPtr.ptr()->store();
@@ -290,7 +282,7 @@ int main(int argc, char *argv[])
 
 
     Info<< "Writing surface." << nl << endl;
-    surfMesh.searchableSurface::write();
+    surfMesh.objectRegistry::write();
 
     Info<< "End\n" << endl;
 
diff --git a/applications/utilities/surface/surfaceSubset/surfaceSubset.C b/applications/utilities/surface/surfaceSubset/surfaceSubset.C
index c0ae6bd241..c2a168ba68 100644
--- a/applications/utilities/surface/surfaceSubset/surfaceSubset.C
+++ b/applications/utilities/surface/surfaceSubset/surfaceSubset.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -387,6 +387,8 @@ int main(int argc, char *argv[])
 
     surf2.write(outFileName);
 
+    Info<< "\nEnd\n" << endl;
+
     return 0;
 }
 
diff --git a/src/edgeMesh/edgeMesh.C b/src/edgeMesh/edgeMesh.C
index 0e8cb2acad..947472af78 100644
--- a/src/edgeMesh/edgeMesh.C
+++ b/src/edgeMesh/edgeMesh.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -286,7 +286,11 @@ void Foam::edgeMesh::scalePoints(const scalar scaleFactor)
 }
 
 
-void Foam::edgeMesh::mergePoints(const scalar mergeDist)
+void Foam::edgeMesh::mergePoints
+(
+    const scalar mergeDist,
+    labelList& reversePointMap
+)
 {
     pointField newPoints;
     labelList pointMap;
@@ -307,6 +311,9 @@ void Foam::edgeMesh::mergePoints(const scalar mergeDist)
 
         points_.transfer(newPoints);
 
+        // connectivity changed
+        pointEdgesPtr_.clear();
+
         // Renumber and make sure e[0] < e[1] (not really necessary)
         forAll(edges_, edgeI)
         {
@@ -383,6 +390,9 @@ void Foam::edgeMesh::mergeEdges()
     {
         edges_[iter()] = iter.key();
     }
+
+    // connectivity changed
+    pointEdgesPtr_.clear();
 }
 
 
diff --git a/src/edgeMesh/edgeMesh.H b/src/edgeMesh/edgeMesh.H
index d3397a290b..9962ca9176 100644
--- a/src/edgeMesh/edgeMesh.H
+++ b/src/edgeMesh/edgeMesh.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -247,11 +247,12 @@ public:
         //- Scale points. A non-positive factor is ignored
         virtual void scalePoints(const scalar);
 
-        //- Merge common points (points within mergeDist)
-        void mergePoints(const scalar mergeDist);
+        //- Merge common points (points within mergeDist). Return map from
+        //  old to new points.
+        virtual void mergePoints(const scalar mergeDist, labelList&);
 
-        //- Merge similar edges
-        void mergeEdges();
+        //- Merge duplicate edges
+        virtual void mergeEdges();
 
 
     // Write
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
index d7bfa4b2f6..f2e379b493 100644
--- a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,6 +32,8 @@ License
 #include "DynamicField.H"
 #include "edgeMeshFormatsCore.H"
 #include "IOmanip.H"
+#include "searchableSurface.H"
+#include "triSurfaceMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -118,7 +120,6 @@ Foam::wordHashSet Foam::extendedEdgeMesh::writeTypes()
 }
 
 
-
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
 bool Foam::extendedEdgeMesh::canReadType
@@ -176,7 +177,7 @@ Foam::extendedEdgeMesh::classifyFeaturePoint
     label ptI
 ) const
 {
-    labelList ptEds(pointEdges()[ptI]);
+    const labelList& ptEds(pointEdges()[ptI]);
 
     label nPtEds = ptEds.size();
     label nExternal = 0;
@@ -217,46 +218,182 @@ Foam::extendedEdgeMesh::classifyFeaturePoint
 }
 
 
-Foam::extendedEdgeMesh::edgeStatus
-Foam::extendedEdgeMesh::classifyEdge
+void Foam::extendedEdgeMesh::cut
 (
-    const List<vector>& norms,
-    const labelList& edNorms,
-    const vector& fC0tofC1
+    const searchableSurface& surf,
+
+    labelList& pointMap,
+    labelList& edgeMap,
+    labelList& pointsFromEdge,
+    labelList& oldEdge,
+    labelList& surfTri
 )
 {
-    label nEdNorms = edNorms.size();
+    const edgeList& edges = this->edges();
+    const pointField& points = this->points();
 
-    if (nEdNorms == 1)
+
+    List<List<pointIndexHit> > edgeHits(edges.size());
     {
-        return OPEN;
+        pointField start(edges.size());
+        pointField end(edges.size());
+        forAll(edges, edgeI)
+        {
+            const edge& e = edges[edgeI];
+            start[edgeI] = points[e[0]];
+            end[edgeI] = points[e[1]];
+        }
+        surf.findLineAll(start, end, edgeHits);
     }
-    else if (nEdNorms == 2)
-    {
-        const vector& n0(norms[edNorms[0]]);
-        const vector& n1(norms[edNorms[1]]);
 
-        if ((n0 & n1) > cosNormalAngleTol_)
+    // Count number of hits
+    label nHits = 0;
+    forAll(edgeHits, edgeI)
+    {
+        nHits += edgeHits[edgeI].size();
+    }
+
+    DynamicField<point> newPoints(points);
+    DynamicList<label> newToOldPoint(identity(points.size()));
+
+    newPoints.setCapacity(newPoints.size()+nHits);
+    newToOldPoint.setCapacity(newPoints.capacity());
+
+    DynamicList<edge> newEdges(edges);
+    DynamicList<label> newToOldEdge(identity(edges.size()));
+
+    newEdges.setCapacity(newEdges.size()+nHits);
+    newToOldEdge.setCapacity(newEdges.capacity());
+
+    // Information on additional points
+    DynamicList<label> dynPointsFromEdge(nHits);
+    DynamicList<label> dynOldEdge(nHits);
+    DynamicList<label> dynSurfTri(nHits);
+
+    forAll(edgeHits, edgeI)
+    {
+        const List<pointIndexHit>& eHits = edgeHits[edgeI];
+
+        if (eHits.size())
         {
-            return FLAT;
-        }
-        else if ((fC0tofC1 & n0) > 0.0)
-        {
-            return INTERNAL;
-        }
-        else
-        {
-            return EXTERNAL;
+            label prevPtI = edges[edgeI][0];
+            forAll(eHits, eHitI)
+            {
+                label newPtI = newPoints.size();
+
+                newPoints.append(eHits[eHitI].hitPoint());
+                newToOldPoint.append(edges[edgeI][0]);  // map from start point
+                dynPointsFromEdge.append(newPtI);
+                dynOldEdge.append(edgeI);
+                dynSurfTri.append(eHits[eHitI].index());
+
+                if (eHitI == 0)
+                {
+                    newEdges[edgeI] = edge(prevPtI, newPtI);
+                }
+                else
+                {
+                    newEdges.append(edge(prevPtI, newPtI));
+                    newToOldEdge.append(edgeI);
+                }
+                prevPtI = newPtI;
+            }
+            newEdges.append(edge(prevPtI, edges[edgeI][1]));
+            newToOldEdge.append(edgeI);
         }
     }
-    else if (nEdNorms > 2)
+
+    pointField allPoints;
+    allPoints.transfer(newPoints);
+    pointMap.transfer(newToOldPoint);
+
+    edgeList allEdges;
+    allEdges.transfer(newEdges);
+    edgeMap.transfer(newToOldEdge);
+
+    pointsFromEdge.transfer(dynPointsFromEdge);
+    oldEdge.transfer(dynOldEdge);
+    surfTri.transfer(dynSurfTri);
+
+    // Update local information
+    autoMap(allPoints, allEdges, pointMap, edgeMap);
+}
+
+
+void Foam::extendedEdgeMesh::select
+(
+    const searchableSurface& surf,
+    const volumeType volType,   // side to keep
+    labelList& pointMap,        // sub to old points
+    labelList& edgeMap          // sub to old edges
+)
+{
+    const edgeList& edges = this->edges();
+    const pointField& points = this->points();
+
+    // Test edge centres for inside/outside
+    if (volType == volumeType::INSIDE || volType == volumeType::OUTSIDE)
     {
-        return MULTIPLE;
+        pointField edgeCentres(edges.size());
+        forAll(edgeCentres, edgeI)
+        {
+            const edge& e = edges[edgeI];
+            edgeCentres[edgeI] = e.centre(points);
+        }
+        List<volumeType> volTypes;
+        surf.getVolumeType(edgeCentres, volTypes);
+
+        // Extract edges on correct side
+        edgeMap.setSize(edges.size());
+        label compactEdgeI = 0;
+
+        forAll(volTypes, edgeI)
+        {
+            if (volTypes[edgeI] == volType)
+            {
+                edgeMap[compactEdgeI++] = edgeI;
+            }
+        }
+        edgeMap.setSize(compactEdgeI);
+
+        // Extract used points
+        labelList pointToCompact(points.size(), -1);
+        forAll(edgeMap, i)
+        {
+            const edge& e = edges[edgeMap[i]];
+            pointToCompact[e[0]] = labelMax;       // tag with a value
+            pointToCompact[e[1]] = labelMax;
+        }
+
+        pointMap.setSize(points.size());
+        label compactPointI = 0;
+        forAll(pointToCompact, pointI)
+        {
+            if (pointToCompact[pointI] != -1)
+            {
+                pointToCompact[pointI] = compactPointI;
+                pointMap[compactPointI++] = pointI;
+            }
+        }
+        pointMap.setSize(compactPointI);
+        pointField subPoints(points, pointMap);
+
+        // Renumber edges
+        edgeList subEdges(edgeMap.size());
+        forAll(edgeMap, i)
+        {
+            const edge& e = edges[edgeMap[i]];
+            subEdges[i][0] = pointToCompact[e[0]];
+            subEdges[i][1] = pointToCompact[e[1]];
+        }
+
+        // Reset primitives and map other quantities
+        autoMap(subPoints, subEdges, pointMap, edgeMap);
     }
     else
     {
-        // There is a problem - the edge has no normals
-        return NONE;
+        pointMap = identity(points.size());
+        edgeMap = identity(edges.size());
     }
 }
 
@@ -287,10 +424,7 @@ Foam::extendedEdgeMesh::extendedEdgeMesh()
 {}
 
 
-Foam::extendedEdgeMesh::extendedEdgeMesh
-(
-    const extendedEdgeMesh& fem
-)
+Foam::extendedEdgeMesh::extendedEdgeMesh(const extendedEdgeMesh& fem)
 :
     edgeMesh(fem),
     concaveStart_(fem.concaveStart()),
@@ -1332,126 +1466,649 @@ void Foam::extendedEdgeMesh::flipNormals()
 }
 
 
-void Foam::extendedEdgeMesh::writeObj
+void Foam::extendedEdgeMesh::autoMap
 (
-    const fileName& prefix
-) const
+    const pointField& subPoints,
+    const edgeList& subEdges,
+    const labelList& pointMap,
+    const labelList& edgeMap
+)
+{
+    // Determine slicing for subEdges
+    label subIntStart = edgeMap.size();
+    label subFlatStart = edgeMap.size();
+    label subOpenStart = edgeMap.size();
+    label subMultipleStart = edgeMap.size();
+
+    forAll(edgeMap, subEdgeI)
+    {
+        label edgeI = edgeMap[subEdgeI];
+        if (edgeI >= internalStart() && subIntStart == edgeMap.size())
+        {
+            subIntStart = subEdgeI;
+        }
+        if (edgeI >= flatStart() && subFlatStart == edgeMap.size())
+        {
+            subFlatStart = subEdgeI;
+        }
+        if (edgeI >= openStart() && subOpenStart == edgeMap.size())
+        {
+            subOpenStart = subEdgeI;
+        }
+        if (edgeI >= multipleStart() && subMultipleStart == edgeMap.size())
+        {
+            subMultipleStart = subEdgeI;
+        }
+    }
+
+
+    // Determine slicing for subPoints
+
+    label subConcaveStart = pointMap.size();
+    label subMixedStart = pointMap.size();
+    label subNonFeatStart = pointMap.size();
+
+    forAll(pointMap, subPointI)
+    {
+        label pointI = pointMap[subPointI];
+        if (pointI >= concaveStart() && subConcaveStart == pointMap.size())
+        {
+            subConcaveStart = subPointI;
+        }
+        if (pointI >= mixedStart() && subMixedStart == pointMap.size())
+        {
+            subMixedStart = subPointI;
+        }
+        if
+        (
+            pointI >= nonFeatureStart()
+         && subNonFeatStart == pointMap.size()
+        )
+        {
+            subNonFeatStart = subPointI;
+        }
+    }
+
+
+
+    // Compact region edges
+    labelList subRegionEdges;
+    {
+        PackedBoolList isRegionEdge(edges().size());
+        forAll(regionEdges(), i)
+        {
+            label edgeI = regionEdges()[i];
+            isRegionEdge[edgeI] = true;
+        }
+
+        DynamicList<label> newRegionEdges(regionEdges().size());
+        forAll(edgeMap, subEdgeI)
+        {
+            label edgeI = edgeMap[subEdgeI];
+            if (isRegionEdge[edgeI])
+            {
+                newRegionEdges.append(subEdgeI);
+            }
+        }
+        subRegionEdges.transfer(newRegionEdges);
+    }
+
+
+    labelListList subFeaturePointEdges;
+    if (featurePointEdges().size())
+    {
+        subFeaturePointEdges.setSize(subNonFeatStart);
+        for (label subPointI = 0; subPointI < subNonFeatStart; subPointI++)
+        {
+            label pointI = pointMap[subPointI];
+            const labelList& pEdges = featurePointEdges()[pointI];
+
+            labelList& subPEdges = subFeaturePointEdges[subPointI];
+            subPEdges.setSize(pEdges.size());
+
+            if (pEdges.size())
+            {
+                forAll(pEdges, i)
+                {
+                    subPEdges[i] = edgeMap[pEdges[i]];
+                }
+            }
+        }
+    }
+
+
+    vectorField subEdgeDirections(edgeDirections(), edgeMap);
+
+    // Find used normals
+    labelList reverseNormalMap(normals().size(), -1);
+    DynamicList<label> normalMap(normals().size());
+
+    {
+        PackedBoolList isSubNormal(normals().size());
+        for (label subPointI = 0; subPointI < subNonFeatStart; subPointI++)
+        {
+            label pointI = pointMap[subPointI];
+            const labelList& pNormals = featurePointNormals()[pointI];
+
+            forAll(pNormals, i)
+            {
+                isSubNormal[pNormals[i]] = true;
+            }
+        }
+        forAll(edgeMap, subEdgeI)
+        {
+            label edgeI = edgeMap[subEdgeI];
+            const labelList& eNormals = edgeNormals()[edgeI];
+
+            forAll(eNormals, i)
+            {
+                isSubNormal[eNormals[i]] = true;
+            }
+        }
+
+        forAll(isSubNormal, normalI)
+        {
+            if (isSubNormal[normalI])
+            {
+                label subNormalI = normalMap.size();
+                reverseNormalMap[normalI] = subNormalI;
+                normalMap.append(subNormalI);
+            }
+        }
+    }
+
+
+    // Use compaction map on data referencing normals
+    labelListList subNormalDirections;
+
+    if (normalDirections().size())
+    {
+        subNormalDirections.setSize(edgeMap.size());
+
+        forAll(edgeMap, subEdgeI)
+        {
+            label edgeI = edgeMap[subEdgeI];
+            const labelList& eNormals = normalDirections()[edgeI];
+
+            labelList& subNormals = subNormalDirections[subEdgeI];
+            subNormals.setSize(eNormals.size());
+            forAll(eNormals, i)
+            {
+                if (eNormals[i] >= 0)
+                {
+                    subNormals[i] = reverseNormalMap[eNormals[i]];
+                }
+                else
+                {
+                    subNormals[i] = -1;
+                }
+            }
+        }
+    }
+
+    labelListList subEdgeNormals(edgeMap.size());
+    forAll(edgeMap, subEdgeI)
+    {
+        label edgeI = edgeMap[subEdgeI];
+        const labelList& eNormals = edgeNormals()[edgeI];
+        labelList& subNormals = subEdgeNormals[subEdgeI];
+
+        subNormals = UIndirectList<label>(reverseNormalMap, eNormals);
+    }
+
+    labelListList subPointNormals(pointMap.size());
+    for (label subPointI = 0; subPointI < subNonFeatStart; subPointI++)
+    {
+        label pointI = pointMap[subPointI];
+        const labelList& pNormals = featurePointNormals()[pointI];
+        labelList& subNormals = subPointNormals[subPointI];
+
+        subNormals = UIndirectList<label>(reverseNormalMap, pNormals);
+    }
+
+    // Use compaction map to compact normal data
+    vectorField subNormals(normals(), normalMap);
+
+    List<extendedEdgeMesh::sideVolumeType> subNormalVolumeTypes;
+    if (normalVolumeTypes().size())
+    {
+        subNormalVolumeTypes =
+            UIndirectList<extendedEdgeMesh::sideVolumeType>
+            (
+                normalVolumeTypes(),
+                normalMap
+            );
+    }
+
+    extendedEdgeMesh subMesh
+    (
+        subPoints,
+        subEdges,
+
+        // Feature points slices
+        subConcaveStart,
+        subMixedStart,
+        subNonFeatStart,
+
+        // Feature edges slices
+        subIntStart,
+        subFlatStart,
+        subOpenStart,
+        subMultipleStart,
+
+        // All normals
+        subNormals,
+        subNormalVolumeTypes,
+
+        // Per edge edge vector
+        subEdgeDirections,
+
+        // Per edge list of normal indices
+        subNormalDirections,
+        // Per edge list of normal indices
+        subEdgeNormals,
+
+        // Per point list of normal indices
+        subPointNormals,
+        subFeaturePointEdges,
+        subRegionEdges
+    );
+
+    transfer(subMesh);
+}
+
+
+void Foam::extendedEdgeMesh::trim
+(
+    const searchableSurface& surf,
+    const volumeType volType,
+    labelList& pointMap,
+    labelList& edgeMap
+)
+{
+    // Cut edges with the other surfaces
+
+    labelList allPointMap;      // from all to original point
+    labelList allEdgeMap;       // from all to original edge
+
+    labelList pointsFromEdge;   // list of new points created by cutting
+    labelList oldEdge;          // for each of these points the orginal edge
+    labelList surfTri;          // for each of these points the surface triangle
+    cut
+    (
+        surf,
+        allPointMap,
+        allEdgeMap,
+        pointsFromEdge,
+        oldEdge,
+        surfTri
+    );
+
+    const label nOldPoints = points().size();
+
+    // Remove outside edges and compact
+
+    labelList subPointMap;  // sub to old points
+    labelList subEdgeMap;   // sub to old edges
+    select(surf, volType, subPointMap, subEdgeMap);
+
+    // Update overall point maps
+    pointMap = UIndirectList<label>(allPointMap, subPointMap);
+    edgeMap = UIndirectList<label>(allEdgeMap, subEdgeMap);
+
+    // Extract current point and edge status
+    List<edgeStatus> edgeStat(edges().size());
+    List<pointStatus> pointStat(points().size());
+    forAll(edgeStat, edgeI)
+    {
+        edgeStat[edgeI] = getEdgeStatus(edgeI);
+    }
+    forAll(pointStat, pointI)
+    {
+        pointStat[pointI] = getPointStatus(pointI);
+    }
+
+    // Re-classify exposed points (from cutting)
+    labelList oldPointToIndex(nOldPoints, -1);
+    forAll(pointsFromEdge, i)
+    {
+        oldPointToIndex[pointsFromEdge[i]] = i;
+    }
+    forAll(subPointMap, pointI)
+    {
+        label oldPointI = subPointMap[pointI];
+        label index = oldPointToIndex[oldPointI];
+        if (index != -1)
+        {
+            pointStat[pointI] = classifyFeaturePoint(pointI);
+        }
+    }
+
+    // Reset based on new point and edge status
+    labelList sortedToOriginalPoint;
+    labelList sortedToOriginalEdge;
+    setFromStatus
+    (
+        pointStat,
+        edgeStat,
+        sortedToOriginalPoint,
+        sortedToOriginalEdge
+    );
+
+    // Update the overall pointMap, edgeMap
+    pointMap = UIndirectList<label>(pointMap, sortedToOriginalPoint)();
+    edgeMap = UIndirectList<label>(edgeMap, sortedToOriginalEdge)();
+}
+
+
+void Foam::extendedEdgeMesh::setFromStatus
+(
+    const List<extendedEdgeMesh::pointStatus>& pointStat,
+    const List<extendedEdgeMesh::edgeStatus>& edgeStat,
+    labelList& sortedToOriginalPoint,
+    labelList& sortedToOriginalEdge
+)
+{
+    // Use pointStatus and edgeStatus to determine new ordering
+    label pointConcaveStart;
+    label pointMixedStart;
+    label pointNonFeatStart;
+
+    label edgeInternalStart;
+    label edgeFlatStart;
+    label edgeOpenStart;
+    label edgeMultipleStart;
+    sortedOrder
+    (
+        pointStat,
+        edgeStat,
+        sortedToOriginalPoint,
+        sortedToOriginalEdge,
+
+        pointConcaveStart,
+        pointMixedStart,
+        pointNonFeatStart,
+
+        edgeInternalStart,
+        edgeFlatStart,
+        edgeOpenStart,
+        edgeMultipleStart
+    );
+
+
+    // Order points and edges
+    labelList reversePointMap(points().size(), -1);
+    forAll(sortedToOriginalPoint, sortedI)
+    {
+        reversePointMap[sortedToOriginalPoint[sortedI]] = sortedI;
+    }
+
+    edgeList sortedEdges(UIndirectList<edge>(edges(), sortedToOriginalEdge)());
+    forAll(sortedEdges, sortedI)
+    {
+        inplaceRenumber(reversePointMap, sortedEdges[sortedI]);
+    }
+
+    // Update local data
+    autoMap
+    (
+        pointField(points(), sortedToOriginalPoint),
+        sortedEdges,
+        sortedToOriginalPoint,
+        sortedToOriginalEdge
+    );
+
+    // Reset the slice starts
+    concaveStart_ = pointConcaveStart;
+    mixedStart_ = pointMixedStart;
+    nonFeatureStart_ = pointNonFeatStart;
+    internalStart_ = edgeInternalStart;
+    flatStart_ = edgeFlatStart;
+    openStart_ = edgeOpenStart;
+    multipleStart_ = edgeMultipleStart;
+}
+
+
+bool Foam::extendedEdgeMesh::mergePointsAndSort
+(
+    const scalar mergeDist,
+    labelList& pointMap,
+    labelList& edgeMap
+)
+{
+    const label nOldPoints = points().size();
+
+    // Detect and merge collocated feature points
+    labelList oldToMerged;
+    label nNewPoints = ::Foam::mergePoints
+    (
+        points(),
+        SMALL,
+        false,
+        oldToMerged
+    );
+
+    pointMap.setSize(nNewPoints);
+    pointMap = -1;
+    forAll(oldToMerged, oldI)
+    {
+        label newI = oldToMerged[oldI];
+        if (pointMap[newI] == -1)
+        {
+            pointMap[newI] = oldI;
+        }
+    }
+
+    // Renumber edges
+    edgeList newEdges(edges().size());
+    forAll(edges(), edgeI)
+    {
+        const edge& oldE = edges()[edgeI];
+        newEdges[edgeI] = edge(oldToMerged[oldE[0]], oldToMerged[oldE[1]]);
+    }
+
+    // Shuffle basic information (reorders point data)
+    autoMap
+    (
+        pointField(points(), pointMap),
+        newEdges,
+        pointMap,
+        identity(newEdges.size())
+    );
+
+    // Re-classify the merged points
+    List<edgeStatus> edgeStat(edges().size());
+    forAll(edgeStat, edgeI)
+    {
+        edgeStat[edgeI] = getEdgeStatus(edgeI);
+    }
+
+    List<pointStatus> pointStat(points().size());
+    forAll(pointStat, pointI)
+    {
+        pointStat[pointI] = getPointStatus(pointI);
+    }
+
+    // Re-classify merged points
+    labelList nPoints(nNewPoints, 0);
+    forAll(oldToMerged, oldPointI)
+    {
+        nPoints[oldToMerged[oldPointI]]++;
+    }
+
+    forAll(nPoints, pointI)
+    {
+        if (nPoints[pointI] != 1)
+        {
+            pointStat[pointI] = classifyFeaturePoint(pointI);
+        }
+    }
+
+    labelList sortedToOriginalPoint;
+    setFromStatus
+    (
+        pointStat,
+        edgeStat,
+        sortedToOriginalPoint,
+        edgeMap             // point merging above did not affect edge order
+    );
+    pointMap = UIndirectList<label>(pointMap, sortedToOriginalPoint)();
+
+    return nNewPoints != nOldPoints;
+}
+
+
+void Foam::extendedEdgeMesh::writeObj(const fileName& prefix) const
 {
     Info<< nl << "Writing extendedEdgeMesh components to " << prefix
         << endl;
 
     edgeMesh::write(prefix + "_edgeMesh.obj");
 
-    OBJstream convexFtPtStr(prefix + "_convexFeaturePts.obj");
-    Info<< "Writing convex feature points to " << convexFtPtStr.name() << endl;
-
-    for(label i = 0; i < concaveStart_; i++)
     {
-        convexFtPtStr.write(points()[i]);
-    }
+        OBJstream convexFtPtStr(prefix + "_convexFeaturePts.obj");
+        Info<< "Writing " << concaveStart_
+            << " convex feature points to " << convexFtPtStr.name() << endl;
 
-    OBJstream concaveFtPtStr(prefix + "_concaveFeaturePts.obj");
-    Info<< "Writing concave feature points to "
-        << concaveFtPtStr.name() << endl;
-
-    for(label i = concaveStart_; i < mixedStart_; i++)
-    {
-        convexFtPtStr.write(points()[i]);
-    }
-
-    OBJstream mixedFtPtStr(prefix + "_mixedFeaturePts.obj");
-    Info<< "Writing mixed feature points to " << mixedFtPtStr.name() << endl;
-
-    for(label i = mixedStart_; i < nonFeatureStart_; i++)
-    {
-        mixedFtPtStr.write(points()[i]);
-    }
-
-    OBJstream mixedFtPtStructureStr(prefix + "_mixedFeaturePtsStructure.obj");
-    Info<< "Writing mixed feature point structure to "
-        << mixedFtPtStructureStr.name() << endl;
-
-    for(label i = mixedStart_; i < nonFeatureStart_; i++)
-    {
-        const labelList& ptEds = pointEdges()[i];
-
-        forAll(ptEds, j)
+        for(label i = 0; i < concaveStart_; i++)
         {
-            const edge& e = edges()[ptEds[j]];
-            mixedFtPtStructureStr.write
-            (
-                linePointRef(points()[e[0]],
-                points()[e[1]])
-            );
+            convexFtPtStr.write(points()[i]);
         }
     }
 
-    OBJstream externalStr(prefix + "_externalEdges.obj");
-    Info<< "Writing external edges to " << externalStr.name() << endl;
-
-    for (label i = externalStart_; i < internalStart_; i++)
     {
-        const edge& e = edges()[i];
-        externalStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+        OBJstream concaveFtPtStr(prefix + "_concaveFeaturePts.obj");
+        Info<< "Writing " << mixedStart_-concaveStart_
+            << " concave feature points to "
+            << concaveFtPtStr.name() << endl;
+
+        for(label i = concaveStart_; i < mixedStart_; i++)
+        {
+            concaveFtPtStr.write(points()[i]);
+        }
     }
 
-    OBJstream internalStr(prefix + "_internalEdges.obj");
-    Info<< "Writing internal edges to " << internalStr.name() << endl;
-
-    for (label i = internalStart_; i < flatStart_; i++)
     {
-        const edge& e = edges()[i];
-        internalStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+        OBJstream mixedFtPtStr(prefix + "_mixedFeaturePts.obj");
+        Info<< "Writing " << nonFeatureStart_-mixedStart_
+            << " mixed feature points to " << mixedFtPtStr.name() << endl;
+
+        for(label i = mixedStart_; i < nonFeatureStart_; i++)
+        {
+            mixedFtPtStr.write(points()[i]);
+        }
     }
 
-    OBJstream flatStr(prefix + "_flatEdges.obj");
-    Info<< "Writing flat edges to " << flatStr.name() << endl;
-
-    for (label i = flatStart_; i < openStart_; i++)
     {
-        const edge& e = edges()[i];
-        flatStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+        OBJstream mixedFtPtStructureStr(prefix+"_mixedFeaturePtsStructure.obj");
+        Info<< "Writing "
+            << nonFeatureStart_-mixedStart_
+            << " mixed feature point structure to "
+            << mixedFtPtStructureStr.name() << endl;
+
+        for(label i = mixedStart_; i < nonFeatureStart_; i++)
+        {
+            const labelList& ptEds = pointEdges()[i];
+
+            forAll(ptEds, j)
+            {
+                const edge& e = edges()[ptEds[j]];
+                mixedFtPtStructureStr.write
+                (
+                    linePointRef(points()[e[0]],
+                    points()[e[1]])
+                );
+            }
+        }
     }
 
-    OBJstream openStr(prefix + "_openEdges.obj");
-    Info<< "Writing open edges to " << openStr.name() << endl;
-
-    for (label i = openStart_; i < multipleStart_; i++)
     {
-        const edge& e = edges()[i];
-        openStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+        OBJstream externalStr(prefix + "_externalEdges.obj");
+        Info<< "Writing " << internalStart_-externalStart_
+            << " external edges to " << externalStr.name() << endl;
+
+        for (label i = externalStart_; i < internalStart_; i++)
+        {
+            const edge& e = edges()[i];
+            externalStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+        }
     }
 
-    OBJstream multipleStr(prefix + "_multipleEdges.obj");
-    Info<< "Writing multiple edges to " << multipleStr.name() << endl;
-
-    for (label i = multipleStart_; i < edges().size(); i++)
     {
-        const edge& e = edges()[i];
-        multipleStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+        OBJstream internalStr(prefix + "_internalEdges.obj");
+        Info<< "Writing " << flatStart_-internalStart_
+            << " internal edges to " << internalStr.name() << endl;
+
+        for (label i = internalStart_; i < flatStart_; i++)
+        {
+            const edge& e = edges()[i];
+            internalStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+        }
     }
 
-    OBJstream regionStr(prefix + "_regionEdges.obj");
-    Info<< "Writing region edges to " << regionStr.name() << endl;
-
-    forAll(regionEdges_, i)
     {
-        const edge& e = edges()[regionEdges_[i]];
-        regionStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+        OBJstream flatStr(prefix + "_flatEdges.obj");
+        Info<< "Writing " << openStart_-flatStart_
+            << " flat edges to " << flatStr.name() << endl;
+
+        for (label i = flatStart_; i < openStart_; i++)
+        {
+            const edge& e = edges()[i];
+            flatStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+        }
     }
 
-    OBJstream edgeDirsStr(prefix + "_edgeDirections.obj");
-    Info<< "Writing edge directions to " << edgeDirsStr.name() << endl;
-
-    forAll(edgeDirections_, i)
     {
-        const vector& eVec = edgeDirections_[i];
-        const edge& e = edges()[i];
+        OBJstream openStr(prefix + "_openEdges.obj");
+        Info<< "Writing " << multipleStart_-openStart_
+            << " open edges to " << openStr.name() << endl;
 
-        edgeDirsStr.write
-        (
-            linePointRef(points()[e.start()], eVec + points()[e.start()])
-        );
+        for (label i = openStart_; i < multipleStart_; i++)
+        {
+            const edge& e = edges()[i];
+            openStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+        }
+    }
+
+    {
+        OBJstream multipleStr(prefix + "_multipleEdges.obj");
+        Info<< "Writing " << edges().size()-multipleStart_
+            << " multiple edges to " << multipleStr.name() << endl;
+
+        for (label i = multipleStart_; i < edges().size(); i++)
+        {
+            const edge& e = edges()[i];
+            multipleStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+        }
+    }
+
+    {
+        OBJstream regionStr(prefix + "_regionEdges.obj");
+        Info<< "Writing " << regionEdges_.size()
+            << " region edges to " << regionStr.name() << endl;
+
+        forAll(regionEdges_, i)
+        {
+            const edge& e = edges()[regionEdges_[i]];
+            regionStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+        }
+    }
+
+    {
+        OBJstream edgeDirsStr(prefix + "_edgeDirections.obj");
+        Info<< "Writing " << edgeDirections_.size()
+            << " edge directions to " << edgeDirsStr.name() << endl;
+
+        forAll(edgeDirections_, i)
+        {
+            const vector& eVec = edgeDirections_[i];
+            const edge& e = edges()[i];
+
+            edgeDirsStr.write
+            (
+                linePointRef(points()[e.start()], eVec + points()[e.start()])
+            );
+        }
     }
 }
 
@@ -1506,6 +2163,228 @@ void Foam::extendedEdgeMesh::writeStats(Ostream& os) const
 }
 
 
+Foam::extendedEdgeMesh::edgeStatus
+Foam::extendedEdgeMesh::classifyEdge
+(
+    const List<vector>& norms,
+    const labelList& edNorms,
+    const vector& fC0tofC1
+)
+{
+    label nEdNorms = edNorms.size();
+
+    if (nEdNorms == 1)
+    {
+        return OPEN;
+    }
+    else if (nEdNorms == 2)
+    {
+        const vector& n0(norms[edNorms[0]]);
+        const vector& n1(norms[edNorms[1]]);
+
+        if ((n0 & n1) > cosNormalAngleTol_)
+        {
+            return FLAT;
+        }
+        else if ((fC0tofC1 & n0) > 0.0)
+        {
+            return INTERNAL;
+        }
+        else
+        {
+            return EXTERNAL;
+        }
+    }
+    else if (nEdNorms > 2)
+    {
+        return MULTIPLE;
+    }
+    else
+    {
+        // There is a problem - the edge has no normals
+        return NONE;
+    }
+}
+
+
+void Foam::extendedEdgeMesh::sortedOrder
+(
+    const List<extendedEdgeMesh::pointStatus>& pointStat,
+    const List<extendedEdgeMesh::edgeStatus>& edgeStat,
+    labelList& sortedToOriginalPoint,
+    labelList& sortedToOriginalEdge,
+
+    label& pointConcaveStart,
+    label& pointMixedStart,
+    label& pointNonFeatStart,
+
+    label& edgeInternalStart,
+    label& edgeFlatStart,
+    label& edgeOpenStart,
+    label& edgeMultipleStart
+)
+{
+    sortedToOriginalPoint.setSize(pointStat.size());
+    sortedToOriginalPoint = -1;
+
+    sortedToOriginalEdge.setSize(edgeStat.size());
+    sortedToOriginalEdge = -1;
+
+
+    // Order edges
+    // ~~~~~~~~~~~
+
+    label nConvex = 0;
+    label nConcave = 0;
+    label nMixed = 0;
+    label nNonFeat = 0;
+
+    forAll(pointStat, pointI)
+    {
+        switch (pointStat[pointI])
+        {
+            case extendedEdgeMesh::CONVEX:
+                nConvex++;
+            break;
+
+            case extendedEdgeMesh::CONCAVE:
+                nConcave++;
+            break;
+
+            case extendedEdgeMesh::MIXED:
+                nMixed++;
+            break;
+
+            case extendedEdgeMesh::NONFEATURE:
+                nNonFeat++;
+            break;
+
+            default:
+                FatalErrorIn("order(..)") << "Problem" << exit(FatalError);
+            break;
+        }
+    }
+
+    label convexStart = 0;
+    label concaveStart = nConvex;
+    label mixedStart = concaveStart+nConcave;
+    label nonFeatStart = mixedStart+nMixed;
+
+
+    // Copy to parameters
+    pointConcaveStart = concaveStart;
+    pointMixedStart = mixedStart;
+    pointNonFeatStart = nonFeatStart;
+
+    forAll(pointStat, pointI)
+    {
+        switch (pointStat[pointI])
+        {
+            case extendedEdgeMesh::CONVEX:
+                sortedToOriginalPoint[convexStart++] = pointI;
+            break;
+
+            case extendedEdgeMesh::CONCAVE:
+                sortedToOriginalPoint[concaveStart++] = pointI;
+            break;
+
+            case extendedEdgeMesh::MIXED:
+                sortedToOriginalPoint[mixedStart++] = pointI;
+            break;
+
+            case extendedEdgeMesh::NONFEATURE:
+                sortedToOriginalPoint[nonFeatStart++] = pointI;
+            break;
+        }
+    }
+
+
+    // Order edges
+    // ~~~~~~~~~~~
+
+    label nExternal = 0;
+    label nInternal = 0;
+    label nFlat = 0;
+    label nOpen = 0;
+    label nMultiple = 0;
+
+    forAll(edgeStat, edgeI)
+    {
+        switch (edgeStat[edgeI])
+        {
+            case extendedEdgeMesh::EXTERNAL:
+                nExternal++;
+            break;
+
+            case extendedEdgeMesh::INTERNAL:
+                nInternal++;
+            break;
+
+            case extendedEdgeMesh::FLAT:
+                nFlat++;
+            break;
+
+            case extendedEdgeMesh::OPEN:
+                nOpen++;
+            break;
+
+            case extendedEdgeMesh::MULTIPLE:
+                nMultiple++;
+            break;
+
+            case extendedEdgeMesh::NONE:
+            default:
+                FatalErrorIn("order(..)") << "Problem" << exit(FatalError);
+            break;
+        }
+    }
+
+    label externalStart = 0;
+    label internalStart = nExternal;
+    label flatStart = internalStart + nInternal;
+    label openStart = flatStart + nFlat;
+    label multipleStart = openStart + nOpen;
+
+
+    // Copy to parameters
+    edgeInternalStart = internalStart;
+    edgeFlatStart = flatStart;
+    edgeOpenStart = openStart;
+    edgeMultipleStart = multipleStart;
+
+    forAll(edgeStat, edgeI)
+    {
+        switch (edgeStat[edgeI])
+        {
+            case extendedEdgeMesh::EXTERNAL:
+                sortedToOriginalEdge[externalStart++] = edgeI;
+            break;
+
+            case extendedEdgeMesh::INTERNAL:
+                sortedToOriginalEdge[internalStart++] = edgeI;
+            break;
+
+            case extendedEdgeMesh::FLAT:
+                sortedToOriginalEdge[flatStart++] = edgeI;
+            break;
+
+            case extendedEdgeMesh::OPEN:
+                sortedToOriginalEdge[openStart++] = edgeI;
+            break;
+
+            case extendedEdgeMesh::MULTIPLE:
+                sortedToOriginalEdge[multipleStart++] = edgeI;
+            break;
+
+            case extendedEdgeMesh::NONE:
+            default:
+                FatalErrorIn("order(..)") << "Problem" << exit(FatalError);
+            break;
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 Foam::Istream& Foam::operator>>
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H
index 58e2ba8cab..1dfb921e48 100644
--- a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -68,7 +68,7 @@ namespace Foam
 {
 
 class surfaceFeatures;
-class objectRegistry;
+class searchableSurface;
 
 // Forward declaration of friend functions and operators
 class extendedEdgeMesh;
@@ -138,7 +138,7 @@ protected:
         static label externalStart_;
 
 
-    // Private data
+    // Protected data
 
         //- Index of the start of the concave feature points
         label concaveStart_;
@@ -200,12 +200,33 @@ protected:
         mutable PtrList<indexedOctree<treeDataEdge> > edgeTreesByType_;
 
 
-    // Private Member Functions
+    // Protected Member Functions
 
-        //- Classify the type of feature point.  Requires valid stored member
+        //- Classify the type of feature point. Requires valid stored member
         //  data for edges and normals.
         pointStatus classifyFeaturePoint(label ptI) const;
 
+        //- Cut edges with surface. Return map from cut points&edges back
+        //  to original
+        void cut
+        (
+            const searchableSurface&,
+            labelList& pMap,
+            labelList& eMap,
+            labelList& pointsFromEdge,  // new points created by cutting
+            labelList& oldEdge,         // the orginal edge
+            labelList& surfTri          // the surface triangle index
+        );
+
+        //- Remove outside/inside edges. volType denotes which side to keep
+        void select
+        (
+            const searchableSurface& surf,
+            const volumeType volType,
+            labelList& pMap,
+            labelList& eMap
+        );
+
         template<class Patch>
         void sortPointsAndEdges
         (
@@ -262,9 +283,8 @@ public:
             const Xfer<edgeList>&
         );
 
-        //- Construct given a surface with selected edges,point
-        //  (surfaceFeatures), an objectRegistry and a
-        //  fileName to write to.
+        //- Construct given a surface with selected edges,points
+        //  (surfaceFeatures)
         //  Extracts, classifies and reorders the data from surfaceFeatures.
         extendedEdgeMesh
         (
@@ -501,6 +521,43 @@ public:
             //  etc.
             void flipNormals();
 
+            //- Update with derived geometry
+            void autoMap
+            (
+                const pointField&,
+                const edgeList&,
+                const labelList& pointMap,
+                const labelList& edgeMap
+            );
+
+            //- Trim to surface. Keep volType side. Return map from current back
+            //  to original points (-1 for newly introduced points), edges
+            void trim
+            (
+                const searchableSurface& surf,
+                const volumeType volType,
+                labelList& pointMap,
+                labelList& edgeMap
+            );
+
+            //- Order according to point and edge status
+            void setFromStatus
+            (
+                const List<extendedEdgeMesh::pointStatus>& pointStat,
+                const List<extendedEdgeMesh::edgeStatus>& edgeStat,
+                labelList& sortedToOriginalPoint,
+                labelList& sortedToOriginalEdge
+            );
+
+            //- Geometric merge points. Returns true if any points merged.
+            //  Return maps from new back to original points/edges.
+            bool mergePointsAndSort
+            (
+                const scalar mergeDist,
+                labelList& pointMap,
+                labelList& edgeMap
+            );
+
 
         // Read
 
@@ -532,6 +589,24 @@ public:
             const vector& fC0tofC1
         );
 
+        //- Determine the ordering
+        static void sortedOrder
+        (
+            const List<extendedEdgeMesh::pointStatus>& pointStat,
+            const List<extendedEdgeMesh::edgeStatus>& edgeStat,
+            labelList& sortedToOriginalPoint,
+            labelList& sortedToOriginalEdge,
+
+            label& pointConcaveStart,
+            label& pointMixedStart,
+            label& pointNonFeatStart,
+
+            label& edgeInternalStart,
+            label& edgeFlatStart,
+            label& edgeOpenStart,
+            label& edgeMultipleStart
+        );
+
 
         // Ostream Operator