diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index 5dd7671b948c6286157ef84b8f138b1fde1c6c2a..0d7c929cf68d42d826604db95a1fb87ff7cea8dd 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -472,6 +472,33 @@ void unmarkBaffles
 }
 
 
+void writeStats(const extendedFeatureEdgeMesh& fem, Ostream& os)
+{
+    os  << "    points : " << fem.points().size() << nl
+        << "    of which" << nl
+        << "        convex             : "
+        << fem.concaveStart() << nl
+        << "        concave            : "
+        << (fem.mixedStart()-fem.concaveStart()) << nl
+        << "        mixed              : "
+        << (fem.nonFeatureStart()-fem.mixedStart()) << nl
+        << "        non-feature        : "
+        << (fem.points().size()-fem.nonFeatureStart()) << nl
+        << "    edges  : " << fem.edges().size() << nl
+        << "    of which" << nl
+        << "        external edges     : "
+        << fem.internalStart() << nl
+        << "        internal edges     : "
+        << (fem.flatStart()- fem.internalStart()) << nl
+        << "        flat edges         : "
+        << (fem.openStart()- fem.flatStart()) << nl
+        << "        open edges         : "
+        << (fem.multipleStart()- fem.openStart()) << nl
+        << "        multiply connected : "
+        << (fem.edges().size()- fem.multipleStart()) << nl;
+}
+
+
 // Main program:
 
 int main(int argc, char *argv[])
@@ -746,16 +773,6 @@ int main(int argc, char *argv[])
         //    newSet.writeObj("final");
         //}
 
-        Info<< nl
-            << "Final feature set after trimming and subsetting:" << nl
-            << "    feature points : " << newSet.featurePoints().size() << nl
-            << "    feature edges  : " << newSet.featureEdges().size() << nl
-            << "    of which" << nl
-            << "        region edges   : " << newSet.nRegionEdges() << nl
-            << "        external edges : " << newSet.nExternalEdges() << nl
-            << "        internal edges : " << newSet.nInternalEdges() << nl
-            << endl;
-
         // Extracting and writing a extendedFeatureEdgeMesh
         extendedFeatureEdgeMesh feMesh
         (
@@ -764,6 +781,36 @@ int main(int argc, char *argv[])
             sFeatFileName + ".extendedFeatureEdgeMesh"
         );
 
+
+        if (surfaceDict.isDict("addFeatures"))
+        {
+            const word addFeName = surfaceDict.subDict("addFeatures")["name"];
+            Info<< "Adding (without merging) features from " << addFeName
+                << nl << endl;
+
+            extendedFeatureEdgeMesh addFeMesh
+            (
+                IOobject
+                (
+                    addFeName,
+                    runTime.time().constant(),
+                    "extendedFeatureEdgeMesh",
+                    runTime.time(),
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                )
+            );
+            Info<< "Read " << addFeMesh.name() << nl;
+            writeStats(addFeMesh, Info);
+
+            feMesh.add(addFeMesh);
+        }
+
+
+        Info<< nl
+            << "Final feature set:" << nl;
+        writeStats(feMesh, Info);
+
         Info<< nl << "Writing extendedFeatureEdgeMesh to "
             << feMesh.objectPath() << endl;
 
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
index e9741543bd59330d94fa79d1ee37db2169cabd1e..5e87e98a2b76731a1f950b04f158eb8f5cba8fcb 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
@@ -70,6 +70,12 @@ surface2.nas
         manifoldEdges       no;
     }
 
+    addFeatures
+    {
+        // Add (without merging) another extendedFeatureEdgeMesh
+        name                axZ.extendedFeatureEdgeMesh;
+    }
+
     // Output the curvature of the surface
     curvature               no;
 
diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
index aa046d6ec6028aefd0b407f07824491ac823e33e..7899c7743a35a44978058bb9ef91f33bba5ceeca 100644
--- a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
+++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
@@ -33,6 +33,7 @@ License
 #include "OFstream.H"
 #include "IFstream.H"
 #include "unitConversion.H"
+#include "DynamicField.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1059,12 +1060,264 @@ Foam::extendedFeatureEdgeMesh::edgeTreesByType() const
 }
 
 
