diff --git a/etc/caseDicts/postProcessing/catalyst/default.cfg b/etc/caseDicts/insitu/catalyst/catalyst.cfg
similarity index 91%
rename from etc/caseDicts/postProcessing/catalyst/default.cfg
rename to etc/caseDicts/insitu/catalyst/catalyst.cfg
index e9ace6b706c1c55b456b01faf4a564c82d13bce9..f3defd4b902e0daf1fdc77bfd3f2fad4eefbce1f 100644
--- a/etc/caseDicts/postProcessing/catalyst/default.cfg
+++ b/etc/caseDicts/insitu/catalyst/catalyst.cfg
@@ -5,7 +5,7 @@
 |   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
-// Insitu processing of finiteVolume fields with ParaView Catalyst
+// Insitu processing with ParaView Catalyst
 
 type            catalyst;
 libs            ("libcatalystFoam.so");
diff --git a/etc/caseDicts/insitu/catalyst/printChannels.py b/etc/caseDicts/insitu/catalyst/printChannels.py
new file mode 100644
index 0000000000000000000000000000000000000000..b9f4e61753dcb550cad35bcf625bf954e53ba73a
--- /dev/null
+++ b/etc/caseDicts/insitu/catalyst/printChannels.py
@@ -0,0 +1,68 @@
+from paraview.simple import *
+from paraview import coprocessing
+
+# The frequency to output everything
+outputfrequency = 1
+
+# Simply print out all channel names that our function object is producing
+
+# ----------------------- CoProcessor definition -----------------------
+
+def CreateCoProcessor():
+  def _CreatePipeline(coprocessor, datadescription):
+    class Pipeline:
+      for i in range(datadescription.GetNumberOfInputDescriptions()):
+        name = datadescription.GetInputDescriptionName(i)
+        input = coprocessor.CreateProducer(datadescription, name)
+        grid = input.GetClientSideObject().GetOutputDataObject(0)
+        print "Channel <" + name + "> is", grid.GetClassName()
+
+    return Pipeline()
+
+  class CoProcessor(coprocessing.CoProcessor):
+    def CreatePipeline(self, datadescription):
+      self.Pipeline = _CreatePipeline(self, datadescription)
+
+  return CoProcessor()
+
+#--------------------------------------------------------------
+# Global variables that will hold the pipeline for each timestep
+# Creating the CoProcessor object, doesn't actually create the ParaView pipeline.
+# It will be automatically setup when coprocessor.UpdateProducers() is called the
+# first time.
+coprocessor = CreateCoProcessor()
+
+#--------------------------------------------------------------
+# Enable Live-Visualizaton with ParaView
+coprocessor.EnableLiveVisualization(False)
+
+
+# ---------------------- Data Selection method ----------------------
+
+def RequestDataDescription(datadescription):
+    "Callback to populate the request for current timestep"
+    global coprocessor
+    if datadescription.GetForceOutput() == True or datadescription.GetTimeStep() % outputfrequency == 0:
+        # We are just going to request all fields and meshes from the simulation
+        # code/adaptor.
+        for i in range(datadescription.GetNumberOfInputDescriptions()):
+            datadescription.GetInputDescription(i).AllFieldsOn()
+            datadescription.GetInputDescription(i).GenerateMeshOn()
+        return
+
+    # setup requests for all inputs based on the requirements of the
+    # pipeline.
+    coprocessor.LoadRequestedData(datadescription)
+
+# ------------------------ Processing method ------------------------
+
+def DoCoProcessing(datadescription):
+    "Callback to do co-processing for current timestep"
+    global coprocessor
+
+    # Update the coprocessor by providing it the newly generated simulation data.
+    # If the pipeline hasn't been setup yet, this will setup the pipeline.
+    coprocessor.UpdateProducers(datadescription)
+
+    # Write output data, if appropriate.
+    coprocessor.WriteData(datadescription);
diff --git a/etc/caseDicts/insitu/catalyst/writeAll.py b/etc/caseDicts/insitu/catalyst/writeAll.py
new file mode 100644
index 0000000000000000000000000000000000000000..658ae6e5545829cec5f2a397eac74046b2bc9595
--- /dev/null
+++ b/etc/caseDicts/insitu/catalyst/writeAll.py
@@ -0,0 +1,77 @@
+from paraview.simple import *
+from paraview import coprocessing
+
+# The frequency to output everything
+outputfrequency = 5
+
+# This is largely identical to the Catalyst allinputsgridwriter.py example
+# but only handle vtkMultiBlockDataSet, since that is what we generate
+
+# ----------------------- CoProcessor definition -----------------------
+
+def CreateCoProcessor():
+  def _CreatePipeline(coprocessor, datadescription):
+    class Pipeline:
+      for i in range(datadescription.GetNumberOfInputDescriptions()):
+        name  = datadescription.GetInputDescriptionName(i)
+        input = coprocessor.CreateProducer(datadescription, name)
+        grid  = input.GetClientSideObject().GetOutputDataObject(0)
+        if grid.IsA('vtkMultiBlockDataSet'):
+          writer = servermanager.writers.XMLMultiBlockDataWriter(Input=input)
+          coprocessor.RegisterWriter(writer, filename=name+'_%t.vtm', freq=outputfrequency)
+
+    return Pipeline()
+
+  class CoProcessor(coprocessing.CoProcessor):
+    def CreatePipeline(self, datadescription):
+      self.Pipeline = _CreatePipeline(self, datadescription)
+
+  return CoProcessor()
+
+#--------------------------------------------------------------
+# Global variables that will hold the pipeline for each timestep
+# Creating the CoProcessor object, doesn't actually create the ParaView pipeline.
+# It will be automatically setup when coprocessor.UpdateProducers() is called the
+# first time.
+coprocessor = CreateCoProcessor()
+
+#--------------------------------------------------------------
+# Enable Live-Visualizaton with ParaView
+coprocessor.EnableLiveVisualization(False)
+
+
+# ---------------------- Data Selection method ----------------------
+
+def RequestDataDescription(datadescription):
+    "Callback to populate the request for current timestep"
+    global coprocessor
+    if datadescription.GetForceOutput() == True or datadescription.GetTimeStep() % outputfrequency == 0:
+        # We are just going to request all fields and meshes from the simulation
+        # code/adaptor.
+        for i in range(datadescription.GetNumberOfInputDescriptions()):
+            datadescription.GetInputDescription(i).AllFieldsOn()
+            datadescription.GetInputDescription(i).GenerateMeshOn()
+        return
+
+    # setup requests for all inputs based on the requirements of the
+    # pipeline.
+    coprocessor.LoadRequestedData(datadescription)
+
+# ------------------------ Processing method ------------------------
+
+def DoCoProcessing(datadescription):
+    "Callback to do co-processing for current timestep"
+    global coprocessor
+
+    # Update the coprocessor by providing it the newly generated simulation data.
+    # If the pipeline hasn't been setup yet, this will setup the pipeline.
+    coprocessor.UpdateProducers(datadescription)
+
+    # Write output data, if appropriate.
+    coprocessor.WriteData(datadescription);
+
+    # Write image capture (Last arg: rescale lookup table), if appropriate.
+    coprocessor.WriteImages(datadescription, rescale_lookuptable=False)
+
+    # Live Visualization, if enabled.
+    coprocessor.DoLiveVisualization(datadescription, "localhost", 22222)
diff --git a/etc/caseDicts/insitu/catalyst/writeMesh.py b/etc/caseDicts/insitu/catalyst/writeMesh.py
new file mode 100644
index 0000000000000000000000000000000000000000..a8517e26bbc80bb6696a8353b28a1d61e2e74704
--- /dev/null
+++ b/etc/caseDicts/insitu/catalyst/writeMesh.py
@@ -0,0 +1,78 @@
+from paraview.simple import *
+from paraview import coprocessing
+
+# The frequency to output everything
+outputfrequency = 5
+
+# As per writeAll, but only for "/mesh" sub-channels.
+
+# ----------------------- CoProcessor definition -----------------------
+
+def CreateCoProcessor():
+  def _CreatePipeline(coprocessor, datadescription):
+    class Pipeline:
+      for i in range(datadescription.GetNumberOfInputDescriptions()):
+        name = datadescription.GetInputDescriptionName(i)
+        if not name.endswith('/mesh'):
+          continue
+        input = coprocessor.CreateProducer(datadescription, name)
+        grid  = input.GetClientSideObject().GetOutputDataObject(0)
+        if grid.IsA('vtkMultiBlockDataSet'):
+          writer = servermanager.writers.XMLMultiBlockDataWriter(Input=input)
+          coprocessor.RegisterWriter(writer, filename=name+'_%t.vtm', freq=outputfrequency)
+
+    return Pipeline()
+
+  class CoProcessor(coprocessing.CoProcessor):
+    def CreatePipeline(self, datadescription):
+      self.Pipeline = _CreatePipeline(self, datadescription)
+
+  return CoProcessor()
+
+#--------------------------------------------------------------
+# Global variables that will hold the pipeline for each timestep
+# Creating the CoProcessor object, doesn't actually create the ParaView pipeline.
+# It will be automatically setup when coprocessor.UpdateProducers() is called the
+# first time.
+coprocessor = CreateCoProcessor()
+
+#--------------------------------------------------------------
+# Enable Live-Visualizaton with ParaView
+coprocessor.EnableLiveVisualization(False)
+
+
+# ---------------------- Data Selection method ----------------------
+
+def RequestDataDescription(datadescription):
+    "Callback to populate the request for current timestep"
+    global coprocessor
+    if datadescription.GetForceOutput() == True or datadescription.GetTimeStep() % outputfrequency == 0:
+        # We are just going to request all fields and meshes from the simulation
+        # code/adaptor.
+        for i in range(datadescription.GetNumberOfInputDescriptions()):
+            datadescription.GetInputDescription(i).AllFieldsOn()
+            datadescription.GetInputDescription(i).GenerateMeshOn()
+        return
+
+    # setup requests for all inputs based on the requirements of the
+    # pipeline.
+    coprocessor.LoadRequestedData(datadescription)
+
+# ------------------------ Processing method ------------------------
+
+def DoCoProcessing(datadescription):
+    "Callback to do co-processing for current timestep"
+    global coprocessor
+
+    # Update the coprocessor by providing it the newly generated simulation data.
+    # If the pipeline hasn't been setup yet, this will setup the pipeline.
+    coprocessor.UpdateProducers(datadescription)
+
+    # Write output data, if appropriate.
+    coprocessor.WriteData(datadescription);
+
+    # Write image capture (Last arg: rescale lookup table), if appropriate.
+    coprocessor.WriteImages(datadescription, rescale_lookuptable=False)
+
+    # Live Visualization, if enabled.
+    coprocessor.DoLiveVisualization(datadescription, "localhost", 22222)
diff --git a/etc/caseDicts/insitu/catalyst/writePatches.py b/etc/caseDicts/insitu/catalyst/writePatches.py
new file mode 100644
index 0000000000000000000000000000000000000000..404267de18e0d1a427e6793dbf20589c0eb55d4b
--- /dev/null
+++ b/etc/caseDicts/insitu/catalyst/writePatches.py
@@ -0,0 +1,78 @@
+from paraview.simple import *
+from paraview import coprocessing
+
+# The frequency to output everything
+outputfrequency = 5
+
+# As per writeAll, but only for "/patches" sub-channels.
+
+# ----------------------- CoProcessor definition -----------------------
+
+def CreateCoProcessor():
+  def _CreatePipeline(coprocessor, datadescription):
+    class Pipeline:
+      for i in range(datadescription.GetNumberOfInputDescriptions()):
+        name = datadescription.GetInputDescriptionName(i)
+        if not name.endswith('/patches'):
+          continue
+        input = coprocessor.CreateProducer(datadescription, name)
+        grid  = input.GetClientSideObject().GetOutputDataObject(0)
+        if grid.IsA('vtkMultiBlockDataSet'):
+          writer = servermanager.writers.XMLMultiBlockDataWriter(Input=input)
+          coprocessor.RegisterWriter(writer, filename=name+'_%t.vtm', freq=outputfrequency)
+
+    return Pipeline()
+
+  class CoProcessor(coprocessing.CoProcessor):
+    def CreatePipeline(self, datadescription):
+      self.Pipeline = _CreatePipeline(self, datadescription)
+
+  return CoProcessor()
+
+#--------------------------------------------------------------
+# Global variables that will hold the pipeline for each timestep
+# Creating the CoProcessor object, doesn't actually create the ParaView pipeline.
+# It will be automatically setup when coprocessor.UpdateProducers() is called the
+# first time.
+coprocessor = CreateCoProcessor()
+
+#--------------------------------------------------------------
+# Enable Live-Visualizaton with ParaView
+coprocessor.EnableLiveVisualization(False)
+
+
+# ---------------------- Data Selection method ----------------------
+
+def RequestDataDescription(datadescription):
+    "Callback to populate the request for current timestep"
+    global coprocessor
+    if datadescription.GetForceOutput() == True or datadescription.GetTimeStep() % outputfrequency == 0:
+        # We are just going to request all fields and meshes from the simulation
+        # code/adaptor.
+        for i in range(datadescription.GetNumberOfInputDescriptions()):
+            datadescription.GetInputDescription(i).AllFieldsOn()
+            datadescription.GetInputDescription(i).GenerateMeshOn()
+        return
+
+    # setup requests for all inputs based on the requirements of the
+    # pipeline.
+    coprocessor.LoadRequestedData(datadescription)
+
+# ------------------------ Processing method ------------------------
+
+def DoCoProcessing(datadescription):
+    "Callback to do co-processing for current timestep"
+    global coprocessor
+
+    # Update the coprocessor by providing it the newly generated simulation data.
+    # If the pipeline hasn't been setup yet, this will setup the pipeline.
+    coprocessor.UpdateProducers(datadescription)
+
+    # Write output data, if appropriate.
+    coprocessor.WriteData(datadescription);
+
+    # Write image capture (Last arg: rescale lookup table), if appropriate.
+    coprocessor.WriteImages(datadescription, rescale_lookuptable=False)
+
+    # Live Visualization, if enabled.
+    coprocessor.DoLiveVisualization(datadescription, "localhost", 22222)
diff --git a/etc/caseDicts/postProcessing/catalyst/area.cfg b/etc/caseDicts/postProcessing/catalyst/area.cfg
deleted file mode 100644
index 2f6850f99531553441118e215c8bd3eca5c86142..0000000000000000000000000000000000000000
--- a/etc/caseDicts/postProcessing/catalyst/area.cfg
+++ /dev/null
@@ -1,16 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  plus                                  |
-|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-// Insitu processing of finiteArea fields with ParaView Catalyst
-
-type            catalyst::area;
-libs            ("libcatalystFoam.so");
-
-executeControl  timeStep;
-writeControl    none;
-
-// ************************************************************************* //
diff --git a/etc/caseDicts/postProcessing/catalyst/cloud.cfg b/etc/caseDicts/postProcessing/catalyst/cloud.cfg
deleted file mode 100644
index d860e901a9003aa5716b796653b2c0e9ae1e19f1..0000000000000000000000000000000000000000
--- a/etc/caseDicts/postProcessing/catalyst/cloud.cfg
+++ /dev/null
@@ -1,16 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  plus                                  |
-|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-// Insitu processing of lagrangian clouds with ParaView Catalyst
-
-type            catalyst::cloud;
-libs            ("libcatalystFoam.so");
-
-executeControl  timeStep;
-writeControl    none;
-
-// ************************************************************************* //
diff --git a/etc/config.csh/example/paraview b/etc/config.csh/example/paraview
index 98d7ad3cd8350f82eb3979b88da498712ab5cf55..d2c0ba0b14e0b2cccf753a69dfd36d7ae3d06911 100644
--- a/etc/config.csh/example/paraview
+++ b/etc/config.csh/example/paraview
@@ -13,8 +13,8 @@
 #     config.csh/example/paraview
 #
 # Description
-#     Example of defining a different ParaView_VERSION but retaining
-#     the standard config.csh/paraview mechanism
+#     Example of defining a different ParaView_VERSION but retaining the
+#     standard config.csh/paraview mechanism
 #
 # Note
 #     This file could be copied to a user or site location, but should never
@@ -25,7 +25,8 @@
 
 set pv=5.5.0
 set pv=5.5.0-mpipy
+set qt=qt-5.9.0
 
-eval `foamEtcFile -csh -config -mode=o paraview -- ParaView_VERSION=$pv`
+eval `foamEtcFile -csh -config -mode=o paraview -- ParaView_VERSION=$pv ParaView_QT=$qt`
 
 #------------------------------------------------------------------------------
diff --git a/etc/config.csh/paraview b/etc/config.csh/paraview
index c3710b59fe09e9665a35fb6b7299e07c6048b458..19239e6980b67c478511e72294cab35cf39cf597 100644
--- a/etc/config.csh/paraview
+++ b/etc/config.csh/paraview
@@ -120,8 +120,16 @@ if ( $?ParaView_VERSION ) then
         #OBSOLETE? endif
 
         # QT libraries as required
+        # Set Qt5_DIR to root directory.
+        # Another possibility: "qtpaths --qt-version"
+
         set qtDir="$archDir/$ParaView_QT"
         if ( -d "$qtDir" ) then
+            switch ($ParaView_QT)
+            case *-qt*:
+                setenv Qt5_DIR $qtDir
+                breaksw
+            endsw
             foreach qtLibDir ("$qtDir/lib$WM_COMPILER_LIB_ARCH" "$qtDir/lib")
                 if ( -d "$qtLibDir" ) then
                     setenv LD_LIBRARY_PATH "${qtLibDir}:${LD_LIBRARY_PATH}"
diff --git a/etc/config.csh/unset b/etc/config.csh/unset
index 9d4a22b6b6da94ba56f49aec59c90cdce41d292e..92cef08d00d9bdb1b35c325aa71f590181fae781 100644
--- a/etc/config.csh/unset
+++ b/etc/config.csh/unset
@@ -111,6 +111,7 @@ unsetenv ParaView_INCLUDE_DIR
 unsetenv ParaView_VERSION
 unsetenv PV_PLUGIN_PATH
 unsetenv VTK_DIR
+unsetenv Qt5_DIR  # Perhaps only unset if it is in WM_THIRD_PARTY_DIR?
 
 #------------------------------------------------------------------------------
 # Unset other ThirdParty environment variables
diff --git a/etc/config.sh/example/paraview b/etc/config.sh/example/paraview
index c1526e8c545ac6760ab222bc2f48fd2ea04b3d92..fb91cc1699a7e6824f929419531d899e256947d8 100644
--- a/etc/config.sh/example/paraview
+++ b/etc/config.sh/example/paraview
@@ -13,8 +13,8 @@
 #     config.sh/example/paraview
 #
 # Description
-#     Example of defining a different ParaView_VERSION but retaining
-#     the standard config.sh/paraview mechanism
+#     Example of defining a different ParaView_VERSION but retaining the
+#     standard config.sh/paraview mechanism
 #
 # Note
 #     This file could be copied to a user or site location, but should never
@@ -25,7 +25,8 @@
 
 pv=5.5.0
 pv=5.5.0-mpipy
+qt=qt-5.9.0
 
-eval $(foamEtcFile -sh -config -mode=o paraview -- ParaView_VERSION=$pv)
+eval $(foamEtcFile -sh -config -mode=o paraview -- ParaView_VERSION=$pv ParaView_QT=$qt)
 
 #------------------------------------------------------------------------------
diff --git a/etc/config.sh/paraview b/etc/config.sh/paraview
index 553378c472949552a9e959f64a3c520d3eb0b47d..33a8b327557d4fde3ad8254d0ca32bdbaf448ae1 100644
--- a/etc/config.sh/paraview
+++ b/etc/config.sh/paraview
@@ -127,10 +127,16 @@ then
         #OBSOLETE?    export PYTHONPATH=$PYTHONPATH:${PYTHONPATH:+:}$pvPython:$pvLibDir
         #OBSOLETE? fi
 
-        # QT libraries as required
+        # QT libraries as required, and Qt5_DIR for the root directory.
+        # Another possibility: "qtpaths --qt-version"
         qtDir="$archDir/$ParaView_QT"
         if [ -d "$qtDir" ]
         then
+            case "$ParaView_QT" in
+            *-5*)
+                export Qt5_DIR=$qtDir
+                ;;
+            esac
             for qtLibDir in $qtDir/lib$WM_COMPILER_LIB_ARCH $qtDir/lib
             do
                 if [ -d "$qtLibDir" ]
diff --git a/etc/config.sh/unset b/etc/config.sh/unset
index 024ba6316074f71661b6b87cd82cb02c8aa9616c..09aa5855793706bc6cb8628eac43b1605857df86 100644
--- a/etc/config.sh/unset
+++ b/etc/config.sh/unset
@@ -97,7 +97,6 @@ then
     unset OPAL_PREFIX
 fi
 
-
 #------------------------------------------------------------------------------
 # Unset Ensight/ParaView-related environment variables
 
@@ -108,6 +107,12 @@ unset ParaView_VERSION
 unset PV_PLUGIN_PATH
 unset VTK_DIR
 
+# Undefine Qt5_DIR if set to one of the paths on foamOldDirs
+if [ -z "$($foamClean -env=Qt5_DIR "$foamOldDirs")" ]
+then
+    unset Qt5_DIR
+fi
+
 #------------------------------------------------------------------------------
 # Unset other ThirdParty environment variables
 
diff --git a/etc/controlDict b/etc/controlDict
index 98a281c5d3f3f6439fbc20500ac15b5dfb10ee66..704374a5fa7e6a4e17c230bc53d1c6067e96646f 100644
--- a/etc/controlDict
+++ b/etc/controlDict
@@ -39,8 +39,8 @@ InfoSwitches
     writeDictionaries 0;
     writeOptionalEntries 0;
 
-    // Write lagrangian "positions" file in v1706 format (at earlier)
-    writeLagrangianPositions 0;
+    // Write lagrangian "positions" file in v1706 format (and earlier)
+    writeLagrangianPositions 1;
 
     // Report hosts used (parallel)
     // - 0 = none
diff --git a/modules/avalanche b/modules/avalanche
index 99e68179fa946476f0be4f957f8d6a2b78887925..6e6e105844897d4bf780bbc8d14031bc827e4b04 160000
--- a/modules/avalanche
+++ b/modules/avalanche
@@ -1 +1 @@
-Subproject commit 99e68179fa946476f0be4f957f8d6a2b78887925
+Subproject commit 6e6e105844897d4bf780bbc8d14031bc827e4b04
diff --git a/modules/catalyst b/modules/catalyst
index 3c0a2e7959755a84f48e25bbe3436ec6437f7cf6..5828d4510816948b034aa9afdf0b7b05659a9f59 160000
--- a/modules/catalyst
+++ b/modules/catalyst
@@ -1 +1 @@
-Subproject commit 3c0a2e7959755a84f48e25bbe3436ec6437f7cf6
+Subproject commit 5828d4510816948b034aa9afdf0b7b05659a9f59
diff --git a/modules/cfmesh b/modules/cfmesh
index 8f6e65ae7c71a95948b53679e57a73899b1dc999..288f05e08f07e693d4222e7b84ea12430947e5bf 160000
--- a/modules/cfmesh
+++ b/modules/cfmesh
@@ -1 +1 @@
-Subproject commit 8f6e65ae7c71a95948b53679e57a73899b1dc999
+Subproject commit 288f05e08f07e693d4222e7b84ea12430947e5bf
diff --git a/src/OpenFOAM/db/error/error.C b/src/OpenFOAM/db/error/error.C
index 89e83c622348240f9dda64359d73a05f6da9f48f..6f618860d9b7a3c2d8339c1f031b5754d9c119ee 100644
--- a/src/OpenFOAM/db/error/error.C
+++ b/src/OpenFOAM/db/error/error.C
@@ -236,7 +236,7 @@ void Foam::error::exit(const int errNo)
     {
         Perr<< endl << *this << endl
             << "\nFOAM exiting\n" << endl;
-        ::exit(1);
+        ::exit(errNo);
     }
 }
 
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
index 4113b1d752b68430395af1e3d0231123cd9eb69a..ff44c9ab3761519fadb8665001008a86afb8698f 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.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.
@@ -86,6 +86,7 @@ void pow
     gf.oriented() = pow(gf1.oriented(), r);
 }
 
+
 template
 <
     class Type,
@@ -181,6 +182,7 @@ void sqr
     gf.oriented() = sqr(gf1.oriented());
 }
 
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 tmp
 <
@@ -217,6 +219,7 @@ sqr(const GeometricField<Type, PatchField, GeoMesh>& gf)
     return tSqr;
 }
 
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 tmp
 <
@@ -270,6 +273,7 @@ void magSqr
     gsf.oriented() = magSqr(gf.oriented());
 }
 
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
 (
@@ -343,6 +347,7 @@ void mag
     gsf.oriented() = mag(gf.oriented());
 }
 
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 tmp<GeometricField<scalar, PatchField, GeoMesh>> mag
 (
@@ -458,6 +463,7 @@ cmptAv(const GeometricField<Type, PatchField, GeoMesh>& gf)
     return CmptAv;
 }
 
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 tmp
 <
@@ -500,7 +506,7 @@ cmptAv(const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf)
 }
 
 
-#define UNARY_REDUCTION_FUNCTION_WITH_BOUNDARY(returnType, func, gFunc)        \
+#define UNARY_REDUCTION_FUNCTION_WITH_BOUNDARY(returnType, func, binaryOp)     \
                                                                                \
 template<class Type, template<class> class PatchField, class GeoMesh>          \
 dimensioned<returnType> func                                                   \
