diff --git a/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C b/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C
index aade46b54c..e48f7b48e3 100644
--- a/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C
+++ b/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C
@@ -88,10 +88,9 @@ int main(int argc, char *argv[])
 
     forAllConstIter(dictionary, agglomDict, iter)
     {
-        labelList patchids = boundary.findIndices(iter().keyword());
-        forAll(patchids, i)
+        labelList patchids = boundary.indices(iter().keyword());
+        for (const label patchi : patchids)
         {
-            label patchi =  patchids[i];
             const polyPatch& pp = boundary[patchi];
 
             if (!pp.coupled())
diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index 30ec3a8c8d..7d575efa82 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -338,12 +338,11 @@ int main(int argc, char *argv[])
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
     const polyBoundaryMesh& coarsePatches = coarseMesh.boundaryMesh();
 
-    labelList viewFactorsPatches(patches.findIndices(viewFactorWall));
-    forAll(viewFactorsPatches, i)
+    labelList viewFactorsPatches(patches.indices(viewFactorWall));
+    for (const label patchi : viewFactorsPatches)
     {
-        label patchI = viewFactorsPatches[i];
-        nCoarseFaces += coarsePatches[patchI].size();
-        nFineFaces += patches[patchI].size();
+        nCoarseFaces += coarsePatches[patchi].size();
+        nFineFaces += patches[patchi].size();
     }
 
     // total number of coarse faces
@@ -370,10 +369,8 @@ int main(int argc, char *argv[])
     DynamicList<label> localAgg(nCoarseFaces);
     labelHashSet includePatches;
 
-    forAll(viewFactorsPatches, i)
+    for (const label patchID : viewFactorsPatches)
     {
-        const label patchID = viewFactorsPatches[i];
-
         const polyPatch& pp = patches[patchID];
         const labelList& agglom = finalAgglom[patchID];
 
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C
index e949298eaa..53faf5c046 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,7 +55,7 @@ readField
     {
         if (iter().isDict() && !iter().keyword().isPattern())
         {
-            label patchi = bmesh_.findPatchID(iter().keyword());
+            const label patchi = bmesh_.findPatchID(iter().keyword());
 
             if (patchi != -1)
             {
@@ -98,16 +98,11 @@ readField
 
             if (e.isDict() && !e.keyword().isPattern())
             {
-                const labelList patchIDs = bmesh_.findIndices
-                (
-                    e.keyword(),
-                    true                    // use patchGroups
-                );
+                const labelList patchIds =
+                    bmesh_.indices(e.keyword(), true); // use patchGroups
 
-                forAll(patchIDs, i)
+                for (const label patchi : patchIds)
                 {
-                    label patchi = patchIDs[i];
-
                     if (!this->set(patchi))
                     {
                         this->set
diff --git a/src/OpenFOAM/meshes/Identifiers/DynamicID/DynamicID.H b/src/OpenFOAM/meshes/Identifiers/DynamicID/DynamicID.H
index a0d82fb779..c3e341296d 100644
--- a/src/OpenFOAM/meshes/Identifiers/DynamicID/DynamicID.H
+++ b/src/OpenFOAM/meshes/Identifiers/DynamicID/DynamicID.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,9 +44,8 @@ Description
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
+// Forward declarations
 template<class> class DynamicID;
-
 template<class ObjectType>
 Ostream& operator<<(Ostream&, const DynamicID<ObjectType>&);
 
@@ -74,18 +73,19 @@ public:
         DynamicID(const keyType& key, const ObjectType& obj)
         :
             key_(key),
-            indices_(obj.findIndices(key_))
+            indices_(obj.indices(key_))
         {}
 
         //- Construct from Istream
         DynamicID(Istream& is, const ObjectType& obj)
         :
             key_(is),
-            indices_(obj.findIndices(key_))
+            indices_(obj.indices(key_))
         {}
 
 
-    // Destructor - default
+    //- Destructor
+    ~DynamicID() = default;
 
 
     // Member Functions
@@ -107,7 +107,7 @@ public:
             //- Return index of first matching zone
             label index() const
             {
-                return indices_.empty() ? -1 : indices_[0];
+                return indices_.empty() ? -1 : indices_.first();
             }
 
             //- Has the zone been found
@@ -122,7 +122,7 @@ public:
             //- Update
             void update(const ObjectType& obj)
             {
-                indices_ = obj.findIndices(key_);
+                indices_ = obj.indices(key_);
             }
 
 
diff --git a/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C b/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C
index b466bbc30f..1475cf6c8a 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C
+++ b/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -58,19 +58,20 @@ Foam::pointBoundaryMesh::pointBoundaryMesh
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::label Foam::pointBoundaryMesh::findPatchID(const word& patchName) const
+
+Foam::labelList Foam::pointBoundaryMesh::indices
+(
+    const keyType& key,
+    const bool useGroups
+) const
 {
-    return mesh()().boundaryMesh().findPatchID(patchName);
+    return mesh()().boundaryMesh().indices(key, useGroups);
 }
 
 
-Foam::labelList Foam::pointBoundaryMesh::findIndices
-(
-    const keyType& key,
-    const bool usePatchGroups
-) const
+Foam::label Foam::pointBoundaryMesh::findPatchID(const word& patchName) const
 {
-    return mesh()().boundaryMesh().findIndices(key, usePatchGroups);
+    return mesh()().boundaryMesh().findPatchID(patchName);
 }
 
 
diff --git a/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.H b/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.H
index f0abdc3beb..71ff76f839 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.H
+++ b/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.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) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,7 +42,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
+// Forward declarations
 class pointMesh;
 class polyBoundaryMesh;
 
@@ -62,7 +62,7 @@ class pointBoundaryMesh
 
     // Private Member Functions
 
-        //- Calculate the geometry for the patches (transformation tensors etc.)
+        //- Calculate geometry for the patches (transformation tensors etc.)
         void calcGeometry();
 
         //- No copy construct
@@ -88,7 +88,7 @@ public:
         );
 
 
-    // Member functions
+    // Member Functions
 
         //- Return the mesh reference
         const pointMesh& mesh() const
@@ -96,17 +96,27 @@ public:
             return mesh_;
         }
 
+        //- Find patch indices given a name
+        labelList indices(const keyType& key, const bool useGroups) const;
+
         //- Find patch index given a name
         label findPatchID(const word& patchName) const;
 
-        //- Find patch indices given a name
-        labelList findIndices(const keyType&, const bool useGroups) const;
-
         //- Correct polyBoundaryMesh after moving points
         void movePoints(const pointField&);
 
         //- Correct polyBoundaryMesh after topology update
         void updateMesh();
+
+
+    // Housekeeping
+
+        //- Identical to the indices() method (AUG-2018)
+        labelList findIndices(const keyType& key, const bool useGroups) const
+        {
+            return indices(key, useGroups);
+        }
+
 };
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
index cfe173a063..42140d7332 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
@@ -40,6 +40,81 @@ namespace Foam
     defineTypeNameAndDebug(polyBoundaryMesh, 0);
 }
 
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    // Templated implementation for types(), names(), etc - file-scope
+    template<class ListType, class GetOp>
+    static ListType getMethodImpl
+    (
+        const polyPatchList& list,
+        const GetOp& getop
+    )
+    {
+        const label len = list.size();
+
+        ListType output(len);
+
+        for (label i = 0; i < len; ++i)
+        {
+            output[i] = getop(list[i]);
+        }
+
+        return output;
+    }
+
+
+    // Templated implementation for indices() - file-scope
+    template<class UnaryMatchPredicate>
+    static labelList indicesImpl
+    (
+        const polyPatchList& list,
+        const UnaryMatchPredicate& matcher
+    )
+    {
+        const label len = list.size();
+
+        labelList output(len);
+
+        label count = 0;
+        for (label i = 0; i < len; ++i)
+        {
+            if (matcher(list[i].name()))
+            {
+                output[count++] = i;
+            }
+        }
+
+        output.resize(count);
+
+        return output;
+    }
+
+
+    // Templated implementation for findIndex() - file-scope
+    template<class UnaryMatchPredicate>
+    label findIndexImpl
+    (
+        const polyPatchList& list,
+        const UnaryMatchPredicate& matcher
+    )
+    {
+        const label len = list.size();
+
+        for (label i = 0; i < len; ++i)
+        {
+            if (matcher(list[i].name()))
+            {
+                return i;
+            }
+        }
+
+        return -1;
+    }
+
+} // End namespace Foam
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -169,9 +244,11 @@ Foam::polyBoundaryMesh::polyBoundaryMesh
 
 void Foam::polyBoundaryMesh::clearGeom()
 {
-    forAll(*this, patchi)
+    polyPatchList& patches = *this;
+
+    for (polyPatch& p : patches)
     {
-        operator[](patchi).clearGeom();
+        p.clearGeom();
     }
 }
 
@@ -182,9 +259,11 @@ void Foam::polyBoundaryMesh::clearAddressing()
     patchIDPtr_.clear();
     groupPatchIDsPtr_.clear();
 
-    forAll(*this, patchi)
+    polyPatchList& patches = *this;
+
+    for (polyPatch& p : patches)
     {
-        operator[](patchi).clearAddressing();
+        p.clearAddressing();
     }
 }
 
@@ -402,11 +481,11 @@ Foam::polyBoundaryMesh::groupPatchIDs() const
         groupPatchIDsPtr_.reset(new HashTable<labelList>(16));
         auto& groupPatchIDs = *groupPatchIDsPtr_;
 
-        const polyBoundaryMesh& bm = *this;
+        const polyBoundaryMesh& patches = *this;
 
-        forAll(bm, patchi)
+        forAll(patches, patchi)
         {
-            const wordList& groups = bm[patchi].inGroups();
+            const wordList& groups = patches[patchi].inGroups();
 
             for (const word& groupName : groups)
             {
@@ -424,12 +503,12 @@ Foam::polyBoundaryMesh::groupPatchIDs() const
         }
 
         // Remove patch names from patchGroups
-        forAll(bm, patchi)
+        forAll(patches, patchi)
         {
-            if (groupPatchIDs.erase(bm[patchi].name()))
+            if (groupPatchIDs.erase(patches[patchi].name()))
             {
                 WarningInFunction
-                    << "Removing patchGroup '" << bm[patchi].name()
+                    << "Removed patchGroup '" << patches[patchi].name()
                     << "' which clashes with patch " << patchi
                     << " of the same name."
                     << endl;
@@ -513,46 +592,24 @@ Foam::label Foam::polyBoundaryMesh::nNonProcessor() const
 
 Foam::wordList Foam::polyBoundaryMesh::names() const
 {
-    const polyPatchList& patches = *this;
-
-    wordList list(patches.size());
-
-    forAll(patches, patchi)
-    {
-        list[patchi] = patches[patchi].name();
-    }
-
-    return list;
+    return getMethodImpl<wordList>(*this, getNameOp<polyPatch>());
 }
 
 
 Foam::wordList Foam::polyBoundaryMesh::types() const
 {
-    const polyPatchList& patches = *this;
-
-    wordList list(patches.size());
-
-    forAll(patches, patchi)
-    {
-        list[patchi] = patches[patchi].type();
-    }
-
-    return list;
+    return getMethodImpl<wordList>(*this, getTypeOp<polyPatch>());
 }
 
 
 Foam::wordList Foam::polyBoundaryMesh::physicalTypes() const
 {
-    const polyPatchList& patches = *this;
-
-    wordList list(patches.size());
-
-    forAll(patches, patchi)
-    {
-        list[patchi] = patches[patchi].physicalType();
-    }
-
-    return list;
+    return
+        getMethodImpl<wordList>
+        (
+            *this,
+            [](const polyPatch& p) { return p.physicalType(); }
+        );
 }
 
 
@@ -578,84 +635,79 @@ Foam::labelRange Foam::polyBoundaryMesh::range(const label patchi) const
 }
 
 
-Foam::labelList Foam::polyBoundaryMesh::findIndices
+Foam::labelList Foam::polyBoundaryMesh::indices
 (
     const keyType& key,
-    const bool usePatchGroups
+    const bool useGroups
 ) const
 {
-    DynamicList<label> indices;
-
     if (key.empty())
     {
-        // no-op
+        return labelList();
     }
-    else if (key.isPattern())
+
+    DynamicList<label> patchIndices;
+
+    if (key.isPattern())
     {
-        const regExp keyRe(key);
+        const regExp matcher(key);
 
-        indices = findStrings(keyRe, this->names());
+        patchIndices = indicesImpl(*this, matcher);
 
-        if (usePatchGroups && groupPatchIDs().size())
+        // Only examine patch groups if requested and when they exist.
+        if (useGroups && !groupPatchIDs().empty())
         {
-            labelHashSet indexSet(indices);
+            const wordList groupNames
+            (
+                groupPatchIDs().tocKeys(matcher)
+            );
 
-            const wordList allGroupNames = groupPatchIDs().toc();
-            labelList groupIDs = findStrings(keyRe, allGroupNames);
-
-            for (const label groupi : groupIDs)
+            if (groupNames.size())
             {
-                const word& grpName = allGroupNames[groupi];
+                labelHashSet groupIndices;
 
-                const labelList& patchIDs = groupPatchIDs()[grpName];
-
-                for (const label patchi : patchIDs)
+                for (const word& grpName : groupNames)
                 {
-                    if (indexSet.insert(patchi))
-                    {
-                        indices.append(patchi);
-                    }
+                    // Hash the patch ids for the group
+                    groupIndices.insert( groupPatchIDs()[grpName] );
                 }
+
+                groupIndices.erase(patchIndices);  // Skip existing
+                patchIndices.append(groupIndices.sortedToc());
             }
         }
     }
     else
     {
-        // Literal string. Special version of above to avoid
-        // unnecessary memory allocations
+        // Literal string.
+        // Special version of above for reduced memory footprint
 
-        indices.setCapacity(1);
-        forAll(*this, i)
+        const word& matcher = key;
+
+        const label patchId = findIndexImpl(*this, matcher);
+
+        if (patchId >= 0)
         {
-            if (key == operator[](i).name())
-            {
-                indices.append(i);
-                break;
-            }
+            patchIndices.append(patchId);
         }
 
-        if (usePatchGroups && groupPatchIDs().size())
+        // Only examine patch groups if requested and when they exist.
+        if (useGroups && !groupPatchIDs().empty())
         {
             const auto iter = groupPatchIDs().cfind(key);
 
             if (iter.found())
             {
-                labelHashSet indexSet(indices);
+                // Hash the patch ids for the group
+                labelHashSet groupIndices(iter.object());
 
-                const labelList& patchIDs = *iter;
-
-                for (const label patchi : patchIDs)
-                {
-                    if (indexSet.insert(patchi))
-                    {
-                        indices.append(patchi);
-                    }
-                }
+                groupIndices.erase(patchIndices);  // Skip existing
+                patchIndices.append(groupIndices.sortedToc());
             }
         }
     }
 
-    return indices;
+    return patchIndices;
 }
 
 
@@ -663,31 +715,20 @@ Foam::label Foam::polyBoundaryMesh::findIndex(const keyType& key) const
 {
     if (key.empty())
     {
-        // no-op
+        return -1;
     }
     else if (key.isPattern())
     {
-        labelList indices = this->findIndices(key);
-
-        // Return the first element
-        if (!indices.empty())
-        {
-            return indices.first();
-        }
+        // Find as regex
+        regExp matcher(key);
+        return findIndexImpl(*this, matcher);
     }
     else
     {
-        forAll(*this, i)
-        {
-            if (key == operator[](i).name())
-            {
-                return i;
-            }
-        }
+        // Find as literal string
+        const word& matcher = key;
+        return findIndexImpl(*this, matcher);
     }
-
-    // Not found, return -1
-    return -1;
 }
 
 
@@ -697,14 +738,11 @@ Foam::label Foam::polyBoundaryMesh::findPatchID
     bool allowNotFound
 ) const
 {
-    const polyPatchList& patches = *this;
+    const label patchId = findIndexImpl(*this, patchName);
 
-    forAll(patches, patchi)
+    if (patchId >= 0)
     {
-        if (patches[patchi].name() == patchName)
-        {
-            return patchi;
-        }
+        return patchId;
     }
 
     if (!allowNotFound)
@@ -789,7 +827,7 @@ Foam::labelHashSet Foam::polyBoundaryMesh::patchSet
 (
     const UList<wordRe>& patchNames,
     const bool warnNotFound,
-    const bool usePatchGroups
+    const bool useGroups
 ) const
 {
     const wordList allPatchNames(this->names());
@@ -799,28 +837,26 @@ Foam::labelHashSet Foam::polyBoundaryMesh::patchSet
     {
         // Treat the given patch names as wild-cards and search the set
         // of all patch names for matches
-        labelList patchIDs = findStrings(patchName, allPatchNames);
-        ids.insert(patchIDs);
+        labelList patchIndices = findStrings(patchName, allPatchNames);
+        ids.insert(patchIndices);
 
-        if (patchIDs.empty())
+        if (patchIndices.empty())
         {
-            if (usePatchGroups)
+            if (useGroups)
             {
-                const wordList allGroupNames = groupPatchIDs().toc();
+                // Treat as group name or regex for group name
 
-                // Regard as group name
-                labelList groupIDs = findStrings(patchName, allGroupNames);
+                const wordList groupNames
+                (
+                    groupPatchIDs().tocKeys(patchName)
+                );
 
-                for (const label groupi : groupIDs)
+                for (const word& grpName : groupNames)
                 {
-                    const word& groupName = allGroupNames[groupi];
-                    ids.insert
-                    (
-                        groupPatchIDs()[groupName]
-                    );
+                    ids.insert( groupPatchIDs()[grpName] );
                 }
 
-                if (groupIDs.empty() && warnNotFound)
+                if (groupNames.empty() && warnNotFound)
                 {
                     WarningInFunction
                         << "Cannot find any patch or group names matching "
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
index 397257ee0a..0947df807d 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
@@ -175,13 +175,15 @@ public:
         labelRange range(const label patchi) const;
 
 
-        //- Return patch indices for all matches. Optionally matches patchGroups
-        labelList findIndices
+        //- Return patch indices for all matches.
+        //  Optionally matches patchGroups
+        labelList indices
         (
             const keyType& key,
-            const bool usePatchGroups = true
+            const bool useGroups = true
         ) const;
 
+
         //- Return patch index for the first match, return -1 if not found
         label findIndex(const keyType& key) const;
 
@@ -215,7 +217,7 @@ public:
         (
             const UList<wordRe>& patchNames,
             const bool warnNotFound = true,
-            const bool usePatchGroups = true
+            const bool useGroups = true
         ) const;
 
         //- Match the patches to groups.
@@ -276,6 +278,19 @@ public:
     // Ostream operator
 
         friend Ostream& operator<<(Ostream& os, const polyBoundaryMesh& pbm);
+
+
+    // Housekeeping
+
+        //- Identical to the indices() method (AUG-2018)
+        labelList findIndices
+        (
+            const keyType& key,
+            const bool useGroups = true
+        ) const
+        {
+            return this->indices(key, useGroups);
+        }
 };
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMeshTemplates.C b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMeshTemplates.C
index 76be9e370f..05548ec49c 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMeshTemplates.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMeshTemplates.C
@@ -28,13 +28,13 @@ License
 template<class Type>
 Foam::labelHashSet Foam::polyBoundaryMesh::findPatchIDs() const
 {
-    const polyBoundaryMesh& bm = *this;
+    const polyBoundaryMesh& patches = *this;
 
-    labelHashSet patchIDs(bm.size());
+    labelHashSet patchIDs(patches.size());
 
-    forAll(bm, patchi)
+    forAll(patches, patchi)
     {
-        if (isA<Type>(bm[patchi]))
+        if (isA<Type>(patches[patchi]))
         {
             patchIDs.insert(patchi);
         }
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C
index 9163b5c19c..308496b555 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C
@@ -176,9 +176,13 @@ Foam::labelList Foam::processorCyclicPolyPatch::patchIDs
     const polyBoundaryMesh& bm
 )
 {
-    return bm.findIndices
+    return bm.indices
     (
-        string("procBoundary.*to.*through" + cyclicPolyPatchName)
+        keyType
+        (
+            "procBoundary.*to.*through" + cyclicPolyPatchName,
+            true  // isPattern
+        )
     );
 }
 
diff --git a/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C b/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C
index 6fe7f6fac1..ab8073d00d 100644
--- a/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C
@@ -46,8 +46,6 @@ namespace Foam
 template<class ZoneType, class MeshType>
 void Foam::ZoneMesh<ZoneType, MeshType>::calcZoneMap() const
 {
-    // It is an error to attempt to recalculate cellEdges
-    // if the pointer is already set
     if (zoneMapPtr_)
     {
         FatalErrorInFunction
@@ -59,9 +57,11 @@ void Foam::ZoneMesh<ZoneType, MeshType>::calcZoneMap() const
         // Count number of objects in all zones
         label nObjects = 0;
 
-        forAll(*this, zonei)
+        const PtrList<ZoneType>& zones = *this;
+
+        for (const ZoneType& zn : zones)
         {
-            nObjects += this->operator[](zonei).size();
+            nObjects += zn.size();
         }
 
         zoneMapPtr_ = new Map<label>(2*nObjects);
@@ -69,14 +69,19 @@ void Foam::ZoneMesh<ZoneType, MeshType>::calcZoneMap() const
 
         // Fill in objects of all zones into the map.
         // The key is the global object index, value is the zone index
-        forAll(*this, zonei)
-        {
-            const labelList& zoneObjects = this->operator[](zonei);
 
-            for (const label idx : zoneObjects)
+        label zonei = 0;
+
+        for (const ZoneType& zn : zones)
+        {
+            const labelList& labels = zn;
+
+            for (const label idx : labels)
             {
                 zm.insert(idx, zonei);
             }
+
+            ++zonei;
         }
     }
 }
@@ -101,7 +106,7 @@ bool Foam::ZoneMesh<ZoneType, MeshType>::read()
         Istream& is = readStream(typeName);
 
         PtrList<entry> patchEntries(is);
-        zones.setSize(patchEntries.size());
+        zones.resize(patchEntries.size());
 
         forAll(zones, zonei)
         {
@@ -167,6 +172,55 @@ Foam::wordList Foam::ZoneMesh<ZoneType, MeshType>::namesImpl
 }
 
 
+template<class ZoneType, class MeshType>
+template<class UnaryMatchPredicate>
+Foam::labelList Foam::ZoneMesh<ZoneType, MeshType>::indicesImpl
+(
+    const PtrList<ZoneType>& list,
+    const UnaryMatchPredicate& matcher
+)
+{
+    const label len = list.size();
+
+    labelList output(len);
+
+    label count = 0;
+    for (label i = 0; i < len; ++i)
+    {
+        if (matcher(list[i].name()))
+        {
+            output[count++] = i;
+        }
+    }
+
+    output.resize(count);
+
+    return output;
+}
+
+
+template<class ZoneType, class MeshType>
+template<class UnaryMatchPredicate>
+Foam::label Foam::ZoneMesh<ZoneType, MeshType>::findIndexImpl
+(
+    const PtrList<ZoneType>& list,
+    const UnaryMatchPredicate& matcher
+)
+{
+    const label len = list.size();
+
+    for (label i = 0; i < len; ++i)
+    {
+        if (matcher(list[i].name()))
+        {
+            return i;
+        }
+    }
+
+    return -1;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class ZoneType, class MeshType>
@@ -220,7 +274,7 @@ Foam::ZoneMesh<ZoneType, MeshType>::ZoneMesh
     {
         // Nothing read. Use supplied zones
         PtrList<ZoneType>& zones = *this;
-        zones.setSize(pzm.size());
+        zones.resize(pzm.size());
 
         forAll(zones, zonei)
         {
@@ -269,14 +323,14 @@ Foam::wordList Foam::ZoneMesh<ZoneType, MeshType>::types() const
 {
     const PtrList<ZoneType>& zones = *this;
 
-    wordList lst(zones.size());
+    wordList list(zones.size());
 
     forAll(zones, zonei)
     {
-        lst[zonei] = zones[zonei].type();
+        list[zonei] = zones[zonei].type();
     }
 
-    return lst;
+    return list;
 }
 
 
@@ -285,14 +339,14 @@ Foam::wordList Foam::ZoneMesh<ZoneType, MeshType>::names() const
 {
     const PtrList<ZoneType>& zones = *this;
 
-    wordList lst(zones.size());
+    wordList list(zones.size());
 
     forAll(zones, zonei)
     {
-        lst[zonei] = zones[zonei].name();
+        list[zonei] = zones[zonei].name();
     }
 
-    return lst;
+    return list;
 }
 
 
@@ -320,10 +374,10 @@ const
 template<class ZoneType, class MeshType>
 Foam::wordList Foam::ZoneMesh<ZoneType, MeshType>::sortedNames() const
 {
-    wordList sortedLst = this->names();
-    sort(sortedLst);
+    wordList sorted(this->names());
+    Foam::sort(sorted);
 
-    return sortedLst;
+    return sorted;
 }
 
 
@@ -349,36 +403,42 @@ const
 
 
 template<class ZoneType, class MeshType>
-Foam::labelList Foam::ZoneMesh<ZoneType, MeshType>::findIndices
+Foam::labelList Foam::ZoneMesh<ZoneType, MeshType>::indices
 (
     const keyType& key
 ) const
 {
-    labelList indices;
-
     if (key.empty())
     {
-        // no-op
+        return labelList();
     }
     else if (key.isPattern())
     {
-        indices = findStrings(key, this->names());
+        // Match as regex
+        regExp matcher(key);
+        return indicesImpl(*this, matcher);
     }
     else
     {
-        indices.resize(this->size());
-        label count = 0;
-        forAll(*this, i)
-        {
-            if (key == operator[](i).name())
-            {
-                indices[count++] = i;
-            }
-        }
-        indices.resize(count);
+        // Compare as literal string
+        const word& matcher = key;
+        return indicesImpl(*this, matcher);
+    }
+}
+
+
+template<class ZoneType, class MeshType>
+Foam::labelList Foam::ZoneMesh<ZoneType, MeshType>::indices
+(
+    const wordRes& matcher
+) const
+{
+    if (matcher.empty())
+    {
+        return labelList();
     }
 
-    return indices;
+    return indicesImpl(*this, matcher);
 }
 
 
@@ -390,30 +450,20 @@ Foam::label Foam::ZoneMesh<ZoneType, MeshType>::findIndex
 {
     if (key.empty())
     {
-        // no-op
+        return -1;
     }
     else if (key.isPattern())
     {
-        labelList indices = this->findIndices(key);
-
-        if (!indices.empty())
-        {
-            return indices.first();  // first match
-        }
+        // Find as regex
+        regExp matcher(key);
+        return findIndexImpl(*this, matcher);
     }
     else
     {
-        forAll(*this, i)
-        {
-            if (key == operator[](i).name())
-            {
-                return i;
-            }
-        }
+        // Find as literal string
+        const word& matcher = key;
+        return findIndexImpl(*this, matcher);
     }
-
-    // Not found
-    return -1;
 }
 
 
@@ -423,14 +473,11 @@ Foam::label Foam::ZoneMesh<ZoneType, MeshType>::findZoneID
     const word& zoneName
 ) const
 {
-    const PtrList<ZoneType>& zones = *this;
+    label zoneId = findIndexImpl(*this, zoneName);
 
-    forAll(zones, zonei)
+    if (zoneId >= 0)
     {
-        if (zoneName == zones[zonei].name())
-        {
-            return zonei;
-        }
+        return zoneId;
     }
 
     // Zone not found
@@ -444,6 +491,7 @@ Foam::label Foam::ZoneMesh<ZoneType, MeshType>::findZoneID
     if (disallowGenericZones != 0)
     {
         // Create a new ...
+
         Info<< "Creating dummy zone " << zoneName << endl;
         dictionary dict;
         dict.set("type", ZoneType::typeName);
@@ -452,44 +500,37 @@ Foam::label Foam::ZoneMesh<ZoneType, MeshType>::findZoneID
         // flipMap only really applicable for face zones, but should not get
         // in the way for cell- and point-zones...
         dict.set("flipMap", boolList());
-        label newZonei = zones.size();
 
-        ZoneMesh<ZoneType, MeshType>& zm =
-            const_cast<ZoneMesh<ZoneType, MeshType>&>(*this);
+        auto& zm = const_cast<ZoneMesh<ZoneType, MeshType>&>(*this);
+        zoneId = zm.size();
 
-        zm.append
-        (
-            new ZoneType
-            (
-                zoneName,
-                dict,
-                newZonei,
-                zm
-            )
-        );
-
-        return newZonei;
-    }
-    else
-    {
-        // Not found
-        return -1;
+        zm.append(new ZoneType(zoneName, dict, zoneId, zm));
     }
+
+    return zoneId;
 }
 
 
 template<class ZoneType, class MeshType>
-Foam::bitSet Foam::ZoneMesh<ZoneType, MeshType>::findMatching
+Foam::bitSet Foam::ZoneMesh<ZoneType, MeshType>::selection
 (
-    const keyType& key
+    const labelUList& zoneIds
 ) const
 {
     bitSet bitset;
 
-    const labelList indices = this->findIndices(key);
-
-    for (const label zonei : indices)
+    for (const label zonei : zoneIds)
     {
+        #ifdef FULLDEBUG
+        if (zonei < 0 || zonei >= this->size())
+        {
+            FatalErrorInFunction
+                << ZoneType::typeName << " "
+                << zonei << " out of range [0," << this->size() << ")"
+                << abort(FatalError);
+        }
+        #endif
+
         bitset.set
         (
             static_cast<const labelList&>(this->operator[](zonei))
@@ -500,6 +541,26 @@ Foam::bitSet Foam::ZoneMesh<ZoneType, MeshType>::findMatching
 }
 
 
+template<class ZoneType, class MeshType>
+Foam::bitSet Foam::ZoneMesh<ZoneType, MeshType>::selection
+(
+    const keyType& key
+) const
+{
+    return this->selection(this->indices(key));
+}
+
+
+template<class ZoneType, class MeshType>
+Foam::bitSet Foam::ZoneMesh<ZoneType, MeshType>::selection
+(
+    const wordRes& matcher
+) const
+{
+    return this->selection(this->indices(matcher));
+}
+
+
 template<class ZoneType, class MeshType>
 const ZoneType* Foam::ZoneMesh<ZoneType, MeshType>::zonePtr
 (
@@ -551,9 +612,9 @@ void Foam::ZoneMesh<ZoneType, MeshType>::clearAddressing()
 
     PtrList<ZoneType>& zones = *this;
 
-    forAll(zones, zonei)
+    for (ZoneType& zn : zones)
     {
-        zones[zonei].clearAddressing();
+        zn.clearAddressing();
     }
 }
 
@@ -576,10 +637,11 @@ bool Foam::ZoneMesh<ZoneType, MeshType>::checkDefinition
 
     const PtrList<ZoneType>& zones = *this;
 
-    forAll(zones, zonei)
+    for (const ZoneType& zn : zones)
     {
-        inError |= zones[zonei].checkDefinition(report);
+        inError |= zn.checkDefinition(report);
     }
+
     return inError;
 }
 
@@ -595,7 +657,6 @@ bool Foam::ZoneMesh<ZoneType, MeshType>::checkParallelSync
         return false;
     }
 
-
     const PtrList<ZoneType>& zones = *this;
 
     bool hasError = false;
@@ -613,7 +674,7 @@ bool Foam::ZoneMesh<ZoneType, MeshType>::checkParallelSync
 
     // Have every processor check but only master print error.
 
-    for (label proci = 1; proci < allNames.size(); proci++)
+    for (label proci = 1; proci < allNames.size(); ++proci)
     {
         if
         (
@@ -639,16 +700,16 @@ bool Foam::ZoneMesh<ZoneType, MeshType>::checkParallelSync
     // Check contents
     if (!hasError)
     {
-        forAll(zones, zonei)
+        for (const ZoneType& zn : zones)
         {
-            if (zones[zonei].checkParallelSync(false))
+            if (zn.checkParallelSync(false))
             {
                 hasError = true;
 
                 if (debug || (report && Pstream::master()))
                 {
-                    Info<< " ***Zone " << zones[zonei].name()
-                        << " of type " << zones[zonei].type()
+                    Info<< " ***Zone " << zn.name()
+                        << " of type " << zn.type()
                         << " is not correctly synchronised"
                         << " across coupled boundaries."
                         << " (coupled faces are either not both"
@@ -667,9 +728,9 @@ void Foam::ZoneMesh<ZoneType, MeshType>::movePoints(const pointField& pts)
 {
     PtrList<ZoneType>& zones = *this;
 
-    forAll(zones, zonei)
+    for (ZoneType& zn : zones)
     {
-        zones[zonei].movePoints(pts);
+        zn.movePoints(pts);
     }
 }
 
@@ -738,8 +799,7 @@ ZoneType& Foam::ZoneMesh<ZoneType, MeshType>::operator()
     if (zoneId < 0)
     {
         zoneId = zones.size();
-
-        zones.setSize(zoneId+1);
+        zones.resize(zoneId+1);
         zones.set(zoneId, new ZoneType(zoneName, zoneId, *this));
 
         if (verbose)
@@ -778,7 +838,7 @@ Foam::Ostream& Foam::operator<<
     {
         os  << sz << nl << token::BEGIN_LIST;
 
-        for (label i=0; i<sz; ++i)
+        for (label i=0; i < sz; ++i)
         {
             zones[i].writeDict(os);
         }
diff --git a/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.H b/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.H
index 2389c8f5b4..e469ae1b00 100644
--- a/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.H
@@ -47,8 +47,6 @@ SourceFiles
 namespace Foam
 {
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // Forward declarations
 
 template<class ZoneType, class MeshType> class ZoneMesh;
@@ -93,6 +91,22 @@ class ZoneMesh
             const bool doSort
         );
 
+        //- Templated implementation for indices()
+        template<class UnaryMatchPredicate>
+        static labelList indicesImpl
+        (
+            const PtrList<ZoneType>& list,
+            const UnaryMatchPredicate& matcher
+        );
+
+        //- Templated implementation for findIndex()
+        template<class UnaryMatchPredicate>
+        static label findIndexImpl
+        (
+            const PtrList<ZoneType>& list,
+            const UnaryMatchPredicate& matcher
+        );
+
 
         //- No copy construct
         ZoneMesh(const ZoneMesh&) = delete;
@@ -159,33 +173,53 @@ public:
         //- A list of the zone names
         wordList names() const;
 
-        //- A list of the zone names satisfying the input matcher
+        //- A list of zone names satisfying the input matcher
         wordList names(const wordRe& matcher) const;
 
-        //- A list of the zone names satisfying the input matchers
+        //- A list of zone names satisfying the input matchers
         wordList names(const wordRes& matcher) const;
 
         //- Sorted list of the zone names
         wordList sortedNames() const;
 
-        //- Sorted list of the zone names satisfying the input matcher
+        //- Sorted list of zone names satisfying the input matcher
         wordList sortedNames(const wordRe& matcher) const;
 
-        //- Sorted list of the zone names satisfying the input matchers
+        //- Sorted list of zone names satisfying the input matchers
         wordList sortedNames(const wordRes& matcher) const;
 
 
-        //- Return zone indices for all matches
-        labelList findIndices(const keyType& key) const;
-
         //- Return zone index for the first match, return -1 if not found
         label findIndex(const keyType& key) const;
 
         //- Find zone index given a name, return -1 if not found
         label findZoneID(const word& zoneName) const;
 
-        //- Mark items (cells, faces, points) that match the zone specification
-        bitSet findMatching(const keyType& key) const;
+
+        //- Return zone indices for all matches
+        labelList indices(const keyType& key) const;
+
+        //- Return zone indices for all matches
+        labelList indices(const wordRes& matcher) const;
+
+
+        //- Return all elements (cells, faces, points) contained in the
+        //- listed zones.
+        //  The bitSet is empty (zero-size) if there are no elements matched
+        //  anywhere.
+        bitSet selection(const labelUList& zoneIds) const;
+
+        //- Return all elements (cells, faces, points) that match the zone
+        //- specification as a bitSet.
+        //  The bitSet is empty (zero-size) if there are no elements matched
+        //  anywhere.
+        bitSet selection(const keyType& key) const;
+
+        //- Return all elements (cells, faces, points) that match the zone
+        //- specification as a bitSet.
+        //  The bitSet is empty (zero-size) if there are no elements matched
+        //  anywhere.
+        bitSet selection(const wordRes& matcher) const;
 
 
         //- Lookup zone by name and return const pointer, nullptr on error.
@@ -254,6 +288,15 @@ public:
             Ostream& os,
             const ZoneMesh<ZoneType, MeshType>& zones
         );
+
+
+    // Housekeeping
+
+        //- Identical to the indices() method (AUG-2018)
+        labelList findIndices(const keyType& key) const
+        {
+            return indices(key);
+        }
 };
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/zones/cellZone/cellZone.H b/src/OpenFOAM/meshes/polyMesh/zones/cellZone/cellZone.H
index 2de745493c..796ccde0c4 100644
--- a/src/OpenFOAM/meshes/polyMesh/zones/cellZone/cellZone.H
+++ b/src/OpenFOAM/meshes/polyMesh/zones/cellZone/cellZone.H
@@ -50,11 +50,9 @@ namespace Foam
 {
 
 // Forward declarations
-
 class cellZone;
 Ostream& operator<<(Ostream& os, const cellZone& zn);
 
-
 /*---------------------------------------------------------------------------*\
                            Class cellZone Declaration
 \*---------------------------------------------------------------------------*/
@@ -63,13 +61,6 @@ class cellZone
 :
     public zone
 {
-
-    // Private Member Functions
-
-        //- No copy construct
-        cellZone(const cellZone&) = delete;
-
-
 protected:
 
     // Protected data
@@ -78,9 +69,15 @@ protected:
         const cellZoneMesh& zoneMesh_;
 
 
+    // Protected Member Functions
+
+        //- No copy construct
+        cellZone(const cellZone&) = delete;
+
+
 public:
 
-    // Static data members
+    // Static Data Members
 
         //- The name associated with the zone-labels dictionary entry
         static const char * const labelsName;
diff --git a/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.H b/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.H
index 384b55777c..45ca5c8cf6 100644
--- a/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.H
+++ b/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.H
@@ -50,10 +50,8 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declarations
 class mapPolyMesh;
-
-// Forward declaration of friend functions and operators
-
 class faceZone;
 Ostream& operator<<(Ostream& os, const faceZone& zn);
 
@@ -80,7 +78,7 @@ class faceZone
 
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Flip map for all faces in the zone.
         //  Use true if the face needs flipping for the correct orientation.
@@ -119,7 +117,7 @@ protected:
 
 public:
 
-    // Static data members
+    // Static Data Members
 
         //- The name associated with the zone-labels dictionary entry
         static const char * const labelsName;
diff --git a/src/OpenFOAM/meshes/polyMesh/zones/pointZone/pointZone.H b/src/OpenFOAM/meshes/polyMesh/zones/pointZone/pointZone.H
index d2b5691cf6..01142dc1e8 100644
--- a/src/OpenFOAM/meshes/polyMesh/zones/pointZone/pointZone.H
+++ b/src/OpenFOAM/meshes/polyMesh/zones/pointZone/pointZone.H
@@ -51,8 +51,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
-
+// Forward declarations
 class pointZone;
 Ostream& operator<<(Ostream& os, const pointZone& zn);
 
@@ -65,23 +64,23 @@ class pointZone
 :
     public zone
 {
-    // Private Member Functions
-
-        //- No copy construct
-        pointZone(const pointZone&) = delete;
-
-
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Reference to zone list
         const pointZoneMesh& zoneMesh_;
 
 
+    // Protected Member Functions
+
+        //- No copy construct
+        pointZone(const pointZone&) = delete;
+
+
 public:
 
-    // Static data members
+    // Static Data Members
 
         //- The name associated with the zone-labels dictionary entry
         static const char * const labelsName;
diff --git a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
index 8000d31334..bc904e8e27 100644
--- a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
+++ b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2016-2017 Wikki Ltd
@@ -23,8 +23,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "faBoundaryMesh.H"
@@ -38,6 +36,81 @@ namespace Foam
     defineTypeNameAndDebug(faBoundaryMesh, 0);
 }
 
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    // Templated implementation for types(), names(), etc - file-scope
+    template<class ListType, class GetOp>
+    static ListType getMethodImpl
+    (
+        const faPatchList& list,
+        const GetOp& getop
+    )
+    {
+        const label len = list.size();
+
+        ListType output(len);
+
+        for (label i = 0; i < len; ++i)
+        {
+            output[i] = getop(list[i]);
+        }
+
+        return output;
+    }
+
+
+    // Templated implementation for indices() - file-scope
+    template<class UnaryMatchPredicate>
+    static labelList indicesImpl
+    (
+        const faPatchList& list,
+        const UnaryMatchPredicate& matcher
+    )
+    {
+        const label len = list.size();
+
+        labelList output(len);
+
+        label count = 0;
+        for (label i = 0; i < len; ++i)
+        {
+            if (matcher(list[i].name()))
+            {
+                output[count++] = i;
+            }
+        }
+
+        output.resize(count);
+
+        return output;
+    }
+
+
+    // Templated implementation for findIndex() - file-scope
+    template<class UnaryMatchPredicate>
+    label findIndexImpl
+    (
+        const faPatchList& list,
+        const UnaryMatchPredicate& matcher
+    )
+    {
+        const label len = list.size();
+
+        for (label i = 0; i < len; ++i)
+        {
+            if (matcher(list[i].name()))
+            {
+                return i;
+            }
+        }
+
+        return -1;
+    }
+
+} // End namespace Foam
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -61,22 +134,21 @@ Foam::faBoundaryMesh::faBoundaryMesh
         PtrList<entry> patchEntries(is);
         patches.setSize(patchEntries.size());
 
-        forAll(patches, patchI)
+        forAll(patches, patchi)
         {
             patches.set
             (
-                patchI,
+                patchi,
                 faPatch::New
                 (
-                    patchEntries[patchI].keyword(),
-                    patchEntries[patchI].dict(),
-                    patchI,
+                    patchEntries[patchi].keyword(),
+                    patchEntries[patchi].dict(),
+                    patchi,
                     *this
                 )
             );
         }
 
-        // Check state of IOstream
         is.check(FUNCTION_NAME);
 
         close();
@@ -139,86 +211,78 @@ Foam::lduInterfacePtrsList Foam::faBoundaryMesh::interfaces() const
 }
 
 
-Foam::wordList Foam::faBoundaryMesh::types() const
+Foam::wordList Foam::faBoundaryMesh::names() const
 {
-    const faPatchList& patches = *this;
-
-    wordList t(patches.size());
-
-    forAll(patches, patchI)
-    {
-        t[patchI] = patches[patchI].type();
-    }
-
-    return t;
+    return getMethodImpl<wordList>(*this, getNameOp<faPatch>());
 }
 
 
-Foam::wordList Foam::faBoundaryMesh::names() const
+Foam::wordList Foam::faBoundaryMesh::types() const
 {
-    const faPatchList& patches = *this;
+    return getMethodImpl<wordList>(*this, getTypeOp<faPatch>());
+}
 
-    wordList t(patches.size());
 
-    forAll(patches, patchI)
+Foam::labelList Foam::faBoundaryMesh::indices
+(
+    const keyType& key,
+    const bool useGroups  // ignored
+) const
+{
+    if (key.empty())
     {
-        t[patchI] = patches[patchI].name();
+        return labelList();
     }
 
-    return t;
+    if (key.isPattern())
+    {
+        const regExp matcher(key);
+
+        return indicesImpl(*this, matcher);
+    }
+    else
+    {
+        // Literal string.
+        // Special version of above for reduced memory footprint
+
+        const word& matcher = key;
+
+        const label patchId = findIndexImpl(*this, matcher);
+
+        if (patchId >= 0)
+        {
+            return labelList(one(), patchId);
+        }
+    }
+
+    return labelList();
+}
+
+
+Foam::label Foam::faBoundaryMesh::findIndex(const keyType& key) const
+{
+    if (key.empty())
+    {
+        return -1;
+    }
+    else if (key.isPattern())
+    {
+        // Find as regex
+        regExp matcher(key);
+        return findIndexImpl(*this, matcher);
+    }
+    else
+    {
+        // Find as literal string
+        const word& matcher = key;
+        return findIndexImpl(*this, matcher);
+    }
 }
 
 
 Foam::label Foam::faBoundaryMesh::findPatchID(const word& patchName) const
 {
-    const faPatchList& patches = *this;
-
-    forAll(patches, patchI)
-    {
-        if (patches[patchI].name() == patchName)
-        {
-            return patchI;
-        }
-    }
-
-    // Patch not found
-    return -1;
-}
-
-
-Foam::labelList Foam::faBoundaryMesh::findIndices
-(
-    const keyType& key,
-    const bool usePatchGroups
-) const
-{
-    DynamicList<label> indices;
-
-    if (key.empty())
-    {
-        // no-op
-    }
-    else if (key.isPattern())
-    {
-        indices = findStrings(key, this->names());
-    }
-    else
-    {
-        // Literal string. Special version of above to avoid
-        // unnecessary memory allocations
-
-        indices.setCapacity(1);
-        forAll(*this, i)
-        {
-            if (key == operator[](i).name())
-            {
-                indices.append(i);
-                break;
-            }
-        }
-    }
-
-    return indices;
+    return findIndexImpl(*this, patchName);
 }
 
 
@@ -228,22 +292,22 @@ Foam::label Foam::faBoundaryMesh::whichPatch(const label edgeIndex) const
     // with patch start labels.
     // If the face is internal, return -1;
     // if it is off the end of the list, abort
-    if (edgeIndex >= mesh().nEdges())
-    {
-        FatalErrorInFunction
-            << "given label greater than the number of edges"
-            << abort(FatalError);
-    }
-
     if (edgeIndex < mesh().nInternalEdges())
     {
         return -1;
     }
-
-    forAll(*this, patchI)
+    else if (edgeIndex >= mesh().nEdges())
     {
-        label start = mesh_.patchStarts()[patchI];
-        label size = operator[](patchI).faPatch::size();
+        FatalErrorInFunction
+            << "Edge " << edgeIndex
+            << " out of bounds. Number of geometric edges " << mesh().nEdges()
+            << abort(FatalError);
+    }
+
+    forAll(*this, patchi)
+    {
+        label start = mesh_.patchStarts()[patchi];
+        label size = operator[](patchi).faPatch::size();
 
         if
         (
@@ -251,7 +315,7 @@ Foam::label Foam::faBoundaryMesh::whichPatch(const label edgeIndex) const
          && edgeIndex < start + size
         )
         {
-            return patchI;
+            return patchi;
         }
     }
 
@@ -269,25 +333,31 @@ bool Foam::faBoundaryMesh::checkDefinition(const bool report) const
     label nextPatchStart = mesh().nInternalEdges();
     const faBoundaryMesh& bm = *this;
 
-    bool boundaryError = false;
+    bool hasError = false;
 
-    forAll(bm, patchI)
+    forAll(bm, patchi)
     {
-        if (bm[patchI].start() != nextPatchStart)
+        if (bm[patchi].start() != nextPatchStart && !hasError)
         {
-            boundaryError = true;
+            hasError = true;
 
             InfoInFunction
-                << "Problem with boundary patch " << patchI
-                << ".\nThe patch should start on face no " << nextPatchStart
-                << " and the boundary file specifies " << bm[patchI].start()
-                << "." << nl << endl;
+                << " ****Problem with boundary patch " << patchi
+                << " named " << bm[patchi].name()
+                << " of type " <<  bm[patchi].type()
+                << ". The patch should start on face no " << nextPatchStart
+                << " and the patch specifies " << bm[patchi].start()
+                << "." << endl
+                << "Possibly consecutive patches have this same problem."
+                << " Suppressing future warnings." << endl;
         }
 
-        nextPatchStart += bm[patchI].faPatch::size();
+        // Warn about duplicate boundary patches?
+
+        nextPatchStart += bm[patchi].faPatch::size();
     }
 
-    if (boundaryError)
+    if (hasError)
     {
         SeriousErrorInFunction
             << "This mesh is not valid: boundary definition is in error."
@@ -301,7 +371,7 @@ bool Foam::faBoundaryMesh::checkDefinition(const bool report) const
         }
     }
 
-    return boundaryError;
+    return hasError;
 }
 
 
@@ -309,14 +379,14 @@ void Foam::faBoundaryMesh::movePoints(const pointField& p)
 {
     faPatchList& patches = *this;
 
-    forAll(patches, patchI)
+    forAll(patches, patchi)
     {
-        patches[patchI].initMovePoints(p);
+        patches[patchi].initMovePoints(p);
     }
 
-    forAll(patches, patchI)
+    forAll(patches, patchi)
     {
-        patches[patchI].movePoints(p);
+        patches[patchi].movePoints(p);
     }
 }
 
@@ -343,19 +413,16 @@ bool Foam::faBoundaryMesh::writeData(Ostream& os) const
 
     os  << patches.size() << nl << token::BEGIN_LIST << incrIndent << nl;
 
-    forAll(patches, patchi)
+    for (const faPatch& p : patches)
     {
-        os  << indent << patches[patchi].name() << nl
-            << indent << token::BEGIN_BLOCK << nl
-            << incrIndent << patches[patchi] << decrIndent
-            << indent << token::END_BLOCK << endl;
+        os.beginBlock(p.name());
+        os << p;
+        os.endBlock();
     }
 
     os  << decrIndent << token::END_LIST;
 
-    // Check state of IOstream
     os.check(FUNCTION_NAME);
-
     return os.good();
 }
 
diff --git a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H
index f5a7c51b80..6030a2425e 100644
--- a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H
+++ b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2016-2017 Wikki Ltd
@@ -54,6 +54,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Forward declarations
 class faMesh;
 class faBoundaryMesh;
 Ostream& operator<<(Ostream&, const faBoundaryMesh&);
@@ -67,7 +68,7 @@ class faBoundaryMesh
     public faPatchList,
     public regIOobject
 {
-    // private data
+    // Private Data
 
         //- Reference to mesh
         const faMesh& mesh_;
@@ -103,10 +104,11 @@ public:
         );
 
 
-    // Destructor - default
+    //- Destructor
+    ~faBoundaryMesh() = default;
 
 
-    // Member functions
+    // Member Functions
 
         // Access
 
@@ -121,23 +123,26 @@ public:
             //  with only those pointing to interfaces being set
             lduInterfacePtrsList interfaces() const;
 
-            //- Return a list of patch types
-            wordList types() const;
-
             //- Return a list of patch names
             wordList names() const;
 
-            //- Find patch index given a name
-            label findPatchID(const word& patchName) const;
+            //- Return a list of patch types
+            wordList types() const;
 
-            //- Find patch indices given a name
-            //  Compatibility change  HJ, 12/Aug/2017
-            labelList findIndices
+            //- Return patch indices for all matches.
+            //  Matching patchGroups currently not supported
+            labelList indices
             (
-                const keyType&,
+                const keyType& key,
                 const bool useGroups = false
             ) const;
 
+            //- Return patch index for the first match, return -1 if not found
+            label findIndex(const keyType& key) const;
+
+            //- Find patch index given a name, return -1 if not found
+            label findPatchID(const word& patchName) const;
+
             //- Return patch index for a given edge label
             label whichPatch(const label edgeIndex) const;
 
@@ -160,6 +165,20 @@ public:
     // Ostream operator
 
         friend Ostream& operator<<(Ostream&, const faBoundaryMesh&);
+
+
+    // Housekeeping
+
+        //- Identical to the indices() method (AUG-2018)
+        labelList findIndices
+        (
+            const keyType& key,
+            const bool useGroups = false
+        ) const
+        {
+            return indices(key, useGroups);
+        }
+
 };
 
 
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
index d2d4305376..34d8717dad 100644
--- a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
@@ -109,7 +109,7 @@ Foam::porosityModel::porosityModel
         dict_.lookup("cellZone") >> zoneName_;
     }
 
-    cellZoneIDs_ = mesh_.cellZones().findIndices(zoneName_);
+    cellZoneIDs_ = mesh_.cellZones().indices(zoneName_);
 
     Info<< "    creating porous zone: " << zoneName_ << endl;
 
@@ -262,8 +262,8 @@ bool Foam::porosityModel::read(const dictionary& dict)
 
     coeffs_ = dict.optionalSubDict(type() + "Coeffs");
 
-    dict.lookup("cellZone") >> zoneName_;
-    cellZoneIDs_ = mesh_.cellZones().findIndices(zoneName_);
+    dict.readEntry("cellZone", zoneName_);
+    cellZoneIDs_ = mesh_.cellZones().indices(zoneName_);
 
     return true;
 }
diff --git a/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.C b/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.C
index dd33d2368b..e79535c3a1 100644
--- a/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.C
+++ b/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.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  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,7 +50,7 @@ Foam::fvBoundaryMesh::fvBoundaryMesh
     const fvMesh& m
 )
 :
-    fvPatchList(0),
+    fvPatchList(),
     mesh_(m)
 {}
 
@@ -70,6 +70,16 @@ Foam::fvBoundaryMesh::fvBoundaryMesh
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::labelList Foam::fvBoundaryMesh::indices
+(
+    const keyType& key,
+    const bool useGroups
+) const
+{
+    return mesh().boundaryMesh().indices(key, useGroups);
+}
+
+
 Foam::label Foam::fvBoundaryMesh::findPatchID(const word& patchName) const
 {
     const fvPatchList& patches = *this;
@@ -87,26 +97,18 @@ Foam::label Foam::fvBoundaryMesh::findPatchID(const word& patchName) const
 }
 
 
-Foam::labelList Foam::fvBoundaryMesh::findIndices
-(
-    const keyType& key,
-    const bool usePatchGroups
-) const
-{
-    return mesh().boundaryMesh().findIndices(key, usePatchGroups);
-}
-
-
 void Foam::fvBoundaryMesh::movePoints()
 {
-    forAll(*this, patchi)
+    fvPatchList& patches = *this;
+
+    for (fvPatch& p : patches)
     {
-        operator[](patchi).initMovePoints();
+        p.initMovePoints();
     }
 
-    forAll(*this, patchi)
+    for (fvPatch& p : patches)
     {
-        operator[](patchi).movePoints();
+        p.movePoints();
     }
 }
 
diff --git a/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H b/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H
index 2fe892bd64..052f435531 100644
--- a/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H
+++ b/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,6 +43,8 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declarations
+
 class fvMesh;
 class polyBoundaryMesh;
 
@@ -62,15 +64,16 @@ class fvBoundaryMesh
 
     // Private Member Functions
 
+        //- Add fvPatches corresponding to the given polyBoundaryMesh
+        void addPatches(const polyBoundaryMesh&);
+
+
         //- No copy construct
         fvBoundaryMesh(const fvBoundaryMesh&) = delete;
 
         //- No copy assignment
         void operator=(const fvBoundaryMesh&) = delete;
 
-        //- Add fvPatches corresponding to the given polyBoundaryMesh
-        void addPatches(const polyBoundaryMesh&);
-
 
 protected:
 
@@ -94,24 +97,21 @@ public:
 
     // Member Functions
 
-        // Access
+        //- Return the mesh reference
+        const fvMesh& mesh() const
+        {
+            return mesh_;
+        }
 
-            //- Return the mesh reference
-            const fvMesh& mesh() const
-            {
-                return mesh_;
-            }
+        //- Return a list of pointers for each patch
+        //  with only those pointing to interfaces being set
+        lduInterfacePtrsList interfaces() const;
 
-            //- Return a list of pointers for each patch
-            //  with only those pointing to interfaces being set
-            lduInterfacePtrsList interfaces() const;
-
-            //- Find patch index given a name
-            label findPatchID(const word& patchName) const;
-
-            //- Find patch indices given a name
-            labelList findIndices(const keyType&, const bool useGroups) const;
+        //- Return patch indices for all matches.
+        labelList indices(const keyType& key, const bool useGroups) const;
 
+        //- Find patch index given a name
+        label findPatchID(const word& patchName) const;
 
         //- Correct patches after moving points
         void movePoints();
@@ -127,6 +127,15 @@ public:
 
         //- Return reference to fvPatch by name.
         fvPatch& operator[](const word&);
+
+
+    // Housekeeping
+
+        //- Identical to the indices() method (AUG-2018)
+        labelList findIndices(const keyType& key, const bool useGroups) const
+        {
+            return indices(key, useGroups);
+        }
 };
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/patchInteractionDataList.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/patchInteractionDataList.C
index 79b0f11fe5..a6e9c63f1b 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/patchInteractionDataList.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/patchInteractionDataList.C
@@ -49,7 +49,7 @@ Foam::patchInteractionDataList::patchInteractionDataList
     forAllReverse(items, i)
     {
         const keyType& patchName = items[i].patchName();
-        labelList ids = bMesh.findIndices(patchName);
+        labelList ids = bMesh.indices(patchName);
 
         if (ids.empty())
         {
diff --git a/src/meshTools/coordinateSystems/coordinateSystems.C b/src/meshTools/coordinateSystems/coordinateSystems.C
index 08902cdb67..5a91c86fd1 100644
--- a/src/meshTools/coordinateSystems/coordinateSystems.C
+++ b/src/meshTools/coordinateSystems/coordinateSystems.C
@@ -45,7 +45,7 @@ namespace Foam
     template<class UnaryMatchPredicate>
     static wordList namesImpl
     (
-        const IOPtrList<coordinateSystem>& list,
+        const PtrList<coordinateSystem>& list,
         const UnaryMatchPredicate& matcher,
         const bool doSort
     )
@@ -75,6 +75,55 @@ namespace Foam
         return output;
     }
 
+
+    // Templated implementation for indices() - file-scope
+    template<class UnaryMatchPredicate>
+    static labelList indicesImpl
+    (
+        const PtrList<coordinateSystem>& list,
+        const UnaryMatchPredicate& matcher
+    )
+    {
+        const label len = list.size();
+
+        labelList output(len);
+
+        label count = 0;
+        for (label i = 0; i < len; ++i)
+        {
+            if (matcher(list[i].name()))
+            {
+                output[count++] = i;
+            }
+        }
+
+        output.resize(count);
+
+        return output;
+    }
+
+
+    // Templated implementation for findIndex() - file-scope
+    template<class UnaryMatchPredicate>
+    label findIndexImpl
+    (
+        const PtrList<coordinateSystem>& list,
+        const UnaryMatchPredicate& matcher
+    )
+    {
+        const label len = list.size();
+
+        for (label i = 0; i < len; ++i)
+        {
+            if (matcher(list[i].name()))
+            {
+                return i;
+            }
+        }
+
+        return -1;
+    }
+
 } // End namespace Foam
 
 
@@ -143,32 +192,30 @@ const Foam::coordinateSystems& Foam::coordinateSystems::New
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::labelList Foam::coordinateSystems::findIndices(const keyType& key) const
+Foam::labelList Foam::coordinateSystems::indices(const keyType& key) const
 {
-    labelList indices;
     if (key.empty())
     {
-        // no-op
+        return labelList();
     }
     else if (key.isPattern())
     {
-        indices = findStrings(key, this->names());
+        // Match as regex
+        regExp matcher(key);
+        return indicesImpl(*this, matcher);
     }
     else
     {
-        indices.setSize(this->size());
-        label count = 0;
-        forAll(*this, i)
-        {
-            if (key == operator[](i).name())
-            {
-                indices[count++] = i;
-            }
-        }
-        indices.setSize(count);
+        // Compare as literal string
+        const word& matcher = key;
+        return indicesImpl(*this, matcher);
     }
+}
 
-    return indices;
+
+Foam::labelList Foam::coordinateSystems::indices(const wordRes& matcher) const
+{
+    return indicesImpl(*this, matcher);
 }
 
 
@@ -176,29 +223,27 @@ Foam::label Foam::coordinateSystems::findIndex(const keyType& key) const
 {
     if (key.empty())
     {
-        // no-op
+        return -1;
     }
-    else if (key.isPattern())
+
+    if (key.isPattern())
     {
-        labelList indices = findIndices(key);
-        if (!indices.empty())
-        {
-            return indices.first();  // first match
-        }
+        // Find as regex
+        regExp matcher(key);
+        return findIndexImpl(*this, matcher);
     }
     else
     {
-        forAll(*this, i)
-        {
-            if (key == operator[](i).name())
-            {
-                return i;
-            }
-        }
+        // Find as literal string
+        const word& matcher = key;
+        return findIndexImpl(*this, matcher);
     }
+}
 
-    // Not found
-    return -1;
+
+Foam::label Foam::coordinateSystems::findIndex(const wordRes& matcher) const
+{
+    return findIndexImpl(*this, matcher);
 }
 
 
@@ -210,25 +255,37 @@ bool Foam::coordinateSystems::found(const keyType& key) const
 
 Foam::wordList Foam::coordinateSystems::names() const
 {
-    wordList output(size());
+    const PtrList<coordinateSystem>& systems = *this;
 
-    forAll(*this, i)
+    wordList list(systems.size());
+
+    forAll(systems, i)
     {
-        output[i] = operator[](i).name();
+        list[i] = systems[i].name();
     }
 
-    return output;
+    return list;
 }
 
 
-Foam::wordList Foam::coordinateSystems::names(const keyType& matcher) const
+Foam::wordList Foam::coordinateSystems::names(const keyType& key) const
 {
-    return
-    (
-        matcher.isPattern()
-      ? namesImpl(*this, regExp(matcher), false)
-      : namesImpl(*this, matcher, false)
-    );
+    if (key.empty())
+    {
+        return wordList();
+    }
+    else if (key.isPattern())
+    {
+        // Find as regex
+        regExp matcher(key);
+        return namesImpl(*this, matcher, false);
+    }
+    else
+    {
+        // Find as literal string
+        const word& matcher = key;
+        return namesImpl(*this, matcher, false);
+    }
 }
 
 
diff --git a/src/meshTools/coordinateSystems/coordinateSystems.H b/src/meshTools/coordinateSystems/coordinateSystems.H
index 1e429b97c2..a982f3c30c 100644
--- a/src/meshTools/coordinateSystems/coordinateSystems.H
+++ b/src/meshTools/coordinateSystems/coordinateSystems.H
@@ -117,19 +117,25 @@ public:
     // Member Functions
 
         //- Find and return indices for all matches
-        labelList findIndices(const keyType& key) const;
+        labelList indices(const keyType& key) const;
+
+        //- Find and return indices for all matches
+        labelList indices(const wordRes& matcher) const;
 
         //- Find and return index for the first match, return -1 if not found
         label findIndex(const keyType& key) const;
 
-        //- Search for given key
+        //- Find and return index for the first match, return -1 if not found
+        label findIndex(const wordRes& matcher) const;
+
+        //- Search if given key exists
         bool found(const keyType& key) const;
 
         //- A list of the coordinate-system names
         wordList names() const;
 
         //- A list of the coordinate-system names satisfying the input matcher
-        wordList names(const keyType& matcher) const;
+        wordList names(const keyType& key) const;
 
         //- A list of the coordinate-system names satisfying the input matcher
         wordList names(const wordRe& matcher) const;
@@ -145,6 +151,15 @@ public:
 
         //- Write data
         bool writeData(Ostream&) const;
+
+
+    // Housekeeping
+
+        //- Identical to the indices() method (AUG-2018)
+        labelList findIndices(const keyType& key) const
+        {
+            return this->indices(key);
+        }
 };
 
 
diff --git a/src/sampling/sampledSurface/sampledPlane/sampledPlane.C b/src/sampling/sampledSurface/sampledPlane/sampledPlane.C
index 1e7799d7a8..d89c3e50ce 100644
--- a/src/sampling/sampledSurface/sampledPlane/sampledPlane.C
+++ b/src/sampling/sampledSurface/sampledPlane/sampledPlane.C
@@ -177,7 +177,10 @@ bool Foam::sampledPlane::update()
     }
 
 
-    labelList selectedCells(mesh().cellZones().findMatching(zoneKey_).toc());
+    labelList selectedCells
+    (
+        mesh().cellZones().selection(zoneKey_).sortedToc()
+    );
 
     bool fullMesh = returnReduce(selectedCells.empty(), andOp<bool>());
 
diff --git a/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.C b/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.C
index bb82d2b89b..e1b8369dff 100644
--- a/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.C
+++ b/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.C
@@ -175,7 +175,7 @@ bool Foam::surfMeshSamplePlane::update()
 
     labelList selectedCells
     (
-        mesh().cellZones().findMatching(zoneKey_).sortedToc()
+        mesh().cellZones().selection(zoneKey_).sortedToc()
     );
 
     bool fullMesh = returnReduce(selectedCells.empty(), andOp<bool>());
diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C
index 5e72e57950..85740809ae 100644
--- a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C
+++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C
@@ -306,11 +306,12 @@ void Foam::radiation::solarLoad::initialise(const dictionary& coeffs)
         verticalDir_ = (-g/mag(g)).value();
     }
 
-    includePatches_ = mesh_.boundaryMesh().findIndices(viewFactorWalls);
+    includePatches_ = mesh_.boundaryMesh().indices(viewFactorWalls);
 
     coeffs.lookup("useVFbeamToDiffuse") >> useVFbeamToDiffuse_;
 
     coeffs.lookup("spectralDistribution") >> spectralDistribution_;
+
     spectralDistribution_ =
         spectralDistribution_/sum(spectralDistribution_);
 
diff --git a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
index 891c759db9..d3ed841a75 100644
--- a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
+++ b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
@@ -53,11 +53,11 @@ void Foam::radiation::viewFactor::initialise()
 {
     const polyBoundaryMesh& coarsePatches = coarseMesh_.boundaryMesh();
 
-    selectedPatches_ = mesh_.boundaryMesh().findIndices(viewFactorWalls);
-    forAll(selectedPatches_, i)
+    selectedPatches_ = mesh_.boundaryMesh().indices(viewFactorWalls);
+
+    for (const label patchi : selectedPatches_)
     {
-        const label patchI = selectedPatches_[i];
-        nLocalCoarseFaces_ += coarsePatches[patchI].size();
+        nLocalCoarseFaces_ += coarsePatches[patchi].size();
     }
 
     if (debug)
diff --git a/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.C b/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.C
index c18aa6f5c9..798e923fa0 100644
--- a/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.C
+++ b/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.C
@@ -228,7 +228,8 @@ bool Foam::fv::multiphaseMangrovesSource::read(const dictionary& dict)
             const dictionary& modelDict = regionsDict.subDict(regionName);
 
             const word zoneName(modelDict.get<word>("cellZone"));
-            zoneIDs_[i] = mesh_.cellZones().findIndices(zoneName);
+
+            zoneIDs_[i] = mesh_.cellZones().indices(zoneName);
             if (zoneIDs_[i].empty())
             {
                 FatalErrorInFunction
diff --git a/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.C b/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.C
index 975537a226..b60a977d28 100644
--- a/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.C
+++ b/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.C
@@ -260,7 +260,7 @@ bool Foam::fv::multiphaseMangrovesTurbulenceModel::read(const dictionary& dict)
 
             const word zoneName(modelDict.get<word>("cellZone"));
 
-            zoneIDs_[i] = mesh_.cellZones().findIndices(zoneName);
+            zoneIDs_[i] = mesh_.cellZones().indices(zoneName);
             if (zoneIDs_[i].empty())
             {
                 FatalErrorInFunction