+void Foam::extendedFeatureEdgeMesh::add(const extendedFeatureEdgeMesh& fem)
+{
+    // Points
+    // ~~~~~~
+
+    // From current points into combined points
+    labelList reversePointMap(points().size());
+    labelList reverseFemPointMap(fem.points().size());
+
+    label newPointI = 0;
+    for (label i = 0; i < concaveStart(); i++)
+    {
+        reversePointMap[i] = newPointI++;
+    }
+    for (label i = 0; i < fem.concaveStart(); i++)
+    {
+        reverseFemPointMap[i] = newPointI++;
+    }
+
+    // Concave
+    label newConcaveStart = newPointI;
+    for (label i = concaveStart(); i < mixedStart(); i++)
+    {
+        reversePointMap[i] = newPointI++;
+    }
+    for (label i = fem.concaveStart(); i < fem.mixedStart(); i++)
+    {
+        reverseFemPointMap[i] = newPointI++;
+    }
+
+    // Mixed
+    label newMixedStart = newPointI;
+    for (label i = mixedStart(); i < nonFeatureStart(); i++)
+    {
+        reversePointMap[i] = newPointI++;
+    }
+    for (label i = fem.mixedStart(); i < fem.nonFeatureStart(); i++)
+    {
+        reverseFemPointMap[i] = newPointI++;
+    }
+
+    // Non-feature
+    label newNonFeatureStart = newPointI;
+    for (label i = nonFeatureStart(); i < points().size(); i++)
+    {
+        reversePointMap[i] = newPointI++;
+    }
+    for (label i = fem.nonFeatureStart(); i < fem.points().size(); i++)
+    {
+        reverseFemPointMap[i] = newPointI++;
+    }
+
+
+//Pout<< "points:" << points().size() << endl;
+//Pout<< "reversePointMap:" << reversePointMap << endl;
+//
+//Pout<< "fem.points:" << fem.points().size() << endl;
+//Pout<< "reverseFemPointMap:" << reverseFemPointMap << endl;
+
+
+
+    pointField newPoints(newPointI);
+    newPoints.rmap(points(), reversePointMap);
+    newPoints.rmap(fem.points(), reverseFemPointMap);
+
+
+    // Edges
+    // ~~~~~
+
+    // From current edges into combined edges
+    labelList reverseEdgeMap(edges().size());
+    labelList reverseFemEdgeMap(fem.edges().size());
+
+    // External
+    label newEdgeI = 0;
+    for (label i = 0; i < internalStart(); i++)
+    {
+        reverseEdgeMap[i] = newEdgeI++;
+    }
+    for (label i = 0; i < fem.internalStart(); i++)
+    {
+        reverseFemEdgeMap[i] = newEdgeI++;
+    }
+
+    // Internal
+    label newInternalStart = newEdgeI;
+    for (label i = internalStart(); i < flatStart(); i++)
+    {
+        reverseEdgeMap[i] = newEdgeI++;
+    }
+    for (label i = fem.internalStart(); i < fem.flatStart(); i++)
+    {
+        reverseFemEdgeMap[i] = newEdgeI++;
+    }
+
+    // Flat
+    label newFlatStart = newEdgeI;
+    for (label i = flatStart(); i < openStart(); i++)
+    {
+        reverseEdgeMap[i] = newEdgeI++;
+    }
+    for (label i = fem.flatStart(); i < fem.openStart(); i++)
+    {
+        reverseFemEdgeMap[i] = newEdgeI++;
+    }
+
+    // Open
+    label newOpenStart = newEdgeI;
+    for (label i = openStart(); i < multipleStart(); i++)
+    {
+        reverseEdgeMap[i] = newEdgeI++;
+    }
+    for (label i = fem.openStart(); i < fem.multipleStart(); i++)
+    {
+        reverseFemEdgeMap[i] = newEdgeI++;
+    }
+
+    // Multiple
+    label newMultipleStart = newEdgeI;
+    for (label i = multipleStart(); i < edges().size(); i++)
+    {
+        reverseEdgeMap[i] = newEdgeI++;
+    }
+    for (label i = fem.multipleStart(); i < fem.edges().size(); i++)
+    {
+        reverseFemEdgeMap[i] = newEdgeI++;
+    }
+
+    edgeList newEdges(newEdgeI);
+    UIndirectList<edge>(newEdges, reverseEdgeMap) = edges();
+    UIndirectList<edge>(newEdges, reverseFemEdgeMap) = fem.edges();
+
+    pointField newEdgeDirections(newEdgeI);
+    newEdgeDirections.rmap(edgeDirections(), reverseEdgeMap);
+    newEdgeDirections.rmap(fem.edgeDirections(), reverseFemEdgeMap);
+
+
+
+
+    // Normals
+    // ~~~~~~~
+
+    // Combine normals
+    DynamicField<point> newNormals(normals().size()+fem.normals().size());
+    newNormals.append(normals());
+    newNormals.append(fem.normals());
+
+
+    // Combine and re-index into newNormals
+    labelListList newEdgeNormals(edgeNormals().size()+fem.edgeNormals().size());
+    UIndirectList<labelList>(newEdgeNormals, reverseEdgeMap) =
+        edgeNormals();
+    UIndirectList<labelList>(newEdgeNormals, reverseFemEdgeMap) =
+        fem.edgeNormals();
+    forAll(reverseFemEdgeMap, i)
+    {
+        label mapI = reverseFemEdgeMap[i];
+        labelList& en = newEdgeNormals[mapI];
+        forAll(en, j)
+        {
+            en[j] += edgeNormals().size();
+        }
+    }
+
+//Pout<< "fem.featurePointNormals().size():" << fem.featurePointNormals().size()
+//    << endl;
+//Pout<< "fem.nonFeatureStart():" << fem.nonFeatureStart() << endl;
+//Pout<< "reverseFemPointMap:" << reverseFemPointMap.size() << endl;
+
+
+    // Combine and re-index into newFeaturePointNormals
+    labelListList newFeaturePointNormals
+    (
+       featurePointNormals().size()
+     + fem.featurePointNormals().size()
+    );
+
+//Pout<< "newFeaturePointNormals:" << newFeaturePointNormals.size() << endl;
+//Pout<< "Doing featurePointNormals." << endl;
+
+    // Note: featurePointNormals only go up to nonFeatureStart
+    UIndirectList<labelList>
+    (
+        newFeaturePointNormals,
+        SubList<label>(reversePointMap, nonFeatureStart())
+    ) = featurePointNormals();
+
+
+//Pout<< "Doing fem.featurePointNormals()." << endl;
+    UIndirectList<labelList>
+    (
+        newFeaturePointNormals,
+        SubList<label>(reverseFemPointMap, fem.nonFeatureStart())
+    ) = fem.featurePointNormals();
+    forAll(fem.featurePointNormals(), i)
+    {
+        label mapI = reverseFemPointMap[i];
+        labelList& fn = newFeaturePointNormals[mapI];
+        forAll(fn, j)
+        {
+            fn[j] += featurePointNormals().size();
+        }
+    }
+
+
+    // Combine regionEdges
+    DynamicList<label> newRegionEdges
+    (
+        regionEdges().size()
+      + fem.regionEdges().size()
+    );
+    forAll(regionEdges(), i)
+    {
+        newRegionEdges.append(reverseEdgeMap[regionEdges()[i]]);
+    }
+    forAll(fem.regionEdges(), i)
+    {
+        newRegionEdges.append(reverseFemEdgeMap[fem.regionEdges()[i]]);
+    }
+
+
+    // Assign
+    // ~~~~~~
+
+    // Transfer
+    concaveStart_ = newConcaveStart;
+    mixedStart_ = newMixedStart;
+    nonFeatureStart_ = newNonFeatureStart;
+
+    // Reset points and edges
+    reset(xferMove(newPoints), newEdges.xfer());
+
+    // Transfer
+    internalStart_ = newInternalStart;
+    flatStart_ = newFlatStart;
+    openStart_ = newOpenStart;
+    multipleStart_ = newMultipleStart;
+
+    edgeDirections_.transfer(newEdgeDirections);
+
+    normals_.transfer(newNormals);
+    edgeNormals_.transfer(newEdgeNormals);
+    featurePointNormals_.transfer(newFeaturePointNormals);
+
+    regionEdges_.transfer(newRegionEdges);
+
+    pointTree_.clear();
+    edgeTree_.clear();
+    edgeTreesByType_.clear();
+}
+
+
 void Foam::extendedFeatureEdgeMesh::writeObj
 (
     const fileName& prefix
 ) const
 {
-    Pout<< nl << "Writing extendedFeatureEdgeMesh components to " << prefix
+    Info<< nl << "Writing extendedFeatureEdgeMesh components to " << prefix
         << endl;
 
     label verti = 0;
@@ -1072,7 +1325,7 @@ void Foam::extendedFeatureEdgeMesh::writeObj
     edgeMesh::write(prefix + "_edgeMesh.obj");
 
     OFstream convexFtPtStr(prefix + "_convexFeaturePts.obj");
-    Pout<< "Writing convex feature points to " << convexFtPtStr.name() << endl;
+    Info<< "Writing convex feature points to " << convexFtPtStr.name() << endl;
 
     for(label i = 0; i < concaveStart_; i++)
     {
@@ -1080,7 +1333,7 @@ void Foam::extendedFeatureEdgeMesh::writeObj
     }
 
     OFstream concaveFtPtStr(prefix + "_concaveFeaturePts.obj");
-    Pout<< "Writing concave feature points to "
+    Info<< "Writing concave feature points to "
         << concaveFtPtStr.name() << endl;
 
     for(label i = concaveStart_; i < mixedStart_; i++)
@@ -1089,7 +1342,7 @@ void Foam::extendedFeatureEdgeMesh::writeObj
     }
 
     OFstream mixedFtPtStr(prefix + "_mixedFeaturePts.obj");