@@ -512,7 +518,15 @@ dimensioned<returnType> func                                                   \
     (                                                                          \
         #func "(" + gf.name() + ')',                                           \
         gf.dimensions(),                                                       \
-        Foam::func(gFunc(gf.primitiveField()), gFunc(gf.boundaryField()))      \
+        returnReduce                                                           \
+        (                                                                      \
+            Foam::func                                                         \
+            (                                                                  \
+                Foam::func(gf.primitiveField()),                               \
+                Foam::func(gf.boundaryField())                                 \
+            ),                                                                 \
+            binaryOp<Type>()                                                   \
+        )                                                                      \
     );                                                                         \
 }                                                                              \
                                                                                \
@@ -527,8 +541,8 @@ dimensioned<returnType> func                                                   \
     return res;                                                                \
 }
 
-UNARY_REDUCTION_FUNCTION_WITH_BOUNDARY(Type, max, gMax)
-UNARY_REDUCTION_FUNCTION_WITH_BOUNDARY(Type, min, gMin)
+UNARY_REDUCTION_FUNCTION_WITH_BOUNDARY(Type, max, maxOp)
+UNARY_REDUCTION_FUNCTION_WITH_BOUNDARY(Type, min, minOp)
 
 #undef UNARY_REDUCTION_FUNCTION_WITH_BOUNDARY
 
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.H
index a98f4a608340a11f99079f1d24f9424686f21fa2..63cdbdd8a9591929184a537ab2bd400252eeae29 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.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.
@@ -212,7 +212,7 @@ tmp
 cmptAv(const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf);
 
 
-#define UNARY_REDUCTION_FUNCTION_WITH_BOUNDARY(returnType, func, gFunc)        \
+#define UNARY_REDUCTION_FUNCTION_WITH_BOUNDARY(returnType, func, binaryOp)     \
                                                                                \
 template<class Type, template<class> class PatchField, class GeoMesh>          \
 dimensioned<returnType> func                                                   \
@@ -226,8 +226,8 @@ dimensioned<returnType> func                                                   \
     const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf1                 \
 );
 
-UNARY_REDUCTION_FUNCTION_WITH_BOUNDARY(Type, max, gMax)
-UNARY_REDUCTION_FUNCTION_WITH_BOUNDARY(Type, min, gMin)
+UNARY_REDUCTION_FUNCTION_WITH_BOUNDARY(Type, max, maxOp)
+UNARY_REDUCTION_FUNCTION_WITH_BOUNDARY(Type, min, minOp)
 
 #undef UNARY_REDUCTION_FUNCTION_WITH_BOUNDARY
 
diff --git a/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C b/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C
index cd2fac95cbef48dc00571352c4aca41ebb2900bd..8bf7480557cbe5b9ba5b790a894ee0f2b9d2bc32 100644
--- a/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C
+++ b/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C
@@ -266,7 +266,7 @@ Foam::faceZone::faceZone
 (
     const word& name,
     const labelUList& addr,
-    const boolList& fm,
+    const boolUList& fm,
     const label index,
     const faceZoneMesh& zm
 )
@@ -328,7 +328,7 @@ Foam::faceZone::faceZone
 (
     const faceZone& origZone,
     const labelUList& addr,
-    const boolList& fm,
+    const boolUList& fm,
     const label index,
     const faceZoneMesh& zm
 )
@@ -468,7 +468,7 @@ void Foam::faceZone::resetAddressing
 void Foam::faceZone::resetAddressing
 (
     const labelUList& addr,
-    const boolList& flipMap
+    const boolUList& flipMap
 )
 {
     clearAddressing();
diff --git a/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.H b/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.H
index 50a06cc512da6581efd5ba558822c4f98a63ac26..770e018b108d0b85229d6e43ca467ec9e0bfc3ef 100644
--- a/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.H
+++ b/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.H
@@ -177,7 +177,7 @@ public:
         (
             const word& name,
             const labelUList& addr,
-            const boolList& fm,
+            const boolUList& fm,
             const label index,
             const faceZoneMesh& zm
         );
@@ -208,7 +208,7 @@ public:
         (
             const faceZone& origZone,
             const labelUList& addr,
-            const boolList& fm,
+            const boolUList& fm,
             const label index,
             const faceZoneMesh& zm
         );
@@ -236,7 +236,7 @@ public:
         virtual autoPtr<faceZone> clone
         (
             const labelUList& addr,
-            const boolList& fm,
+            const boolUList& fm,
             const label index,
             const faceZoneMesh& zm
         ) const
@@ -309,7 +309,7 @@ public:
         virtual void resetAddressing
         (
             const labelUList& addr,
-            const boolList& flipMap
+            const boolUList& flipMap
         );
 
         //- Move reset addressing - use uniform flip map value
diff --git a/src/Pstream/mpi/Make/options b/src/Pstream/mpi/Make/options
index 43a48312c42c76a6f500f502d306a5c7e1016c94..8fcb7016f81a290e860556bfe163cac4d36838bd 100644
--- a/src/Pstream/mpi/Make/options
+++ b/src/Pstream/mpi/Make/options
@@ -1,5 +1,5 @@
 sinclude $(GENERAL_RULES)/mplib$(WM_MPLIB)
-sinclude $(RULES)/mplib$(WM_MPLIB)
+sinclude $(DEFAULT_RULES)/mplib$(WM_MPLIB)
 
 EXE_INC  = $(PFLAGS) $(PINC) $(c++LESSWARN)
 LIB_LIBS = $(PLIBS)
diff --git a/src/conversion/vtk/part/foamVtuSizingTemplates.C b/src/conversion/vtk/part/foamVtuSizingTemplates.C
index 4bfda81ac636f81ef536340979c8e2cb7407429e..b070ab67da4b002dbd5a52ced2d6cede56917907 100644
--- a/src/conversion/vtk/part/foamVtuSizingTemplates.C
+++ b/src/conversion/vtk/part/foamVtuSizingTemplates.C
@@ -571,7 +571,8 @@ void Foam::vtk::vtuSizing::populateArrays
             if (output == contentType::LEGACY)
             {
                 // Update size for legacy face stream
-                faceOutput[startLabel] = (faceIndexer - startLabel);
+                // (subtract 1 to avoid counting the storage location)
+                faceOutput[startLabel] = (faceIndexer - 1 - startLabel);
             }
             else
             {
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
index 4206ca325b54de63034dd01e268e0de48fb17289..d232a5455975709384c98e952d44b4d1e28eb88f 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -42,6 +42,7 @@ License
 #include "fvMeshTools.H"
 #include "labelPairHashes.H"
 #include "ListOps.H"
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -688,25 +689,54 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::repatch
 // merge those points.
 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
 (
+    const labelList& pointToGlobalMaster,
     labelListList& constructPointMap
 )
 {
     // Find out which sets of points get merged and create a map from
     // mesh point to unique point.
-    Map<label> pointToMaster
-    (
-        fvMeshAdder::findSharedPoints
-        (
-            mesh_,
-            mergeTol_
-        )
-    );
+
+    label nShared = 0;
+    forAll(pointToGlobalMaster, pointi)
+    {
+        if (pointToGlobalMaster[pointi] != -1)
+        {
+            nShared++;
+        }
+    }
+
+    Map<label> globalMasterToLocalMaster(2*nShared);
+    Map<label> pointToMaster(2*nShared);
+
+    forAll(pointToGlobalMaster, pointi)
+    {
+        label globali = pointToGlobalMaster[pointi];
+        if (globali != -1)
+        {
+            Map<label>::const_iterator iter = globalMasterToLocalMaster.find
+            (
+                globali
+            );
+
+            if (iter == globalMasterToLocalMaster.end())
+            {
+                // Found first point. Designate as master
+                globalMasterToLocalMaster.insert(globali, pointi);
+                pointToMaster.insert(pointi, pointi);
+            }
+            else
+            {
+                pointToMaster.insert(pointi, iter());
+            }
+        }
+    }
 
     if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
     {
         return nullptr;
     }
 
+
     polyTopoChange meshMod(mesh_);
 
     fvMeshAdder::mergePoints(mesh_, pointToMaster, meshMod);
@@ -750,16 +780,19 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
 }
 
 
-// Construct the local environment of all boundary faces.
-void Foam::fvMeshDistribute::getNeighbourData
+void Foam::fvMeshDistribute::getCouplingData
 (
     const labelList& distribution,
     labelList& sourceFace,
     labelList& sourceProc,
     labelList& sourcePatch,
-    labelList& sourceNewNbrProc
+    labelList& sourceNewNbrProc,
+    labelList& sourcePointMaster
 ) const
 {
+    // Construct the coupling information for all (boundary) faces and
+    // points
+
     label nBnd = mesh_.nFaces() - mesh_.nInternalFaces();
     sourceFace.setSize(nBnd);
     sourceProc.setSize(nBnd);
@@ -890,13 +923,62 @@ void Foam::fvMeshDistribute::getNeighbourData
             }
         }
     }
+
+
+    // Collect coupled (collocated) points
+    sourcePointMaster.setSize(mesh_.nPoints());
+    sourcePointMaster = -1;
+    {
+        // Assign global master point
+        const globalIndex globalPoints(mesh_.nPoints());
+
+        const globalMeshData& gmd = mesh_.globalData();
+        const indirectPrimitivePatch& cpp = gmd.coupledPatch();
+        const labelList& meshPoints = cpp.meshPoints();
+        const mapDistribute& slavesMap = gmd.globalCoPointSlavesMap();
+        const labelListList& slaves = gmd.globalCoPointSlaves();
+
+        labelList elems(slavesMap.constructSize(), -1);
+        forAll(meshPoints, pointi)
+        {
+            const labelList& slots = slaves[pointi];
+
+            if (slots.size())
+            {
+                // pointi is a master. Assign a unique label.
+
+                label globalPointi = globalPoints.toGlobal(meshPoints[pointi]);
+                elems[pointi] = globalPointi;
+                forAll(slots, i)
+                {
+                    label sloti = slots[i];
+                    if (sloti >= meshPoints.size())
+                    {
+                        // Filter out local collocated points. We don't want
+                        // to merge these
+                        elems[slots[i]] = globalPointi;
+                    }
+                }
+            }
+        }
+
+        // Push slave-slot data back to slaves
+        slavesMap.reverseDistribute(elems.size(), elems, false);
+
+        // Extract back onto mesh
+        forAll(meshPoints, pointi)
+        {
+            sourcePointMaster[meshPoints[pointi]] = elems[pointi];
+        }
+    }
 }
 
 
 // Subset the neighbourCell/neighbourProc fields
-void Foam::fvMeshDistribute::subsetBoundaryData
+void Foam::fvMeshDistribute::subsetCouplingData
 (
     const fvMesh& mesh,
+    const labelList& pointMap,
     const labelList& faceMap,
     const labelList& cellMap,
 
@@ -909,11 +991,13 @@ void Foam::fvMeshDistribute::subsetBoundaryData
     const labelList& sourceProc,
     const labelList& sourcePatch,
     const labelList& sourceNewNbrProc,
+    const labelList& sourcePointMaster,
 
     labelList& subFace,
     labelList& subProc,
     labelList& subPatch,
-    labelList& subNewNbrProc
+    labelList& subNewNbrProc,
+    labelList& subPointMaster
 )
 {
     subFace.setSize(mesh.nFaces() - mesh.nInternalFaces());
@@ -959,6 +1043,9 @@ void Foam::fvMeshDistribute::subsetBoundaryData
             subNewNbrProc[newBFacei] = sourceNewNbrProc[oldBFacei];
         }
     }
+
+
+    subPointMaster = UIndirectList<label>(sourcePointMaster, pointMap);
 }
 
 
@@ -1042,9 +1129,9 @@ Foam::labelList Foam::fvMeshDistribute::mapBoundaryData
 (
     const primitiveMesh& mesh,      // mesh after adding
     const mapAddedPolyMesh& map,
-    const labelList& boundaryData0, // mesh before adding
+    const labelList& boundaryData0, // on mesh before adding
     const label nInternalFaces1,
-    const labelList& boundaryData1  // added mesh
+    const labelList& boundaryData1  // on added mesh
 )
 {
     labelList newBoundaryData(mesh.nFaces() - mesh.nInternalFaces());
@@ -1076,6 +1163,41 @@ Foam::labelList Foam::fvMeshDistribute::mapBoundaryData
 }
 
 
+Foam::labelList Foam::fvMeshDistribute::mapPointData
+(
+    const primitiveMesh& mesh,      // mesh after adding
+    const mapAddedPolyMesh& map,
+    const labelList& boundaryData0, // on mesh before adding
+    const labelList& boundaryData1  // on added mesh
+)
+{
+    labelList newBoundaryData(mesh.nPoints());
+
+    forAll(boundaryData0, oldPointi)
+    {
+        label newPointi = map.oldPointMap()[oldPointi];
+
+        // Point still exists (is necessary?)
+        if (newPointi >= 0)
+        {
+            newBoundaryData[newPointi] = boundaryData0[oldPointi];
+        }
+    }
+
+    forAll(boundaryData1, addedPointi)
+    {
+        label newPointi = map.addedPointMap()[addedPointi];
+
+        if (newPointi >= 0)
+        {
+            newBoundaryData[newPointi] = boundaryData1[addedPointi];
+        }
+    }
+
+    return newBoundaryData;
+}
+
+
 // Remove cells. Add all exposed faces to patch oldInternalPatchi
 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::doRemoveCells
 (
@@ -1297,6 +1419,7 @@ void Foam::fvMeshDistribute::sendMesh
     const labelList& sourceProc,
     const labelList& sourcePatch,
     const labelList& sourceNewNbrProc,
+    const labelList& sourcePointMaster,
     Ostream& toDomain
 )
 {
@@ -1433,7 +1556,8 @@ void Foam::fvMeshDistribute::sendMesh
         << sourceFace
         << sourceProc
         << sourcePatch
-        << sourceNewNbrProc;
+        << sourceNewNbrProc
+        << sourcePointMaster;
 
 
     if (debug)
@@ -1456,6 +1580,7 @@ Foam::autoPtr<Foam::fvMesh> Foam::fvMeshDistribute::receiveMesh
     labelList& domainSourceProc,
     labelList& domainSourcePatch,
     labelList& domainSourceNewNbrProc,
+    labelList& domainSourcePointMaster,
     Istream& fromNbr
 )
 {
@@ -1474,7 +1599,8 @@ Foam::autoPtr<Foam::fvMesh> Foam::fvMeshDistribute::receiveMesh
         >> domainSourceFace
         >> domainSourceProc
         >> domainSourcePatch
-        >> domainSourceNewNbrProc;
+        >> domainSourceNewNbrProc
+        >> domainSourcePointMaster;
 
     // Construct fvMesh
     auto domainMeshPtr = autoPtr<fvMesh>::New
@@ -1709,13 +1835,15 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
     labelList sourceFace;
     labelList sourceProc;
     labelList sourceNewNbrProc;
-    getNeighbourData
+    labelList sourcePointMaster;
+    getCouplingData
     (
         distribution,
         sourceFace,
         sourceProc,
         sourcePatch,
-        sourceNewNbrProc
+        sourceNewNbrProc,
+        sourcePointMaster
     );
 
 
@@ -1924,11 +2052,13 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
             labelList procSourceProc;
             labelList procSourcePatch;
             labelList procSourceNewNbrProc;
+            labelList procSourcePointMaster;
 
-            subsetBoundaryData
+            subsetCouplingData
             (
                 subsetter.subMesh(),
-                subsetter.faceMap(),        // from subMesh to mesh
+                subsetter.pointMap(),       // from subMesh to mesh
+                subsetter.faceMap(),        //      ,,      ,,
                 subsetter.cellMap(),        //      ,,      ,,
 
                 distribution,               // old mesh distribution
@@ -1940,15 +2070,16 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
                 sourceProc,
                 sourcePatch,
                 sourceNewNbrProc,
+                sourcePointMaster,
 
                 procSourceFace,
                 procSourceProc,
                 procSourcePatch,
-                procSourceNewNbrProc
+                procSourceNewNbrProc,
+                procSourcePointMaster
             );
 
 
-
             // Send to neighbour
             sendMesh
             (
@@ -1963,6 +2094,8 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
                 procSourceProc,
                 procSourcePatch,
                 procSourceNewNbrProc,
+                procSourcePointMaster,
+
                 str
             );
 
@@ -2121,10 +2254,12 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
         labelList domainSourceProc;
         labelList domainSourcePatch;
         labelList domainSourceNewNbrProc;
+        labelList domainSourcePointMaster;
 
-        subsetBoundaryData
+        subsetCouplingData
         (
             mesh_,                          // new mesh
+            subMap().pointMap(),            // from new to original mesh
             subMap().faceMap(),             // from new to original mesh
             subMap().cellMap(),
 
@@ -2137,17 +2272,20 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
             sourceProc,
             sourcePatch,
             sourceNewNbrProc,
+            sourcePointMaster,
 
             domainSourceFace,
             domainSourceProc,
             domainSourcePatch,
-            domainSourceNewNbrProc
+            domainSourceNewNbrProc,
+            domainSourcePointMaster
         );
 
         sourceFace.transfer(domainSourceFace);
         sourceProc.transfer(domainSourceProc);
         sourcePatch.transfer(domainSourcePatch);
         sourceNewNbrProc.transfer(domainSourceNewNbrProc);
+        sourcePointMaster.transfer(domainSourcePointMaster);
     }
 
 
@@ -2205,6 +2343,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
             labelList domainSourceProc;
             labelList domainSourcePatch;
             labelList domainSourceNewNbrProc;
+            labelList domainSourcePointMaster;
 
             autoPtr<fvMesh> domainMeshPtr;
 
@@ -2241,6 +2380,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
                     domainSourceProc,
                     domainSourcePatch,
                     domainSourceNewNbrProc,
+                    domainSourcePointMaster,
                     str
                 );
                 fvMesh& domainMesh = domainMeshPtr();
@@ -2508,6 +2648,15 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
                 domainMesh.nInternalFaces(),
                 domainSourceNewNbrProc
             );
+            // Update pointMaster data
+            sourcePointMaster = mapPointData
+            (
+                mesh_,
+                map(),
+                sourcePointMaster,
+                domainSourcePointMaster
+            );
+
 
             // Update all addressing so xxProcAddressing points to correct
             // item in masterMesh.
@@ -2621,6 +2770,10 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
     }
 
 
+    // See if any originally shared points need to be merged. Note: does
+    // parallel comms. After this points and edges should again be consistent.
+    mergeSharedPoints(sourcePointMaster, constructPointMap);
+
 
     // Add processorPatches
     // ~~~~~~~~~~~~~~~~~~~~
@@ -2650,16 +2803,8 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
 
     // Change patches. Since this might change ordering of coupled faces
     // we also need to adapt our constructMaps.
-    // NOTE: there is one very particular problem with this structure.
-    // We first create the processor patches and use these to merge out
-    // shared points (see mergeSharedPoints below). So temporarily points
-    // and edges do not match!
     repatch(newPatchID, constructFaceMap);
 
-    // See if any geometrically shared points need to be merged. Note: does
-    // parallel comms. After this points and edges should again be consistent.
-    mergeSharedPoints(constructPointMap);
-
     // Bit of hack: processorFvPatchField does not get reset since created
     // from nothing so explicitly reset.
     initPatchFields<volScalarField, processorFvPatchField<scalar>>
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
index 752dab30459ab807893466ec9f5c77358838ba40..e68e7a534a83902bf97c9d03e78572389d29bcbf 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -163,11 +163,11 @@ class fvMeshDistribute
                 labelListList& constructFaceMap
             );
 
-            //- Merge any shared points that are geometrically shared. Needs
-            //  parallel valid mesh - uses globalMeshData.
+            //- Merge any local points that were remotely coupled.
             //  constructPointMap is adapted for the new point labels.
             autoPtr<mapPolyMesh> mergeSharedPoints
             (
+                const labelList& pointToGlobalMaster,
                 labelListList& constructPointMap
             );
 
@@ -175,19 +175,21 @@ class fvMeshDistribute
         // Coupling information
 
             //- Construct the local environment of all boundary faces.
-            void getNeighbourData
+            void getCouplingData
             (
                 const labelList& distribution,
                 labelList& sourceFace,
                 labelList& sourceProc,
                 labelList& sourcePatch,
-                labelList& sourceNewProc
+                labelList& sourceNewProc,
+                labelList& sourcePointMaster
             ) const;
 
             // Subset the neighbourCell/neighbourProc fields
-            static void subsetBoundaryData
+            static void subsetCouplingData
             (
                 const fvMesh& mesh,
+                const labelList& pointMap,
                 const labelList& faceMap,
                 const labelList& cellMap,
 
@@ -200,11 +202,13 @@ class fvMeshDistribute
                 const labelList& sourceProc,
                 const labelList& sourcePatch,
                 const labelList& sourceNewProc,
+                const labelList& sourcePointMaster,
 
                 labelList& subFace,
                 labelList& subProc,
                 labelList& subPatch,
-                labelList& subNewProc
+                labelList& subNewProc,
+                labelList& subPointMaster
             );
 
             //- Find cells on mesh whose faceID/procID match the neighbour
@@ -237,6 +241,14 @@ class fvMeshDistribute
                 const labelList& boundaryData1  // added mesh
             );
 
+            //- Map point data to new mesh (resulting from adding two meshes)
+            static labelList mapPointData
+            (
+                const primitiveMesh& mesh,      // mesh after adding
+                const mapAddedPolyMesh& map,
+                const labelList& boundaryData0, // on mesh before adding
+                const labelList& boundaryData1  // on added mesh
+            );
 
         // Other
 
@@ -276,6 +288,7 @@ class fvMeshDistribute
                 const labelList& sourceProc,
                 const labelList& sourcePatch,
                 const labelList& sourceNewProc,
+                const labelList& sourcePointMaster,
                 Ostream& toDomain
             );
             //- Send subset of fields
@@ -300,6 +313,7 @@ class fvMeshDistribute
                 labelList& domainSourceProc,
                 labelList& domainSourcePatch,
                 labelList& domainSourceNewProc,
+                labelList& domainSourcePointMaster,
                 Istream& fromNbr
             );
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/swirlFlowRateInletVelocity/swirlFlowRateInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/swirlFlowRateInletVelocity/swirlFlowRateInletVelocityFvPatchVectorField.C
index ebe4ebdd58a4ba90db0ba11e609bb918683008c7..89c1ce000847bfcdd5549e95b22eb99e0c382bec 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/swirlFlowRateInletVelocity/swirlFlowRateInletVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/swirlFlowRateInletVelocity/swirlFlowRateInletVelocityFvPatchVectorField.C
@@ -75,7 +75,7 @@ swirlFlowRateInletVelocityFvPatchVectorField
         dict.lookupOrDefault
         (
             "axis",
-            patch().size()
+            returnReduce(patch().size(), maxOp<label>())
           ? -gSum(patch().Sf())/gSum(patch().magSf())
           : Zero
         )
@@ -145,48 +145,53 @@ void Foam::swirlFlowRateInletVelocityFvPatchVectorField::updateCoeffs()
     {
         return;
     }
-
-    const scalar t = this->db().time().timeOutputValue();
-    const scalar flowRate = flowRate_->value(t);
-    const scalar rpm = rpm_->value(t);
-
     const scalar totArea = gSum(patch().magSf());
-    const scalar avgU = -flowRate/totArea;
-
-    const vector axisHat = axis_/mag(axis_);
 
-    // Update angular velocity - convert [rpm] to [rad/s]
-    tmp<vectorField> tangentialVelocity
-    (
-        axisHat ^ (rpm*constant::mathematical::pi/30.0)*(patch().Cf() - origin_)
-    );
+    if (totArea > ROOTVSMALL && axis_ != vector(Zero))
+    {
+        const scalar t = this->db().time().timeOutputValue();
+        const scalar flowRate = flowRate_->value(t);
+        const scalar rpm = rpm_->value(t);
 
-    tmp<vectorField> n = patch().nf();
+        const scalar avgU = -flowRate/totArea;
 
-    const surfaceScalarField& phi =
-        db().lookupObject<surfaceScalarField>(phiName_);
+        const vector axisHat = axis_/mag(axis_);
 
-    if (phi.dimensions() == dimVelocity*dimArea)
-    {
-        // volumetric flow-rate
-        operator==(tangentialVelocity + n*avgU);
-    }
-    else if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
-    {
-        const fvPatchField<scalar>& rhop =
-            patch().lookupPatchField<volScalarField, scalar>(rhoName_);
-
-        // mass flow-rate
-        operator==(tangentialVelocity + n*avgU/rhop);
-    }
-    else
-    {
-        FatalErrorInFunction
-            << "dimensions of " << phiName_ << " are incorrect" << nl
-            << "    on patch " << this->patch().name()
-            << " of field " << this->internalField().name()
-            << " in file " << this->internalField().objectPath()
-            << nl << exit(FatalError);
+        // Update angular velocity - convert [rpm] to [rad/s]
+        tmp<vectorField> tangentialVelocity
+        (
+            axisHat
+           ^(rpm*constant::mathematical::pi/30.0)
+           *(patch().Cf() - origin_)
+        );
+
+        tmp<vectorField> n = patch().nf();
+
+        const surfaceScalarField& phi =
+            db().lookupObject<surfaceScalarField>(phiName_);
+
+        if (phi.dimensions() == dimVelocity*dimArea)
+        {
+            // volumetric flow-rate
+            operator==(tangentialVelocity + n*avgU);
+        }
+        else if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
+        {
+            const fvPatchField<scalar>& rhop =
+                patch().lookupPatchField<volScalarField, scalar>(rhoName_);
+
+            // mass flow-rate
+            operator==(tangentialVelocity + n*avgU/rhop);
+        }
+        else
+        {
+            FatalErrorInFunction
+                << "dimensions of " << phiName_ << " are incorrect" << nl
+                << "    on patch " << this->patch().name()
+                << " of field " << this->internalField().name()
+                << " in file " << this->internalField().objectPath()
+                << nl << exit(FatalError);
+        }
     }
 
     fixedValueFvPatchField<vector>::updateCoeffs();
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
index 1c328d3f47bf0c64f4ee32dee96da7b58d3f6c7a..20a045815cfe9fde91774ee6c4ae08cc482e9cde 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
@@ -261,7 +261,7 @@ EulerDdtScheme<Type>::fvcDdt
             new GeometricField<Type, fvPatchField, volMesh>
             (
                 ddtIOobject,
-                rDeltaT.value()*
+                rDeltaT*
                 (
                     alpha()
                    *rho()
diff --git a/src/functionObjects/field/turbulenceFields/turbulenceFields.C b/src/functionObjects/field/turbulenceFields/turbulenceFields.C
index 054ad322cbea005b24396ab82ca189c1913a0501..0fc33a74e74da5e1d795baf556e001a8a3e7f8b6 100644
--- a/src/functionObjects/field/turbulenceFields/turbulenceFields.C
+++ b/src/functionObjects/field/turbulenceFields/turbulenceFields.C
@@ -61,7 +61,8 @@ Foam::functionObjects::turbulenceFields::compressibleFieldNames_
     { compressibleField::cfAlphaEff, "alphaEff" },
     { compressibleField::cfR, "R" },
     { compressibleField::cfDevRhoReff, "devRhoReff" },
-    { compressibleField::cfL, "L" }
+    { compressibleField::cfL, "L" },
+    { compressibleField::cfI, "I" }
 };
 
 
@@ -80,6 +81,7 @@ Foam::functionObjects::turbulenceFields::incompressibleFieldNames_
     { incompressibleField::ifR, "R" },
     { incompressibleField::ifDevReff, "devReff" },
     { incompressibleField::ifL, "L" },
+    { incompressibleField::ifI, "I" }
 };
 
 
@@ -236,6 +238,11 @@ bool Foam::functionObjects::turbulenceFields::execute()
                     processField<scalar>(f, L(model));
                     break;
                 }
