ENH: multiple surfaces, self-intersection in surfaceFeatureExtract (issue #450)

- If the dictionary is named 'surfaces', a 'surfaces' entry is mandatory.
  This is a list of wordRe, which is used to load multiple surfaces from
  constant/triSurface directory.

- Other dictionaries may contain a 'surfaces' entry.
  In which case the behaviour is as above (loading multiple surfaces).
  The dictionary name will *NOT* be taken as a surface name itself.

- Regardless of how the surfaces are loaded or features extracted,
  an additional selfIntersection test may be used.

  Eg,

    surfaces
    {
        extractionMethod    extractFromSurface;

        surfaces            (surface1.stl surface2.nas);

        // Generate features from self-intersect
        selfIntersection    true;

        // Base output name (optiona)
        output              surfaces;

        // Tolerance for self-intersect
        planarTolerance     1e-3;

        extractFromSurfaceCoeffs
        {
            includedAngle   120;

            // Do not mark region edges
            geometricTestOnly       yes;
        }
    }

applications/utilities/surface/surfaceFeatureExtract/Make/options

@@ -1,9 +1,11 @@
 EXE_INC = \
+    -IextractionMethod/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \
+    -lsurfaceFeatureExtract \
     -lmeshTools \
     -lsampling

applications/utilities/surface/surfaceFeatureExtract/extractionMethod/Make/files

@@ -0,0 +1,6 @@
+method = .
+$(method)/surfaceFeaturesExtraction.C
+$(method)/extractFromSurface.C
+$(method)/extractFromFile.C
+
+LIB = $(FOAM_LIBBIN)/libsurfaceFeatureExtract

applications/utilities/surface/surfaceFeatureExtract/extractionMethod/Make/options

@@ -0,0 +1,11 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/edgeMesh/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude
+
+LIB_LIBS = \
+    -lmeshTools \
+    -ledgeMesh \
+    -ltriSurface

applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromFile.C

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |  Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "extractFromFile.H"
+#include "ListOps.H"
+#include "edgeMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceFeaturesExtraction
+{
+    addNamedToRunTimeSelectionTable
+    (
+        method,
+        extractFromFile,
+        dictionary,
+        extractFromFile
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceFeaturesExtraction::extractFromFile::extractFromFile
+(
+    const dictionary& dict
+)
+:
+    method()
+{
+    const dictionary& coeffDict = dict.subDict("extractFromFileCoeffs");
+
+    coeffDict.lookup("featureEdgeFile") >> featureEdgeFile_;
+    coeffDict.readIfPresent("geometricTestOnly", geometricTestOnly_);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::surfaceFeaturesExtraction::extractFromFile::~extractFromFile()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+Foam::autoPtr<Foam::surfaceFeatures>
+Foam::surfaceFeaturesExtraction::extractFromFile::features
+(
+    const triSurface& surf
+) const
+{
+    edgeMesh eMesh(featureEdgeFile_);
+
+    // Sometimes duplicate edges are present. Remove them.
+    eMesh.mergeEdges();
+
+    Info<< nl << "Reading existing feature edges from file "
+        << featureEdgeFile_ << nl
+        << "Selecting edges based purely on geometric tests: "
+        << geometricTestOnly().asText() << endl;
+
+    return autoPtr<surfaceFeatures>
+    (
+        new surfaceFeatures
+        (
+            surf,
+            eMesh.points(),
+            eMesh.edges(),
+            1e-6,  // mergeTol
+            geometricTestOnly()
+        )
+    );
+}
+
+
+// ************************************************************************* //

applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromFile.H

@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     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/>.
+
+Class
+    Foam::surfaceFeaturesExtraction::extractFromFile
+
+Description
+    Run-time selectable surface feature extraction.
+
+SourceFiles
+    extractionMethod.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfaceFeaturesExtraction_extractFromFile_H
+#define surfaceFeaturesExtraction_extractFromFile_H
+
+#include "surfaceFeaturesExtraction.H"
+#include "dictionary.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceFeaturesExtraction
+{
+
+/*---------------------------------------------------------------------------*\
+         Class surfaceFeaturesExtraction::extractFromFile Declaration
+\*---------------------------------------------------------------------------*/
+
+class extractFromFile
+:
+    public method
+{
+    fileName featureEdgeFile_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from dictionary
+        extractFromFile(const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~extractFromFile();
+
+
+    //- Extracted features from surface
+    virtual autoPtr<surfaceFeatures> features
+    (
+        const triSurface& surf
+    ) const override;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFeaturesExtraction
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromSurface.C

@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |  Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "extractFromSurface.H"
+#include "addToRunTimeSelectionTable.H"
+
+#include "triSurface.H"
+#include "triSurfaceSearch.H"
+#include "scalarField.H"
+#include "edgeIntersections.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceFeaturesExtraction
+{
+    addNamedToRunTimeSelectionTable
+    (
+        method,
+        extractFromSurface,
+        dictionary,
+        extractFromSurface
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceFeaturesExtraction::extractFromSurface::extractFromSurface
+(
+    const dictionary& dict
+)
+:
+    method()
+{
+    const dictionary& coeffDict = dict.subDict("extractFromSurfaceCoeffs");
+
+    coeffDict.lookup("includedAngle") >> includedAngle_;
+    coeffDict.readIfPresent("geometricTestOnly", geometricTestOnly_);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::surfaceFeaturesExtraction::extractFromSurface::~extractFromSurface()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+Foam::autoPtr<Foam::surfaceFeatures>
+Foam::surfaceFeaturesExtraction::extractFromSurface::features
+(
+    const triSurface& surf
+) const
+{
+    Info<< nl << "Constructing feature set from included angle "
+        << includedAngle() << nl
+        << "Selecting edges based purely on geometric tests: "
+        << geometricTestOnly().asText() << endl;
+
+    return autoPtr<surfaceFeatures>
+    (
+        new surfaceFeatures
+        (
+            surf,
+            includedAngle(),
+            0,  // minLen
+            0,  // minElems
+            geometricTestOnly()
+        )
+    );
+}
+
+
+// ************************************************************************* //

applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromSurface.H

@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     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/>.
+
+Class
+    Foam::surfaceFeaturesExtraction::extractFromSurface
+
+Description
+    Run-time selectable surface feature extraction.
+
+SourceFiles
+    extractionMethod.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfaceFeaturesExtraction_extractFromSurface_H
+#define surfaceFeaturesExtraction_extractFromSurface_H
+
+#include "surfaceFeaturesExtraction.H"
+#include "dictionary.H"
+#include "Switch.H"
+#include "triSurface.H"
+#include "edgeIntersections.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceFeaturesExtraction
+{
+
+/*---------------------------------------------------------------------------*\
+         Class surfaceFeaturesExtraction::extractFromSurface Declaration
+\*---------------------------------------------------------------------------*/
+
+class extractFromSurface
+:
+    public method
+{
+
+    bool selfIntersection_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from dictionary
+        extractFromSurface(const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~extractFromSurface();
+
+
+    //- Extracted features from surface
+    virtual autoPtr<surfaceFeatures> features
+    (
+        const triSurface& surf
+    ) const override;
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFeaturesExtraction
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/surface/surfaceFeatureExtract/extractionMethod/surfaceFeaturesExtraction.C

@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |  Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceFeaturesExtraction.H"
+#include "dictionary.H"
+#include "ListOps.H"
+#include "error.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceFeaturesExtraction
+{
+    defineTypeName(method);
+    defineRunTimeSelectionTable
+    (
+        method,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceFeaturesExtraction::method::method()
+:
+    includedAngle_(0),
+    geometricTestOnly_(Switch::NO)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::surfaceFeaturesExtraction::method::~method()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::surfaceFeaturesExtraction::method>
+Foam::surfaceFeaturesExtraction::method::New
+(
+    const dictionary& dict
+)
+{
+    const word methodName = dict.lookup("extractionMethod");
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(methodName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalIOErrorInFunction
+        (
+            dict
+        )   << "Unknown extractionMethod "
+            << methodName << nl << nl
+            << "Valid extraction methods: :" << nl
+            << flatOutput(dictionaryConstructorTablePtr_->sortedToc())
+            << exit(FatalIOError);
+    }
+
+    return autoPtr<method>(cstrIter()(dict));
+}
+
+
+// ************************************************************************* //

applications/utilities/surface/surfaceFeatureExtract/extractionMethod/surfaceFeaturesExtraction.H

@@ -0,0 +1,136 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     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/>.
+
+Class
+    Foam::surfaceFeaturesExtraction::method
+
+Description
+    Run-time selectable surface feature extraction methods.
+
+SourceFiles
+    surfaceFeaturesExtraction.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfaceFeaturesExtraction_method_H
+#define surfaceFeaturesExtraction_method_H
+
+#include "surfaceFeatures.H"
+#include "dictionary.H"
+#include "Switch.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceFeaturesExtraction
+{
+
+/*---------------------------------------------------------------------------*\
+              Class surfaceFeaturesExtraction::method Declaration
+\*---------------------------------------------------------------------------*/
+
+class method
+{
+protected:
+
+    scalar includedAngle_;
+    Switch geometricTestOnly_;
+
+    //- Construct null
+    method();
+
+
+public:
+
+    //- Runtime type information
+    ClassNameNoDebug("method");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        method(const dictionary& dict);
+
+
+    // Declare run-time constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        method,
+        dictionary,
+        (
+            const dictionary& dict
+        ),
+        (dict)
+    );
+
+
+    // Selectors
+
+        //- Select constructed from dictionary
+        static autoPtr<method> New
+        (
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~method();
+
+
+    // Member Functions
+
+        //- The included angle, if set
+        inline scalar includedAngle() const
+        {
+            return includedAngle_;
+        }
+
+        //- Use geometric test only
+        inline Switch geometricTestOnly() const
+        {
+            return geometricTestOnly_;
+        }
+
+        //- Extracted features
+        virtual autoPtr<surfaceFeatures> features
+        (
+            const triSurface& surf
+        ) const = 0;
+
+};
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFeaturesExtraction
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,7 +41,8 @@ Description
 #include "argList.H"
 #include "Time.H"
 #include "triSurface.H"
-#include "surfaceFeatures.H"
+#include "surfaceFeaturesExtraction.H"
+#include "surfaceIntersection.H"
 #include "featureEdgeMesh.H"
 #include "extendedFeatureEdgeMesh.H"
 #include "treeBoundBox.H"
@@ -865,6 +866,61 @@ void writeStats(const extendedFeatureEdgeMesh& fem, Ostream& os)
 }
 
 
+// Read and combine all surfaces into a single one
+
+autoPtr<triSurface> loadSurfaces(Time& runTime, const wordList& surfNames)
+{
+    List<labelledTri> faces;
+    pointField points;
+    label regoff = 0; // region offset
+
+    forAll(surfNames, surfi)
+    {
+        const word& surfName = surfNames[surfi];
+
+        triSurface addsurf(runTime.constantPath()/"triSurface"/surfName);
+
+        List<labelledTri> addfaces(addsurf.xferFaces());
+        List<point> addpoints(addsurf.xferPoints());
+
+        if (surfi)
+        {
+            const label ptoff = points.size(); // point offset
+
+            forAll(addfaces, facei)
+            {
+                labelledTri& f = addfaces[facei];
+                forAll(f, fi)
+                {
+                    f[fi] += ptoff;
+                }
+                f.region() += regoff;
+            }
+
+            faces.append(addfaces);
+            points.append(addpoints);
+        }
+        else
+        {
+            faces.transfer(addfaces);
+            points.transfer(addpoints);
+        }
+
+        regoff += addsurf.patches().size();
+    }
+
+    return autoPtr<triSurface>
+    (
+        new triSurface
+        (
+            faces,
+            geometricSurfacePatchList(),
+            points,
+            true
+        )
+    );
+}
+
 
 int main(int argc, char *argv[])
 {
@@ -873,6 +929,7 @@ int main(int argc, char *argv[])
         "extract and write surface features to file"
     );
     argList::noParallel();
+    argList::noFunctionObjects();
 
     #include "addDictOption.H"
 
@@ -882,6 +939,21 @@ int main(int argc, char *argv[])
     const word dictName("surfaceFeatureExtractDict");
     #include "setSystemRunTimeDictionaryIO.H"
 
+    // Will be using triSurface, so filter according to what is supported
+    fileNameList validSurfaceFiles = readDir
+    (
+        runTime.path()/runTime.constant()/"triSurface",
+        fileName::FILE
+    );
+    inplaceSubsetList
+    (
+        validSurfaceFiles,
+        [](const fileName& f){ return triSurface::canRead(f); }
+    );
+
+    // sort and eliminate duplicates (eg, files with/without .gz)
+    inplaceUniqueSort(validSurfaceFiles);
+
     Info<< "Reading " << dictName << nl << endl;
 
     const IOdictionary dict(dictIO);
@@ -900,132 +972,132 @@ int main(int argc, char *argv[])
             continue;
         }
 
-        const word extractionMethod = surfaceDict.lookup("extractionMethod");
+        autoPtr<surfaceFeaturesExtraction::method> extractPtr =
+            surfaceFeaturesExtraction::method::New
+            (
+                surfaceDict
+            );
 
-        const fileName surfFileName = iter().keyword();
-        const fileName sFeatFileName = surfFileName.lessExt().name();
+        const surfaceFeaturesExtraction::method& extract = extractPtr();
 
-        Info<< "Surface            : " << surfFileName << nl << endl;
+        // Output name, cleansed of extensions
+        const word sFeatFileName =
+            fileName
+            (
+                surfaceDict.lookupOrDefault<word>("output", iter().keyword())
+            ).lessExt();
 
-        const Switch writeVTK =
-            surfaceDict.lookupOrDefault<Switch>("writeVTK", "off");
-        const Switch writeObj =
-            surfaceDict.lookupOrDefault<Switch>("writeObj", "off");
+        wordList surfFileNames;
 
-        const Switch curvature =
-            surfaceDict.lookupOrDefault<Switch>("curvature", "off");
-        const Switch featureProximity =
-            surfaceDict.lookupOrDefault<Switch>("featureProximity", "off");
-        const Switch closeness =
-            surfaceDict.lookupOrDefault<Switch>("closeness", "off");
+        if
+        (
+            iter().keyword() == "surfaces"   // mandatory
+         || surfaceDict.found("surfaces")    // or optional
+        )
+        {
+            wordReList regexs(surfaceDict.lookup("surfaces"));
 
+            labelList selected = findStrings(regexs, validSurfaceFiles);
+            surfFileNames.setSize(selected.size());
+            forAll(selected, i)
+            {
+                surfFileNames[i] = validSurfaceFiles[selected[i]].name();
+            }
+
+            if (surfFileNames.empty())
+            {
+                FatalErrorInFunction
+                    << "No surfaces specified/found for entry: "
+                    << iter().keyword()
+                    << exit(FatalError);
+            }
+        }
+        else
+        {
+            surfFileNames.setSize(1);
+            surfFileNames[0] = iter().keyword();
+
+            const fileName file
+            (
+                runTime.constantPath()/"triSurface"/surfFileNames[0]
+            );
+
+            if (!isFile(file))
+            {
+                FatalErrorInFunction
+                    << "No surface: " << file.name()
+                    << exit(FatalError);
+            }
+        }
+
+        // DebugVar(surfFileNames);
+        // DebugVar(sFeatFileName);
+
+        Info<< "Surfaces           : ";
+        if (surfFileNames.size() == 1)
+        {
+            Info<< surfFileNames[0] << nl;
+        }
+        else
+        {
+            Info<< flatOutput(surfFileNames) << nl;
+        }
+
+        Info<< "Output             : " << sFeatFileName << nl;
+
+        const Switch writeVTK = surfaceDict.lookupOrDefault<Switch>
+        (
+            "writeVTK", Switch::OFF
+        );
+        const Switch writeObj = surfaceDict.lookupOrDefault<Switch>
+        (
+            "writeObj", Switch::OFF
+        );
+
+        const Switch curvature = surfaceDict.lookupOrDefault<Switch>
+        (
+            "curvature", Switch::OFF
+        );
+        const Switch featureProximity = surfaceDict.lookupOrDefault<Switch>
+        (
+            "featureProximity", Switch::OFF
+        );
+        const Switch closeness = surfaceDict.lookupOrDefault<Switch>
+        (
+            "closeness", Switch::OFF
+        );
+
+        Info<< "write VTK: " <<  writeVTK << nl;
 
         Info<< nl << "Feature line extraction is only valid on closed manifold "
             << "surfaces." << endl;
 
-        // Read
-        // ~~~~
 
-        triSurface surf(runTime.constantPath()/"triSurface"/surfFileName);
+        // Read and combine all surfaces into a single one
+
+        autoPtr<triSurface> surfPtr = loadSurfaces(runTime, surfFileNames);
+        triSurface surf = surfPtr();
 
         Info<< "Statistics:" << endl;
         surf.writeStats(Info);
         Info<< endl;
 
-        faceList faces(surf.size());
-
-        forAll(surf, fI)
+        // need plain faces if outputting VTK format
+        faceList faces;
+        if (writeVTK)
         {
-            faces[fI] = surf[fI].triFaceFace();
+            faces.setSize(surf.size());
+            forAll(surf, fi)
+            {
+                faces[fi] = surf[fi].triFaceFace();
+            }
         }
 
 
         // Either construct features from surface & featureAngle or read set.
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        autoPtr<surfaceFeatures> set;
-
-        scalar includedAngle = 0.0;
-
-        if (extractionMethod == "extractFromFile")
-        {
-            const dictionary& extractFromFileDict =
-                surfaceDict.subDict("extractFromFileCoeffs");
-
-            const fileName featureEdgeFile =
-                extractFromFileDict.lookup("featureEdgeFile");
-
-            const Switch geometricTestOnly =
-                extractFromFileDict.lookupOrDefault<Switch>
-                (
-                    "geometricTestOnly",
-                    "no"
-                );
-
-            edgeMesh eMesh(featureEdgeFile);
-
-            // Sometimes duplicate edges are present. Remove them.
-            eMesh.mergeEdges();
-
-            Info<< nl << "Reading existing feature edges from file "
-                << featureEdgeFile << nl
-                << "Selecting edges purely based on geometric tests: "
-                << geometricTestOnly.asText() << endl;
-
-            set.set
-            (
-                new surfaceFeatures
-                (
-                    surf,
-                    eMesh.points(),
-                    eMesh.edges(),
-                    1e-6,
-                    geometricTestOnly
-                )
-            );
-        }
-        else if (extractionMethod == "extractFromSurface")
-        {
-            const dictionary& extractFromSurfaceDict =
-                surfaceDict.subDict("extractFromSurfaceCoeffs");
-
-            includedAngle =
-                readScalar(extractFromSurfaceDict.lookup("includedAngle"));
-
-            const Switch geometricTestOnly =
-                extractFromSurfaceDict.lookupOrDefault<Switch>
-                (
-                    "geometricTestOnly",
-                    "no"
-                );
-
-            Info<< nl
-                << "Constructing feature set from included angle "
-                << includedAngle << nl
-                << "Selecting edges purely based on geometric tests: "
-                << geometricTestOnly.asText() << endl;
-
-            set.set
-            (
-                new surfaceFeatures
-                (
-                    surf,
-                    includedAngle,
-                    0,
-                    0,
-                    geometricTestOnly
-                )
-            );
-        }
-        else
-        {
-            FatalErrorInFunction
-                << "No initial feature set. Provide either one"
-                << " of extractFromFile (to read existing set)" << nl
-                << " or extractFromSurface (to construct new set from angle)"
-                << exit(FatalError);
-        }
-
+        autoPtr<surfaceFeatures> set = extract.features(surf);
 
         // Trim set
         // ~~~~~~~~
@@ -1033,21 +1105,27 @@ int main(int argc, char *argv[])
         if (surfaceDict.isDict("trimFeatures"))
         {
             dictionary trimDict = surfaceDict.subDict("trimFeatures");
-
-            scalar minLen =
+            const scalar minLen =
                 trimDict.lookupOrAddDefault<scalar>("minLen", -GREAT);
 
-            label minElem = trimDict.lookupOrAddDefault<label>("minElem", 0);
+            const label minElem =
+                trimDict.lookupOrAddDefault<label>("minElem", 0);
 
             // Trim away small groups of features
-            if (minElem > 0 || minLen > 0)
+            if (minLen > 0 || minElem > 0)
             {
-                Info<< "Removing features of length < "
-                    << minLen << endl;
-                Info<< "Removing features with number of edges < "
-                    << minElem << endl;
+                if (minLen > 0)
+                {
+                    Info<< "Removing features of length < "
+                        << minLen << endl;
+                }
+                if (minElem > 0)
+                {
+                    Info<< "Removing features with number of edges < "
+                        << minElem << endl;
+                }
 
-                set().trimFeatures(minLen, minElem, includedAngle);
+                set().trimFeatures(minLen, minElem, extract.includedAngle());
             }
         }
 
@@ -1108,7 +1186,7 @@ int main(int argc, char *argv[])
                         (
                             surf,
                             1e-5,   //tol,
-                            includedAngle,
+                            extract.includedAngle(),
                             edgeI
                         );
                     }
@@ -1147,7 +1225,7 @@ int main(int argc, char *argv[])
 
 
         surfaceFeatures newSet(surf);
-        newSet.setFromStatus(edgeStat, includedAngle);
+        newSet.setFromStatus(edgeStat, extract.includedAngle());
 
         Info<< nl
             << "Initial feature set:" << nl
@@ -1215,6 +1293,39 @@ int main(int argc, char *argv[])
         }
 
 
+        if (surfaceDict.lookupOrDefault<bool>("selfIntersection", false))
+        {
+            // TODO: perturb tolerance
+
+            triSurfaceSearch query(surf);
+            surfaceIntersection intersect
+            (
+                query,
+                surfaceDict.lookupOrDefault<scalar>
+                (
+                    "planarTolerance",
+                    surfaceIntersection::defaultTolerance
+                )
+            );
+
+            // surf.write("selfIntersection-input.obj");
+
+            Info<<"self-intersection "
+                << intersect.cutEdges().size() << " edges "
+                << intersect.cutPoints().size() << " points" << nl;
+
+            if (intersect.cutEdges().size())
+            {
+                extendedEdgeMesh addMesh
+                (
+                    xferCopy<pointField>(intersect.cutPoints()),
+                    xferCopy<edgeList>(intersect.cutEdges())
+                );
+
+                feMesh.add(addMesh);
+            }
+        }
+
         Info<< nl
             << "Final feature set:" << nl;
         writeStats(feMesh, Info);
@@ -1226,7 +1337,7 @@ int main(int argc, char *argv[])
 
         if (writeObj)
         {
-            feMesh.writeObj(feMesh.path()/surfFileName.lessExt().name());
+            feMesh.writeObj(feMesh.path()/sFeatFileName);
         }
 
         feMesh.write();
@@ -1236,10 +1347,10 @@ int main(int argc, char *argv[])
         (
             IOobject
             (
-                surfFileName.lessExt().name() + ".eMesh",   // name
-                runTime.constant(),                         // instance
+                sFeatFileName + ".eMesh",   // name
+                runTime.constant(),         // instance
                 "triSurface",
-                runTime,                                    // registry
+                runTime,                    // registry
                 IOobject::NO_READ,
                 IOobject::AUTO_WRITE,
                 false

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict

@@ -114,4 +114,41 @@ surface2.nas
 }
 
 
+// Handle multiple surfaces
+//
+// - If the dictionary is named 'surfaces', it must also contain a 'surfaces'
+//   entry (wordRe list).
+//
+// - If other dictionaries may contain a 'surfaces' entry, it will be taken
+//   for the input.
+//
+surfaces
+{
+    extractionMethod    extractFromSurface;
+
+    surfaces            (surface1.stl surface2.nas);
+
+    // Base output name (optional)
+    // output              surfaces;
+
+    // Generate features from self-intersect
+    selfIntersection    true;
+
+    // Tolerance for self-intersect
+    planarTolerance     1e-3;
+
+    extractFromSurfaceCoeffs
+    {
+        includedAngle   120;
+
+        // Do not mark region edges
+        geometricTestOnly       yes;
+    }
+
+    // Write options
+
+        // Write features to obj format for postprocessing
+        writeObj                yes;
+}
+
 // ************************************************************************* //

src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C

@@ -829,7 +829,7 @@ Foam::labelList Foam::surfaceFeatures::trimFeatures
 }
 
 
-void Foam::surfaceFeatures::writeDict(Ostream& writeFile) const
+void Foam::surfaceFeatures::writeDict(Ostream& os) const
 {
 
     dictionary featInfoDict;
@@ -838,15 +838,14 @@ void Foam::surfaceFeatures::writeDict(Ostream& writeFile) const
     featInfoDict.add("featureEdges", featureEdges_);
     featInfoDict.add("featurePoints", featurePoints_);
 
-    featInfoDict.write(writeFile);
+    featInfoDict.write(os);
 }
 
 
 void Foam::surfaceFeatures::write(const fileName& fName) const
 {
-    OFstream str(fName);
-
-    writeDict(str);
+    OFstream os(fName);
+    writeDict(os);
 }
 
 

src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H

@@ -177,12 +177,12 @@ public:
     // Constructors
 
         //- Construct from surface
-        surfaceFeatures(const triSurface&);
+        surfaceFeatures(const triSurface& surf);
 
         //- Construct from components
         surfaceFeatures
         (
-            const triSurface&,
+            const triSurface& surf,
             const labelList& featurePoints,
             const labelList& featureEdges,
             const label externalStart,
@@ -195,7 +195,7 @@ public:
         //  geometric criteria
         surfaceFeatures
         (
-            const triSurface&,
+            const triSurface& surf,
             const scalar includedAngle,
             const scalar minLen = 0,
             const label minElems = 0,
@@ -203,15 +203,15 @@ public:
         );
 
         //- Construct from dictionary
-        surfaceFeatures(const triSurface&, const dictionary& dict);
+        surfaceFeatures(const triSurface& surf, const dictionary& dict);
 
         //- Construct from file
-        surfaceFeatures(const triSurface&, const fileName& fName);
+        surfaceFeatures(const triSurface& surf, const fileName& fName);
 
         //- Construct from pointField and edgeList (edgeMesh)
         surfaceFeatures
         (
-            const triSurface&,
+            const triSurface& surf,
             const pointField& points,
             const edgeList& edges,
             const scalar mergeTol = 1e-6,
@@ -219,7 +219,7 @@ public:
         );
 
         //- Construct as copy
-        surfaceFeatures(const surfaceFeatures&);
+        surfaceFeatures(const surfaceFeatures& sf);
 
 
     // Member Functions
@@ -405,7 +405,7 @@ public:
         // Write
 
             //- Write as dictionary
-            void writeDict(Ostream&) const;
+            void writeDict(Ostream& os) const;
 
             //- Write as dictionary to file
             void write(const fileName& fName) const;
@@ -418,7 +418,7 @@ public:
 
     // Member Operators
 
-        void operator=(const surfaceFeatures&);
+        void operator=(const surfaceFeatures& rhs);
 
 
 };