-    Pout<< "Writing mixed feature points to " << mixedFtPtStr.name() << endl;
+    Info<< "Writing mixed feature points to " << mixedFtPtStr.name() << endl;
 
     for(label i = mixedStart_; i < nonFeatureStart_; i++)
     {
@@ -1097,7 +1350,7 @@ void Foam::extendedFeatureEdgeMesh::writeObj
     }
 
     OFstream mixedFtPtStructureStr(prefix + "_mixedFeaturePtsStructure.obj");
-    Pout<< "Writing mixed feature point structure to "
+    Info<< "Writing mixed feature point structure to "
         << mixedFtPtStructureStr.name() << endl;
 
     verti = 0;
@@ -1116,7 +1369,7 @@ void Foam::extendedFeatureEdgeMesh::writeObj
     }
 
     OFstream externalStr(prefix + "_externalEdges.obj");
-    Pout<< "Writing external edges to " << externalStr.name() << endl;
+    Info<< "Writing external edges to " << externalStr.name() << endl;
 
     verti = 0;
     for (label i = externalStart_; i < internalStart_; i++)
@@ -1129,7 +1382,7 @@ void Foam::extendedFeatureEdgeMesh::writeObj
     }
 
     OFstream internalStr(prefix + "_internalEdges.obj");
-    Pout<< "Writing internal edges to " << internalStr.name() << endl;
+    Info<< "Writing internal edges to " << internalStr.name() << endl;
 
     verti = 0;
     for (label i = internalStart_; i < flatStart_; i++)