+                case cfI:
+                {
+                    processField<scalar>(f, I(model));
+                    break;
+                }
                 default:
                 {
                     FatalErrorInFunction
@@ -298,6 +305,11 @@ bool Foam::functionObjects::turbulenceFields::execute()
                     processField<scalar>(f, L(model));
                     break;
                 }
+                case ifI:
+                {
+                    processField<scalar>(f, I(model));
+                    break;
+                }
                 default:
                 {
                     FatalErrorInFunction
diff --git a/src/functionObjects/field/turbulenceFields/turbulenceFields.H b/src/functionObjects/field/turbulenceFields/turbulenceFields.H
index a4ae31d19970328e82d8f0788effc3671d2ce1b2..881a772f3bfc52bb393b44709f013f8594b6de51 100644
--- a/src/functionObjects/field/turbulenceFields/turbulenceFields.H
+++ b/src/functionObjects/field/turbulenceFields/turbulenceFields.H
@@ -76,6 +76,7 @@ Usage
         devReff     | Deviatoric part of the effective Reynolds stress
         devRhoReff  | Divergence of the Reynolds stress
         L           | turbulence length scale
+        I           | turbulence intensity
     \endplaintable
 
 See also
@@ -125,7 +126,8 @@ public:
         cfAlphaEff,
         cfR,
         cfDevRhoReff,
-        cfL
+        cfL,
+        cfI
     };
     static const Enum<compressibleField> compressibleFieldNames_;
 
@@ -139,7 +141,8 @@ public:
         ifNuEff,
         ifR,
         ifDevReff,
-        ifL
+        ifL,
+        ifI
     };
     static const Enum<incompressibleField> incompressibleFieldNames_;
 
@@ -179,6 +182,10 @@ protected:
         template<class Model>
         tmp<volScalarField> L(const Model& model) const;
 
+        //- Return I calculated from k and U
+        template<class Model>
+        tmp<volScalarField> I(const Model& model) const;
+
 
 private:
 
diff --git a/src/functionObjects/field/turbulenceFields/turbulenceFieldsTemplates.C b/src/functionObjects/field/turbulenceFields/turbulenceFieldsTemplates.C
index ec6b70d92dfe6320d6205d28adb774d0f0e80589..4ded70c6be41bcf266ffe605634bce6026d784f7 100644
--- a/src/functionObjects/field/turbulenceFields/turbulenceFieldsTemplates.C
+++ b/src/functionObjects/field/turbulenceFields/turbulenceFieldsTemplates.C
@@ -85,19 +85,16 @@ Foam::functionObjects::turbulenceFields::omega
     const volScalarField k(model.k());
     const volScalarField epsilon(model.epsilon());
 
-    return tmp<volScalarField>
+    return tmp<volScalarField>::New
     (
-        new volScalarField
+        IOobject
         (
-            IOobject
-            (
-                "omega.tmp",
-                k.mesh().time().timeName(),
-                k.mesh()
-            ),
-            epsilon/(Cmu*k),
-            epsilon.boundaryField().types()
-        )
+            "omega.tmp",
+            k.mesh().time().timeName(),
+            k.mesh()
+        ),
+        epsilon/(Cmu*k),
+        epsilon.boundaryField().types()
     );
 }
 
@@ -109,9 +106,10 @@ Foam::functionObjects::turbulenceFields::nuTilda
     const Model& model
 ) const
 {
-    return tmp<volScalarField>
+    return tmp<volScalarField>::New
     (
-        new volScalarField("nuTilda.tmp", model.k()/omega(model))
+        "nuTilda.tmp",
+        model.k()/omega(model)
     );
 }
 
@@ -130,15 +128,30 @@ Foam::functionObjects::turbulenceFields::L
     const volScalarField epsilon(model.epsilon());
     const dimensionedScalar eps0("eps0", epsilon.dimensions(), SMALL);
 
-    return tmp<volScalarField>
+    return tmp<volScalarField>::New
     (
-        new volScalarField
-        (
-            "L.tmp",
-            pow(Cmu, 0.75)*pow(k, 1.5)/(epsilon + eps0)
-        )
+        "L.tmp",
+        pow(Cmu, 0.75)*pow(k, 1.5)/(epsilon + eps0)
     );
 }
 
 
+template<class Model>
+Foam::tmp<Foam::volScalarField>
+Foam::functionObjects::turbulenceFields::I
+(
+    const Model& model
+) const
+{
+    // Assume k is available
+    const volScalarField uPrime(sqrt((2.0/3.0)*model.k()));
+    const dimensionedScalar U0("U0", dimVelocity, SMALL);
+
+    return tmp<volScalarField>::New
+    (
+        "I.tmp",
+        uPrime/max(max(uPrime, mag(model.U())), U0)
+    );
+}
+
 // ************************************************************************* //
diff --git a/src/functionObjects/forces/forces/forces.C b/src/functionObjects/forces/forces/forces.C
index f1a945b302af34da830047d5f1c7e8450f68ebf0..84015def879d2fa54d0e94f30419d99b178232f6 100644
--- a/src/functionObjects/forces/forces/forces.C
+++ b/src/functionObjects/forces/forces/forces.C
@@ -353,7 +353,7 @@ Foam::functionObjects::forces::devRhoReff() const
         const dictionary& transportProperties =
             lookupObject<dictionary>("transportProperties");
 
-        dimensionedScalar nu(transportProperties.lookup("nu"));
+        dimensionedScalar nu("nu", dimViscosity, transportProperties);
 
         const volVectorField& U = lookupObject<volVectorField>(UName_);
 
diff --git a/src/functionObjects/utilities/vtkWrite/vtkWrite.C b/src/functionObjects/utilities/vtkWrite/vtkWrite.C
index f5d28883ca3a533d478660a011722f2ffd210e3f..c0d66f29785dd5372bd7eb3879c93de5562e10c1 100644
--- a/src/functionObjects/utilities/vtkWrite/vtkWrite.C
+++ b/src/functionObjects/utilities/vtkWrite/vtkWrite.C
@@ -108,6 +108,7 @@ bool Foam::functionObjects::vtkWrite::read(const dictionary& dict)
     //
     dict.readIfPresent("directory", dirName_);
 
+    decompose_ = dict.lookupOrDefault("decompose", false);
     writeIds_ = dict.lookupOrDefault("writeIds", false);
 
 
@@ -185,7 +186,7 @@ bool Foam::functionObjects::vtkWrite::write()
         (
             mesh_,
             writeOpts_,
-            true  // decompose
+            decompose_
         );
 
         // Write mesh
diff --git a/src/functionObjects/utilities/vtkWrite/vtkWrite.H b/src/functionObjects/utilities/vtkWrite/vtkWrite.H
index cf69d91492adfb9d89f4d1753e79ad9401f3a847..fca1d3c605c6b30d9836af6d5f056bf1a03a0ee0 100644
--- a/src/functionObjects/utilities/vtkWrite/vtkWrite.H
+++ b/src/functionObjects/utilities/vtkWrite/vtkWrite.H
@@ -44,6 +44,7 @@ Description
         writeInterval   1;
         format          binary;
         legacy          false;
+        decompose       false;
         ...
         fields          (U p);
     }
@@ -51,14 +52,15 @@ Description
 
 Usage
     \table
-        Property     | Description             | Required    | Default value
-        type         | Type name: vtkWrite     | yes         |
-        fields       | Fields to output        | yes         |
-        writeControl | Output control          | recommended | timeStep
-        directory    | The output directory name | no        | "VTK"
-        format       | ASCII or binary format  | no          | binary
-        legacy       | Legacy VTK output       | no          | false
-        writeIds     | Write cell ids as field | no          | true
+        Property     | Description                      | Required    | Default
+        type         | Type name: vtkWrite              | yes         |
+        fields       | Fields to output                 | yes         |
+        writeControl | Output control                   | recommended | timeStep
+        directory    | The output directory name        | no          | "VTK"
+        format       | ASCII or binary format           | no          | binary
+        legacy       | Legacy VTK output                | no          | false
+        decompose    | decompose polyhedra              | no          | false
+        writeIds     | Write cell ids as field          | no          | true
     \endtable
 
 See also
@@ -106,6 +108,9 @@ class vtkWrite
         //- Output directory name
         fileName dirName_;
 
+        //- Decompose polyhedra
+        bool decompose_;
+
         //- Write cell ids field
         bool writeIds_;
 
@@ -119,7 +124,11 @@ class vtkWrite
 
         //- Write selected fields for GeoField type.
         template<class GeoField>
-        label writeFields(vtk::internalWriter& writer, bool verbose=true) const;
+        label writeFields
+        (
+            vtk::internalWriter& writer,
+            bool verbose=true
+        ) const;
 
 
         //- Write selected fields for GeoField type.
@@ -131,10 +140,10 @@ class vtkWrite
         ) const;
 
 
-        //- Disallow default bitwise copy construct
+        //- No copy construct
         vtkWrite(const vtkWrite&) = delete;
 
-        //- Disallow default bitwise assignment
+        //- No copy assignment
         void operator=(const vtkWrite&) = delete;
 
 
@@ -150,7 +159,7 @@ public:
         vtkWrite
         (
             const word& name,
-            const Time& t,
+            const Time& runTime,
             const dictionary& dict
         );
 
diff --git a/src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.H b/src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.H
index e56970999cf007c8c48b59a7de98ad56862ec63e..69618baabf2bec00b207a0b21611cc933f2812be 100644
--- a/src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.H
+++ b/src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,7 @@ License
 #include "ParticleErosion.H"
 #include "ParticleTracks.H"
 #include "ParticleTrap.H"
+#include "PatchCollisionDensity.H"
 #include "PatchPostProcessing.H"
 #include "VoidFraction.H"
 
@@ -48,6 +49,7 @@ License
     makeCloudFunctionObjectType(ParticleErosion, CloudType);                   \
     makeCloudFunctionObjectType(ParticleTracks, CloudType);                    \
     makeCloudFunctionObjectType(ParticleTrap, CloudType);                      \
