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 e9ace6b706..f3defd4b90 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 0000000000..b9f4e61753
--- /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 0000000000..658ae6e554
--- /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 0000000000..a8517e26bb
--- /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 0000000000..404267de18
--- /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 2f6850f995..0000000000
--- 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 d860e901a9..0000000000
--- 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 98d7ad3cd8..d2c0ba0b14 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 c3710b59fe..19239e6980 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 9d4a22b6b6..92cef08d00 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 c1526e8c54..fb91cc1699 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 553378c472..33a8b32755 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 024ba63160..09aa585579 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 98a281c5d3..704374a5fa 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 99e68179fa..6e6e105844 160000
--- a/modules/avalanche
+++ b/modules/avalanche
@@ -1 +1 @@
-Subproject commit 99e68179fa946476f0be4f957f8d6a2b78887925
+Subproject commit 6e6e105844897d4bf780bbc8d14031bc827e4b04
diff --git a/modules/catalyst b/modules/catalyst
index 3c0a2e7959..5828d45108 160000
--- a/modules/catalyst
+++ b/modules/catalyst
@@ -1 +1 @@
-Subproject commit 3c0a2e7959755a84f48e25bbe3436ec6437f7cf6
+Subproject commit 5828d4510816948b034aa9afdf0b7b05659a9f59
diff --git a/modules/cfmesh b/modules/cfmesh
index 8f6e65ae7c..288f05e08f 160000
--- a/modules/cfmesh
+++ b/modules/cfmesh
@@ -1 +1 @@
-Subproject commit 8f6e65ae7c71a95948b53679e57a73899b1dc999
+Subproject commit 288f05e08f07e693d4222e7b84ea12430947e5bf
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
index 4113b1d752..ff44c9ab37 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 a98f4a6083..63cdbdd8a9 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 cd2fac95cb..8bf7480557 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 50a06cc512..770e018b10 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 43a48312c4..8fcb7016f8 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 4bfda81ac6..b070ab67da 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/finiteVolume/fields/fvPatchFields/derived/swirlFlowRateInletVelocity/swirlFlowRateInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/swirlFlowRateInletVelocity/swirlFlowRateInletVelocityFvPatchVectorField.C
index ebe4ebdd58..89c1ce0008 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_)
-    );
-
-    tmp<vectorField> n = patch().nf();
-
-    const surfaceScalarField& phi =
-        db().lookupObject<surfaceScalarField>(phiName_);
-
-    if (phi.dimensions() == dimVelocity*dimArea)
+    if (totArea > ROOTVSMALL && axis_ != vector(Zero))
     {
-        // volumetric flow-rate
-        operator==(tangentialVelocity + n*avgU);
-    }
-    else if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
-    {
-        const fvPatchField<scalar>& rhop =
-            patch().lookupPatchField<volScalarField, scalar>(rhoName_);
+        const scalar t = this->db().time().timeOutputValue();
+        const scalar flowRate = flowRate_->value(t);
+        const scalar rpm = rpm_->value(t);
 
-        // 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);
+        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_)
+        );
+
+        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/functionObjects/utilities/vtkWrite/vtkWrite.C b/src/functionObjects/utilities/vtkWrite/vtkWrite.C
index f5d28883ca..c0d66f2978 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 cf69d91492..fca1d3c605 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/lumpedPointMotion/lumpedPointDisplacementPointPatchVectorField.C b/src/lumpedPointMotion/lumpedPointDisplacementPointPatchVectorField.C
index 3368d6f59f..1ac55bd8b8 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 4e74f34a2f..00c34c7209 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;
-    }
-
-    const fileName output(coupler().resolveFile(outputName_));
-    OFstream os(output); // ASCII
-
-    if (outputFormat_ == outputFormatType::PLAIN)
-    {
-        os  <<"# output from OpenFOAM" << nl
-            <<"# N, points, forces" << nl
-            << this->size() << nl;
-
-        const char* zeroVector = "0 0 0";
-
-        forAll(locations_, i)
+        os  <<"########" << nl;
+        if (timeinfo)
         {
-            const vector pos = locations_[i] * axis_;
+            os <<"# ";
+            putTime(os, *timeinfo) << nl;
+        }
+        os  <<"# size=" << this->size() << nl
+            <<"# columns (points) (forces)";
 
-            os  << pos.x() << ' '
-                << pos.y() << ' '
-                << pos.z() << ' ';
+        if (writeMoments)
+        {
+            os << " (moments)";
+        }
 
-            if (i < forces.size())
+        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)
             {
-                os  << forces[i].x() << ' '
-                    << forces[i].y() << ' '
-                    << forces[i].z();
+                report = true;
             }
             else
             {
-                os << zeroVector;
+                scaleMoment = 1.0;
+            }
+        }
+
+        if (report)
+        {
+            os  <<"# scaling points=" << scaleLength
+                <<" forces=" << scaleForce;
+
+            if (writeMoments)
+            {
+                os  <<" moments=" << scaleMoment;
             }
 
-            os  << nl;
+            os << nl;
+        }
+
+        os <<"########" << nl;
+
+        forAll(locations_, i)
+        {
+            const vector pos = scaleLength * (locations_[i] * axis_);
+
+            putPlain(os, pos) << ' ';
+
+            if (i < forces.size())
+            {
+                const vector val(scaleForce * forces[i]);
+                putPlain(os, val);
+            }
+            else
+            {
+                putPlain(os, vector::zero);
+            }
+
+            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 fileName output(coupler().resolveFile(outputName_));
+        OFstream os(output, IOstream::ASCII);
 
-        const char* zeroVector = "0 0 0";
-
-        forAll(locations_, i)
-        {
-            const vector pos = locations_[i] * axis_;
-
-            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 dc0f1fb189..7587ec8378 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 0000000000..5223b7d253
--- /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 43bcdb3094..ef73c4ddea 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 a84e963fb6..8bddc73d8e 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 0d32fa8927..2d69c8c1b2 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/coordinateSystems/coordinateRotation/cylindrical.C b/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.C
index de8b54f402..17303cf7a5 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 1c3cd012cf..9f7a2463e7 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/regionToCell/regionToCell.C b/src/meshTools/sets/cellSources/regionToCell/regionToCell.C
index 804808d97e..e0843a19ab 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.
@@ -86,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;
             }
diff --git a/src/parallel/decompose/ptscotchDecomp/Make/options b/src/parallel/decompose/ptscotchDecomp/Make/options
index 72c6816566..0fcafd2559 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 ebf3fa21b4..d98ed6ce43 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 79dc4f2bc8..2d1fd2e105 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 e957d6fccc..90e934b728 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;
+            vectorField norms;
+            surfPtr_().getNormal(nearest, norms);
 
-            surfPtr_().getVolumeType(cc, volType);
-
-            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;
+                vectorField norms;
+                surfPtr_().getNormal(nearest, norms);
 
-                surfPtr_().getVolumeType(cc, volType);
-
-                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)
-                    {
-                        // 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);
-                    }
+                    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 adaa146555..1a2f36b07e 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]];
+        const bool fpLower = (pVals[f[fp]] < iso_);
 