@@ -1142,7 +1395,7 @@ void Foam::extendedFeatureEdgeMesh::writeObj
     }
 
     OFstream flatStr(prefix + "_flatEdges.obj");
-    Pout<< "Writing flat edges to " << flatStr.name() << endl;
+    Info<< "Writing flat edges to " << flatStr.name() << endl;
 
     verti = 0;
     for (label i = flatStart_; i < openStart_; i++)
@@ -1155,7 +1408,7 @@ void Foam::extendedFeatureEdgeMesh::writeObj
     }
 
     OFstream openStr(prefix + "_openEdges.obj");
-    Pout<< "Writing open edges to " << openStr.name() << endl;
+    Info<< "Writing open edges to " << openStr.name() << endl;
 
     verti = 0;
     for (label i = openStart_; i < multipleStart_; i++)
@@ -1168,7 +1421,7 @@ void Foam::extendedFeatureEdgeMesh::writeObj
     }
 
     OFstream multipleStr(prefix + "_multipleEdges.obj");
-    Pout<< "Writing multiple edges to " << multipleStr.name() << endl;
+    Info<< "Writing multiple edges to " << multipleStr.name() << endl;
 
     verti = 0;
     for (label i = multipleStart_; i < edges().size(); i++)
@@ -1181,7 +1434,7 @@ void Foam::extendedFeatureEdgeMesh::writeObj
     }
 
     OFstream regionStr(prefix + "_regionEdges.obj");