+    makeCloudFunctionObjectType(PatchCollisionDensity, CloudType);             \
     makeCloudFunctionObjectType(PatchPostProcessing, CloudType);               \
     makeCloudFunctionObjectType(VoidFraction, CloudType);
 
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchCollisionDensity/PatchCollisionDensity.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchCollisionDensity/PatchCollisionDensity.C
new file mode 100644
index 0000000000000000000000000000000000000000..111ffcbedf4e19ae3d0897dc1041af7d4814d6bf
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchCollisionDensity/PatchCollisionDensity.C
@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PatchCollisionDensity.H"
+#include "Pstream.H"
+#include "stringListOps.H"
+#include "ListOps.H"
+#include "ListListOps.H"
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::PatchCollisionDensity<CloudType>::write()
+{
+    const scalarField z(this->owner().mesh().nCells(), 0);
+
+    volScalarField
+    (
+        IOobject
+        (
+            this->owner().name() + ":collisionDensity",
+            this->owner().mesh().time().timeName(),
+            this->owner().mesh()
+        ),
+        this->owner().mesh(),
+        dimless/dimArea,
+        z,
+        collisionDensity_
+    )
+   .write();
+
+    volScalarField
+    (
+        IOobject
+        (
+            this->owner().name() + ":collisionDensityRate",
+            this->owner().mesh().time().timeName(),
+            this->owner().mesh()
+        ),
+        this->owner().mesh(),
+        dimless/dimArea/dimTime,
+        z,
+        (collisionDensity_ - collisionDensity0_)
+       /(this->owner().mesh().time().value() - time0_)
+    )
+   .write();
+
+    collisionDensity0_ == collisionDensity_;
+    time0_ = this->owner().mesh().time().value();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::PatchCollisionDensity<CloudType>::PatchCollisionDensity
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& modelName
+)
+:
+    CloudFunctionObject<CloudType>(dict, owner, modelName, typeName),
+    minSpeed_(dict.lookupOrDefault<scalar>("minSpeed", -1)),
+    collisionDensity_
+    (
+        this->owner().mesh().boundary(),
+        volScalarField::Internal::null(),
+        calculatedFvPatchField<scalar>::typeName
+    ),
+    collisionDensity0_
+    (
+        this->owner().mesh().boundary(),
+        volScalarField::Internal::null(),
+        calculatedFvPatchField<scalar>::typeName
+    ),
+    time0_(this->owner().mesh().time().value())
+{
+    collisionDensity_ == 0;
+    collisionDensity0_ == 0;
+
+    IOobject io
+    (
+        this->owner().name() + ":collisionDensity",
+        this->owner().mesh().time().timeName(),
+        this->owner().mesh(),
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (io.typeHeaderOk<volScalarField>())
+    {
+        const volScalarField collisionDensity(io, this->owner().mesh());
+        collisionDensity_ == collisionDensity.boundaryField();
+        collisionDensity0_ == collisionDensity.boundaryField();
+    }
+}
+
+
+template<class CloudType>
+Foam::PatchCollisionDensity<CloudType>::PatchCollisionDensity
+(
+    const PatchCollisionDensity<CloudType>& ppm
+)
+:
+    CloudFunctionObject<CloudType>(ppm),
+    minSpeed_(ppm.minSpeed_),
+    collisionDensity_(ppm.collisionDensity_),
+    collisionDensity0_(ppm.collisionDensity0_),
+    time0_(ppm.time0_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::PatchCollisionDensity<CloudType>::~PatchCollisionDensity()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::PatchCollisionDensity<CloudType>::postPatch
+(
+    const parcelType& p,
+    const polyPatch& pp,
+    bool&
+)
+{
+    const label patchi = pp.index();
+    const label patchFacei = p.face() - pp.start();
+
+    vector nw, Up;
+    this->owner().patchData(p, pp, nw, Up);
+
+    const scalar speed = (p.U() - Up) & nw;
+    if (speed > minSpeed_)
+    {
+        collisionDensity_[patchi][patchFacei] +=
+            1/this->owner().mesh().magSf().boundaryField()[patchi][patchFacei];
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchCollisionDensity/PatchCollisionDensity.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchCollisionDensity/PatchCollisionDensity.H
new file mode 100644
index 0000000000000000000000000000000000000000..8cb4d76bb49cca5ebffc3b5bef9ce64b60e1c4d7
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchCollisionDensity/PatchCollisionDensity.H
@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::PatchCollisionDensity
+
+Description
+    Function object which generates fields of the number and rate of collisions
+    per unit area on all patches. Can optionally take a minimum speed below
+    which a collision is not counted.
+
+    Example usage:
+    \verbatim
+    patchCollisionDensity1
+    {
+        type        patchCollisionDensity;
+        minSpeed    1e-3;
+    }
+    \endverbatim
+
+SourceFiles
+    PatchCollisionDensity.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PatchCollisionDensity_H
+#define PatchCollisionDensity_H
+
+#include "CloudFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class PatchCollisionDensity Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class PatchCollisionDensity
+:
+    public CloudFunctionObject<CloudType>
+{
+    // Private data
+
+        typedef typename CloudType::particleType parcelType;
+
+        //- The threshold for a collision
+        const scalar minSpeed_;
+
+        //- The field of the number of collisions per unit area
+        volScalarField::Boundary collisionDensity_;
+
+        //- The field of the number of collisions per unit area at the last
+        //  output
+        volScalarField::Boundary collisionDensity0_;
+
+        //- The time at the last output
+        scalar time0_;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Write post-processing info
+        void write();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("patchCollisionDensity");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        PatchCollisionDensity
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& modelName
+        );
+
+        //- Construct copy
+        PatchCollisionDensity(const PatchCollisionDensity<CloudType>& ppm);
+
+        //- Construct and return a clone
+        virtual autoPtr<CloudFunctionObject<CloudType>> clone() const
+        {
+            return autoPtr<CloudFunctionObject<CloudType>>
+            (
+                new PatchCollisionDensity<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~PatchCollisionDensity();
+
+
+    // Member Functions
+
+        // Evaluation
+
+            //- Post-patch hook
+            virtual void postPatch
+            (
+                const parcelType& p,
+                const polyPatch& pp,
+                bool& keepParticle
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "PatchCollisionDensity.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lumpedPointMotion/lumpedPointDisplacementPointPatchVectorField.C b/src/lumpedPointMotion/lumpedPointDisplacementPointPatchVectorField.C
index 3368d6f59f58e35cd69833be51921811d313e343..1ac55bd8b86a539de5eead5c201ec3e06e56a50f 100644
--- a/src/lumpedPointMotion/lumpedPointDisplacementPointPatchVectorField.C
+++ b/src/lumpedPointMotion/lumpedPointDisplacementPointPatchVectorField.C
@@ -250,7 +250,7 @@ void Foam::lumpedPointDisplacementPointPatchVectorField::updateCoeffs()
 
             if (Pstream::master())
             {
-                movement().writeData(forces, moments);
+                movement().writeData(forces, moments, &(db().time()));
 
                 // Signal external source to execute
                 movement().coupler().useSlave();
diff --git a/src/lumpedPointMotion/lumpedPointMovement.C b/src/lumpedPointMotion/lumpedPointMovement.C
index 4e74f34a2f733893753895fefa18e1347aaf0160..00c34c7209875ecb7b6c13c032389af6539f644a 100644
--- a/src/lumpedPointMotion/lumpedPointMovement.C
+++ b/src/lumpedPointMotion/lumpedPointMovement.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,29 +48,65 @@ Foam::lumpedPointMovement::formatNames
 };
 
 
+const Foam::Enum
+<
+    Foam::lumpedPointMovement::scalingType
+>
+Foam::lumpedPointMovement::scalingNames
+{
+    { scalingType::LENGTH, "plain" },
+    { scalingType::FORCE, "force" },
+    { scalingType::MOMENT, "moment" }
+};
+
+
 const Foam::word
 Foam::lumpedPointMovement::dictionaryName("lumpedPointMovement");
 
 
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
 namespace Foam
 {
+
+//! \cond fileScope
+//- Space-separated vector value (ASCII)
+static inline Ostream& putPlain(Ostream& os, const vector& val)
+{
+    os  << val.x() << ' ' << val.y() << ' ' << val.z();
+    return os;
+}
+
+
+//! \cond fileScope
+//- Space-separated vector value (ASCII)
+static inline Ostream& putTime(Ostream& os, const Time& t)
+{
+    os  <<"Time index=" << t.timeIndex()
+        << " value=" << t.timeOutputValue();
+
+    return os;
+}
+
+
 //! \cond fileScope
 //- Write list content with size, bracket, content, bracket one-per-line.
 //  This makes for consistent for parsing, regardless of the list length.
 template <class T>
-static void writeList(Ostream& os, const string& header, const UList<T>& lst)
+static void writeList(Ostream& os, const string& header, const UList<T>& list)
 {
+    const label len = list.size();
+
     // Header string
     os  << header.c_str() << nl;
 
     // Write size and start delimiter
-    os  << lst.size() << nl
-        << token::BEGIN_LIST << nl;
+    os  << len << nl << token::BEGIN_LIST << nl;
 
     // Write contents
-    forAll(lst, i)
+    for (label i=0; i < len; ++i)
     {
-        os << lst[i] << nl;
+        os << list[i] << nl;
     }
 
     // Write end delimiter
@@ -165,8 +201,11 @@ Foam::lumpedPointMovement::lumpedPointMovement()
     coupler_(),
     inputName_("positions.in"),
     outputName_("forces.out"),
+    logName_("movement.log"),
     inputFormat_(lumpedPointState::inputFormatType::DICTIONARY),
     outputFormat_(outputFormatType::DICTIONARY),
+    scaleInput_(-1.0),
+    scaleOutput_(-1.0),
     state0_(),
     state_(),
     thresholdPtr_(0),
@@ -198,6 +237,11 @@ Foam::lumpedPointMovement::lumpedPointMovement
     autoCentre_(true),
     forcesDict_(),
     coupler_(),
+    inputName_("positions.in"),
+    outputName_("forces.out"),
+    logName_("movement.log"),
+    scaleInput_(-1.0),
+    scaleOutput_(-1.0),
     state0_(),
     state_(),
     thresholdPtr_(0),
@@ -262,6 +306,7 @@ void Foam::lumpedPointMovement::readDict(const dictionary& dict)
 
     commDict.lookup("inputName")  >> inputName_;
     commDict.lookup("outputName") >> outputName_;
+    commDict.readIfPresent("logName", logName_);
 
     inputFormat_ = lumpedPointState::formatNames.lookup
     (
@@ -274,6 +319,47 @@ void Foam::lumpedPointMovement::readDict(const dictionary& dict)
         "outputFormat",
         commDict
     );
+
+    scaleInput_  = -1;
+    scaleOutput_ = -1;
+
+    const dictionary* scaleDict = nullptr;
+
+    if ((scaleDict = commDict.subDictPtr("scaleInput")))
+    {
+        for (int i=0; i < scaleInput_.size(); ++i)
+        {
+            const word& key = scalingNames[scalingType(i)];
+
+            if
+            (
+                scaleDict->readIfPresent(key, scaleInput_[i])
+             && scaleInput_[i] > 0
+            )
+            {
+                Info<<"Using input " << key << " multiplier: "
+                    << scaleInput_[i] << nl;
+            }
+        }
+    }
+
+    if ((scaleDict = commDict.subDictPtr("scaleOutput")))
+    {
+        for (int i=0; i < scaleOutput_.size(); ++i)
+        {
+            const word& key = scalingNames[scalingType(i)];
+
+            if
+            (
+                scaleDict->readIfPresent(key, scaleOutput_[i])
+             && scaleOutput_[i] > 0
+            )
+            {
+                Info<<"Using output " << key << " multiplier: "
+                    << scaleOutput_[i] << nl;
+            }
+        }
+    }
 }
 
 
@@ -638,6 +724,8 @@ bool Foam::lumpedPointMovement::readState()
         coupler().resolveFile(inputName_)
     );
 
+    state_.scalePoints(scaleInput_[scalingType::LENGTH]);
+
     state_.relax(relax_, prev);
 
     return status;
@@ -646,45 +734,114 @@ bool Foam::lumpedPointMovement::readState()
 
 bool Foam::lumpedPointMovement::writeData
 (
-    const UList<vector>& forces
+    Ostream& os,
+    const UList<vector>& forces,
+    const UList<vector>& moments,
+    const outputFormatType fmt,
+    const Time* timeinfo
 ) const
 {
-    if (!Pstream::master())
+    const bool writeMoments = (moments.size() == forces.size());
+
+    if (fmt == outputFormatType::PLAIN)
     {
-        return false;
-    }
+        os  <<"########" << nl;
+        if (timeinfo)
+        {
+            os <<"# ";
+            putTime(os, *timeinfo) << nl;
+        }
+        os  <<"# size=" << this->size() << nl
+            <<"# columns (points) (forces)";
 
-    const fileName output(coupler().resolveFile(outputName_));
-    OFstream os(output); // ASCII
+        if (writeMoments)
+        {
+            os << " (moments)";
+        }
 
-    if (outputFormat_ == outputFormatType::PLAIN)
-    {
-        os  <<"# output from OpenFOAM" << nl
-            <<"# N, points, forces" << nl
-            << this->size() << nl;
+        os << nl;
+
+        bool report = false;
+        scalar scaleLength = scaleOutput_[scalingType::LENGTH];
+        scalar scaleForce  = scaleOutput_[scalingType::FORCE];
+        scalar scaleMoment = scaleOutput_[scalingType::MOMENT];
+
+        if (scaleLength > 0)
+        {
+            report = true;
+        }
+        else
+        {
+            scaleLength = 1.0;
+        }
+
+        if (scaleForce > 0)
+        {
+            report = true;
+        }
+        else
+        {
+            scaleForce = 1.0;
+        }
+
+        if (writeMoments)
+        {
+            if (scaleMoment > 0)
+            {
+                report = true;
+            }
+            else
+            {
+                scaleMoment = 1.0;
+            }
+        }
+
+        if (report)
+        {
+            os  <<"# scaling points=" << scaleLength
+                <<" forces=" << scaleForce;
+
+            if (writeMoments)
+            {
+                os  <<" moments=" << scaleMoment;
+            }
+
+            os << nl;
+        }
 
-        const char* zeroVector = "0 0 0";
+        os <<"########" << nl;
 
         forAll(locations_, i)
         {
-            const vector pos = locations_[i] * axis_;
+            const vector pos = scaleLength * (locations_[i] * axis_);
 
-            os  << pos.x() << ' '
-                << pos.y() << ' '
-                << pos.z() << ' ';
+            putPlain(os, pos) << ' ';
 
             if (i < forces.size())
             {
-                os  << forces[i].x() << ' '
-                    << forces[i].y() << ' '
-                    << forces[i].z();
+                const vector val(scaleForce * forces[i]);
+                putPlain(os, val);
             }
             else
             {
-                os << zeroVector;
+                putPlain(os, vector::zero);
             }
 
-            os  << nl;
+            if (writeMoments)
+            {
+                os << ' ';
+                if (i < moments.size())
+                {
+                    const vector val(scaleMoment * moments[i]);
+                    putPlain(os, val);
+                }
+                else
+                {
+                    putPlain(os, vector::zero);
+                }
+            }
+
+            os << nl;
         }
     }
     else
@@ -693,10 +850,21 @@ bool Foam::lumpedPointMovement::writeData
         // - exclude the usual OpenFOAM 'FoamFile' header
         // - ensure lists have consistent format
 
-        os  <<"// output from OpenFOAM" << nl << nl;
+        os  <<"////////" << nl;
+        if (timeinfo)
+        {
+            os <<"// ";
+            putTime(os, *timeinfo) << nl;
+        }
+        os  << nl;
 
         writeList(os, "points", (locations_*axis_)());
         writeList(os, "forces", forces);
+
+        if (writeMoments)
+        {
+            writeList(os, "moments", moments);
+        }
     }
 
     return true;
@@ -706,7 +874,8 @@ bool Foam::lumpedPointMovement::writeData
 bool Foam::lumpedPointMovement::writeData
 (
     const UList<vector>& forces,
-    const UList<vector>& moments
+    const UList<vector>& moments,
+    const Time* timeinfo
 ) const
 {
     if (!Pstream::master())
@@ -714,60 +883,28 @@ bool Foam::lumpedPointMovement::writeData
         return false;
     }
 
-    const fileName output(coupler().resolveFile(outputName_));
-    OFstream os(output); // ASCII
-
-    if (outputFormat_ == outputFormatType::PLAIN)
+    // Regular output
     {
-        os  <<"# output from OpenFOAM" << nl
-            <<"# N, points, forces, moments" << nl
-            << this->size() << nl;
-
-        const char* zeroVector = "0 0 0";
-
-        forAll(locations_, i)
-        {
-            const vector pos = locations_[i] * axis_;
+        const fileName output(coupler().resolveFile(outputName_));
+        OFstream os(output, IOstream::ASCII);
 
-            os  << pos.x() << ' '
-                << pos.y() << ' '
-                << pos.z() << ' ';
-
-            if (i < forces.size())
-            {
-                os  << forces[i].x() << ' '
-                    << forces[i].y() << ' '
-                    << forces[i].z() << ' ';
-            }
-            else
-            {
-                os << zeroVector << ' ';
-            }
-
-            if (i < moments.size())
-            {
-                os  << moments[i].x() << ' '
-                    << moments[i].y() << ' '
-                    << moments[i].z();
-            }
-            else
-            {
-                os  << zeroVector;
-            }
-            os  << nl;
-        }
+        writeData(os, forces, moments, outputFormat_, timeinfo);
     }
-    else
+
+    // Log output
     {
-        // Make it easier for external programs to parse
-        // - exclude the usual OpenFOAM 'FoamFile' header
-        // - ensure lists have consistent format
+        const fileName output(coupler().resolveFile(logName_));
 
-        os  <<"// output from OpenFOAM" << nl << nl;
+        OFstream os
+        (
+            output,
+            IOstream::ASCII,
+            IOstream::currentVersion,
+            IOstream::UNCOMPRESSED,
+            true // append mode
+        );
 
-        writeList(os, "points", (locations_*axis_)());
-        writeList(os, "forces", forces);
-        writeList(os, "moments", moments);
+        writeData(os, forces, moments, outputFormatType::PLAIN, timeinfo);
     }
 
     return true;
diff --git a/src/lumpedPointMotion/lumpedPointMovement.H b/src/lumpedPointMotion/lumpedPointMovement.H
index dc0f1fb189676a10bb4a46925c13634b0776e35d..7587ec8378cf8d2e743b25e4d6c4f3a98bc0b070 100644
--- a/src/lumpedPointMotion/lumpedPointMovement.H
+++ b/src/lumpedPointMotion/lumpedPointMovement.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -60,8 +60,10 @@ SourceFiles
 
 namespace Foam
 {
+
 // Forward declarations
 class polyMesh;
+class Time;
 
 /*---------------------------------------------------------------------------*\
                      Class lumpedPointMovement Declaration
@@ -78,11 +80,22 @@ public:
         DICTIONARY
     };
 
+    //- Output format types
+    enum scalingType
+    {
+        LENGTH = 0,
+        FORCE,
+        MOMENT
+    };
+
     // Static data
 
     //- Names for the output format types
     static const Enum<outputFormatType> formatNames;
 
+    //- Names for the scaling types
+    static const Enum<scalingType> scalingNames;
+
 
 private:
 
@@ -125,9 +138,15 @@ private:
         //- File io
         word inputName_;
         word outputName_;
+        word logName_;
+
         lumpedPointState::inputFormatType inputFormat_;
         outputFormatType outputFormat_;
 
+        //- Optional scale factors for input/output files
+        FixedList<scalar, 1> scaleInput_;
+        FixedList<scalar, 3> scaleOutput_;
+
 
     // Demand-driven private data
 
@@ -246,6 +265,9 @@ public:
         //- The output (forces) file name
         inline const word& outputName() const;
 
+        //- The log file name
+        inline const word& logName() const;
+
         //- The input (state) file format
         inline lumpedPointState::inputFormatType inputFormat() const;
 
@@ -324,21 +346,24 @@ public:
         //- Write axis, locations, division as a dictionary
         void writeDict(Ostream& os) const;
 
-
-        //- Write points, forces
+        //- Write points, forces, moments. Only call from the master process
         bool writeData
         (
-            const UList<vector>& forces
+            Ostream& os,
+            const UList<vector>& forces,
+            const UList<vector>& moments,
+            const outputFormatType fmt = outputFormatType::PLAIN,
+            const Time* timeinfo = nullptr
         ) const;
 
         //- Write points, forces, moments
         bool writeData
         (
             const UList<vector>& forces,
-            const UList<vector>& moments
+            const UList<vector>& moments = List<vector>(),
+            const Time* timeinfo = nullptr
         ) const;
 
-
         //- Read state from file, applying relaxation as requested
         bool readState();
 
diff --git a/src/lumpedPointMotion/lumpedPointMovement.dict b/src/lumpedPointMotion/lumpedPointMovement.dict
new file mode 100644
index 0000000000000000000000000000000000000000..5223b7d253a14e5ed2dc7ec3460cee9fb27d1a61
--- /dev/null
+++ b/src/lumpedPointMotion/lumpedPointMovement.dict
@@ -0,0 +1,100 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      lumpedPointMovement;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Reference axis for the locations
+axis            (0 0 1);
+
+// Locations of the lumped points
+locations       11(0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5);
+
+// Division for pressure forces (0-1)
+division        0.5;
+
+//- If present, the offset of patch points compared to the locations
+//  Otherwise determined from the bounding box
+// centre          (0 0 0);
+
+//- The interpolation scheme
+interpolationScheme linear;
+
+//- Relaxation/scaling factor when updating positions
+relax           1.0;
+
+
+forces
+{
+    //- The pressure name (default: p)
+    p           p;
+
+    //- Reference pressure [Pa] (default: 0)
+    pRef        0;
+
+    //- Reference density for incompressible calculations (default: 1)
+    rhoRef      1;
+}
+
+
+communication
+{
+    commsDir        "comms";
+
+    log             on;
+
+    waitInterval    1;
+
+    timeOut         100;
+
+    initByExternal  false;
+
+    // Input file of positions/rotation, written by external application
+    inputName       positions.in;
+
+    // Output file of forces, written by OpenFOAM
+    outputName      forces.out;
+
+    // Log of points/forces/moments during the simulation
+    logName         movement.log;
+
+    inputFormat     dictionary;
+    outputFormat    dictionary;
+
+    debugTable      "$FOAM_CASE/output.txt";
+
+
+    // Scaling applied to values read from 'inputName'
+    scaleInput
+    {
+        //- Length multiplier (to metres). Eg 0.001 for [mm] -> [m]
+        length      1;
+    }
+
+    // Scaling applied to values written to 'outputName'
+    scaleOutput
+    {
+        //- Length multiplier (from metres). Eg 1000 for [m] -> [mm]
+        length      1;
+
+        //- Force units multiplier (from Pa)
+        force       1;
+
+        //- Moment units multiplier (from N.m)
+        moment      1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lumpedPointMotion/lumpedPointMovementI.H b/src/lumpedPointMotion/lumpedPointMovementI.H
index 43bcdb30947a63cbcb501d6110debecb811ff056..ef73c4ddea837b7522710dee5671faba906103a5 100644
--- a/src/lumpedPointMotion/lumpedPointMovementI.H
+++ b/src/lumpedPointMotion/lumpedPointMovementI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2017-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -133,6 +133,12 @@ inline const Foam::word& Foam::lumpedPointMovement::outputName() const
 }
 
 
+inline const Foam::word& Foam::lumpedPointMovement::logName() const
+{
+    return logName_;
+}
+
+
 inline Foam::lumpedPointState::inputFormatType
 Foam::lumpedPointMovement::inputFormat() const
 {
diff --git a/src/lumpedPointMotion/lumpedPointState.C b/src/lumpedPointMotion/lumpedPointState.C
index a84e963fb66ae959cdd619595ee2b011fe2b19be..8bddc73d8e71cac798340691028cb6c0d6e4df4d 100644
--- a/src/lumpedPointMotion/lumpedPointState.C
+++ b/src/lumpedPointMotion/lumpedPointState.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -94,8 +94,8 @@ void Foam::lumpedPointState::readDict(const dictionary& dict)
 
 Foam::lumpedPointState::lumpedPointState()
 :
-    points_(0),
-    angles_(0),
+    points_(),
+    angles_(),
     degrees_(false),
     rotationPtr_(nullptr)
 {}
@@ -110,10 +110,7 @@ Foam::lumpedPointState::lumpedPointState(const lumpedPointState& rhs)
 {}
 
 
-Foam::lumpedPointState::lumpedPointState
-(
-    const pointField& pts
-)
+Foam::lumpedPointState::lumpedPointState(const pointField& pts)
 :
     points_(pts),
     angles_(points_.size(), Zero),
@@ -122,10 +119,7 @@ Foam::lumpedPointState::lumpedPointState
 {}
 
 
-Foam::lumpedPointState::lumpedPointState
-(
-    tmp<pointField>& pts
-)
+Foam::lumpedPointState::lumpedPointState(tmp<pointField>& pts)
 :
     points_(pts),
     angles_(points_.size(), Zero),
@@ -134,13 +128,10 @@ Foam::lumpedPointState::lumpedPointState
 {}
 
 
-Foam::lumpedPointState::lumpedPointState
-(
-    const dictionary& dict
-)
+Foam::lumpedPointState::lumpedPointState(const dictionary& dict)
 :
-    points_(0),
-    angles_(0),
+    points_(),
+    angles_(),
     degrees_(false),
     rotationPtr_(nullptr)
 {
@@ -168,6 +159,15 @@ void Foam::lumpedPointState::operator=(const lumpedPointState& rhs)
 }
 
 
+void Foam::lumpedPointState::scalePoints(const scalar scaleFactor)
+{
+    if (scaleFactor > 0)
+    {
+        points_ *= scaleFactor;
+    }
+}
+
+
 void Foam::lumpedPointState::relax
 (
     const scalar alpha,
@@ -273,19 +273,17 @@ void Foam::lumpedPointState::writePlain(Ostream& os) const
     {
         const vector& pt = points_[i];
 
-        os  << pt.x() << ' '
-            << pt.y() << ' '
-            << pt.z() << ' ';
+        os  << pt.x() << ' ' << pt.y() << ' ' << pt.z();
 
         if (i < angles_.size())
         {
-            os  << angles_[i].x() << ' '
-                << angles_[i].y() << ' '
-                << angles_[i].z() << '\n';
+            os  << ' ' << angles_[i].x()
+                << ' ' << angles_[i].y()
+                << ' ' << angles_[i].z() << '\n';
         }
         else
         {
-            os  << "0 0 0\n";
+            os  << " 0 0 0\n";
         }
     }
 }
diff --git a/src/lumpedPointMotion/lumpedPointState.H b/src/lumpedPointMotion/lumpedPointState.H
index 0d32fa892770ab638f3105c0f0906ebdcd52f818..2d69c8c1b2772df07d900c150d05954f6ae86210 100644
--- a/src/lumpedPointMotion/lumpedPointState.H
+++ b/src/lumpedPointMotion/lumpedPointState.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -143,13 +143,16 @@ public:
         //- The local-to-global transformation for each point
         inline const tensorField& rotations() const;
 
+        //- Scale points by given factor.
+        //  Zero and negative values are ignored.
+        void scalePoints(const scalar scaleFactor);
+
         //- Relax the state
         //  alpha = 1 : no relaxation
         //  alpha < 1 : relaxation
         //  alpha = 0 : do nothing
         void relax(const scalar alpha, const lumpedPointState& prev);
 
-
         //- Read input as dictionary content
         bool readData(Istream& is);
 
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index d9a73c1254872fa8b4d003c74b21c860dba1e818..45f839d1a0defcd51550b297c566c4d808a09b1f 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -170,6 +170,8 @@ $(faceSources)/boundaryToFace/boundaryToFace.C
 $(faceSources)/zoneToFace/zoneToFace.C
 $(faceSources)/boxToFace/boxToFace.C
 $(faceSources)/regionToFace/regionToFace.C
+$(faceSources)/cylinderToFace/cylinderToFace.C
+$(faceSources)/cylinderAnnulusToFace/cylinderAnnulusToFace.C
 
 pointSources = sets/pointSources
 $(pointSources)/labelToPoint/labelToPoint.C
diff --git a/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.C b/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.C
index de8b54f4024a904ec21ce00a4fc91dd54df5d16d..17303cf7a50cd7dae89d7b03e39eaeba46309178 100644
--- a/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.C
+++ b/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.C
@@ -54,7 +54,7 @@ namespace Foam
 void Foam::cylindrical::init
 (
     const objectRegistry& obr,
-    const List<label>& cells
+    const labelUList& cells
 )
 {
     const polyMesh& mesh = refCast<const polyMesh>(obr);
@@ -196,7 +196,7 @@ void Foam::cylindrical::updateCells
 
     forAll(cells, i)
     {
-        label celli = cells[i];
+        const label celli = cells[i];
         vector dir = cc[celli] - origin_;
         dir /= mag(dir) + VSMALL;
 
diff --git a/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.H b/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.H
index 1c3cd012cf1b946171673edcdeac78c751ab68f5..9f7a2463e724211c93abee9e82a45c0c512f8825 100644
--- a/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.H
+++ b/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.H
@@ -52,6 +52,7 @@ SourceFiles
 
 #include "point.H"
 #include "vector.H"
+#include "ListOps.H"
 #include "coordinateRotation.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -85,7 +86,7 @@ class cylindrical
         void init
         (
             const objectRegistry& obr,
-            const List<label>& cells = List<label>()
+            const labelUList& cells = Foam::emptyLabelList
         );
 
 
diff --git a/src/meshTools/sets/cellSources/boxToCell/boxToCell.C b/src/meshTools/sets/cellSources/boxToCell/boxToCell.C
index bd01674c0a8d34aeb0049c2dc48a21313e7c0bbb..60a7b7436bbfb655c8524c49e8cf131ce4c87e0b 100644
--- a/src/meshTools/sets/cellSources/boxToCell/boxToCell.C
+++ b/src/meshTools/sets/cellSources/boxToCell/boxToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,20 +25,15 @@ License
 
 #include "boxToCell.H"
 #include "polyMesh.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(boxToCell, 0);
-
-addToRunTimeSelectionTable(topoSetSource, boxToCell, word);
-
-addToRunTimeSelectionTable(topoSetSource, boxToCell, istream);
-
+    defineTypeNameAndDebug(boxToCell, 0);
+    addToRunTimeSelectionTable(topoSetSource, boxToCell, word);
+    addToRunTimeSelectionTable(topoSetSource, boxToCell, istream);
 }
 
 
@@ -72,7 +67,6 @@ void Foam::boxToCell::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::boxToCell::boxToCell
 (
     const polyMesh& mesh,
@@ -84,7 +78,6 @@ Foam::boxToCell::boxToCell
 {}
 
 
-// Construct from dictionary
 Foam::boxToCell::boxToCell
 (
     const polyMesh& mesh,
@@ -101,7 +94,6 @@ Foam::boxToCell::boxToCell
 {}
 
 
-// Construct from Istream
 Foam::boxToCell::boxToCell
 (
     const polyMesh& mesh,
@@ -112,6 +104,7 @@ Foam::boxToCell::boxToCell
     bbs_(1, treeBoundBox(checkIs(is)))
 {}
 
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::boxToCell::~boxToCell()
diff --git a/src/meshTools/sets/cellSources/boxToCell/boxToCell.H b/src/meshTools/sets/cellSources/boxToCell/boxToCell.H
index 1f100c77f9efd6df41b2c6cb0d2d6b4f8bda72ea..24a8b482733b71d059367995f193e19a8d17819d 100644
--- a/src/meshTools/sets/cellSources/boxToCell/boxToCell.H
+++ b/src/meshTools/sets/cellSources/boxToCell/boxToCell.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,13 +51,11 @@ class boxToCell
 :
     public topoSetSource
 {
-
     // Private data
 
         //- Add usage string
         static addToUsageTable usage_;
 
-
         //- Bounding box.
         treeBoundBoxList bbs_;
 
@@ -112,7 +110,6 @@ public:
             const topoSetSource::setAction action,
             topoSet&
         ) const;
-
 };
 
 
diff --git a/src/meshTools/sets/cellSources/cellToCell/cellToCell.C b/src/meshTools/sets/cellSources/cellToCell/cellToCell.C
index f7394cf43c3d2da28a3e5258ac83befec0ed256b..4644f755af36014a92c59fff543e966d7c513c89 100644
--- a/src/meshTools/sets/cellSources/cellToCell/cellToCell.C
+++ b/src/meshTools/sets/cellSources/cellToCell/cellToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,20 +26,15 @@ License
 #include "cellToCell.H"
 #include "polyMesh.H"
 #include "cellSet.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(cellToCell, 0);
-
-addToRunTimeSelectionTable(topoSetSource, cellToCell, word);
-
-addToRunTimeSelectionTable(topoSetSource, cellToCell, istream);
-
+    defineTypeNameAndDebug(cellToCell, 0);
+    addToRunTimeSelectionTable(topoSetSource, cellToCell, word);
+    addToRunTimeSelectionTable(topoSetSource, cellToCell, istream);
 }
 
 
@@ -53,7 +48,6 @@ Foam::topoSetSource::addToUsageTable Foam::cellToCell::usage_
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::cellToCell::cellToCell
 (
     const polyMesh& mesh,
@@ -65,7 +59,6 @@ Foam::cellToCell::cellToCell
 {}
 
 
-// Construct from dictionary
 Foam::cellToCell::cellToCell
 (
     const polyMesh& mesh,
@@ -77,7 +70,6 @@ Foam::cellToCell::cellToCell
 {}
 
 
-// Construct from Istream
 Foam::cellToCell::cellToCell
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/cellSources/cylinderAnnulusToCell/cylinderAnnulusToCell.C b/src/meshTools/sets/cellSources/cylinderAnnulusToCell/cylinderAnnulusToCell.C
index 06d89f6e3e8fd8facdc212b5fa99571ba4f2d829..a566caf656cf2a433500707287e8631750323e3a 100644
--- a/src/meshTools/sets/cellSources/cylinderAnnulusToCell/cylinderAnnulusToCell.C
+++ b/src/meshTools/sets/cellSources/cylinderAnnulusToCell/cylinderAnnulusToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,7 +107,6 @@ Foam::cylinderAnnulusToCell::cylinderAnnulusToCell
 {}
 
 
-// Construct from Istream
 Foam::cylinderAnnulusToCell::cylinderAnnulusToCell
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/cellSources/cylinderToCell/cylinderToCell.C b/src/meshTools/sets/cellSources/cylinderToCell/cylinderToCell.C
index 787c73da692a061b50e1846f9d864b2a698b2658..6094b9c02fc699aef9de7d87ef93d773d7d3ca47 100644
--- a/src/meshTools/sets/cellSources/cylinderToCell/cylinderToCell.C
+++ b/src/meshTools/sets/cellSources/cylinderToCell/cylinderToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,7 +102,6 @@ Foam::cylinderToCell::cylinderToCell
 {}
 
 
-// Construct from Istream
 Foam::cylinderToCell::cylinderToCell
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/cellSources/faceToCell/faceToCell.C b/src/meshTools/sets/cellSources/faceToCell/faceToCell.C
index 5a9b3a3f630224affaa56501b5640827d5a2f6eb..5847c624dd29541cedc4808b6841dccfae870a42 100644
--- a/src/meshTools/sets/cellSources/faceToCell/faceToCell.C
+++ b/src/meshTools/sets/cellSources/faceToCell/faceToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,6 @@ License
 #include "faceToCell.H"
 #include "polyMesh.H"
 #include "faceSet.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -148,7 +147,6 @@ void Foam::faceToCell::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::faceToCell::faceToCell
 (
     const polyMesh& mesh,
@@ -162,7 +160,6 @@ Foam::faceToCell::faceToCell
 {}
 
 
-// Construct from dictionary
 Foam::faceToCell::faceToCell
 (
     const polyMesh& mesh,
@@ -175,7 +172,6 @@ Foam::faceToCell::faceToCell
 {}
 
 
-// Construct from Istream
 Foam::faceToCell::faceToCell
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/cellSources/faceZoneToCell/faceZoneToCell.C b/src/meshTools/sets/cellSources/faceZoneToCell/faceZoneToCell.C
index 7640e0519229eaef3e7554e144fc31d02733d8ce..96ddde1489176e6f6b4c0e4ffa6185d472b4f03f 100644
--- a/src/meshTools/sets/cellSources/faceZoneToCell/faceZoneToCell.C
+++ b/src/meshTools/sets/cellSources/faceZoneToCell/faceZoneToCell.C
@@ -25,7 +25,6 @@ License
 
 #include "faceZoneToCell.H"
 #include "polyMesh.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -105,7 +104,6 @@ void Foam::faceZoneToCell::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::faceZoneToCell::faceZoneToCell
 (
     const polyMesh& mesh,
@@ -119,7 +117,6 @@ Foam::faceZoneToCell::faceZoneToCell
 {}
 
 
-// Construct from dictionary
 Foam::faceZoneToCell::faceZoneToCell
 (
     const polyMesh& mesh,
@@ -132,7 +129,6 @@ Foam::faceZoneToCell::faceZoneToCell
 {}
 
 
-// Construct from Istream
 Foam::faceZoneToCell::faceZoneToCell
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/cellSources/fieldToCell/fieldToCell.C b/src/meshTools/sets/cellSources/fieldToCell/fieldToCell.C
index 08ce936c46523f7ef01b02781bdada71dc3af4d0..c3914d5bb9f81013de8bd35ae62be9d6a97c2418 100644
--- a/src/meshTools/sets/cellSources/fieldToCell/fieldToCell.C
+++ b/src/meshTools/sets/cellSources/fieldToCell/fieldToCell.C
@@ -93,7 +93,6 @@ void Foam::fieldToCell::applyToSet
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::fieldToCell::fieldToCell
 (
     const polyMesh& mesh,
@@ -109,7 +108,6 @@ Foam::fieldToCell::fieldToCell
 {}
 
 
-// Construct from dictionary
 Foam::fieldToCell::fieldToCell
 (
     const polyMesh& mesh,
@@ -123,7 +121,6 @@ Foam::fieldToCell::fieldToCell
 {}
 
 
-// Construct from Istream
 Foam::fieldToCell::fieldToCell
 (
     const polyMesh& mesh,
@@ -151,34 +148,6 @@ void Foam::fieldToCell::applyToSet
     topoSet& set
 ) const
 {
-
-//    // Construct temporary fvMesh from polyMesh
-//    fvMesh fMesh
-//    (
-//        mesh(), // IOobject
-//        mesh().points(),
-//        mesh().faces(),
-//        mesh().cells()
-//    );
-//
-//    const polyBoundaryMesh& patches = mesh().boundaryMesh();
-//
-//    List<polyPatch*> newPatches(patches.size());
-//    forAll(patches, patchi)
-//    {
-//        const polyPatch& pp = patches[patchi];
-//
-//        newPatches[patchi] =
-//            patches[patchi].clone
-//            (
-//                fMesh.boundaryMesh(),
-//                patchi,
-//                pp.size(),
-//                pp.start()
-//            ).ptr();
-//    }
-//    fMesh.addFvPatches(newPatches);
-
     // Try to load field
     IOobject fieldObject
     (
diff --git a/src/meshTools/sets/cellSources/labelToCell/labelToCell.C b/src/meshTools/sets/cellSources/labelToCell/labelToCell.C
index f7b350046c8eb9305269c072dd8f810ad3c2859d..2f71548cdeb4189c191005c19a9b21133a869d2d 100644
--- a/src/meshTools/sets/cellSources/labelToCell/labelToCell.C
+++ b/src/meshTools/sets/cellSources/labelToCell/labelToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,20 +25,15 @@ License
 
 #include "labelToCell.H"
 #include "polyMesh.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(labelToCell, 0);
-
-addToRunTimeSelectionTable(topoSetSource, labelToCell, word);
-
-addToRunTimeSelectionTable(topoSetSource, labelToCell, istream);
-
+    defineTypeNameAndDebug(labelToCell, 0);
+    addToRunTimeSelectionTable(topoSetSource, labelToCell, word);
+    addToRunTimeSelectionTable(topoSetSource, labelToCell, istream);
 }
 
 
@@ -63,7 +58,6 @@ void Foam::labelToCell::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::labelToCell::labelToCell
 (
     const polyMesh& mesh,
@@ -75,7 +69,6 @@ Foam::labelToCell::labelToCell
 {}
 
 
-// Construct from dictionary
 Foam::labelToCell::labelToCell
 (
     const polyMesh& mesh,
@@ -87,7 +80,6 @@ Foam::labelToCell::labelToCell
 {}
 
 
-// Construct from Istream
 Foam::labelToCell::labelToCell
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/cellSources/nbrToCell/nbrToCell.C b/src/meshTools/sets/cellSources/nbrToCell/nbrToCell.C
index cb52cb01db3ede9963fd5235c82dfb1a40483a2f..85f5f494b356d295a7f2dd2db09d01e56418cb80 100644
--- a/src/meshTools/sets/cellSources/nbrToCell/nbrToCell.C
+++ b/src/meshTools/sets/cellSources/nbrToCell/nbrToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,20 +25,15 @@ License
 
 #include "nbrToCell.H"
 #include "polyMesh.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(nbrToCell, 0);
-
-addToRunTimeSelectionTable(topoSetSource, nbrToCell, word);
-
-addToRunTimeSelectionTable(topoSetSource, nbrToCell, istream);
-
+    defineTypeNameAndDebug(nbrToCell, 0);
+    addToRunTimeSelectionTable(topoSetSource, nbrToCell, word);
+    addToRunTimeSelectionTable(topoSetSource, nbrToCell, istream);
 }
 
 
@@ -104,7 +99,6 @@ void Foam::nbrToCell::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::nbrToCell::nbrToCell
 (
     const polyMesh& mesh,
@@ -116,7 +110,6 @@ Foam::nbrToCell::nbrToCell
 {}
 
 
-// Construct from dictionary
 Foam::nbrToCell::nbrToCell
 (
     const polyMesh& mesh,
@@ -128,7 +121,6 @@ Foam::nbrToCell::nbrToCell
 {}
 
 
-// Construct from Istream
 Foam::nbrToCell::nbrToCell
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/cellSources/nearestToCell/nearestToCell.C b/src/meshTools/sets/cellSources/nearestToCell/nearestToCell.C
index cb10f3a0492b237974627eacbe8fc80c1567ea9f..d7e2e11f8ff22197c79c5a13993d7609c6603928 100644
--- a/src/meshTools/sets/cellSources/nearestToCell/nearestToCell.C
+++ b/src/meshTools/sets/cellSources/nearestToCell/nearestToCell.C
@@ -32,13 +32,9 @@ License
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(nearestToCell, 0);
-
-addToRunTimeSelectionTable(topoSetSource, nearestToCell, word);
-
-addToRunTimeSelectionTable(topoSetSource, nearestToCell, istream);
-
+    defineTypeNameAndDebug(nearestToCell, 0);
+    addToRunTimeSelectionTable(topoSetSource, nearestToCell, word);
+    addToRunTimeSelectionTable(topoSetSource, nearestToCell, istream);
 }
 
 
@@ -84,7 +80,6 @@ void Foam::nearestToCell::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::nearestToCell::nearestToCell
 (
     const polyMesh& mesh,
@@ -96,7 +91,6 @@ Foam::nearestToCell::nearestToCell
 {}
 
 
-// Construct from dictionary
 Foam::nearestToCell::nearestToCell
 (
     const polyMesh& mesh,
@@ -108,7 +102,6 @@ Foam::nearestToCell::nearestToCell
 {}
 
 
-// Construct from Istream
 Foam::nearestToCell::nearestToCell
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/cellSources/pointToCell/pointToCell.C b/src/meshTools/sets/cellSources/pointToCell/pointToCell.C
index cfafa4d2e7191a614d3049532a9c97d3b21652bb..415969ad8adbeba0ff765eaffe1cfe7628fe484b 100644
--- a/src/meshTools/sets/cellSources/pointToCell/pointToCell.C
+++ b/src/meshTools/sets/cellSources/pointToCell/pointToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,6 @@ License
 #include "pointToCell.H"
 #include "polyMesh.H"
 #include "pointSet.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -105,7 +104,6 @@ void Foam::pointToCell::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::pointToCell::pointToCell
 (
     const polyMesh& mesh,
@@ -119,7 +117,6 @@ Foam::pointToCell::pointToCell
 {}
 
 
-// Construct from dictionary
 Foam::pointToCell::pointToCell
 (
     const polyMesh& mesh,
@@ -132,7 +129,6 @@ Foam::pointToCell::pointToCell
 {}
 
 
-// Construct from Istream
 Foam::pointToCell::pointToCell
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/cellSources/regionToCell/regionToCell.C b/src/meshTools/sets/cellSources/regionToCell/regionToCell.C
index e85fef3cc2ff6d374e35fb05bd204ee16a78ccda..e0843a19ab0e8e3608cbf4e25ad8d7b6219e0f55 100644
--- a/src/meshTools/sets/cellSources/regionToCell/regionToCell.C
+++ b/src/meshTools/sets/cellSources/regionToCell/regionToCell.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 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,13 +34,9 @@ License
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(regionToCell, 0);
-
-addToRunTimeSelectionTable(topoSetSource, regionToCell, word);
-
-addToRunTimeSelectionTable(topoSetSource, regionToCell, istream);
-
+    defineTypeNameAndDebug(regionToCell, 0);
+    addToRunTimeSelectionTable(topoSetSource, regionToCell, word);
+    addToRunTimeSelectionTable(topoSetSource, regionToCell, istream);
 }
 
 
@@ -90,11 +86,7 @@ void Foam::regionToCell::markRegionFaces
         {
             label facei = pp.start()+i;
             label bFacei = facei-mesh_.nInternalFaces();
-            if
-            (
-                selectedCell[faceCells[i]]
-             != selectedCell[nbrSelected[bFacei]]
-            )
+            if (selectedCell[faceCells[i]] != nbrSelected[bFacei])
             {
                 regionFace[facei] = true;
             }
@@ -385,7 +377,6 @@ void Foam::regionToCell::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::regionToCell::regionToCell
 (
     const polyMesh& mesh,
@@ -401,7 +392,6 @@ Foam::regionToCell::regionToCell
 {}
 
 
-// Construct from dictionary
 Foam::regionToCell::regionToCell
 (
     const polyMesh& mesh,
@@ -420,7 +410,6 @@ Foam::regionToCell::regionToCell
 {}
 
 
-// Construct from Istream
 Foam::regionToCell::regionToCell
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.C b/src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.C
index 2e0d9bd223a755ac4fcd48049269dbce29bad138..799d4d5237d4c6155aea6cdb6fca6dc4f6966d23 100644
--- a/src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.C
+++ b/src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,20 +26,15 @@ License
 #include "rotatedBoxToCell.H"
 #include "polyMesh.H"
 #include "cellModel.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(rotatedBoxToCell, 0);
-
-addToRunTimeSelectionTable(topoSetSource, rotatedBoxToCell, word);
-
-addToRunTimeSelectionTable(topoSetSource, rotatedBoxToCell, istream);
-
+    defineTypeNameAndDebug(rotatedBoxToCell, 0);
+    addToRunTimeSelectionTable(topoSetSource, rotatedBoxToCell, word);
+    addToRunTimeSelectionTable(topoSetSource, rotatedBoxToCell, istream);
 }
 
 
@@ -117,7 +112,6 @@ void Foam::rotatedBoxToCell::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::rotatedBoxToCell::rotatedBoxToCell
 (
     const polyMesh& mesh,
@@ -135,7 +129,6 @@ Foam::rotatedBoxToCell::rotatedBoxToCell
 {}
 
 
-// Construct from dictionary
 Foam::rotatedBoxToCell::rotatedBoxToCell
 (
     const polyMesh& mesh,
@@ -150,7 +143,6 @@ Foam::rotatedBoxToCell::rotatedBoxToCell
 {}
 
 
-// Construct from Istream
 Foam::rotatedBoxToCell::rotatedBoxToCell(const polyMesh& mesh, Istream& is)
 :
     topoSetSource(mesh),
diff --git a/src/meshTools/sets/cellSources/shapeToCell/shapeToCell.C b/src/meshTools/sets/cellSources/shapeToCell/shapeToCell.C
index f98f9be0819ee780b9e1dbaef2abc042c2329a9c..961c678708eea231331b0853959570c7962cf3d2 100644
--- a/src/meshTools/sets/cellSources/shapeToCell/shapeToCell.C
+++ b/src/meshTools/sets/cellSources/shapeToCell/shapeToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,20 +28,15 @@ License
 #include "unitConversion.H"
 #include "hexMatcher.H"
 #include "cellFeatures.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(shapeToCell, 0);
-
-addToRunTimeSelectionTable(topoSetSource, shapeToCell, word);
-
-addToRunTimeSelectionTable(topoSetSource, shapeToCell, istream);
-
+    defineTypeNameAndDebug(shapeToCell, 0);
+    addToRunTimeSelectionTable(topoSetSource, shapeToCell, word);
+    addToRunTimeSelectionTable(topoSetSource, shapeToCell, istream);
 }
 
 
diff --git a/src/meshTools/sets/cellSources/sphereToCell/sphereToCell.C b/src/meshTools/sets/cellSources/sphereToCell/sphereToCell.C
index d3c7cfc3ae8bcec3230a7e978ba8f5c2a6e130db..cd356da9057f6a02e0b80d4111dba6b5c24a6ef7 100644
--- a/src/meshTools/sets/cellSources/sphereToCell/sphereToCell.C
+++ b/src/meshTools/sets/cellSources/sphereToCell/sphereToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,7 +91,6 @@ Foam::sphereToCell::sphereToCell
 {}
 
 
-// Construct from Istream
 Foam::sphereToCell::sphereToCell
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.C b/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.C
index b8e633b9f8587d74f940cb2440096619c6e94a59..4ab71b43d060536090b698a0fd1b5836de2abf0a 100644
--- a/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.C
+++ b/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,20 +28,15 @@ License
 #include "globalMeshData.H"
 #include "plane.H"
 #include "cellSet.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(targetVolumeToCell, 0);
-
-addToRunTimeSelectionTable(topoSetSource, targetVolumeToCell, word);
-
-addToRunTimeSelectionTable(topoSetSource, targetVolumeToCell, istream);
-
+    defineTypeNameAndDebug(targetVolumeToCell, 0);
+    addToRunTimeSelectionTable(topoSetSource, targetVolumeToCell, word);
+    addToRunTimeSelectionTable(topoSetSource, targetVolumeToCell, istream);
 }
 
 
@@ -269,7 +264,6 @@ void Foam::targetVolumeToCell::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::targetVolumeToCell::targetVolumeToCell
 (
     const polyMesh& mesh,
@@ -283,7 +277,6 @@ Foam::targetVolumeToCell::targetVolumeToCell
 {}
 
 
-// Construct from dictionary
 Foam::targetVolumeToCell::targetVolumeToCell
 (
     const polyMesh& mesh,
@@ -297,7 +290,6 @@ Foam::targetVolumeToCell::targetVolumeToCell
 {}
 
 
-// Construct from Istream
 Foam::targetVolumeToCell::targetVolumeToCell
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/cellSources/zoneToCell/zoneToCell.C b/src/meshTools/sets/cellSources/zoneToCell/zoneToCell.C
index 56ed36a172660cb0ee20fa5e8e6daccb15a0e1c0..612da1e918833a89355ed58d48de55cf64bca3c4 100644
--- a/src/meshTools/sets/cellSources/zoneToCell/zoneToCell.C
+++ b/src/meshTools/sets/cellSources/zoneToCell/zoneToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,20 +25,15 @@ License
 
 #include "zoneToCell.H"
 #include "polyMesh.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(zoneToCell, 0);
-
-addToRunTimeSelectionTable(topoSetSource, zoneToCell, word);
-
-addToRunTimeSelectionTable(topoSetSource, zoneToCell, istream);
-
+    defineTypeNameAndDebug(zoneToCell, 0);
+    addToRunTimeSelectionTable(topoSetSource, zoneToCell, word);
+    addToRunTimeSelectionTable(topoSetSource, zoneToCell, istream);
 }
 
 
@@ -92,7 +87,6 @@ void Foam::zoneToCell::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::zoneToCell::zoneToCell
 (
     const polyMesh& mesh,
@@ -104,7 +98,6 @@ Foam::zoneToCell::zoneToCell
 {}
 
 
-// Construct from dictionary
 Foam::zoneToCell::zoneToCell
 (
     const polyMesh& mesh,
@@ -116,7 +109,6 @@ Foam::zoneToCell::zoneToCell
 {}
 
 
-// Construct from Istream
 Foam::zoneToCell::zoneToCell
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/cellZoneSources/setToCellZone/setToCellZone.C b/src/meshTools/sets/cellZoneSources/setToCellZone/setToCellZone.C
index afd7f703316ba0f1adfd77e1e7c2e066e66c07ce..db01c781693dd03bbc9b87b384de639e162ab756 100644
--- a/src/meshTools/sets/cellZoneSources/setToCellZone/setToCellZone.C
+++ b/src/meshTools/sets/cellZoneSources/setToCellZone/setToCellZone.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,20 +26,15 @@ License
 #include "setToCellZone.H"
 #include "polyMesh.H"
 #include "cellZoneSet.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(setToCellZone, 0);
-
-addToRunTimeSelectionTable(topoSetSource, setToCellZone, word);
-
-addToRunTimeSelectionTable(topoSetSource, setToCellZone, istream);
-
+    defineTypeNameAndDebug(setToCellZone, 0);
+    addToRunTimeSelectionTable(topoSetSource, setToCellZone, word);
+    addToRunTimeSelectionTable(topoSetSource, setToCellZone, istream);
 }
 
 
@@ -53,7 +48,6 @@ Foam::topoSetSource::addToUsageTable Foam::setToCellZone::usage_
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::setToCellZone::setToCellZone
 (
     const polyMesh& mesh,
@@ -65,7 +59,6 @@ Foam::setToCellZone::setToCellZone
 {}
 
 
-// Construct from dictionary
 Foam::setToCellZone::setToCellZone
 (
     const polyMesh& mesh,
@@ -77,7 +70,6 @@ Foam::setToCellZone::setToCellZone
 {}
 
 
-// Construct from Istream
 Foam::setToCellZone::setToCellZone
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/faceSources/boundaryToFace/boundaryToFace.C b/src/meshTools/sets/faceSources/boundaryToFace/boundaryToFace.C
index 7ae3ab1a6c1ae69fa1ce8f843329ecccc9204d57..a1085fdbf2af09cbae8a9410fe9785d4dd90de1d 100644
--- a/src/meshTools/sets/faceSources/boundaryToFace/boundaryToFace.C
+++ b/src/meshTools/sets/faceSources/boundaryToFace/boundaryToFace.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,20 +25,15 @@ License
 
 #include "boundaryToFace.H"
 #include "polyMesh.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(boundaryToFace, 0);
-
-addToRunTimeSelectionTable(topoSetSource, boundaryToFace, word);
-
-addToRunTimeSelectionTable(topoSetSource, boundaryToFace, istream);
-
+    defineTypeNameAndDebug(boundaryToFace, 0);
+    addToRunTimeSelectionTable(topoSetSource, boundaryToFace, word);
+    addToRunTimeSelectionTable(topoSetSource, boundaryToFace, istream);
 }
 
 
@@ -68,21 +63,18 @@ void Foam::boundaryToFace::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::boundaryToFace::boundaryToFace(const polyMesh& mesh)
 :
     topoSetSource(mesh)
 {}
 
 
-// Construct from dictionary
 Foam::boundaryToFace::boundaryToFace(const polyMesh& mesh, const dictionary&)
 :
     topoSetSource(mesh)
 {}
 
 
-// Construct from Istream
 Foam::boundaryToFace::boundaryToFace
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/faceSources/boxToFace/boxToFace.C b/src/meshTools/sets/faceSources/boxToFace/boxToFace.C
index 5ad4d8e5e1bafc8b68e13842a9377726b3a11e8f..c15b2963ad1dc3a25078b45561de672a1fcfbd01 100644
--- a/src/meshTools/sets/faceSources/boxToFace/boxToFace.C
+++ b/src/meshTools/sets/faceSources/boxToFace/boxToFace.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,20 +25,15 @@ License
 
 #include "boxToFace.H"
 #include "polyMesh.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(boxToFace, 0);
-
-addToRunTimeSelectionTable(topoSetSource, boxToFace, word);
-
-addToRunTimeSelectionTable(topoSetSource, boxToFace, istream);
-
+    defineTypeNameAndDebug(boxToFace, 0);
+    addToRunTimeSelectionTable(topoSetSource, boxToFace, word);
+    addToRunTimeSelectionTable(topoSetSource, boxToFace, istream);
 }
 
 
@@ -72,7 +67,6 @@ void Foam::boxToFace::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::boxToFace::boxToFace
 (
     const polyMesh& mesh,
@@ -84,7 +78,6 @@ Foam::boxToFace::boxToFace
 {}
 
 
-// Construct from dictionary
 Foam::boxToFace::boxToFace
 (
     const polyMesh& mesh,
@@ -101,7 +94,6 @@ Foam::boxToFace::boxToFace
 {}
 
 
-// Construct from Istream
 Foam::boxToFace::boxToFace
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/faceSources/boxToFace/boxToFace.H b/src/meshTools/sets/faceSources/boxToFace/boxToFace.H
index f3b4d4707d4569753e4746fd3f6e9f209dabbb28..2cde01e1f9850132bb6ef772af16e9472d516018 100644
--- a/src/meshTools/sets/faceSources/boxToFace/boxToFace.H
+++ b/src/meshTools/sets/faceSources/boxToFace/boxToFace.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,7 +51,6 @@ class boxToFace
 :
     public topoSetSource
 {
-
     // Private data
 
         //- Add usage string
@@ -111,7 +110,6 @@ public:
             const topoSetSource::setAction action,
             topoSet&
         ) const;
-
 };
 
 
diff --git a/src/meshTools/sets/faceSources/cellToFace/cellToFace.C b/src/meshTools/sets/faceSources/cellToFace/cellToFace.C
index ac8336a6b8a95d81dceccc4882eac958d042c44b..62912ef7f4bfa9fd9a147395177296a81d3a0f4b 100644
--- a/src/meshTools/sets/faceSources/cellToFace/cellToFace.C
+++ b/src/meshTools/sets/faceSources/cellToFace/cellToFace.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/meshTools/sets/faceSources/cylinderAnnulusToFace/cylinderAnnulusToFace.C b/src/meshTools/sets/faceSources/cylinderAnnulusToFace/cylinderAnnulusToFace.C
new file mode 100644
index 0000000000000000000000000000000000000000..0dbbe7c9ff90137c603c2e4fd7fb3c36023507ba
--- /dev/null
+++ b/src/meshTools/sets/faceSources/cylinderAnnulusToFace/cylinderAnnulusToFace.C
@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "cylinderAnnulusToFace.H"
+#include "polyMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(cylinderAnnulusToFace, 0);
+    addToRunTimeSelectionTable(topoSetSource, cylinderAnnulusToFace, word);
+    addToRunTimeSelectionTable(topoSetSource, cylinderAnnulusToFace, istream);
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::cylinderAnnulusToFace::usage_
+(
+    cylinderAnnulusToFace::typeName,
+    "\n    Usage: cylinderAnnulusToFace (p1X p1Y p1Z) (p2X p2Y p2Z)"
+    " outerRadius innerRadius\n\n"
+    "    Select all faces with face centre within bounding cylinder annulus\n\n"
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::cylinderAnnulusToFace::combine(topoSet& set, const bool add) const
+{
+    const vector axis = p2_ - p1_;
+    const scalar orad2 = sqr(outerRadius_);
+    const scalar irad2 = sqr(innerRadius_);
+    const scalar magAxis2 = magSqr(axis);
+
+    const pointField& ctrs = mesh_.faceCentres();
+
+    forAll(ctrs, facei)
+    {
+        vector d = ctrs[facei] - p1_;
+        scalar magD = d & axis;
+
+        if ((magD > 0) && (magD < magAxis2))
+        {
+            scalar d2 = (d & d) - sqr(magD)/magAxis2;
+            if ((d2 < orad2) && (d2 > irad2))
+            {
+                addOrDelete(set, facei, add);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cylinderAnnulusToFace::cylinderAnnulusToFace
+(
+    const polyMesh& mesh,
+    const vector& p1,
+    const vector& p2,
+    const scalar outerRadius,
+    const scalar innerRadius
+)
+:
+    topoSetSource(mesh),
+    p1_(p1),
+    p2_(p2),
+    outerRadius_(outerRadius),
+    innerRadius_(innerRadius)
+{}
+
+
+Foam::cylinderAnnulusToFace::cylinderAnnulusToFace
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    p1_(dict.lookup("p1")),
+    p2_(dict.lookup("p2")),
+    outerRadius_(readScalar(dict.lookup("outerRadius"))),
+    innerRadius_(readScalar(dict.lookup("innerRadius")))
+{}
+
+
+Foam::cylinderAnnulusToFace::cylinderAnnulusToFace
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetSource(mesh),
+    p1_(checkIs(is)),
+    p2_(checkIs(is)),
+    outerRadius_(readScalar(checkIs(is))),
+    innerRadius_(readScalar(checkIs(is)))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::cylinderAnnulusToFace::~cylinderAnnulusToFace()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::cylinderAnnulusToFace::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+    {
+        Info<< "    Adding faces with centre within cylinder annulus,"
+            << " with p1 = "
+            << p1_ << ", p2 = " << p2_ << " and outer radius = " << outerRadius_
+        << " and inner radius = " << innerRadius_
+        << endl;
+
+        combine(set, true);
+    }
+    else if (action == topoSetSource::DELETE)
+    {
+        Info<< "    Removing faces with centre within cylinder, with p1 = "
+            << p1_ << ", p2 = " << p2_ << " and outer radius = " << outerRadius_
+        << " and inner radius " << innerRadius_
+        << endl;
+
+        combine(set, false);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceSources/cylinderAnnulusToFace/cylinderAnnulusToFace.H b/src/meshTools/sets/faceSources/cylinderAnnulusToFace/cylinderAnnulusToFace.H
new file mode 100644
index 0000000000000000000000000000000000000000..3149a60429c4b296354bbfab0d93420e405cca5d
--- /dev/null
+++ b/src/meshTools/sets/faceSources/cylinderAnnulusToFace/cylinderAnnulusToFace.H
@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::cylinderAnnulusToFace
+
+Description
+    A topoSetSource to select faces based on face centres inside a
+    cylinder annulus.
+
+SourceFiles
+    cylinderAnnulusToFace.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cylinderAnnulusToFace_H
+#define cylinderAnnulusToFace_H
+
+#include "topoSetSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class cylinderAnnulusToFace Declaration
+\*---------------------------------------------------------------------------*/
+
+class cylinderAnnulusToFace
+:
+    public topoSetSource
+{
+
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- First point on cylinder axis
+        vector p1_;
+
+        //- Second point on cylinder axis
+        vector p2_;
+
+        //- Outer Radius
+        scalar outerRadius_;
+
+        //- Inner Radius
+        scalar innerRadius_;
+
+
+    // Private Member Functions
+
+        void combine(topoSet& set, const bool add) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("cylinderAnnulusToFace");
+
+
+    // Constructors
+
+        //- Construct from components
+        cylinderAnnulusToFace
+        (
+            const polyMesh& mesh,
+            const vector& p1,
+            const vector& p2,
+            const scalar outerRadius,
+            const scalar innerRadius
+        );
+
+        //- Construct from dictionary
+        cylinderAnnulusToFace
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from Istream
+        cylinderAnnulusToFace
+        (
+            const polyMesh& mesh,
+            Istream&
+        );
+
+
+    // Destructor
+    virtual ~cylinderAnnulusToFace();
+
+    // Member Functions
+
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceSources/cylinderToFace/cylinderToFace.C b/src/meshTools/sets/faceSources/cylinderToFace/cylinderToFace.C
new file mode 100644
index 0000000000000000000000000000000000000000..a360e4f926c229a09af8b21f055f0411437ae4e4
--- /dev/null
+++ b/src/meshTools/sets/faceSources/cylinderToFace/cylinderToFace.C
@@ -0,0 +1,149 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "cylinderToFace.H"
+#include "polyMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(cylinderToFace, 0);
+    addToRunTimeSelectionTable(topoSetSource, cylinderToFace, word);
+    addToRunTimeSelectionTable(topoSetSource, cylinderToFace, istream);
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::cylinderToFace::usage_
+(
+    cylinderToFace::typeName,
+    "\n    Usage: cylinderToFace (p1X p1Y p1Z) (p2X p2Y p2Z) radius\n\n"
+    "    Select all faces with face centre within bounding cylinder\n\n"
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::cylinderToFace::combine(topoSet& set, const bool add) const
+{
+    const vector axis = p2_ - p1_;
+    const scalar rad2 = sqr(radius_);
+    const scalar magAxis2 = magSqr(axis);
+
+    const pointField& ctrs = mesh_.faceCentres();
+
+    forAll(ctrs, facei)
+    {
+        vector d = ctrs[facei] - p1_;
+        scalar magD = d & axis;
+
+        if ((magD > 0) && (magD < magAxis2))
+        {
+            scalar d2 = (d & d) - sqr(magD)/magAxis2;
+            if (d2 < rad2)
+            {
+                addOrDelete(set, facei, add);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cylinderToFace::cylinderToFace
+(
+    const polyMesh& mesh,
+    const vector& p1,
+    const vector& p2,
+    const scalar radius
+)
+:
+    topoSetSource(mesh),
+    p1_(p1),
+    p2_(p2),
+    radius_(radius)
+{}
+
+
+Foam::cylinderToFace::cylinderToFace
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    p1_(dict.lookup("p1")),
+    p2_(dict.lookup("p2")),
+    radius_(readScalar(dict.lookup("radius")))
+{}
+
+
+Foam::cylinderToFace::cylinderToFace
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetSource(mesh),
+    p1_(checkIs(is)),
+    p2_(checkIs(is)),
+    radius_(readScalar(checkIs(is)))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::cylinderToFace::~cylinderToFace()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::cylinderToFace::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+    {
+        Info<< "    Adding faces with centre within cylinder, with p1 = "
+            << p1_ << ", p2 = " << p2_ << " and radius = " << radius_ << endl;
+
+        combine(set, true);
+    }
+    else if (action == topoSetSource::DELETE)
+    {
+        Info<< "    Removing faces with centre within cylinder, with p1 = "
+            << p1_ << ", p2 = " << p2_ << " and radius = " << radius_ << endl;
+
+        combine(set, false);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceSources/cylinderToFace/cylinderToFace.H b/src/meshTools/sets/faceSources/cylinderToFace/cylinderToFace.H
new file mode 100644
index 0000000000000000000000000000000000000000..b51083aebf562594711a61b9ea9d545930073388
--- /dev/null
+++ b/src/meshTools/sets/faceSources/cylinderToFace/cylinderToFace.H
@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::cylinderToFace
+
+Description
+    A topoSetSource to select faces based on face centres inside a cylinder.
+
+SourceFiles
+    cylinderToFace.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cylinderToFace_H
+#define cylinderToFace_H
+
+#include "topoSetSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class cylinderToFace Declaration
+\*---------------------------------------------------------------------------*/
+
+class cylinderToFace
+:
+    public topoSetSource
+{
+
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- First point on cylinder axis
+        vector p1_;
+
+        //- Second point on cylinder axis
+        vector p2_;
+
+        //- Radius
+        scalar radius_;
+
+
+    // Private Member Functions
+
+        void combine(topoSet& set, const bool add) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("cylinderToFace");
+
+
+    // Constructors
+
+        //- Construct from components
+        cylinderToFace
+        (
+            const polyMesh& mesh,
+            const vector& p1,
+            const vector& p2,
+            const scalar radius
+        );
+
+        //- Construct from dictionary
+        cylinderToFace
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from Istream
+        cylinderToFace
+        (
+            const polyMesh& mesh,
+            Istream&
+        );
+
+
+    //- Destructor
+    virtual ~cylinderToFace();
+
+
+    // Member Functions
+
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceSources/faceToFace/faceToFace.C b/src/meshTools/sets/faceSources/faceToFace/faceToFace.C
index 1e05fdefde2402a5cd60ef4ee33577ea3524f3bb..fb78cb52eda1c5bec397f8121bf6eace136346e9 100644
--- a/src/meshTools/sets/faceSources/faceToFace/faceToFace.C
+++ b/src/meshTools/sets/faceSources/faceToFace/faceToFace.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,20 +26,15 @@ License
 #include "faceToFace.H"
 #include "polyMesh.H"
 #include "faceSet.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(faceToFace, 0);
-
-addToRunTimeSelectionTable(topoSetSource, faceToFace, word);
-
-addToRunTimeSelectionTable(topoSetSource, faceToFace, istream);
-
+    defineTypeNameAndDebug(faceToFace, 0);
+    addToRunTimeSelectionTable(topoSetSource, faceToFace, word);
+    addToRunTimeSelectionTable(topoSetSource, faceToFace, istream);
 }
 
 
@@ -53,7 +48,6 @@ Foam::topoSetSource::addToUsageTable Foam::faceToFace::usage_
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::faceToFace::faceToFace
 (
     const polyMesh& mesh,
@@ -65,7 +59,6 @@ Foam::faceToFace::faceToFace
 {}
 
 
-// Construct from dictionary
 Foam::faceToFace::faceToFace
 (
     const polyMesh& mesh,
@@ -77,7 +70,6 @@ Foam::faceToFace::faceToFace
 {}
 
 
-// Construct from Istream
 Foam::faceToFace::faceToFace
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/faceSources/labelToFace/labelToFace.C b/src/meshTools/sets/faceSources/labelToFace/labelToFace.C
index ca61be907170f8d4cf5eb050dc1e21def9c6aaf3..9f7efd81e91a787acdcdef1cced0a96c7b1b9ef9 100644
--- a/src/meshTools/sets/faceSources/labelToFace/labelToFace.C
+++ b/src/meshTools/sets/faceSources/labelToFace/labelToFace.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,20 +25,15 @@ License
 
 #include "labelToFace.H"
 #include "polyMesh.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(labelToFace, 0);
-
-addToRunTimeSelectionTable(topoSetSource, labelToFace, word);
-
-addToRunTimeSelectionTable(topoSetSource, labelToFace, istream);
-
+    defineTypeNameAndDebug(labelToFace, 0);
+    addToRunTimeSelectionTable(topoSetSource, labelToFace, word);
+    addToRunTimeSelectionTable(topoSetSource, labelToFace, istream);
 }
 
 
@@ -63,7 +58,6 @@ void Foam::labelToFace::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::labelToFace::labelToFace
 (
     const polyMesh& mesh,
@@ -75,7 +69,6 @@ Foam::labelToFace::labelToFace
 {}
 
 
-// Construct from dictionary
 Foam::labelToFace::labelToFace
 (
     const polyMesh& mesh,
@@ -87,7 +80,6 @@ Foam::labelToFace::labelToFace
 {}
 
 
-// Construct from Istream
 Foam::labelToFace::labelToFace
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/faceSources/normalToFace/normalToFace.C b/src/meshTools/sets/faceSources/normalToFace/normalToFace.C
index 576aeae449499ff29e94ae22c769e275a7c7cd2c..9cfdc240ee9bc62b73a15ed4ba864cb8028a4057 100644
--- a/src/meshTools/sets/faceSources/normalToFace/normalToFace.C
+++ b/src/meshTools/sets/faceSources/normalToFace/normalToFace.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,20 +26,15 @@ License
 #include "normalToFace.H"
 #include "polyMesh.H"
 #include "faceSet.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(normalToFace, 0);
-
-addToRunTimeSelectionTable(topoSetSource, normalToFace, word);
-
-addToRunTimeSelectionTable(topoSetSource, normalToFace, istream);
-
+    defineTypeNameAndDebug(normalToFace, 0);
+    addToRunTimeSelectionTable(topoSetSource, normalToFace, word);
+    addToRunTimeSelectionTable(topoSetSource, normalToFace, istream);
 }
 
 
@@ -71,7 +66,6 @@ void Foam::normalToFace::setNormal()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::normalToFace::normalToFace
 (
     const polyMesh& mesh,
@@ -87,7 +81,6 @@ Foam::normalToFace::normalToFace
 }
 
 
-// Construct from dictionary
 Foam::normalToFace::normalToFace(const polyMesh& mesh, const dictionary& dict)
 :
     topoSetSource(mesh),
@@ -98,7 +91,6 @@ Foam::normalToFace::normalToFace(const polyMesh& mesh, const dictionary& dict)
 }
 
 
-// Construct from Istream
 Foam::normalToFace::normalToFace(const polyMesh& mesh, Istream& is)
 :
     topoSetSource(mesh),
diff --git a/src/meshTools/sets/faceSources/patchToFace/patchToFace.C b/src/meshTools/sets/faceSources/patchToFace/patchToFace.C
index 26f7d7459c5a640ca2125567a9b5629411e8604b..94d3170143a7cd0fe0b9366f7c0ee56e7fe557ac 100644
--- a/src/meshTools/sets/faceSources/patchToFace/patchToFace.C
+++ b/src/meshTools/sets/faceSources/patchToFace/patchToFace.C
@@ -25,20 +25,15 @@ License
 
 #include "patchToFace.H"
 #include "polyMesh.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(patchToFace, 0);
-
-addToRunTimeSelectionTable(topoSetSource, patchToFace, word);
-
-addToRunTimeSelectionTable(topoSetSource, patchToFace, istream);
-
+    defineTypeNameAndDebug(patchToFace, 0);
+    addToRunTimeSelectionTable(topoSetSource, patchToFace, word);
+    addToRunTimeSelectionTable(topoSetSource, patchToFace, istream);
 }
 
 
@@ -90,7 +85,6 @@ void Foam::patchToFace::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::patchToFace::patchToFace
 (
     const polyMesh& mesh,
@@ -102,7 +96,6 @@ Foam::patchToFace::patchToFace
 {}
 
 
-// Construct from dictionary
 Foam::patchToFace::patchToFace
 (
     const polyMesh& mesh,
@@ -114,7 +107,6 @@ Foam::patchToFace::patchToFace
 {}
 
 
-// Construct from Istream
 Foam::patchToFace::patchToFace
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/faceSources/pointToFace/pointToFace.C b/src/meshTools/sets/faceSources/pointToFace/pointToFace.C
index 31e5322a39e28c21d5df2b69e7c2b41aacd8f2de..2ecc19a91b72a6d3f3d87e5ac7a95e3b52c0496e 100644
--- a/src/meshTools/sets/faceSources/pointToFace/pointToFace.C
+++ b/src/meshTools/sets/faceSources/pointToFace/pointToFace.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,6 @@ License
 #include "pointToFace.H"
 #include "polyMesh.H"
 #include "pointSet.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -146,7 +145,6 @@ void Foam::pointToFace::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::pointToFace::pointToFace
 (
     const polyMesh& mesh,
@@ -160,7 +158,6 @@ Foam::pointToFace::pointToFace
 {}
 
 
-// Construct from dictionary
 Foam::pointToFace::pointToFace
 (
     const polyMesh& mesh,
@@ -173,7 +170,6 @@ Foam::pointToFace::pointToFace
 {}
 
 
-// Construct from Istream
 Foam::pointToFace::pointToFace
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/faceSources/regionToFace/regionToFace.C b/src/meshTools/sets/faceSources/regionToFace/regionToFace.C
index 29548096e737223f5b5cc435bd00acc5471e8153..e9fbf00645eceb2bd5e97b3ab02d59f013c89d7c 100644
--- a/src/meshTools/sets/faceSources/regionToFace/regionToFace.C
+++ b/src/meshTools/sets/faceSources/regionToFace/regionToFace.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,13 +37,9 @@ License
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(regionToFace, 0);
-
-addToRunTimeSelectionTable(topoSetSource, regionToFace, word);
-
-addToRunTimeSelectionTable(topoSetSource, regionToFace, istream);
-
+    defineTypeNameAndDebug(regionToFace, 0);
+    addToRunTimeSelectionTable(topoSetSource, regionToFace, word);
+    addToRunTimeSelectionTable(topoSetSource, regionToFace, istream);
 }
 
 
@@ -176,7 +172,6 @@ void Foam::regionToFace::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::regionToFace::regionToFace
 (
     const polyMesh& mesh,
@@ -190,7 +185,6 @@ Foam::regionToFace::regionToFace
 {}
 
 
-// Construct from dictionary
 Foam::regionToFace::regionToFace
 (
     const polyMesh& mesh,
@@ -203,7 +197,6 @@ Foam::regionToFace::regionToFace
 {}
 
 
-// Construct from Istream
 Foam::regionToFace::regionToFace
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/faceSources/zoneToFace/zoneToFace.C b/src/meshTools/sets/faceSources/zoneToFace/zoneToFace.C
index fa29739e75f011d81d23a42a333196f785292fce..9285cc1e74203211b58b7b82fe2a1878d02993ff 100644
--- a/src/meshTools/sets/faceSources/zoneToFace/zoneToFace.C
+++ b/src/meshTools/sets/faceSources/zoneToFace/zoneToFace.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,20 +25,15 @@ License
 
 #include "zoneToFace.H"
 #include "polyMesh.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(zoneToFace, 0);
-
-addToRunTimeSelectionTable(topoSetSource, zoneToFace, word);
-
-addToRunTimeSelectionTable(topoSetSource, zoneToFace, istream);
-
+    defineTypeNameAndDebug(zoneToFace, 0);
+    addToRunTimeSelectionTable(topoSetSource, zoneToFace, word);
+    addToRunTimeSelectionTable(topoSetSource, zoneToFace, istream);
 }
 
 
@@ -92,7 +87,6 @@ void Foam::zoneToFace::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::zoneToFace::zoneToFace
 (
     const polyMesh& mesh,
@@ -104,7 +98,6 @@ Foam::zoneToFace::zoneToFace
 {}
 
 
-// Construct from dictionary
 Foam::zoneToFace::zoneToFace
 (
     const polyMesh& mesh,
@@ -116,7 +109,6 @@ Foam::zoneToFace::zoneToFace
 {}
 
 
-// Construct from Istream
 Foam::zoneToFace::zoneToFace
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/faceZoneSources/faceZoneToFaceZone/faceZoneToFaceZone.C b/src/meshTools/sets/faceZoneSources/faceZoneToFaceZone/faceZoneToFaceZone.C
index e3f36bb25c151986f8abfd70afa9b8e5c099797f..fe0ee3c5d25e8e36864b08856ec33b7c37b5ff09 100644
--- a/src/meshTools/sets/faceZoneSources/faceZoneToFaceZone/faceZoneToFaceZone.C
+++ b/src/meshTools/sets/faceZoneSources/faceZoneToFaceZone/faceZoneToFaceZone.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,20 +26,15 @@ License
 #include "faceZoneToFaceZone.H"
 #include "polyMesh.H"
 #include "faceZoneSet.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(faceZoneToFaceZone, 0);
-
-addToRunTimeSelectionTable(topoSetSource, faceZoneToFaceZone, word);
-
-addToRunTimeSelectionTable(topoSetSource, faceZoneToFaceZone, istream);
-
+    defineTypeNameAndDebug(faceZoneToFaceZone, 0);
+    addToRunTimeSelectionTable(topoSetSource, faceZoneToFaceZone, word);
+    addToRunTimeSelectionTable(topoSetSource, faceZoneToFaceZone, istream);
 }
 
 
@@ -53,7 +48,6 @@ Foam::topoSetSource::addToUsageTable Foam::faceZoneToFaceZone::usage_
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::faceZoneToFaceZone::faceZoneToFaceZone
 (
     const polyMesh& mesh,
@@ -65,7 +59,6 @@ Foam::faceZoneToFaceZone::faceZoneToFaceZone
 {}
 
 
-// Construct from dictionary
 Foam::faceZoneToFaceZone::faceZoneToFaceZone
 (
     const polyMesh& mesh,
@@ -77,7 +70,6 @@ Foam::faceZoneToFaceZone::faceZoneToFaceZone
 {}
 
 
-// Construct from Istream
 Foam::faceZoneToFaceZone::faceZoneToFaceZone
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/faceZoneSources/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.C b/src/meshTools/sets/faceZoneSources/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.C
index 50b423a4ad72f0372c09fbfd693d31b0a9bd8f72..8739feb0657c6bc52b24503a136cad8cec5bb31e 100644
--- a/src/meshTools/sets/faceZoneSources/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.C
+++ b/src/meshTools/sets/faceZoneSources/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,7 +57,6 @@ Foam::topoSetSource::addToUsageTable Foam::searchableSurfaceToFaceZone::usage_
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from dictionary
 Foam::searchableSurfaceToFaceZone::searchableSurfaceToFaceZone
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/faceZoneSources/setToFaceZone/setToFaceZone.C b/src/meshTools/sets/faceZoneSources/setToFaceZone/setToFaceZone.C
index 821d0777f29a4b19e8342dfb84a9ba2cc084c473..b9b418f2f2598b59cf1d05cabf361f203f893801 100644
--- a/src/meshTools/sets/faceZoneSources/setToFaceZone/setToFaceZone.C
+++ b/src/meshTools/sets/faceZoneSources/setToFaceZone/setToFaceZone.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,20 +26,15 @@ License
 #include "setToFaceZone.H"
 #include "polyMesh.H"
 #include "faceZoneSet.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(setToFaceZone, 0);
-
-addToRunTimeSelectionTable(topoSetSource, setToFaceZone, word);
-
-addToRunTimeSelectionTable(topoSetSource, setToFaceZone, istream);
-
+    defineTypeNameAndDebug(setToFaceZone, 0);
+    addToRunTimeSelectionTable(topoSetSource, setToFaceZone, word);
+    addToRunTimeSelectionTable(topoSetSource, setToFaceZone, istream);
 }
 
 
@@ -54,7 +49,6 @@ Foam::topoSetSource::addToUsageTable Foam::setToFaceZone::usage_
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::setToFaceZone::setToFaceZone
 (
     const polyMesh& mesh,
@@ -66,7 +60,6 @@ Foam::setToFaceZone::setToFaceZone
 {}
 
 
-// Construct from dictionary
 Foam::setToFaceZone::setToFaceZone
 (
     const polyMesh& mesh,
@@ -78,7 +71,6 @@ Foam::setToFaceZone::setToFaceZone
 {}
 
 
-// Construct from Istream
 Foam::setToFaceZone::setToFaceZone
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C b/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C
index 89d666e810cd0d6cecc4341bb1e4f327220f005b..70892ab62f216213794f33f98e505bcaf88895f6 100644
--- a/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C
+++ b/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,7 +51,6 @@ Foam::topoSetSource::addToUsageTable Foam::setsToFaceZone::usage_
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::setsToFaceZone::setsToFaceZone
 (
     const polyMesh& mesh,
@@ -67,7 +66,6 @@ Foam::setsToFaceZone::setsToFaceZone
 {}
 
 
-// Construct from dictionary
 Foam::setsToFaceZone::setsToFaceZone
 (
     const polyMesh& mesh,
@@ -81,7 +79,6 @@ Foam::setsToFaceZone::setsToFaceZone
 {}
 
 
-// Construct from Istream
 Foam::setsToFaceZone::setsToFaceZone
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/pointSources/boxToPoint/boxToPoint.C b/src/meshTools/sets/pointSources/boxToPoint/boxToPoint.C
index c11a49064f97e708284b7ec87257b2bd81b4e20f..7f55b9eb07077335efe2af0301742dc9e04efacf 100644
--- a/src/meshTools/sets/pointSources/boxToPoint/boxToPoint.C
+++ b/src/meshTools/sets/pointSources/boxToPoint/boxToPoint.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,20 +25,15 @@ License
 
 #include "boxToPoint.H"
 #include "polyMesh.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(boxToPoint, 0);
-
-addToRunTimeSelectionTable(topoSetSource, boxToPoint, word);
-
-addToRunTimeSelectionTable(topoSetSource, boxToPoint, istream);
-
+    defineTypeNameAndDebug(boxToPoint, 0);
+    addToRunTimeSelectionTable(topoSetSource, boxToPoint, word);
+    addToRunTimeSelectionTable(topoSetSource, boxToPoint, istream);
 }
 
 
@@ -71,7 +66,6 @@ void Foam::boxToPoint::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::boxToPoint::boxToPoint
 (
     const polyMesh& mesh,
@@ -83,7 +77,6 @@ Foam::boxToPoint::boxToPoint
 {}
 
 
-// Construct from dictionary
 Foam::boxToPoint::boxToPoint
 (
     const polyMesh& mesh,
@@ -100,7 +93,6 @@ Foam::boxToPoint::boxToPoint
 {}
 
 
-// Construct from Istream
 Foam::boxToPoint::boxToPoint
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/pointSources/cellToPoint/cellToPoint.C b/src/meshTools/sets/pointSources/cellToPoint/cellToPoint.C
index 244dc1bba1b53e4041abba8689ee425761bc836e..bd17ab30ebfe059bca68b0a9d3927d5b7bed1d7d 100644
--- a/src/meshTools/sets/pointSources/cellToPoint/cellToPoint.C
+++ b/src/meshTools/sets/pointSources/cellToPoint/cellToPoint.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,6 @@ License
 #include "cellToPoint.H"
 #include "polyMesh.H"
 #include "cellSet.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -84,7 +83,6 @@ void Foam::cellToPoint::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::cellToPoint::cellToPoint
 (
     const polyMesh& mesh,
@@ -98,7 +96,6 @@ Foam::cellToPoint::cellToPoint
 {}
 
 
-// Construct from dictionary
 Foam::cellToPoint::cellToPoint
 (
     const polyMesh& mesh,
@@ -111,7 +108,6 @@ Foam::cellToPoint::cellToPoint
 {}
 
 
-// Construct from Istream
 Foam::cellToPoint::cellToPoint
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/pointSources/faceToPoint/faceToPoint.C b/src/meshTools/sets/pointSources/faceToPoint/faceToPoint.C
index 6e8392f885196cab5adfe4e5a282e09af6ec46ce..1010ae8000c166c6f5efc19ada3caed169e183f2 100644
--- a/src/meshTools/sets/pointSources/faceToPoint/faceToPoint.C
+++ b/src/meshTools/sets/pointSources/faceToPoint/faceToPoint.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,6 @@ License
 #include "faceToPoint.H"
 #include "polyMesh.H"
 #include "faceSet.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -77,7 +76,6 @@ void Foam::faceToPoint::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::faceToPoint::faceToPoint
 (
     const polyMesh& mesh,
@@ -91,7 +89,6 @@ Foam::faceToPoint::faceToPoint
 {}
 
 
-// Construct from dictionary
 Foam::faceToPoint::faceToPoint
 (
     const polyMesh& mesh,
@@ -104,7 +101,6 @@ Foam::faceToPoint::faceToPoint
 {}
 
 
-// Construct from Istream
 Foam::faceToPoint::faceToPoint
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/pointSources/labelToPoint/labelToPoint.C b/src/meshTools/sets/pointSources/labelToPoint/labelToPoint.C
index cd8cd3aec27f070cc2799707bca8d0a51582fbb6..1eed08530889b06ac5f2aee8f7cfd0220ace9198 100644
--- a/src/meshTools/sets/pointSources/labelToPoint/labelToPoint.C
+++ b/src/meshTools/sets/pointSources/labelToPoint/labelToPoint.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,20 +25,15 @@ License
 
 #include "labelToPoint.H"
 #include "polyMesh.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(labelToPoint, 0);
-
-addToRunTimeSelectionTable(topoSetSource, labelToPoint, word);
-
-addToRunTimeSelectionTable(topoSetSource, labelToPoint, istream);
-
+    defineTypeNameAndDebug(labelToPoint, 0);
+    addToRunTimeSelectionTable(topoSetSource, labelToPoint, word);
+    addToRunTimeSelectionTable(topoSetSource, labelToPoint, istream);
 }
 
 
@@ -63,7 +58,6 @@ void Foam::labelToPoint::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::labelToPoint::labelToPoint
 (
     const polyMesh& mesh,
@@ -75,7 +69,6 @@ Foam::labelToPoint::labelToPoint
 {}
 
 
-// Construct from dictionary
 Foam::labelToPoint::labelToPoint
 (
     const polyMesh& mesh,
@@ -87,7 +80,6 @@ Foam::labelToPoint::labelToPoint
 {}
 
 
-// Construct from Istream
 Foam::labelToPoint::labelToPoint
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/pointSources/nearestToPoint/nearestToPoint.C b/src/meshTools/sets/pointSources/nearestToPoint/nearestToPoint.C
index 0f5d134ae9083b74324ae46b297e02773c2832ca..339c9fac7c2d55344babaf6c11c438ffce9416e9 100644
--- a/src/meshTools/sets/pointSources/nearestToPoint/nearestToPoint.C
+++ b/src/meshTools/sets/pointSources/nearestToPoint/nearestToPoint.C
@@ -32,13 +32,9 @@ License
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(nearestToPoint, 0);
-
-addToRunTimeSelectionTable(topoSetSource, nearestToPoint, word);
-
-addToRunTimeSelectionTable(topoSetSource, nearestToPoint, istream);
-
+    defineTypeNameAndDebug(nearestToPoint, 0);
+    addToRunTimeSelectionTable(topoSetSource, nearestToPoint, word);
+    addToRunTimeSelectionTable(topoSetSource, nearestToPoint, istream);
 }
 
 
@@ -104,7 +100,6 @@ void Foam::nearestToPoint::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::nearestToPoint::nearestToPoint
 (
     const polyMesh& mesh,
@@ -116,7 +111,6 @@ Foam::nearestToPoint::nearestToPoint
 {}
 
 
-// Construct from dictionary
 Foam::nearestToPoint::nearestToPoint
 (
     const polyMesh& mesh,
@@ -128,7 +122,6 @@ Foam::nearestToPoint::nearestToPoint
 {}
 
 
-// Construct from Istream
 Foam::nearestToPoint::nearestToPoint
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/pointSources/pointToPoint/pointToPoint.C b/src/meshTools/sets/pointSources/pointToPoint/pointToPoint.C
index ecad9fa04bd885eb79699bb4d505539edec8ddec..3a72bdda6ce97430d5849a268e37756fcadb24c2 100644
--- a/src/meshTools/sets/pointSources/pointToPoint/pointToPoint.C
+++ b/src/meshTools/sets/pointSources/pointToPoint/pointToPoint.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,20 +26,15 @@ License
 #include "pointToPoint.H"
 #include "polyMesh.H"
 #include "pointSet.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(pointToPoint, 0);
-
-addToRunTimeSelectionTable(topoSetSource, pointToPoint, word);
-
-addToRunTimeSelectionTable(topoSetSource, pointToPoint, istream);
-
+    defineTypeNameAndDebug(pointToPoint, 0);
+    addToRunTimeSelectionTable(topoSetSource, pointToPoint, word);
+    addToRunTimeSelectionTable(topoSetSource, pointToPoint, istream);
 }
 
 
@@ -53,7 +48,6 @@ Foam::topoSetSource::addToUsageTable Foam::pointToPoint::usage_
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::pointToPoint::pointToPoint
 (
     const polyMesh& mesh,
@@ -65,7 +59,6 @@ Foam::pointToPoint::pointToPoint
 {}
 
 
-// Construct from dictionary
 Foam::pointToPoint::pointToPoint
 (
     const polyMesh& mesh,
@@ -77,7 +70,6 @@ Foam::pointToPoint::pointToPoint
 {}
 
 
-// Construct from Istream
 Foam::pointToPoint::pointToPoint
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/pointSources/zoneToPoint/zoneToPoint.C b/src/meshTools/sets/pointSources/zoneToPoint/zoneToPoint.C
index 8e9f62097392aed3083d744d01ae8da3548dc6c7..cdb7eb28d9919ea7413c393f36d324af44eefcb0 100644
--- a/src/meshTools/sets/pointSources/zoneToPoint/zoneToPoint.C
+++ b/src/meshTools/sets/pointSources/zoneToPoint/zoneToPoint.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,13 +32,9 @@ License
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(zoneToPoint, 0);
-
-addToRunTimeSelectionTable(topoSetSource, zoneToPoint, word);
-
-addToRunTimeSelectionTable(topoSetSource, zoneToPoint, istream);
-
+    defineTypeNameAndDebug(zoneToPoint, 0);
+    addToRunTimeSelectionTable(topoSetSource, zoneToPoint, word);
+    addToRunTimeSelectionTable(topoSetSource, zoneToPoint, istream);
 }
 
 
@@ -92,7 +88,6 @@ void Foam::zoneToPoint::combine(topoSet& set, const bool add) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::zoneToPoint::zoneToPoint
 (
     const polyMesh& mesh,
@@ -104,7 +99,6 @@ Foam::zoneToPoint::zoneToPoint
 {}
 
 
-// Construct from dictionary
 Foam::zoneToPoint::zoneToPoint
 (
     const polyMesh& mesh,
@@ -116,7 +110,6 @@ Foam::zoneToPoint::zoneToPoint
 {}
 
 
-// Construct from Istream
 Foam::zoneToPoint::zoneToPoint
 (
     const polyMesh& mesh,
diff --git a/src/meshTools/sets/pointZoneSources/setToPointZone/setToPointZone.C b/src/meshTools/sets/pointZoneSources/setToPointZone/setToPointZone.C
index b272752adc46cf7fe63bebe5028c47e6f446e852..b44321f048a349eafccc227c449fa15a7697a675 100644
--- a/src/meshTools/sets/pointZoneSources/setToPointZone/setToPointZone.C
+++ b/src/meshTools/sets/pointZoneSources/setToPointZone/setToPointZone.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,20 +26,15 @@ License
 #include "setToPointZone.H"
 #include "polyMesh.H"
 #include "pointZoneSet.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(setToPointZone, 0);
-
-addToRunTimeSelectionTable(topoSetSource, setToPointZone, word);
-
-addToRunTimeSelectionTable(topoSetSource, setToPointZone, istream);
-
+    defineTypeNameAndDebug(setToPointZone, 0);
+    addToRunTimeSelectionTable(topoSetSource, setToPointZone, word);
+    addToRunTimeSelectionTable(topoSetSource, setToPointZone, istream);
 }
 
 
@@ -53,7 +48,6 @@ Foam::topoSetSource::addToUsageTable Foam::setToPointZone::usage_
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::setToPointZone::setToPointZone
 (
     const polyMesh& mesh,
@@ -65,7 +59,6 @@ Foam::setToPointZone::setToPointZone
 {}
 
 
-// Construct from dictionary
 Foam::setToPointZone::setToPointZone
 (
     const polyMesh& mesh,
@@ -77,7 +70,6 @@ Foam::setToPointZone::setToPointZone
 {}
 
 
-// Construct from Istream
 Foam::setToPointZone::setToPointZone
 (
     const polyMesh& mesh,
diff --git a/src/parallel/decompose/ptscotchDecomp/Make/options b/src/parallel/decompose/ptscotchDecomp/Make/options
index 72c681656641cb76f35515d6f49326244c29f02b..0fcafd255977ebd21454be422e771870e3269072 100644
--- a/src/parallel/decompose/ptscotchDecomp/Make/options
+++ b/src/parallel/decompose/ptscotchDecomp/Make/options
@@ -1,5 +1,5 @@
 sinclude $(GENERAL_RULES)/mplib$(WM_MPLIB)
-sinclude $(RULES)/mplib$(WM_MPLIB)
+sinclude $(DEFAULT_RULES)/mplib$(WM_MPLIB)
 
 EXE_INC = \
     $(PFLAGS) $(PINC) \
diff --git a/src/parallel/decompose/scotchDecomp/Make/options b/src/parallel/decompose/scotchDecomp/Make/options
index ebf3fa21b439dfefba0f3a5ba56039285b333a21..d98ed6ce43f08e67396e8f1170f9cdba9e716ddc 100644
--- a/src/parallel/decompose/scotchDecomp/Make/options
+++ b/src/parallel/decompose/scotchDecomp/Make/options
@@ -3,7 +3,7 @@
  * This is purely to avoid scotch.h including mpicxx.h, which causes problems.
  */
 sinclude $(GENERAL_RULES)/mplib$(WM_MPLIB)
-sinclude $(RULES)/mplib$(WM_MPLIB)
+sinclude $(DEFAULT_RULES)/mplib$(WM_MPLIB)
 
 EXE_INC = \
     $(PFLAGS) $(PINC) \
diff --git a/src/renumber/zoltanRenumber/Make/options b/src/renumber/zoltanRenumber/Make/options
index 79dc4f2bc899c9fe6d9fa890fb3e03075bd7a11f..2d1fd2e105f4e02cc2b11f81d11a74aa131f783c 100644
--- a/src/renumber/zoltanRenumber/Make/options
+++ b/src/renumber/zoltanRenumber/Make/options
@@ -1,5 +1,5 @@
 sinclude $(GENERAL_RULES)/mplib$(WM_MPLIB)
-sinclude $(RULES)/mplib$(WM_MPLIB)
+sinclude $(DEFAULT_RULES)/mplib$(WM_MPLIB)
 
 EXE_INC = \
     $(PFLAGS) $(PINC) \
diff --git a/src/sampling/surface/distanceSurface/distanceSurface.C b/src/sampling/surface/distanceSurface/distanceSurface.C
index e957d6fccc2f52cbccaf616f9786ac6470ae9765..90e934b72809371aa0ed59fac7db42f9d7ac9f71 100644
--- a/src/sampling/surface/distanceSurface/distanceSurface.C
+++ b/src/sampling/surface/distanceSurface/distanceSurface.C
@@ -29,7 +29,6 @@ License
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
-#include "volumeType.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -134,8 +133,6 @@ void Foam::distanceSurface::createGeometry()
 
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh_);
 
-    const labelList& own = fvm.faceOwner();
-
     // Distance to cell centres
     // ~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -173,35 +170,14 @@ void Foam::distanceSurface::createGeometry()
 
         if (signed_)
         {
-            List<volumeType> volType;
-
-            surfPtr_().getVolumeType(cc, volType);
+            vectorField norms;
+            surfPtr_().getNormal(nearest, norms);
 
-            forAll(volType, i)
+            forAll(norms, i)
             {
-                volumeType vT = volType[i];
+                const point diff(cc[i] - nearest[i].hitPoint());
 
-                if (vT == volumeType::OUTSIDE)
-                {
-                    fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
-                }
-                else if (vT == volumeType::INSIDE)
-                {
-                    fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint());
-                }
-                else if (vT == volumeType::UNKNOWN)
-                {
-                    // Treat as very far outside
-                    fld[i] = GREAT;
-                }
-                else
-                {
-                    FatalErrorInFunction
-                        << "getVolumeType failure:"
-                        << " neither INSIDE or OUTSIDE but "
-                        << volumeType::names[vT]
-                        << exit(FatalError);
-                }
+                fld[i] = sign(diff & norms[i]) * Foam::mag(diff);
             }
         }
         else
@@ -223,9 +199,6 @@ void Foam::distanceSurface::createGeometry()
             const pointField& cc = fvm.C().boundaryField()[patchi];
             fvPatchScalarField& fld = cellDistanceBf[patchi];
 
-            const label patchStarti = fvm.boundaryMesh()[patchi].start();
-
-
             List<pointIndexHit> nearest;
             surfPtr_().findNearest
             (
@@ -236,41 +209,14 @@ void Foam::distanceSurface::createGeometry()
 
             if (signed_)
             {
-                List<volumeType> volType;
-
-                surfPtr_().getVolumeType(cc, volType);
+                vectorField norms;
+                surfPtr_().getNormal(nearest, norms);
 
-                forAll(volType, i)
+                forAll(norms, i)
                 {
-                    volumeType vT = volType[i];
-
-                    if (vT == volumeType::OUTSIDE)
-                    {
-                        fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
-                    }
-                    else if (vT == volumeType::INSIDE)
-                    {
-                        fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint());
-                    }
-                    else if (vT == volumeType::UNKNOWN)
-                    {
-                        // Nothing known, so use the cell value.
-                        // - this avoids spurious changes on the boundary
-
-                        // The cell value
-                        const label meshFacei = i+patchStarti;
-                        const scalar& cellVal = cellDistance[own[meshFacei]];
-
-                        fld[i] = cellVal;
-                    }
-                    else
-                    {
-                        FatalErrorInFunction
-                            << "getVolumeType failure:"
-                            << " neither INSIDE or OUTSIDE but "
-                            << volumeType::names[vT]
-                            << exit(FatalError);
-                    }
+                    const point diff(cc[i] - nearest[i].hitPoint());
+
+                    fld[i] = sign(diff & norms[i]) * Foam::mag(diff);
                 }
             }
             else
@@ -303,44 +249,21 @@ void Foam::distanceSurface::createGeometry()
 
         if (signed_)
         {
-            List<volumeType> volType;
+            vectorField norms;
+            surfPtr_().getNormal(nearest, norms);
 
-            surfPtr_().getVolumeType(pts, volType);
-
-            forAll(volType, i)
+            forAll(norms, i)
             {
-                volumeType vT = volType[i];
+                const point diff(pts[i] - nearest[i].hitPoint());
 
-                if (vT == volumeType::OUTSIDE)
-                {
-                    pointDistance_[i] =
-                        Foam::mag(pts[i] - nearest[i].hitPoint());
-                }
-                else if (vT == volumeType::INSIDE)
-                {
-                    pointDistance_[i] =
-                        -Foam::mag(pts[i] - nearest[i].hitPoint());
-                }
-                else if (vT == volumeType::UNKNOWN)
-                {
-                    // Treat as very far outside
-                    pointDistance_[i] = GREAT;
-                }
-                else
-                {
-                    FatalErrorInFunction
-                        << "getVolumeType failure:"
-                        << " neither INSIDE or OUTSIDE but "
-                        << volumeType::names[vT]
-                        << exit(FatalError);
-                }
+                pointDistance_[i] = sign(diff & norms[i]) * Foam::mag(diff);
             }
         }
         else
         {
             forAll(nearest, i)
             {
-                pointDistance_[i] = Foam::mag(pts[i]-nearest[i].hitPoint());
+                pointDistance_[i] = Foam::mag(pts[i] - nearest[i].hitPoint());
             }
         }
     }
diff --git a/src/sampling/surface/isoSurface/isoSurface.C b/src/sampling/surface/isoSurface/isoSurface.C
index adaa14655562038a9ed3a095662ceba00c3f6864..1a2f36b07ecec3b4692139cf33353d1dca2fbeae 100644
--- a/src/sampling/surface/isoSurface/isoSurface.C
+++ b/src/sampling/surface/isoSurface/isoSurface.C
@@ -60,14 +60,6 @@ namespace Foam
         }
     };
 
-
-    // Avoid detecting change if the cells have been marked as GREAT
-    // (ie, ignore them)
-    static inline constexpr bool ignoreValue(const scalar val)
-    {
-        return (val >= 0.5*Foam::GREAT);
-    }
-
 } // End namespace Foam
 
 
@@ -165,7 +157,7 @@ void Foam::isoSurface::syncUnseparatedPoints
 
                 forAll(nbrPts, pointi)
                 {
-                    label nbrPointi = nbrPts[pointi];
+                    const label nbrPointi = nbrPts[pointi];
                     patchInfo[nbrPointi] = pointValues[meshPts[pointi]];
                 }
 
@@ -314,39 +306,19 @@ bool Foam::isoSurface::isEdgeOfFaceCut
     const bool neiLower
 ) const
 {
-    // Could also count number of edges cut and return when they are > 1
-    // but doesn't appear to improve anything
-
     forAll(f, fp)
     {
-        const scalar& pt0Value = pVals[f[fp]];
-
-        if (ignoreValue(pt0Value))
-        {
-            continue;
-        }
-
-        const bool fpLower = (pt0Value < iso_);
+        const bool fpLower = (pVals[f[fp]] < iso_);
 
-        if (fpLower != ownLower || fpLower != neiLower)
+        if
+        (
+            fpLower != ownLower
+         || fpLower != neiLower
+         || fpLower != (pVals[f[f.fcIndex(fp)]] < iso_)
+        )
         {
-            // ++ncut;
             return true;
         }
-        else
-        {
-            const scalar& pt1Value = pVals[f[f.fcIndex(fp)]];
-
-            if (!ignoreValue(pt1Value) && (fpLower != (pt1Value < iso_)))
-            {
-                // ++ncut;
-                return true;
-            }
-        }
-        // if (ncut > 1)
-        // {
-        //     return true;
-        // }
     }
 
     return false;
@@ -401,17 +373,9 @@ void Foam::isoSurface::calcCutTypes
     faceCutType_.setSize(mesh_.nFaces());
     faceCutType_ = NOTCUT;
 
-    // Avoid detecting change if the cells have been marked as GREAT
-    // (ie, ignore them)
-
     for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
     {
         const scalar& ownValue = cVals[own[facei]];
-        if (ignoreValue(ownValue))
-        {
-            continue;
-        }
-
         const bool ownLower = (ownValue < iso_);
 
         scalar nbrValue;
@@ -427,11 +391,6 @@ void Foam::isoSurface::calcCutTypes
             nbrPoint
         );
 
-        if (ignoreValue(nbrValue))
-        {
-            continue;
-        }
-
         const bool neiLower = (nbrValue < iso_);
 
         if (ownLower != neiLower)
@@ -503,7 +462,6 @@ void Foam::isoSurface::calcCutTypes
 
 
     // Propagate internal face cuts into the cells.
-    // For cells marked as ignore (eg, GREAT) - skip this.
 
     for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
     {
@@ -512,20 +470,12 @@ void Foam::isoSurface::calcCutTypes
             continue;
         }
 
-        if
-        (
-            cellCutType_[own[facei]] == NOTCUT
-         && !ignoreValue(cVals[own[facei]])
-        )
+        if (cellCutType_[own[facei]] == NOTCUT)
         {
             cellCutType_[own[facei]] = CUT;
             ++nCutCells_;
         }
-        if
-        (
-            cellCutType_[nei[facei]] == NOTCUT
-         && !ignoreValue(cVals[nei[facei]])
-        )
+        if (cellCutType_[nei[facei]] == NOTCUT)
         {
             cellCutType_[nei[facei]] = CUT;
             ++nCutCells_;
@@ -534,8 +484,6 @@ void Foam::isoSurface::calcCutTypes
 
 
     // Propagate boundary face cuts into the cells.
-    // For cells marked as ignore (eg, GREAT) - skip this and
-    // also suppress the boundary face cut to prevent dangling face cuts.
 
     for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); ++facei)
     {
@@ -544,12 +492,7 @@ void Foam::isoSurface::calcCutTypes
             continue;
         }
 
-        if (ignoreValue(cVals[own[facei]]))
-        {
-            // Suppress dangling boundary face cut
-            faceCutType_[facei] = NOTCUT;
-        }
-        else if (cellCutType_[own[facei]] == NOTCUT)
+        if (cellCutType_[own[facei]] == NOTCUT)
         {
             cellCutType_[own[facei]] = CUT;
             ++nCutCells_;
@@ -774,10 +717,8 @@ void Foam::isoSurface::calcSnappedPoint
 
         bool anyCut = false;
 
-        forAll(pFaces, i)
+        for (const label facei : pFaces)
         {
-            label facei = pFaces[i];
-
             if (faceCutType_[facei] == CUT)
             {
                 anyCut = true;
@@ -795,12 +736,10 @@ void Foam::isoSurface::calcSnappedPoint
         label nOther = 0;
         point otherPointSum = Zero;
 
-        forAll(pFaces, pFacei)
+        for (const label facei : pFaces)
         {
             // Create points for all intersections close to point
             // (i.e. from pyramid edges)
-
-            label facei = pFaces[pFacei];
             const face& f = mesh_.faces()[facei];
             label own = mesh_.faceOwner()[facei];
 
diff --git a/src/sampling/surface/isoSurface/isoSurfaceCell.C b/src/sampling/surface/isoSurface/isoSurfaceCell.C
index e0b1d6bbe1dfea186e11122196ea60acfbafd202..b110d1552ccca45d63ecc4620ee187ed6003aab8 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceCell.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceCell.C
@@ -44,20 +44,6 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    // Avoid detecting change if the cells have been marked as GREAT
-    // (ie, ignore them)
-    static inline constexpr bool ignoreValue(const scalar val)
-    {
-        return (val >= 0.5*Foam::GREAT);
-    }
-
-} // End namespace Foam
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 Foam::scalar Foam::isoSurfaceCell::isoFraction
@@ -99,11 +85,6 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
     const label celli
 ) const
 {
-    if (ignoreValue(cellValues[celli]))
-    {
-        return NOTCUT;
-    }
-
     const cell& cFaces = mesh_.cells()[celli];
 
     if (isTet.test(celli))
@@ -137,11 +118,7 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
         // Check pyramids cut
         for (const label labi : f)
         {
-            if
-            (
-                !ignoreValue(pointValues[labi])
-             && cellLower != (pointValues[labi] < iso_)
-            )
+            if (cellLower != (pointValues[labi] < iso_))
             {
                 edgeCut = true;
                 break;
@@ -187,11 +164,7 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
 
         for (const label pointi : cPoints)
         {
-            if
-            (
-                !ignoreValue(pointValues[pointi])
-             && cellLower != (pointValues[pointi] < iso_)
-            )
+            if (cellLower != (pointValues[pointi] < iso_))
             {
                 ++nCuts;
             }
@@ -201,7 +174,7 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
         {
             return SPHERE;
         }
-        else if (nCuts > 1)
+        else
         {
             return CUT;
         }
diff --git a/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.H b/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.H
index ae3363984a62d9d1e12e42a2ac4cd273028b2c68..6b830147500829f540c478c6fa6a45f165242e2b 100644
--- a/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.H
+++ b/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.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  | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,6 +40,7 @@ SourceFiles
 #include "pointField.H"
 #include "labelledTri.H"
 #include "HashSet.H"
+#include "ListOps.H"
 #include "surfZoneList.H"
 #include "surfaceFormatsCore.H"
 #include "runTimeSelectionTables.H"
@@ -101,7 +102,7 @@ public:
             const pointField& pointLst,
             const UList<Face>& faceLst,
             const UList<surfZone>& zoneLst = List<surfZone>(),
-            const labelUList& faceMap = List<label>()
+            const labelUList& faceMap = Foam::emptyLabelList
         );
 
 
diff --git a/src/surfMesh/meshedSurf/meshedSurfRef.H b/src/surfMesh/meshedSurf/meshedSurfRef.H
index a69cf87c0fa1a6d286b4123fe5fab2bf9f3d1042..a0a15b8df05953eb096b77643fa8f04712be97e7 100644
--- a/src/surfMesh/meshedSurf/meshedSurfRef.H
+++ b/src/surfMesh/meshedSurf/meshedSurfRef.H
@@ -54,10 +54,10 @@ class meshedSurfRef
 
     // Private Member Functions
 
-        //- Disallow construct as copy
+        //- No copy construct
         meshedSurfRef(const meshedSurfRef&) = delete;
 
-        //- Disallow default bitwise assignment
+        //- No copy construct assignment
         void operator=(const meshedSurfRef&) = delete;
 
 public:
@@ -69,7 +69,7 @@ public:
         (
             const pointField& pts,
             const faceList& faces,
-            const labelList& ids = Foam::emptyLabelList
+            const labelUList& ids = Foam::emptyLabelList
         )
         :
             points_(pts),
diff --git a/tutorials/incompressible/lumpedPointMotion/building/steady/system/lumpedPointMovement b/tutorials/incompressible/lumpedPointMotion/building/steady/system/lumpedPointMovement
index f1b33cbbfdfd3ccde4ae5f719bed07ec7e08b20b..5c65da66dabc095f09159cc0039ecb57c20e4bce 100644
--- a/tutorials/incompressible/lumpedPointMotion/building/steady/system/lumpedPointMovement
+++ b/tutorials/incompressible/lumpedPointMotion/building/steady/system/lumpedPointMovement
@@ -66,10 +66,33 @@ communication
     // Output file of forces, written by OpenFOAM
     outputName      forces.out;
 
+    // Log of points/forces/moments during the simulation
+    logName         movement.log;
+
     inputFormat     dictionary;
     outputFormat    dictionary;
 
     debugTable      "<case>/output.txt";
+
+    // Scaling applied to values read from 'inputName'
+    scaleInput
+    {
+        //- Length multiplier (to metres). Eg 0.001 for [mm] -> [m]
+        length      1;
+    }
+
+    // Scaling applied to values written to 'outputName'
+    scaleOutput
+    {
+        //- Length multiplier (from metres). Eg 1000 for [m] -> [mm]
+        length      1;
+
+        //- Force units multiplier (from Pa)
+        force       1;
+
+        //- Moment units multiplier (from N.m)
+        moment      1;
+    }
 }
 
 // ************************************************************************* //
diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/catalyst b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/catalyst
index 6c744b7a9d5b24c48aab6eca8111be1489f04c7d..c70afe7d4e1459440ddb5bbbedcd1aea6ff101d2 100644
--- a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/catalyst
+++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/catalyst
@@ -4,19 +4,23 @@ functions
 {
     catalyst
     {
-        #includeEtc "caseDicts/postProcessing/catalyst/default.cfg"
-
-        mkdir "<case>/insitu";
-
-        // Selected fields (words or regex). Must have cellMask for overset!
-        fields      ( cellMask U p );
+        #includeEtc "caseDicts/insitu/catalyst/catalyst.cfg"
 
         scripts
         (
-            "<system>/scripts/overset.py"
-            "<system>/scripts/writeOverset.py"
+            "<system>/scripts/pressure.py"
+            // "<system>/scripts/vorticity.py"
+            // "<etc>/caseDicts/insitu/catalyst/writeMesh.py"
         );
 
+        inputs
+        {
+            region
+            {
+                // Selected fields (words or regex).
+                fields      ( U p );
+            }
+        }
     }
 }
 
diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/overset.py b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/pressure.py
similarity index 52%
rename from tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/overset.py
rename to tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/pressure.py
index 2d95e8f8bf9cce75b9224dabc54c7dd96bdd5ea7..571bc31d020ef381ac604f17c931d53356dd6a22 100644
--- a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/overset.py
+++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/pressure.py
@@ -10,11 +10,13 @@ from paraview import coprocessing
 imageFileNamePadding=4
 rescale_lookuptable=False
 
+
 # ----------------------- CoProcessor definition -----------------------
 
 def CreateCoProcessor():
   def _CreatePipeline(coprocessor, datadescription):
     class Pipeline:
+      # state file generated using paraview version 5.5.0
 
       # ----------------------------------------------------------------
       # setup views used in the visualization
@@ -27,16 +29,15 @@ def CreateCoProcessor():
 
       # Create a new 'Render View'
       renderView1 = CreateView('RenderView')
-      renderView1.ViewSize = [1077, 763]
+      renderView1.ViewSize = [1091, 766]
       renderView1.AxesGrid = 'GridAxes3DActor'
-      renderView1.CenterOfRotation = [0.00784385809674859, 0.005000000004656613, 0.004999999888241291]
+      renderView1.CenterOfRotation = [0.009999999776482582, 0.004999999888241291, 0.004999999888241291]
       renderView1.StereoType = 0
-      renderView1.CameraPosition = [0.0072242101003740155, 0.0002877833685303474, 0.035060283710920806]
-      renderView1.CameraFocalPoint = [0.00868966107678934, 0.004150999005211765, 0.0049322758242629034]
-      renderView1.CameraViewUp = [0.3542102656908786, 0.9252429122682538, 0.135869941401907]
-      renderView1.CameraParallelScale = 0.00787069031419879
+      renderView1.CameraPosition = [0.009999999776482582, 0.004999999888241291, 0.04819751509880177]
+      renderView1.CameraFocalPoint = [0.009999999776482582, 0.004999999888241291, 0.004999999888241291]
+      renderView1.CameraParallelScale = 0.011180339637598877
       renderView1.CameraParallelProjection = 1
-      renderView1.Background = [0.32, 0.34, 0.43]
+      renderView1.Background = [0, 0, 0]
 
       # init the 'GridAxes3DActor' selected for 'AxesGrid'
       renderView1.AxesGrid.XTitleFontFile = ''
@@ -50,7 +51,7 @@ def CreateCoProcessor():
       # and provide it with information such as the filename to use,
       # how frequently to write the images, etc.
       coprocessor.RegisterView(renderView1,
-          filename='insitu/image_%t.png', freq=1, fittoscreen=0, magnification=1, width=1077, height=763, cinema={})
+          filename='press_%t.png', freq=1, fittoscreen=0, magnification=1, width=1091, height=766, cinema={})
       renderView1.ViewTime = datadescription.GetTime()
 
       # ----------------------------------------------------------------
@@ -62,89 +63,101 @@ def CreateCoProcessor():
       # setup the data processing pipelines
       # ----------------------------------------------------------------
 
-      # a producer from a simulation input
-      input1 = coprocessor.CreateProducer(datadescription, 'mesh')
+      # create a new 'XML MultiBlock Data Reader'
+      # create a producer from a simulation input
+      regionmesh = coprocessor.CreateProducer(datadescription, 'region/mesh')
+
+      # create a new 'Slice'
+      slice1 = Slice(Input=regionmesh)
+      slice1.SliceType = 'Plane'
+      slice1.SliceOffsetValues = [0.0]
 
-      # cellMask [0,1]
-      threshold1 = Threshold(Input=input1)
-      threshold1.Scalars = ['CELLS', 'cellMask']
-      threshold1.ThresholdRange = [0.9, 1.1]
+      # init the 'Plane' selected for 'SliceType'
+      slice1.SliceType.Origin = [0.01, 0.005, 0.005]
+      slice1.SliceType.Normal = [0.0, 0.0, 1.0]
 
       # ----------------------------------------------------------------
       # setup the visualization in view 'renderView1'
       # ----------------------------------------------------------------
 
-      # show data from threshold1
-      threshold1Display = Show(threshold1, renderView1)
-
-      # get color transfer function/color map for 'cellTypes'
-      cellTypesLUT = GetColorTransferFunction('cellTypes')
-      cellTypesLUT.RGBPoints = [0.0, 0.231373, 0.298039, 0.752941, 1.000244140625, 0.865003, 0.865003, 0.865003, 2.00048828125, 0.705882, 0.0156863, 0.14902]
-      cellTypesLUT.ScalarRangeInitialized = 1.0
+      # show data from slice1
+      slice1Display = Show(slice1, renderView1)
 
-      # get opacity transfer function/opacity map for 'cellTypes'
-      cellTypesPWF = GetOpacityTransferFunction('cellTypes')
-      cellTypesPWF.Points = [0.0, 0.0, 0.5, 0.0, 2.00048828125, 1.0, 0.5, 0.0]
-      cellTypesPWF.ScalarRangeInitialized = 1
+      # get color transfer function/color map for 'p'
+      pLUT = GetColorTransferFunction('p')
+      pLUT.RGBPoints = [-0.2227432131767273, 0.231373, 0.298039, 0.752941, 0.0011433586478233337, 0.865003, 0.865003, 0.865003, 0.22502993047237396, 0.705882, 0.0156863, 0.14902]
+      pLUT.ScalarRangeInitialized = 1.0
 
       # trace defaults for the display properties.
-      threshold1Display.Representation = 'Surface With Edges'
-      threshold1Display.ColorArrayName = ['CELLS', 'cellTypes']
-      threshold1Display.LookupTable = cellTypesLUT
-      threshold1Display.OSPRayScaleArray = 'U'
-      threshold1Display.OSPRayScaleFunction = 'PiecewiseFunction'
-      threshold1Display.SelectOrientationVectors = 'None'
-      threshold1Display.ScaleFactor = 0.0019999999552965165
-      threshold1Display.SelectScaleArray = 'None'
-      threshold1Display.GlyphType = 'Arrow'
-      threshold1Display.GlyphTableIndexArray = 'None'
-      threshold1Display.GaussianRadius = 9.999999776482583e-05
-      threshold1Display.SetScaleArray = ['POINTS', 'U']
-      threshold1Display.ScaleTransferFunction = 'PiecewiseFunction'
-      threshold1Display.OpacityArray = ['POINTS', 'U']
-      threshold1Display.OpacityTransferFunction = 'PiecewiseFunction'
-      threshold1Display.DataAxesGrid = 'GridAxesRepresentation'
-      threshold1Display.SelectionCellLabelFontFile = ''
-      threshold1Display.SelectionPointLabelFontFile = ''
-      threshold1Display.PolarAxes = 'PolarAxesRepresentation'
-      threshold1Display.ScalarOpacityFunction = cellTypesPWF
-      threshold1Display.ScalarOpacityUnitDistance = 0.0017065741933059136
+      slice1Display.Representation = 'Surface'
+      slice1Display.ColorArrayName = ['POINTS', 'p']
+      slice1Display.LookupTable = pLUT
+      slice1Display.OSPRayScaleArray = 'U'
+      slice1Display.OSPRayScaleFunction = 'PiecewiseFunction'
+      slice1Display.SelectOrientationVectors = 'None'
+      slice1Display.ScaleFactor = 0.0019999999552965165
+      slice1Display.SelectScaleArray = 'None'
+      slice1Display.GlyphType = 'Arrow'
+      slice1Display.GlyphTableIndexArray = 'None'
+      slice1Display.GaussianRadius = 9.999999776482583e-05
+      slice1Display.SetScaleArray = ['POINTS', 'U']
+      slice1Display.ScaleTransferFunction = 'PiecewiseFunction'
+      slice1Display.OpacityArray = ['POINTS', 'U']
+      slice1Display.OpacityTransferFunction = 'PiecewiseFunction'
+      slice1Display.DataAxesGrid = 'GridAxesRepresentation'
+      slice1Display.SelectionCellLabelFontFile = ''
+      slice1Display.SelectionPointLabelFontFile = ''
+      slice1Display.PolarAxes = 'PolarAxesRepresentation'
 
       # init the 'PiecewiseFunction' selected for 'ScaleTransferFunction'
-      threshold1Display.ScaleTransferFunction.Points = [-0.2505497634410858, 0.0, 0.5, 0.0, 0.3270378112792969, 1.0, 0.5, 0.0]
+      slice1Display.ScaleTransferFunction.Points = [-0.2436095029115677, 0.0, 0.5, 0.0, 0.2753259241580963, 1.0, 0.5, 0.0]
 
       # init the 'PiecewiseFunction' selected for 'OpacityTransferFunction'
-      threshold1Display.OpacityTransferFunction.Points = [-0.2505497634410858, 0.0, 0.5, 0.0, 0.3270378112792969, 1.0, 0.5, 0.0]
+      slice1Display.OpacityTransferFunction.Points = [-0.2436095029115677, 0.0, 0.5, 0.0, 0.2753259241580963, 1.0, 0.5, 0.0]
 
       # init the 'GridAxesRepresentation' selected for 'DataAxesGrid'
-      threshold1Display.DataAxesGrid.XTitleFontFile = ''
-      threshold1Display.DataAxesGrid.YTitleFontFile = ''
-      threshold1Display.DataAxesGrid.ZTitleFontFile = ''
-      threshold1Display.DataAxesGrid.XLabelFontFile = ''
-      threshold1Display.DataAxesGrid.YLabelFontFile = ''
-      threshold1Display.DataAxesGrid.ZLabelFontFile = ''
+      slice1Display.DataAxesGrid.XTitleFontFile = ''
+      slice1Display.DataAxesGrid.YTitleFontFile = ''
+      slice1Display.DataAxesGrid.ZTitleFontFile = ''
+      slice1Display.DataAxesGrid.XLabelFontFile = ''
+      slice1Display.DataAxesGrid.YLabelFontFile = ''
+      slice1Display.DataAxesGrid.ZLabelFontFile = ''
 
       # init the 'PolarAxesRepresentation' selected for 'PolarAxes'
-      threshold1Display.PolarAxes.PolarAxisTitleFontFile = ''
-      threshold1Display.PolarAxes.PolarAxisLabelFontFile = ''
-      threshold1Display.PolarAxes.LastRadialAxisTextFontFile = ''
-      threshold1Display.PolarAxes.SecondaryRadialAxesTextFontFile = ''
+      slice1Display.PolarAxes.PolarAxisTitleFontFile = ''
+      slice1Display.PolarAxes.PolarAxisLabelFontFile = ''
+      slice1Display.PolarAxes.LastRadialAxisTextFontFile = ''
+      slice1Display.PolarAxes.SecondaryRadialAxesTextFontFile = ''
 
       # setup the color legend parameters for each legend in this view
 
-      # get color legend/bar for cellTypesLUT in view renderView1
-      cellTypesLUTColorBar = GetScalarBar(cellTypesLUT, renderView1)
-      cellTypesLUTColorBar.Title = 'cellTypes'
-      cellTypesLUTColorBar.ComponentTitle = ''
-      cellTypesLUTColorBar.TitleFontFile = ''
-      cellTypesLUTColorBar.LabelFontFile = ''
+      # get color legend/bar for pLUT in view renderView1
+      pLUTColorBar = GetScalarBar(pLUT, renderView1)
+      pLUTColorBar.Title = 'p'
+      pLUTColorBar.ComponentTitle = ''
+      pLUTColorBar.TitleFontFile = ''
+      pLUTColorBar.LabelFontFile = ''
 
       # set color bar visibility
-      cellTypesLUTColorBar.Visibility = 1
+      pLUTColorBar.Visibility = 1
 
       # show color legend
-      threshold1Display.SetScalarBarVisibility(renderView1, True)
+      slice1Display.SetScalarBarVisibility(renderView1, True)
 
+      # ----------------------------------------------------------------
+      # setup color maps and opacity mapes used in the visualization
+      # note: the Get..() functions create a new object, if needed
+      # ----------------------------------------------------------------
+
+      # get opacity transfer function/opacity map for 'p'
+      pPWF = GetOpacityTransferFunction('p')
+      pPWF.Points = [-0.2227432131767273, 0.0, 0.5, 0.0, 0.22502993047237396, 1.0, 0.5, 0.0]
+      pPWF.ScalarRangeInitialized = 1
+
+      # ----------------------------------------------------------------
+      # finally, restore active source
+      SetActiveSource(slice1)
+      # ----------------------------------------------------------------
     return Pipeline()
 
   class CoProcessor(coprocessing.CoProcessor):
@@ -152,8 +165,8 @@ def CreateCoProcessor():
       self.Pipeline = _CreatePipeline(self, datadescription)
 
   coprocessor = CoProcessor()
-  # Frequencies at which the coprocessor updates.
-  freqs = {'mesh': [1, 1, 1]}
+  # these are the frequencies at which the coprocessor updates.
+  freqs = {'region/mesh': [1, 1, 1]}
   coprocessor.SetUpdateFrequencies(freqs)
   return coprocessor
 
@@ -167,7 +180,7 @@ coprocessor = CreateCoProcessor()
 
 #--------------------------------------------------------------
 # Enable Live-Visualizaton with ParaView and the update frequency
-coprocessor.EnableLiveVisualization(True, 1)
+coprocessor.EnableLiveVisualization(False, 1)
 
 # ---------------------- Data Selection method ----------------------
 
diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/vorticity.py b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/vorticity.py
new file mode 100644
index 0000000000000000000000000000000000000000..743da93805553713d22c2df569a88b1e019dfe26
--- /dev/null
+++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/vorticity.py
@@ -0,0 +1,233 @@
+
+from paraview.simple import *
+from paraview import coprocessing
+
+
+#--------------------------------------------------------------
+# Code generated from cpstate.py to create the CoProcessor.
+# paraview version 5.5.0
+
+#--------------------------------------------------------------
+# Global screenshot output options
+imageFileNamePadding=4
+rescale_lookuptable=False
+
+
+# ----------------------- CoProcessor definition -----------------------
+
+def CreateCoProcessor():
+  def _CreatePipeline(coprocessor, datadescription):
+    class Pipeline:
+      # state file generated using paraview version 5.5.0
+
+      # ----------------------------------------------------------------
+      # setup views used in the visualization
+      # ----------------------------------------------------------------
+
+      # trace generated using paraview version 5.5.0
+
+      #### disable automatic camera reset on 'Show'
+      paraview.simple._DisableFirstRenderCameraReset()
+
+      # Create a new 'Render View'
+      renderView1 = CreateView('RenderView')
+      renderView1.ViewSize = [1091, 766]
+      renderView1.AxesGrid = 'GridAxes3DActor'
+      renderView1.CenterOfRotation = [0.009999999776482582, 0.004999999888241291, 0.004999999888241291]
+      renderView1.StereoType = 0
+      renderView1.CameraPosition = [0.009999999776482582, 0.004999999888241291, 0.05232050690623429]
+      renderView1.CameraFocalPoint = [0.009999999776482582, 0.004999999888241291, 0.004999999888241291]
+      renderView1.CameraParallelScale = 0.01224744844016408
+      renderView1.CameraParallelProjection = 1
+      renderView1.Background = [0, 0, 0]
+
+      # init the 'GridAxes3DActor' selected for 'AxesGrid'
+      renderView1.AxesGrid.XTitleFontFile = ''
+      renderView1.AxesGrid.YTitleFontFile = ''
+      renderView1.AxesGrid.ZTitleFontFile = ''
+      renderView1.AxesGrid.XLabelFontFile = ''
+      renderView1.AxesGrid.YLabelFontFile = ''
+      renderView1.AxesGrid.ZLabelFontFile = ''
+
+      # register the view with coprocessor
+      # and provide it with information such as the filename to use,
+      # how frequently to write the images, etc.
+      coprocessor.RegisterView(renderView1,
+          filename='vorticity_%t.png', freq=1, fittoscreen=0, magnification=1, width=1091, height=766, cinema={})
+      renderView1.ViewTime = datadescription.GetTime()
+
+      # ----------------------------------------------------------------
+      # restore active view
+      SetActiveView(renderView1)
+      # ----------------------------------------------------------------
+
+      # ----------------------------------------------------------------
+      # setup the data processing pipelines
+      # ----------------------------------------------------------------
+
+      # create a new 'XML MultiBlock Data Reader'
+      # create a producer from a simulation input
+      regionmesh = coprocessor.CreateProducer(datadescription, 'region/mesh')
+
+      # create a new 'Slice'
+      slice1 = Slice(Input=regionmesh)
+      slice1.SliceType = 'Plane'
+      slice1.SliceOffsetValues = [0.0]
+
+      # init the 'Plane' selected for 'SliceType'
+      slice1.SliceType.Origin = [0.01, 0.005, 0.005]
+      slice1.SliceType.Normal = [0.0, 0.0, 1.0]
+
+      # create a new 'Stream Tracer'
+      streamTracer1 = StreamTracer(Input=slice1,
+          SeedType='High Resolution Line Source')
+      streamTracer1.Vectors = ['POINTS', 'U']
+      streamTracer1.MaximumStreamlineLength = 0.019999999552965164
+
+      # init the 'High Resolution Line Source' selected for 'SeedType'
+      streamTracer1.SeedType.Point1 = [0.0, 0.0, 0.005]
+      streamTracer1.SeedType.Point2 = [0.02, 0.01, 0.005]
+
+      # ----------------------------------------------------------------
+      # setup the visualization in view 'renderView1'
+      # ----------------------------------------------------------------
+
+      # show data from streamTracer1
+      streamTracer1Display = Show(streamTracer1, renderView1)
+
+      # get color transfer function/color map for 'Vorticity'
+      vorticityLUT = GetColorTransferFunction('Vorticity')
+      vorticityLUT.RGBPoints = [0.0, 0.229806, 0.298718, 0.753683, 37.5, 0.303869, 0.406535, 0.844959, 75.0, 0.383013, 0.509419, 0.917388, 112.5, 0.466667, 0.604563, 0.968155, 150.0, 0.552953, 0.688929, 0.995376, 187.5, 0.639176, 0.7596, 0.998151, 225.0, 0.722193, 0.813953, 0.976575, 262.5, 0.798692, 0.849786, 0.931689, 300.0, 0.865395, 0.86541, 0.865396, 337.5, 0.924128, 0.827385, 0.774508, 375.0, 0.958853, 0.769768, 0.678008, 412.5, 0.969954, 0.694267, 0.579375, 450.0, 0.958003, 0.602842, 0.481776, 487.50000000000006, 0.923945, 0.497309, 0.38797, 525.0, 0.869187, 0.378313, 0.300267, 562.5, 0.795632, 0.241284, 0.220526, 600.0, 0.705673, 0.0155562, 0.150233]
+      vorticityLUT.ColorSpace = 'Lab'
+      vorticityLUT.ScalarRangeInitialized = 1.0
+
+      # trace defaults for the display properties.
+      streamTracer1Display.Representation = 'Surface'
+      streamTracer1Display.ColorArrayName = ['POINTS', 'Vorticity']
+      streamTracer1Display.LookupTable = vorticityLUT
+      streamTracer1Display.OSPRayScaleArray = 'AngularVelocity'
+      streamTracer1Display.OSPRayScaleFunction = 'PiecewiseFunction'
+      streamTracer1Display.SelectOrientationVectors = 'Normals'
+      streamTracer1Display.ScaleFactor = 0.001999993808567524
+      streamTracer1Display.SelectScaleArray = 'AngularVelocity'
+      streamTracer1Display.GlyphType = 'Arrow'
+      streamTracer1Display.GlyphTableIndexArray = 'AngularVelocity'
+      streamTracer1Display.GaussianRadius = 9.99996904283762e-05
+      streamTracer1Display.SetScaleArray = ['POINTS', 'AngularVelocity']
+      streamTracer1Display.ScaleTransferFunction = 'PiecewiseFunction'
+      streamTracer1Display.OpacityArray = ['POINTS', 'AngularVelocity']
+      streamTracer1Display.OpacityTransferFunction = 'PiecewiseFunction'
+      streamTracer1Display.DataAxesGrid = 'GridAxesRepresentation'
+      streamTracer1Display.SelectionCellLabelFontFile = ''
+      streamTracer1Display.SelectionPointLabelFontFile = ''
+      streamTracer1Display.PolarAxes = 'PolarAxesRepresentation'
+
+      # init the 'PiecewiseFunction' selected for 'ScaleTransferFunction'
+      streamTracer1Display.ScaleTransferFunction.Points = [-1.1626180405813291e-11, 0.0, 0.5, 0.0, 1.7840937690112886e-11, 1.0, 0.5, 0.0]
+
+      # init the 'PiecewiseFunction' selected for 'OpacityTransferFunction'
+      streamTracer1Display.OpacityTransferFunction.Points = [-1.1626180405813291e-11, 0.0, 0.5, 0.0, 1.7840937690112886e-11, 1.0, 0.5, 0.0]
+
+      # init the 'GridAxesRepresentation' selected for 'DataAxesGrid'
+      streamTracer1Display.DataAxesGrid.XTitleFontFile = ''
+      streamTracer1Display.DataAxesGrid.YTitleFontFile = ''
+      streamTracer1Display.DataAxesGrid.ZTitleFontFile = ''
+      streamTracer1Display.DataAxesGrid.XLabelFontFile = ''
+      streamTracer1Display.DataAxesGrid.YLabelFontFile = ''
+      streamTracer1Display.DataAxesGrid.ZLabelFontFile = ''
+
+      # init the 'PolarAxesRepresentation' selected for 'PolarAxes'
+      streamTracer1Display.PolarAxes.PolarAxisTitleFontFile = ''
+      streamTracer1Display.PolarAxes.PolarAxisLabelFontFile = ''
+      streamTracer1Display.PolarAxes.LastRadialAxisTextFontFile = ''
+      streamTracer1Display.PolarAxes.SecondaryRadialAxesTextFontFile = ''
+
+      # setup the color legend parameters for each legend in this view
+
+      # get color legend/bar for vorticityLUT in view renderView1
+      vorticityLUTColorBar = GetScalarBar(vorticityLUT, renderView1)
+      vorticityLUTColorBar.Title = 'Vorticity'
+      vorticityLUTColorBar.ComponentTitle = 'Magnitude'
+      vorticityLUTColorBar.TitleFontFile = ''
+      vorticityLUTColorBar.LabelFontFile = ''
+
+      # set color bar visibility
+      vorticityLUTColorBar.Visibility = 1
+
+      # show color legend
+      streamTracer1Display.SetScalarBarVisibility(renderView1, True)
+
+      # ----------------------------------------------------------------
+      # setup color maps and opacity mapes used in the visualization
+      # note: the Get..() functions create a new object, if needed
+      # ----------------------------------------------------------------
+
+      # get opacity transfer function/opacity map for 'Vorticity'
+      vorticityPWF = GetOpacityTransferFunction('Vorticity')
+      vorticityPWF.Points = [0.0, 0.0, 0.5, 0.0, 600.0, 1.0, 0.5, 0.0]
+      vorticityPWF.ScalarRangeInitialized = 1
+
+      # ----------------------------------------------------------------
+      # finally, restore active source
+      SetActiveSource(streamTracer1)
+      # ----------------------------------------------------------------
+    return Pipeline()
+
+  class CoProcessor(coprocessing.CoProcessor):
+    def CreatePipeline(self, datadescription):
+      self.Pipeline = _CreatePipeline(self, datadescription)
+
+  coprocessor = CoProcessor()
+  # these are the frequencies at which the coprocessor updates.
+  freqs = {'region/mesh': [1, 1, 1]}
+  coprocessor.SetUpdateFrequencies(freqs)
+  return coprocessor
+
+
+#--------------------------------------------------------------
+# Global variable that will hold the pipeline for each timestep
+# Creating the CoProcessor object, doesn't actually create the ParaView pipeline.
+# It will be automatically setup when coprocessor.UpdateProducers() is called the
+# first time.
+coprocessor = CreateCoProcessor()
+
+#--------------------------------------------------------------
+# Enable Live-Visualizaton with ParaView and the update frequency
+coprocessor.EnableLiveVisualization(True, 1)
+
+# ---------------------- Data Selection method ----------------------
+
+def RequestDataDescription(datadescription):
+    "Callback to populate the request for current timestep"
+    global coprocessor
+    if datadescription.GetForceOutput() == True:
+        # We are just going to request all fields and meshes from the simulation
+        # code/adaptor.
+        for i in range(datadescription.GetNumberOfInputDescriptions()):
+            datadescription.GetInputDescription(i).AllFieldsOn()
+            datadescription.GetInputDescription(i).GenerateMeshOn()
+        return
+
+    # setup requests for all inputs based on the requirements of the
+    # pipeline.
+    coprocessor.LoadRequestedData(datadescription)
+
+# ------------------------ Processing method ------------------------
+
+def DoCoProcessing(datadescription):
+    "Callback to do co-processing for current timestep"
+    global coprocessor
+
+    # Update the coprocessor by providing it the newly generated simulation data.
+    # If the pipeline hasn't been setup yet, this will setup the pipeline.
+    coprocessor.UpdateProducers(datadescription)
+
+    # Write output data, if appropriate.
+    coprocessor.WriteData(datadescription);
+
+    # Write image capture (Last arg: rescale lookup table), if appropriate.
+    coprocessor.WriteImages(datadescription, rescale_lookuptable=rescale_lookuptable,
+        image_quality=0, padding_amount=imageFileNamePadding)
+
+    # Live Visualization, if enabled.
+    coprocessor.DoLiveVisualization(datadescription, "localhost", 22222)
diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/writeOverset.py b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/writeMesh.py
similarity index 91%
rename from tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/writeOverset.py
rename to tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/writeMesh.py
index ba083271fcb049c7d1662592e87938f0cbec9ef1..fa5889ddde3adbc39629eeee9b8ce5835873ab39 100644
--- a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/writeOverset.py
+++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/scripts/writeMesh.py
@@ -18,13 +18,8 @@ def CreateCoProcessor():
       # a producer from a simulation input
       input1 = coprocessor.CreateProducer(datadescription, 'mesh')
 
-      # cellMask [0,1]
-      threshold1 = Threshold(Input=input1)
-      threshold1.Scalars = ['CELLS', 'cellMask']
-      threshold1.ThresholdRange = [0.9, 1.1]
-
-      writer1 = servermanager.writers.XMLMultiBlockDataWriter(Input=threshold1)
-      coprocessor.RegisterWriter(writer1, filename='insitu/overset_%t.vtm', freq=1, paddingamount=0)
+      writer1 = servermanager.writers.XMLMultiBlockDataWriter(Input=input1)
+      coprocessor.RegisterWriter(writer1, filename='insitu/mesh_%t.vtm', freq=1, paddingamount=0)
 
     return Pipeline()
 
diff --git a/wmake/rules/General/ADIOS b/wmake/rules/General/ADIOS
index 87bda2c75204b30516a278532120340376ee487e..c2cb3558ae21d83cdda2c88023684b2af08a95c0 100644
--- a/wmake/rules/General/ADIOS
+++ b/wmake/rules/General/ADIOS
@@ -2,7 +2,7 @@
 # ADIOS includes/libraries
 
 sinclude $(GENERAL_RULES)/mplib$(WM_MPLIB)
-sinclude $(RULES)/mplib$(WM_MPLIB)
+sinclude $(DEFAULT_RULES)/mplib$(WM_MPLIB)
 
 # Obtain compile/link flags via adios_config
 ADIOS_INC  := $(shell $(ADIOS_ARCH_PATH)/bin/adios_config -c)
diff --git a/wmake/rules/General/ADIOS2 b/wmake/rules/General/ADIOS2
index f158236942ce7a3883a793b4b94c33143bf26726..c880d0c3b97ee62265070b54d73ba142f059eb5b 100644
--- a/wmake/rules/General/ADIOS2
+++ b/wmake/rules/General/ADIOS2
@@ -2,7 +2,7 @@
 # ADIOS2 includes/libraries
 
 sinclude $(GENERAL_RULES)/mplib$(WM_MPLIB)
-sinclude $(RULES)/mplib$(WM_MPLIB)
+sinclude $(DEFAULT_RULES)/mplib$(WM_MPLIB)
 
 # Obtain prefix and library information via adios2-config
 ADIOS_PREFIX := $(shell $(ADIOS2_ARCH_PATH)/bin/adios2-config --prefix)