-        if (ignoreValue(pt0Value))
+        if
+        (
+            fpLower != ownLower
+         || fpLower != neiLower
+         || fpLower != (pVals[f[f.fcIndex(fp)]] < iso_)
+        )
         {
-            continue;
-        }
-
-        const bool fpLower = (pt0Value < iso_);
-
-        if (fpLower != ownLower || fpLower != neiLower)
-        {
-            // ++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 e0b1d6bbe1..b110d1552c 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 ae3363984a..6b83014750 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 a69cf87c0f..a0a15b8df0 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 f1b33cbbfd..5c65da66da 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 6c744b7a9d..c70afe7d4e 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 2d95e8f8bf..571bc31d02 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')
 
-      # cellMask [0,1]
-      threshold1 = Threshold(Input=input1)
-      threshold1.Scalars = ['CELLS', 'cellMask']
-      threshold1.ThresholdRange = [0.9, 1.1]
+      # 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]
 
       # ----------------------------------------------------------------
       # setup the visualization in view 'renderView1'
       # ----------------------------------------------------------------
 
-      # show data from threshold1
-      threshold1Display = Show(threshold1, renderView1)
+      # show data from slice1
+      slice1Display = Show(slice1, 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
-
-      # 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 0000000000..743da93805
--- /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 ba083271fc..fa5889ddde 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 87bda2c752..c2cb3558ae 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 f158236942..c880d0c3b9 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)