-    Pout<< "Writing region edges to " << regionStr.name() << endl;
+    Info<< "Writing region edges to " << regionStr.name() << endl;
 
     verti = 0;
     forAll(regionEdges_, i)
@@ -1197,23 +1450,26 @@ void Foam::extendedFeatureEdgeMesh::writeObj
 
 bool Foam::extendedFeatureEdgeMesh::writeData(Ostream& os) const
 {
-    os  << "// points, edges, concaveStart, mixedStart, nonFeatureStart, " << nl
-        << "// internalStart, flatStart, openStart, multipleStart, " << nl
-        << "// normals, edgeNormals, featurePointNormals, regionEdges" << nl
-        << endl;
-
-    os  << points() << nl
+    os  << "// points" << nl
+        << points() << nl
+        << "// edges" << nl
         << edges() << nl
+        << "// concaveStart mixedStart nonFeatureStart" << nl
         << concaveStart_ << token::SPACE
         << mixedStart_ << token::SPACE
-        << nonFeatureStart_ << token::SPACE
+        << nonFeatureStart_ << nl
+        << "// internalStart flatStart openStart multipleStart" << nl
         << internalStart_ << token::SPACE
         << flatStart_ << token::SPACE
         << openStart_ << token::SPACE
         << multipleStart_ << nl
+        << "// normals" << nl
         << normals_ << nl
+        << "// edgeNormals" << nl
         << edgeNormals_ << nl
+        << "// featurePointNormals" << nl
         << featurePointNormals_ << nl
+        << "// regionEdges" << nl
         << regionEdges_
         << endl;
 
diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
index 28601a26ddf5906af58925f4725506d64b7e6fea..85dcfcf7dcf76fb537b21dcb5bb0624cde133740 100644
--- a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
+++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,8 +30,8 @@ Description
 
     Feature points are a sorted subset at the start of the overall points list:
         0 .. concaveStart_-1                : convex points (w.r.t normals)
-        concaveStart_-1 .. mixedStart_-1    : concave points
-        mixedStart_ .. nonFeatureStart_     : mixed internal/external points
+        concaveStart_ .. mixedStart_-1      : concave points
+        mixedStart_ .. nonFeatureStart_-1   : mixed internal/external points
         nonFeatureStart_ .. size-1          : non-feature points
 
     Feature edges are the edgeList of the edgeMesh and are sorted:
@@ -151,6 +151,7 @@ private:
         labelListList edgeNormals_;
 
         //- Indices of the normals that are adjacent to the feature points
+        //  (only valid for 0..nonFeatureStart_-1)
         labelListList featurePointNormals_;
 
         //- Feature edges which are on the boundary between regions
@@ -212,9 +213,10 @@ public:
             const Xfer<edgeList>&
         );
 
-        //- Construct (read) given surfaceFeatures, an objectRegistry and a
-        //  fileName to write to.  Extracts, classifies and reorders the data
-        //  from surfaceFeatures.
+        //- Construct given a surface with selected edges,point
+        //  (surfaceFeatures), an objectRegistry and a
+        //  fileName to write to.
+        //  Extracts, classifies and reorders the data from surfaceFeatures.
         extendedFeatureEdgeMesh
         (
             const surfaceFeatures& sFeat,
@@ -381,6 +383,12 @@ public:
             edgeTreesByType() const;
 
 
+        // Edit
+
+            //- Add extendedFeatureEdgeMesh. No filtering of duplicates.
+            void add(const extendedFeatureEdgeMesh&);
+
+
         // Write
 
             //- Write all components of the extendedFeatureEdgeMesh as obj files