diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 0bcfcf48a575f19a1d2aa01748f8e5ba00385db1..d7ea0f00996f6b171414068c5ff2c8463bbe7dd3 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -251,6 +251,8 @@ surfaceInterpolation = interpolation/surfaceInterpolation
 $(surfaceInterpolation)/surfaceInterpolation/surfaceInterpolation.C
 $(surfaceInterpolation)/surfaceInterpolationScheme/surfaceInterpolationSchemes.C
 
+$(surfaceInterpolation)/blendedSchemeBase/blendedSchemeBaseName.C
+
 schemes = $(surfaceInterpolation)/schemes
 $(schemes)/linear/linear.C
 $(schemes)/pointLinear/pointLinear.C
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.C
index c59d8d89ee43617be7a238f83b7a176e39cba1a8..e05c34aa208c3e90f42a13a70961851554010c46 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.C
@@ -199,10 +199,7 @@ void Foam::flowRateInletVelocityFvPatchVectorField::updateCoeffs()
     else
     {
         // mass flow-rate
-        if
-        (
-            patch().boundaryMesh().mesh().foundObject<volScalarField>(rhoName_)
-        )
+        if (db().foundObject<volScalarField>(rhoName_))
         {
             const fvPatchField<scalar>& rhop =
                 patch().lookupPatchField<volScalarField, scalar>(rhoName_);
diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.C b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.C
index c814f220601e2618adc0cc1976e8bed0e15e6529..1c68b9cf12eb497f20b35a12bd84e913a6838f5d 100644
--- a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.C
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,6 +39,14 @@ namespace fv
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+template<class Type>
+const surfaceInterpolationScheme<Type>&
+gaussConvectionScheme<Type>::interpScheme() const
+{
+    return tinterpScheme_();
+}
+
+
 template<class Type>
 tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
 gaussConvectionScheme<Type>::interpolate
diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H
index 61d6c22d37aadd055fa8546a2ddbab561760b6ad..a3ce0d1ec0d02ac8e23db84cd220c707023ae09d 100644
--- a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -134,6 +134,8 @@ public:
 
     // Member Functions
 
+        const surfaceInterpolationScheme<Type>& interpScheme() const;
+
         tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > interpolate
         (
             const surfaceScalarField&,
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcCellReduce.C b/src/finiteVolume/finiteVolume/fvc/fvcCellReduce.C
new file mode 100644
index 0000000000000000000000000000000000000000..63c718375bf950119b8592ddd0966eb3cc796d64
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/fvc/fvcCellReduce.C
@@ -0,0 +1,118 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvcCellReduce.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "zeroGradientFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fvc
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class CombineOp>
+tmp<GeometricField<Type, fvPatchField, volMesh> > cellReduce
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf,
+    const CombineOp& cop
+)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> volFieldType;
+
+    const fvMesh& mesh = ssf.mesh();
+
+    tmp<volFieldType> tresult
+    (
+        new volFieldType
+        (
+            IOobject
+            (
+                "cellReduce(" + ssf.name() + ')',
+                ssf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensioned<Type>("0", ssf.dimensions(), pTraits<Type>::zero),
+            zeroGradientFvPatchField<Type>::typeName
+        )
+    );
+
+    volFieldType& result = tresult();
+
+    const labelUList& own = mesh.owner();
+    const labelUList& nbr = mesh.neighbour();
+
+    forAll(own, i)
+    {
+        label cellI = own[i];
+        cop(result[cellI], ssf[i]);
+    }
+    forAll(nbr, i)
+    {
+        label cellI = nbr[i];
+        cop(result[cellI], ssf[i]);
+    }
+
+    result.correctBoundaryConditions();
+
+    return tresult;
+}
+
+
+template<class Type, class CombineOp>
+tmp<GeometricField<Type, fvPatchField, volMesh> > cellReduce
+(
+    const tmp<GeometricField<Type, fvsPatchField, surfaceMesh>&> tssf,
+    const CombineOp& cop
+)
+{
+    tmp<GeometricField<Type, fvPatchField, volMesh> >
+        tvf(cellReduce(cop, tssf));
+
+   tssf.clear();
+    return tvf;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fvc
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcCellReduce.H b/src/finiteVolume/finiteVolume/fvc/fvcCellReduce.H
new file mode 100644
index 0000000000000000000000000000000000000000..19bd24135af0a4dd75a11ca446192301e0e6443a
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/fvc/fvcCellReduce.H
@@ -0,0 +1,82 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+InNamespace
+    Foam::fvc
+
+Description
+    Construct a volume field from a surface field using a combine operator.
+
+SourceFiles
+    fvcCellReduce.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fvcCellReduce_H
+#define fvcCellReduce_H
+
+#include "volFieldsFwd.H"
+#include "surfaceFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Namespace fvc functions Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace fvc
+{
+    template<class Type, class CombineOp>
+    tmp<GeometricField<Type, fvPatchField, volMesh> > cellReduce
+    (
+        const GeometricField<Type, fvsPatchField, surfaceMesh>&,
+        const CombineOp& cop
+    );
+
+    template<class Type, class CombineOp>
+    tmp<GeometricField<Type, fvPatchField, volMesh> > cellReduce
+    (
+        const tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >&,
+        const CombineOp& cop
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "fvcCellReduce.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/blendedSchemeBase/blendedSchemeBase.H b/src/finiteVolume/interpolation/surfaceInterpolation/blendedSchemeBase/blendedSchemeBase.H
new file mode 100644
index 0000000000000000000000000000000000000000..2c485e01477edcf0740317437b98d62d0fc4d75f
--- /dev/null
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/blendedSchemeBase/blendedSchemeBase.H
@@ -0,0 +1,87 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::blendedSchemeBase
+
+Description
+    Base class for blended schemes to provide access to the blending factor
+    surface field
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef blendedSchemeBase_H
+#define blendedSchemeBase_H
+
+#include "className.H"
+#include "tmp.H"
+#include "surfaceFieldsFwd.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+TemplateName(blendedSchemeBase);
+
+/*---------------------------------------------------------------------------*\
+                      Class blendedSchemeBase Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class blendedSchemeBase
+:
+    public blendedSchemeBaseName
+{
+
+public:
+
+    //- Constructor
+    blendedSchemeBase()
+    {}
+
+    //- Destructor
+    virtual ~blendedSchemeBase()
+    {}
+
+
+    // Memeber Functions
+
+        //- Return the face-based blending factor
+        virtual tmp<surfaceScalarField> blendingFactor
+        (
+             const GeometricField<Type, fvPatchField, volMesh>& vf
+        ) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/blendedSchemeBase/blendedSchemeBaseName.C b/src/finiteVolume/interpolation/surfaceInterpolation/blendedSchemeBase/blendedSchemeBaseName.C
new file mode 100644
index 0000000000000000000000000000000000000000..0a97d7f617b49805ea7502176dc0acca2295be4c
--- /dev/null
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/blendedSchemeBase/blendedSchemeBaseName.C
@@ -0,0 +1,36 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "blendedSchemeBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(blendedSchemeBaseName, 0);
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/blended/blended.H b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/blended/blended.H
index 6d207fcdcde08bfedfa59287d128464109a5b2b0..f58b8c2393b4119a847702eb342d28e72e537b57 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/blended/blended.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/blended/blended.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,7 @@ SourceFiles
 #define blended_H
 
 #include "limitedSurfaceInterpolationScheme.H"
+#include "blendedSchemeBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,7 +50,8 @@ namespace Foam
 template<class Type>
 class blended
 :
-    public limitedSurfaceInterpolationScheme<Type>
+    public limitedSurfaceInterpolationScheme<Type>,
+    public blendedSchemeBase<Type>
 {
     // Private data
 
@@ -111,8 +113,40 @@ public:
         {}
 
 
+    //- Destructor
+    virtual ~blended()
+    {}
+
+
     // Member Functions
 
+        //- Return the face-based blending factor
+        virtual tmp<surfaceScalarField> blendingFactor
+        (
+             const GeometricField<Type, fvPatchField, volMesh>& vf
+        ) const
+        {
+            return tmp<surfaceScalarField>
+            (
+                new surfaceScalarField
+                (
+                    IOobject
+                    (
+                        vf.name() + "BlendingFactor",
+                        this->mesh().time().timeName(),
+                        this->mesh()
+                    ),
+                    this->mesh(),
+                    dimensionedScalar
+                    (
+                        "blendingFactor",
+                        dimless,
+                        blendingFactor_
+                    )
+                )
+            );
+        }
+
         //- Return the interpolation limiter
         virtual tmp<surfaceScalarField> limiter
         (
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CoBlended/CoBlended.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CoBlended/CoBlended.H
index e8ed67a9d836ebc1f7e01ed3f8de9e8608b2947b..3f15520061c879c341d7093a997462daa162c12a 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CoBlended/CoBlended.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CoBlended/CoBlended.H
@@ -63,6 +63,7 @@ SourceFiles
 #define CoBlended_H
 
 #include "surfaceInterpolationScheme.H"
+#include "blendedSchemeBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -76,7 +77,8 @@ namespace Foam
 template<class Type>
 class CoBlended
 :
-    public surfaceInterpolationScheme<Type>
+    public surfaceInterpolationScheme<Type>,
+    public blendedSchemeBase<Type>
 {
     // Private data
 
@@ -135,10 +137,7 @@ public:
             ),
             faceFlux_
             (
-                mesh.lookupObject<surfaceScalarField>
-                (
-                    word(is)
-                )
+                mesh.lookupObject<surfaceScalarField>(word(is))
             )
         {
             if (Co1_ < 0 || Co2_ < 0 || Co1_ >= Co2_)
@@ -170,7 +169,7 @@ public:
             (
                 surfaceInterpolationScheme<Type>::New(mesh, faceFlux, is)
             ),
-             faceFlux_(faceFlux)
+            faceFlux_(faceFlux)
         {
             if (Co1_ < 0 || Co2_ < 0 || Co1_ >= Co2_)
             {
@@ -182,17 +181,39 @@ public:
         }
 
 
+    //- Destructor
+    virtual ~CoBlended()
+    {}
+
+
     // Member Functions
 
-        //- Return the face-based Courant number blending factor
-        tmp<surfaceScalarField> blendingFactor() const
+        //- Return the face-based blending factor
+        virtual tmp<surfaceScalarField> blendingFactor
+        (
+             const GeometricField<Type, fvPatchField, volMesh>& vf
+        ) const
         {
-            const fvMesh& mesh = faceFlux_.mesh();
+            const fvMesh& mesh = this->mesh();
 
-            return
+            tmp<surfaceScalarField> tbf
             (
-                scalar(1) -
-                max
+                new surfaceScalarField
+                (
+                    IOobject
+                    (
+                        vf.name() + "BlendingFactor",
+                        mesh.time().timeName(),
+                        mesh
+                    ),
+                    mesh,
+                    dimensionedScalar("blendingFactor", dimless, 0.0)
+                )
+            );
+
+            tbf() =
+                scalar(1)
+              - max
                 (
                     min
                     (
@@ -204,8 +225,9 @@ public:
                         scalar(1)
                     ),
                     scalar(0)
-                )
-            );
+                );
+
+            return tbf;
         }
 
 
@@ -216,9 +238,10 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& vf
         ) const
         {
-            surfaceScalarField bf(blendingFactor());
+            surfaceScalarField bf(blendingFactor(vf));
 
-            Info<< "weights " << max(bf) << " " << min(bf) << endl;
+            Info<< "weights " << max(bf).value() << " " << min(bf).value()
+                << endl;
 
             return
                 bf*tScheme1_().weights(vf)
@@ -234,7 +257,7 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& vf
         ) const
         {
-            surfaceScalarField bf(blendingFactor());
+            surfaceScalarField bf(blendingFactor(vf));
 
             return
                 bf*tScheme1_().interpolate(vf)
@@ -257,7 +280,7 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& vf
         ) const
         {
-            surfaceScalarField bf(blendingFactor());
+            surfaceScalarField bf(blendingFactor(vf));
 
             if (tScheme1_().corrected())
             {
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/localBlended/localBlended.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/localBlended/localBlended.H
index 5f49be98fd2348f75b2586ccf225c602fdaa13b3..0e85f0ea7620ff97734e3ec8012a91e420194664 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/localBlended/localBlended.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/localBlended/localBlended.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,7 @@ SourceFiles
 #define localBlended_H
 
 #include "surfaceInterpolationScheme.H"
+#include "blendedSchemeBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,7 +50,8 @@ namespace Foam
 template<class Type>
 class localBlended
 :
-    public surfaceInterpolationScheme<Type>
+    public surfaceInterpolationScheme<Type>,
+    public blendedSchemeBase<Type>
 {
     // Private Member Functions
 
@@ -115,8 +117,27 @@ public:
         {}
 
 
+    //- Destructor
+    virtual ~localBlended()
+    {}
+
+
     // Member Functions
 
+        //- Return the face-based blending factor
+        virtual tmp<surfaceScalarField> blendingFactor
+        (
+             const GeometricField<Type, fvPatchField, volMesh>& vf
+        ) const
+        {
+            return
+                this->mesh().objectRegistry::template
+                    lookupObject<const surfaceScalarField>
+                    (
+                        word(vf.name() + "BlendingFactor")
+                    );
+        }
+
         //- Return the interpolation weighting factors
         tmp<surfaceScalarField> weights
         (
@@ -125,10 +146,10 @@ public:
         {
             const surfaceScalarField& blendingFactor =
                 this->mesh().objectRegistry::template
-                lookupObject<const surfaceScalarField>
-                (
-                    word(vf.name() + "BlendingFactor")
-                );
+                    lookupObject<const surfaceScalarField>
+                    (
+                        word(vf.name() + "BlendingFactor")
+                    );
 
             return
                 blendingFactor*tScheme1_().weights(vf)
diff --git a/src/fvOptions/fvOptions/fvIOoptionList.C b/src/fvOptions/fvOptions/fvIOoptionList.C
index 287e4171be455dd9856a37eeaf0530260324ffe4..0a7ab4e6fc9f50b891253d378699cbc804b0101b 100644
--- a/src/fvOptions/fvOptions/fvIOoptionList.C
+++ b/src/fvOptions/fvOptions/fvIOoptionList.C
@@ -45,7 +45,7 @@ Foam::IOobject Foam::fv::IOoptionList::createIOobject
 
     if (io.headerOk())
     {
-        Info<< "Creating fintite volume options from " << io.name() << nl
+        Info<< "Creating finite volume options from " << io.name() << nl
             << endl;
 
         io.readOpt() = IOobject::MUST_READ_IF_MODIFIED;
diff --git a/src/postProcessing/functionObjects/utilities/Make/files b/src/postProcessing/functionObjects/utilities/Make/files
index f711fc0880616f7edc11475d876a79dac7d76ec1..6a905f28a1e0245d752a0275f83a15469dfdbaa4 100644
--- a/src/postProcessing/functionObjects/utilities/Make/files
+++ b/src/postProcessing/functionObjects/utilities/Make/files
@@ -12,6 +12,9 @@ Peclet/PecletFunctionObject.C
 Q/Q.C
 Q/QFunctionObject.C
 
+blendingFactor/blendingFactor.C
+blendingFactor/blendingFactorFunctionObject.C
+
 DESModelRegions/DESModelRegions.C
 DESModelRegions/DESModelRegionsFunctionObject.C
 
diff --git a/src/postProcessing/functionObjects/utilities/Peclet/Peclet.C b/src/postProcessing/functionObjects/utilities/Peclet/Peclet.C
index 07c7e436b55a4f32b2347162076342c8b0cd25c8..0c77608b13ffb1877685e1884749adfa6b739c0b 100644
--- a/src/postProcessing/functionObjects/utilities/Peclet/Peclet.C
+++ b/src/postProcessing/functionObjects/utilities/Peclet/Peclet.C
@@ -197,12 +197,12 @@ void Foam::Peclet::execute()
 
 void Foam::Peclet::end()
 {
-    // Do nothing - only valid on write
+    // Do nothing
 }
 
 void Foam::Peclet::timeSet()
 {
-    // Do nothing - only valid on write
+    // Do nothing
 }
 
 
diff --git a/src/postProcessing/functionObjects/utilities/blendingFactor/IOblendingFactor.H b/src/postProcessing/functionObjects/utilities/blendingFactor/IOblendingFactor.H
new file mode 100644
index 0000000000000000000000000000000000000000..2d1ae3dc1e4885eeed8272ef93d3331f45a1a17c
--- /dev/null
+++ b/src/postProcessing/functionObjects/utilities/blendingFactor/IOblendingFactor.H
@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Typedef
+    Foam::IOblendingFactor
+
+Description
+    Instance of the generic IOOutputFilter for blendingFactor.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IOblendingFactor_H
+#define IOblendingFactor_H
+
+#include "blendingFactor.H"
+#include "IOOutputFilter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef IOOutputFilter<blendingFactor> IOblendingFactor;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactor.C b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactor.C
new file mode 100644
index 0000000000000000000000000000000000000000..351c48bfb4b6f4994753dcc7221b1764fed87199
--- /dev/null
+++ b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactor.C
@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "blendingFactor.H"
+#include "dictionary.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(blendingFactor, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::blendingFactor::blendingFactor
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    const bool loadFromFiles
+)
+:
+    name_(name),
+    obr_(obr),
+    active_(true),
+    phiName_("unknown-phiName"),
+    fieldName_("unknown-fieldName")
+{
+    // Check if the available mesh is an fvMesh, otherwise deactivate
+    if (!isA<fvMesh>(obr_))
+    {
+        active_ = false;
+        WarningIn
+        (
+            "blendingFactor::blendingFactor"
+            "("
+                "const word&, "
+                "const objectRegistry&, "
+                "const dictionary&, "
+                "const bool"
+            ")"
+        )   << "No fvMesh available, deactivating " << name_ << nl
+            << endl;
+    }
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::blendingFactor::~blendingFactor()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::blendingFactor::read(const dictionary& dict)
+{
+    if (active_)
+    {
+        phiName_ = dict.lookupOrDefault<word>("phiName", "phi");
+        dict.lookup("fieldName") >> fieldName_;
+    }
+}
+
+
+void Foam::blendingFactor::execute()
+{
+    if (active_)
+    {
+        calc<scalar>();
+        calc<vector>();
+    }
+}
+
+
+void Foam::blendingFactor::end()
+{
+    // Do nothing
+}
+
+void Foam::blendingFactor::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::blendingFactor::write()
+{
+    if (active_)
+    {
+        const word fieldName = "blendingFactor:" + fieldName_;
+
+        const volScalarField& blendingFactor =
+            obr_.lookupObject<volScalarField>(fieldName);
+
+        Info<< type() << " " << name_ << " output:" << nl
+            << "    writing field " << blendingFactor.name() << nl
+            << endl;
+
+        blendingFactor.write();
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactor.H b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactor.H
new file mode 100644
index 0000000000000000000000000000000000000000..8b86a25cea60707f2cd5dc7208c6af2344f064b3
--- /dev/null
+++ b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactor.H
@@ -0,0 +1,175 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::blendingFactor
+
+Group
+    grpUtilitiesFunctionObjects
+
+Description
+    This function object calculates and outputs the blendingFactor as used by
+    the bended convection schemes.  The output is a volume field (cells) whose
+    value is calculated via the maximum blending factor for any cell face.
+
+
+SourceFiles
+    blendingFactor.C
+    IOblendingFactor.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef blendingFactor_H
+#define blendingFactor_H
+
+#include "volFieldsFwd.H"
+#include "surfaceFieldsFwd.H"
+#include "OFstream.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class objectRegistry;
+class dictionary;
+class polyMesh;
+class mapPolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                          Class blendingFactor Declaration
+\*---------------------------------------------------------------------------*/
+
+class blendingFactor
+{
+    // Private data
+
+        //- Name of this set of blendingFactor objects
+        word name_;
+
+        //- Reference to the database
+        const objectRegistry& obr_;
+
+        //- On/off switch
+        bool active_;
+
+        //- Name of flux field, default is "phi"
+        word phiName_;
+
+        //- Field name
+        word fieldName_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        blendingFactor(const blendingFactor&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const blendingFactor&);
+
+        //- Return the blending factor field from the database
+        template<class Type>
+        volScalarField& factor
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& field
+        );
+
+        //- Calculate the blending factor
+        template<class Type>
+        void calc();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("blendingFactor");
+
+
+    // Constructors
+
+        //- Construct for given objectRegistry and dictionary.
+        //  Allow the possibility to load fields from files
+        blendingFactor
+        (
+            const word& name,
+            const objectRegistry&,
+            const dictionary&,
+            const bool loadFromFiles = false
+        );
+
+
+    //- Destructor
+    virtual ~blendingFactor();
+
+
+    // Member Functions
+
+        //- Return name of the set of blendingFactor
+        virtual const word& name() const
+        {
+            return name_;
+        }
+
+        //- Read the blendingFactor data
+        virtual void read(const dictionary&);
+
+        //- Execute, currently does nothing
+        virtual void execute();
+
+        //- Execute at the final time-loop, currently does nothing
+        virtual void end();
+
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
+        //- Calculate the blendingFactor and write
+        virtual void write();
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh&)
+        {}
+
+        //- Update for changes of mesh
+        virtual void movePoints(const polyMesh&)
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "blendingFactorTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorFunctionObject.C b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorFunctionObject.C
new file mode 100644
index 0000000000000000000000000000000000000000..6f32ea9c6f5d800583c1f4ef9ac9a387980e25bf
--- /dev/null
+++ b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorFunctionObject.C
@@ -0,0 +1,42 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "blendingFactorFunctionObject.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineNamedTemplateTypeNameAndDebug(blendingFactorFunctionObject, 0);
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        blendingFactorFunctionObject,
+        dictionary
+    );
+}
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorFunctionObject.H b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorFunctionObject.H
new file mode 100644
index 0000000000000000000000000000000000000000..b2e2369ebe191138b2ee919473590a3e917b2096
--- /dev/null
+++ b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorFunctionObject.H
@@ -0,0 +1,54 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Typedef
+    Foam::blendingFactorFunctionObject
+
+Description
+    FunctionObject wrapper around blendingFactor to allow it to be created
+    via the functions entry within controlDict.
+
+SourceFiles
+    blendingFactorFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef blendingFactorFunctionObject_H
+#define blendingFactorFunctionObject_H
+
+#include "blendingFactor.H"
+#include "OutputFilterFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef OutputFilterFunctionObject<blendingFactor>
+        blendingFactorFunctionObject;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorTemplates.C b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..3a284b8c731b0d83025cb989674f7c9c93662135
--- /dev/null
+++ b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorTemplates.C
@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "gaussConvectionScheme.H"
+#include "blendedSchemeBase.H"
+#include "fvcCellReduce.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::volScalarField& Foam::blendingFactor::factor
+(
+    const GeometricField<Type, fvPatchField, volMesh>& field
+)
+{
+    const word fieldName = "blendingFactor:" + field.name();
+
+    if (!obr_.foundObject<volScalarField>(fieldName))
+    {
+        const fvMesh& mesh = refCast<const fvMesh>(obr_);
+
+        volScalarField* factorPtr =
+            new volScalarField
+            (
+                IOobject
+                (
+                    fieldName,
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedScalar("0", dimless, 0.0),
+                zeroGradientFvPatchScalarField::typeName
+            );
+
+        obr_.store(factorPtr);
+    }
+
+    return
+        const_cast<volScalarField&>
+        (
+            obr_.lookupObject<volScalarField>(fieldName)
+        );
+}
+
+
+template<class Type>
+void Foam::blendingFactor::calc()
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if (!obr_.foundObject<fieldType>(fieldName_))
+    {
+        return;
+    }
+
+    const fvMesh& mesh = refCast<const fvMesh>(obr_);
+
+    const fieldType& field = mesh.lookupObject<fieldType>(fieldName_);
+
+    const word divScheme("div(" + phiName_ + ',' + fieldName_ + ')');
+    ITstream& its = mesh.divScheme(divScheme);
+
+    const surfaceScalarField& phi =
+        mesh.lookupObject<surfaceScalarField>(phiName_);
+
+    tmp<fv::convectionScheme<Type> > cs =
+        fv::convectionScheme<Type>::New(mesh, phi, its);
+
+    const fv::gaussConvectionScheme<Type>& gcs =
+        refCast<const fv::gaussConvectionScheme<Type> >(cs());
+
+    const surfaceInterpolationScheme<Type>& interpScheme =
+        gcs.interpScheme();
+
+    if (!isA<blendedSchemeBase<Type> >(interpScheme))
+    {
+        FatalErrorIn("void Foam::blendingFactor::execute()")
+            << interpScheme.typeName << " is not a blended scheme"
+            << exit(FatalError);
+    }
+
+    // retrieve the face-based blending factor
+    const blendedSchemeBase<Type>& blendedScheme =
+        refCast<const blendedSchemeBase<Type> >(interpScheme);
+    const surfaceScalarField factorf(blendedScheme.blendingFactor(field));
+
+    // convert into vol field whose values represent the local face maxima
+    volScalarField& factor = this->factor(field);
+    factor = fvc::cellReduce(factorf, maxEqOp<scalar>());
+    factor.correctBoundaryConditions();
+}
+
+
+// ************************************************************************* //
diff --git a/src/sixDoFRigidBodyMotion/Make/files b/src/sixDoFRigidBodyMotion/Make/files
index 022baa43df6e0ebe4c8f36a2193cf8ca7c56418c..75435920ea441c167af884370570c97ad5eae050 100644
--- a/src/sixDoFRigidBodyMotion/Make/files
+++ b/src/sixDoFRigidBodyMotion/Make/files
@@ -3,7 +3,7 @@ sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
 sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.C
 sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateIO.C
 
-restraints = sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint
+restraints = sixDoFRigidBodyMotion/restraints
 
 $(restraints)/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C
 $(restraints)/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraintNew.C
@@ -12,7 +12,7 @@ $(restraints)/linearSpring/linearSpring.C
 $(restraints)/sphericalAngularSpring/sphericalAngularSpring.C
 $(restraints)/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C
 
-constraints = sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint
+constraints = sixDoFRigidBodyMotion/constraints
 
 $(constraints)/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C
 $(constraints)/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraintNew.C
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedAxis/fixedAxis.C
similarity index 60%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.C
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedAxis/fixedAxis.C
index 0e2737842535eba35bbd5a64b08e37c9a1a2b701..db7b66a6ff993e33abb9c248b43a83e00e774c8d 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedAxis/fixedAxis.C
@@ -49,10 +49,11 @@ namespace sixDoFRigidBodyMotionConstraints
 
 Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::fixedAxis
 (
+    const word& name,
     const dictionary& sDoFRBMCDict
 )
 :
-    sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
+    sixDoFRigidBodyMotionConstraint(name, sDoFRBMCDict),
     fixedAxis_()
 {
     read(sDoFRBMCDict);
@@ -67,79 +68,19 @@ Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::~fixedAxis()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-bool Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::constrain
+void Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::constrainTranslation
 (
-    const sixDoFRigidBodyMotion& motion,
-    const vector& existingConstraintForce,
-    const vector& existingConstraintMoment,
-    scalar deltaT,
-    vector& constraintPosition,
-    vector& constraintForceIncrement,
-    vector& constraintMomentIncrement
+    pointConstraint& pc
 ) const
-{
-    constraintMomentIncrement = vector::zero;
-
-    vector predictedDir = motion.predictedOrientation
-    (
-        fixedAxis_,
-        existingConstraintMoment,
-        deltaT
-    );
-
-    scalar theta = acos(min(predictedDir & fixedAxis_, 1.0));
-
-    vector rotationAxis = fixedAxis_ ^ predictedDir;
-
-    scalar magRotationAxis = mag(rotationAxis);
-
-    if (magRotationAxis > VSMALL)
-    {
-        rotationAxis /= magRotationAxis;
-
-        const tensor& Q = motion.orientation();
-
-        // Transform rotationAxis to body local system
-        rotationAxis = Q.T() & rotationAxis;
-
-        constraintMomentIncrement =
-           -relaxationFactor_
-           *(motion.momentOfInertia() & rotationAxis)
-           *theta/sqr(deltaT);
-
-        // Transform moment increment to global system
-        constraintMomentIncrement = Q & constraintMomentIncrement;
-
-        // Remove any moment that is around the fixedAxis_
-        constraintMomentIncrement -=
-            (constraintMomentIncrement & fixedAxis_)*fixedAxis_;
-    }
-
-    constraintPosition = motion.centreOfMass();
-
-    constraintForceIncrement = vector::zero;
-
-    bool converged(mag(theta) < tolerance_);
+{}
 
-    if (sixDoFRigidBodyMotionConstraint::debug)
-    {
-        Info<< " angle " << theta
-            << " force " << constraintForceIncrement
-            << " moment " << constraintMomentIncrement;
-
-        if (converged)
-        {
-            Info<< " converged";
-        }
-        else
-        {
-            Info<< " not converged";
-        }
-
-        Info<< endl;
-    }
 
-    return converged;
+void Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::constrainRotation
+(
+    pointConstraint& pc
+) const
+{
+    pc.combine(pointConstraint(Tuple2<label, vector>(2, vector(0,1,0))));
 }
 
 
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedAxis/fixedAxis.H
similarity index 83%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.H
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedAxis/fixedAxis.H
index 7772833f4d53fccb4f81be74b404cc92c9433e2e..e2a2cc26da0ebff96788e1188869933228072678 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedAxis/fixedAxis.H
@@ -37,8 +37,6 @@ SourceFiles
 #define fixedAxis_H
 
 #include "sixDoFRigidBodyMotionConstraint.H"
-#include "point.H"
-#include "tensor.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -74,6 +72,7 @@ public:
         //- Construct from components
         fixedAxis
         (
+            const word& name,
             const dictionary& sDoFRBMCDict
         );
 
@@ -93,19 +92,11 @@ public:
 
     // Member Functions
 
-        //- Calculate the constraint position, force and moment.
-        //  Global reference frame vectors.  Returns boolean stating
-        //  whether the constraint been converged to tolerance.
-        virtual bool constrain
-        (
-            const sixDoFRigidBodyMotion& motion,
-            const vector& existingConstraintForce,
-            const vector& existingConstraintMoment,
-            scalar deltaT,
-            vector& constraintPosition,
-            vector& constraintForceIncrement,
-            vector& constraintMomentIncrement
-        ) const;
+        //- Apply and accumulate translational constraints
+        virtual void constrainTranslation(pointConstraint&) const;
+
+        //- Apply and accumulate rotational constraints
+        virtual void constrainRotation(pointConstraint&) const;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& sDoFRBMCCoeff);
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedLine/fixedLine.C
similarity index 63%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.C
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedLine/fixedLine.C
index f77afa2b10ade1dc45492a1995ab8dd67f2fd1da..5aa4d669daa7d1ef5400ee1ac884bdc09e831320 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedLine/fixedLine.C
@@ -49,11 +49,11 @@ namespace sixDoFRigidBodyMotionConstraints
 
 Foam::sixDoFRigidBodyMotionConstraints::fixedLine::fixedLine
 (
+    const word& name,
     const dictionary& sDoFRBMCDict
 )
 :
-    sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
-    refPt_(),
+    sixDoFRigidBodyMotionConstraint(name, sDoFRBMCDict),
     dir_()
 {
     read(sDoFRBMCDict);
@@ -68,65 +68,20 @@ Foam::sixDoFRigidBodyMotionConstraints::fixedLine::~fixedLine()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-bool Foam::sixDoFRigidBodyMotionConstraints::fixedLine::constrain
+void Foam::sixDoFRigidBodyMotionConstraints::fixedLine::constrainTranslation
 (
-    const sixDoFRigidBodyMotion& motion,
-    const vector& existingConstraintForce,
-    const vector& existingConstraintMoment,
-    scalar deltaT,
-    vector& constraintPosition,
-    vector& constraintForceIncrement,
-    vector& constraintMomentIncrement
+    pointConstraint& pc
 ) const
 {
-    point predictedPosition = motion.predictedPosition
-    (
-        refPt_,
-        existingConstraintForce,
-        existingConstraintMoment,
-        deltaT
-    );
-
-    constraintPosition = motion.currentPosition(refPt_);
-
-    // Info<< "current position " << constraintPosition << nl
-    //     << "next predictedPosition " << predictedPosition
-    //     << endl;
-
-    // Vector from reference point to predicted point
-    vector rC = predictedPosition - refPt_;
-
-    vector error = rC - ((rC) & dir_)*dir_;
-
-    // Info<< "error " << error << endl;
-
-    constraintForceIncrement =
-        -relaxationFactor_*error*motion.mass()/sqr(deltaT);
-
-    constraintMomentIncrement = vector::zero;
-
-    bool converged(mag(error) < tolerance_);
+    pc.combine(pointConstraint(Tuple2<label, vector>(2, dir_)));
+}
 
-    if (sixDoFRigidBodyMotionConstraint::debug)
-    {
-        Info<< " error " << error
-            << " force " << constraintForceIncrement
-            << " moment " << constraintMomentIncrement;
-
-        if (converged)
-        {
-            Info<< " converged";
-        }
-        else
-        {
-            Info<< " not converged";
-        }
-
-        Info<< endl;
-    }
 
-    return converged;
-}
+void Foam::sixDoFRigidBodyMotionConstraints::fixedLine::constrainRotation
+(
+    pointConstraint& pc
+) const
+{}
 
 
 bool Foam::sixDoFRigidBodyMotionConstraints::fixedLine::read
@@ -136,8 +91,6 @@ bool Foam::sixDoFRigidBodyMotionConstraints::fixedLine::read
 {
     sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict);
 
-    sDoFRBMCCoeffs_.lookup("refPoint") >> refPt_;
-
     sDoFRBMCCoeffs_.lookup("direction") >> dir_;
 
     scalar magDir(mag(dir_));
@@ -168,9 +121,6 @@ void Foam::sixDoFRigidBodyMotionConstraints::fixedLine::write
     Ostream& os
 ) const
 {
-    os.writeKeyword("refPoint")
-        << refPt_ << token::END_STATEMENT << nl;
-
     os.writeKeyword("direction")
         << dir_ << token::END_STATEMENT << nl;
 }
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedLine/fixedLine.H
similarity index 81%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.H
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedLine/fixedLine.H
index 8d0710e0908780deb5364a99e63c4f5ab695e388..2e48f564a74a068d1642b308cea38d4b900f9d52 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedLine/fixedLine.H
@@ -37,7 +37,6 @@ SourceFiles
 #define fixedLine_H
 
 #include "sixDoFRigidBodyMotionConstraint.H"
-#include "point.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -57,9 +56,6 @@ class fixedLine
 {
     // Private data
 
-        //- Reference point for the constraining line
-        point refPt_;
-
         //- Direction of the constraining line
         vector dir_;
 
@@ -75,6 +71,7 @@ public:
         //- Construct from components
         fixedLine
         (
+            const word& name,
             const dictionary& sDoFRBMCDict
         );
 
@@ -94,19 +91,11 @@ public:
 
     // Member Functions
 
-        //- Calculate the constraint position, force and moment.
-        //  Global reference frame vectors.  Returns boolean stating
-        //  whether the constraint been converged to tolerance.
-        virtual bool constrain
-        (
-            const sixDoFRigidBodyMotion& motion,
-            const vector& existingConstraintForce,
-            const vector& existingConstraintMoment,
-            scalar deltaT,
-            vector& constraintPosition,
-            vector& constraintForceIncrement,
-            vector& constraintMomentIncrement
-        ) const;
+        //- Apply and accumulate translational constraints
+        virtual void constrainTranslation(pointConstraint&) const;
+
+        //- Apply and accumulate rotational constraints
+        virtual void constrainRotation(pointConstraint&) const;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& sDoFRBMCCoeff);
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedOrientation/fixedOrientation.C
similarity index 51%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.C
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedOrientation/fixedOrientation.C
index c66bdbb3623383e897bf77d359c2ca110314030e..5d4cc8cff045eb32764069b25d880556dc7bacc1 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedOrientation/fixedOrientation.C
@@ -49,10 +49,11 @@ namespace sixDoFRigidBodyMotionConstraints
 
 Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::fixedOrientation
 (
+    const word& name,
     const dictionary& sDoFRBMCDict
 )
 :
-    sixDoFRigidBodyMotionConstraint(sDoFRBMCDict)
+    sixDoFRigidBodyMotionConstraint(name, sDoFRBMCDict)
 {
     read(sDoFRBMCDict);
 }
@@ -66,102 +67,21 @@ Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::~fixedOrientation()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-bool Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::constrain
+void Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::
+constrainTranslation
 (
-    const sixDoFRigidBodyMotion& motion,
-    const vector& existingConstraintForce,
-    const vector& existingConstraintMoment,
-    scalar deltaT,
-    vector& constraintPosition,
-    vector& constraintForceIncrement,
-    vector& constraintMomentIncrement
+    pointConstraint& pc
 ) const
-{
-    constraintMomentIncrement = vector::zero;
-
-    scalar maxTheta = -SMALL;
-
-    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
-    {
-        vector axis = vector::zero;
-
-        axis[cmpt] = 1;
-
-        vector refDir = vector::zero;
-
-        refDir[(cmpt + 1) % 3] = 1;
-
-        vector predictedDir = motion.predictedOrientation
-        (
-            refDir,
-            existingConstraintMoment,
-            deltaT
-        );
-
-        // Removing any axis component from predictedDir
-        predictedDir -= (axis & predictedDir)*axis;
-
-        scalar theta = GREAT;
-
-        scalar magPredictedDir = mag(predictedDir);
-
-        if (magPredictedDir > VSMALL)
-        {
-            predictedDir /= magPredictedDir;
-
-            theta = acos(min(predictedDir & refDir, 1.0));
-
-            // Recalculating axis to give correct sign to angle
-            axis = (refDir ^ predictedDir);
-
-            scalar magAxis = mag(axis);
-
-            if (magAxis > VSMALL)
-            {
-                axis /= magAxis;
-            }
-            else
-            {
-                axis = vector::zero;
-            }
-        }
-
-        if (theta > maxTheta)
-        {
-            maxTheta = theta;
-        }
-
-        constraintMomentIncrement +=
-           -relaxationFactor_
-           *theta*axis
-           *motion.momentOfInertia()[cmpt]/sqr(deltaT);
-    }
-
-    constraintPosition = motion.centreOfMass();
-
-    constraintForceIncrement = vector::zero;
-
-    bool converged(mag(maxTheta) < tolerance_);
-
-    if (sixDoFRigidBodyMotionConstraint::debug)
-    {
-        Info<< " max angle " << maxTheta
-            << " force " << constraintForceIncrement
-            << " moment " << constraintMomentIncrement;
-
-        if (converged)
-        {
-            Info<< " converged";
-        }
-        else
-        {
-            Info<< " not converged";
-        }
+{}
 
-        Info<< endl;
-    }
 
-    return converged;
+void Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::
+constrainRotation
+(
+    pointConstraint& pc
+) const
+{
+    pc.combine(pointConstraint(Tuple2<label, vector>(3, vector::zero)));
 }
 
 
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedOrientation/fixedOrientation.H
similarity index 79%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.H
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedOrientation/fixedOrientation.H
index fcf7a62167ac24d29205cc07a68cd133c0d4a89a..b41945837a82eaa85779db2d9df568e27d06df28 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedOrientation/fixedOrientation.H
@@ -25,9 +25,7 @@ Class
     Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation
 
 Description
-    sixDoFRigidBodyMotionConstraint.  Orientation of body fixed global
-    space. Only valid where the predicted deviation from alignment is
-    < 90 degrees.
+    Fix orientation of body in global space.
 
 SourceFiles
     fixedOrientation.C
@@ -38,8 +36,6 @@ SourceFiles
 #define fixedOrientation_H
 
 #include "sixDoFRigidBodyMotionConstraint.H"
-#include "point.H"
-#include "tensor.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -69,6 +65,7 @@ public:
         //- Construct from components
         fixedOrientation
         (
+            const word& name,
             const dictionary& sDoFRBMCDict
         );
 
@@ -88,19 +85,11 @@ public:
 
     // Member Functions
 
-        //- Calculate the constraint position, force and moment.
-        //  Global reference frame vectors.  Returns boolean stating
-        //  whether the constraint been converged to tolerance.
-        virtual bool constrain
-        (
-            const sixDoFRigidBodyMotion& motion,
-            const vector& existingConstraintForce,
-            const vector& existingConstraintMoment,
-            scalar deltaT,
-            vector& constraintPosition,
-            vector& constraintForceIncrement,
-            vector& constraintMomentIncrement
-        ) const;
+        //- Apply and accumulate translational constraints
+        virtual void constrainTranslation(pointConstraint&) const;
+
+        //- Apply and accumulate rotational constraints
+        virtual void constrainRotation(pointConstraint&) const;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& sDoFRBMCCoeff);
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedPlane/fixedPlane.C
similarity index 56%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.C
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedPlane/fixedPlane.C
index 263694ee3f2ee1a8c449055844088dcb1632e88f..edcde4d07c22a35b40ee01d6b2cc2dc00e3526d0 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedPlane/fixedPlane.C
@@ -49,11 +49,12 @@ namespace sixDoFRigidBodyMotionConstraints
 
 Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::fixedPlane
 (
+    const word& name,
     const dictionary& sDoFRBMCDict
 )
 :
-    sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
-    fixedPlane_(vector::one)
+    sixDoFRigidBodyMotionConstraint(name, sDoFRBMCDict),
+    normal_(vector::zero)
 {
     read(sDoFRBMCDict);
 }
@@ -67,66 +68,20 @@ Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::~fixedPlane()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-bool Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::constrain
+void Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::constrainTranslation
 (
-    const sixDoFRigidBodyMotion& motion,
-    const vector& existingConstraintForce,
-    const vector& existingConstraintMoment,
-    scalar deltaT,
-    vector& constraintPosition,
-    vector& constraintForceIncrement,
-    vector& constraintMomentIncrement
+    pointConstraint& pc
 ) const
 {
-    const point& refPt = fixedPlane_.refPoint();
-
-    const vector& n = fixedPlane_.normal();
-
-    point predictedPosition = motion.predictedPosition
-    (
-        refPt,
-        existingConstraintForce,
-        existingConstraintMoment,
-        deltaT
-    );
-
-    constraintPosition = motion.currentPosition(refPt);
-
-    // Info<< "current position " << constraintPosition << nl
-    //     << "next predictedPosition " << predictedPosition
-    //     << endl;
-
-    vector error = ((predictedPosition - refPt) & n)*n;
-
-    // Info<< "error " << error << endl;
-
-    constraintForceIncrement =
-        -relaxationFactor_*error*motion.mass()/sqr(deltaT);
-
-    constraintMomentIncrement = vector::zero;
-
-    bool converged(mag(error) < tolerance_);
-
-    if (sixDoFRigidBodyMotionConstraint::debug)
-    {
-        Info<< " error " << error
-            << " force " << constraintForceIncrement
-            << " moment " << constraintMomentIncrement;
-
-        if (converged)
-        {
-            Info<< " converged";
-        }
-        else
-        {
-            Info<< " not converged";
-        }
+    pc.applyConstraint(normal_);
+}
 
-        Info<< endl;
-    }
 
-    return converged;
-}
+void Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::constrainRotation
+(
+    pointConstraint& pc
+) const
+{}
 
 
 bool Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::read
@@ -136,11 +91,7 @@ bool Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::read
 {
     sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict);
 
-    point refPt = sDoFRBMCCoeffs_.lookup("refPoint");
-
-    vector normal = sDoFRBMCCoeffs_.lookup("normal");
-
-    fixedPlane_ = plane(refPt, normal);
+    normal_ = sDoFRBMCCoeffs_.lookup("normal");
 
     return true;
 }
@@ -151,11 +102,8 @@ void Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::write
     Ostream& os
 ) const
 {
-    os.writeKeyword("refPoint")
-        << fixedPlane_.refPoint() << token::END_STATEMENT << nl;
-
     os.writeKeyword("normal")
-        << fixedPlane_.normal() << token::END_STATEMENT << nl;
+        << normal_ << token::END_STATEMENT << nl;
 }
 
 // ************************************************************************* //
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedPlane/fixedPlane.H
similarity index 80%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.H
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedPlane/fixedPlane.H
index b7f100ac5eae4e7dfd6d0d556dbe702fe8cdf055..82ed9fab9e714dfd32c32556711fb95d6cad256b 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedPlane/fixedPlane.H
@@ -37,7 +37,6 @@ SourceFiles
 #define fixedPlane_H
 
 #include "sixDoFRigidBodyMotionConstraint.H"
-#include "plane.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -57,9 +56,8 @@ class fixedPlane
 {
     // Private data
 
-        //- Plane which the transformed reference point is constrained
-        //  to move along
-        plane fixedPlane_;
+        //- Normal to plane
+        vector normal_;
 
 
 public:
@@ -73,6 +71,7 @@ public:
         //- Construct from components
         fixedPlane
         (
+            const word& name,
             const dictionary& sDoFRBMCDict
         );
 
@@ -92,19 +91,11 @@ public:
 
     // Member Functions
 
-        //- Calculate the constraint position, force and moment.
-        //  Global reference frame vectors.  Returns boolean stating
-        //  whether the constraint been converged to tolerance.
-        virtual bool constrain
-        (
-            const sixDoFRigidBodyMotion& motion,
-            const vector& existingConstraintForce,
-            const vector& existingConstraintMoment,
-            scalar deltaT,
-            vector& constraintPosition,
-            vector& constraintForceIncrement,
-            vector& constraintMomentIncrement
-        ) const;
+        //- Apply and accumulate translational constraints
+        virtual void constrainTranslation(pointConstraint&) const;
+
+        //- Apply and accumulate rotational constraints
+        virtual void constrainRotation(pointConstraint&) const;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& sDoFRBMCCoeff);
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedPoint/fixedPoint.C
similarity index 56%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.C
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedPoint/fixedPoint.C
index 084334868dbee3ab2cc9a8f20b9625e29c22666c..c0f4bbf311665af81d24f8a07832bebf62c3e562 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedPoint/fixedPoint.C
@@ -49,10 +49,11 @@ namespace sixDoFRigidBodyMotionConstraints
 
 Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::fixedPoint
 (
+    const word& name,
     const dictionary& sDoFRBMCDict
 )
 :
-    sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
+    sixDoFRigidBodyMotionConstraint(name, sDoFRBMCDict),
     fixedPoint_()
 {
     read(sDoFRBMCDict);
@@ -67,75 +68,20 @@ Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::~fixedPoint()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-bool Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::constrain
+void Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::constrainTranslation
 (
-    const sixDoFRigidBodyMotion& motion,
-    const vector& existingConstraintForce,
-    const vector& existingConstraintMoment,
-    scalar deltaT,
-    vector& constraintPosition,
-    vector& constraintForceIncrement,
-    vector& constraintMomentIncrement
+    pointConstraint& pc
 ) const
 {
-    point predictedPosition = motion.predictedPosition
-    (
-        fixedPoint_,
-        existingConstraintForce,
-        existingConstraintMoment,
-        deltaT
-    );
-
-    constraintPosition = motion.currentPosition(fixedPoint_);
-
-    // Info<< "current position " << constraintPosition << nl
-    //     << "next predictedPosition " << predictedPosition
-    //     << endl;
-
-    vector error = predictedPosition - fixedPoint_;
-
-    // Info<< "error " << error << endl;
-
-    // Correction force derived from Lagrange multiplier:
-    //     G = -lambda*grad(sigma)
-    // where
-    //     sigma = mag(error) = 0
-    // so
-    //     grad(sigma) = error/mag(error)
-    // Solving for lambda using the SHAKE methodology gives
-    //     lambda = mass*mag(error)/sqr(deltaT)
-    // This is only strictly applicable (i.e. will converge in one
-    // iteration) to constraints at the centre of mass.  Everything
-    // else will need to iterate, and may need under-relaxed to be
-    // stable.
-
-    constraintForceIncrement =
-        -relaxationFactor_*error*motion.mass()/sqr(deltaT);
-
-    constraintMomentIncrement = vector::zero;
-
-    bool converged(mag(error) < tolerance_);
-
-    if (sixDoFRigidBodyMotionConstraint::debug)
-    {
-        Info<< " error " << error
-            << " force " << constraintForceIncrement
-            << " moment " << constraintMomentIncrement;
-
-        if (converged)
-        {
-            Info<< " converged";
-        }
-        else
-        {
-            Info<< " not converged";
-        }
+    pc.combine(pointConstraint(Tuple2<label, vector>(3, vector::zero)));
+}
 
-        Info<< endl;
-    }
 
-    return converged;
-}
+void Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::constrainRotation
+(
+    pointConstraint& pc
+) const
+{}
 
 
 bool Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::read
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedPoint/fixedPoint.H
similarity index 83%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.H
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedPoint/fixedPoint.H
index c361d3a8b3d0ad730ec293f2fdc39b02c50f379a..4745adc8a9e6af28fcb9d5c8778ff54bd7fb23b6 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/fixedPoint/fixedPoint.H
@@ -25,7 +25,7 @@ Class
     Foam::sixDoFRigidBodyMotionConstraints::fixedPoint
 
 Description
-    sixDoFRigidBodyMotionConstraint.  Point fixed in space.
+    Point fixed in space.
 
 SourceFiles
     fixedPoint.C
@@ -74,6 +74,7 @@ public:
         //- Construct from components
         fixedPoint
         (
+            const word& name,
             const dictionary& sDoFRBMCDict
         );
 
@@ -93,19 +94,11 @@ public:
 
     // Member Functions
 
-        //- Calculate the constraint position, force and moment.
-        //  Global reference frame vectors.  Returns boolean stating
-        //  whether the constraint been converged to tolerance.
-        virtual bool constrain
-        (
-            const sixDoFRigidBodyMotion& motion,
-            const vector& existingConstraintForce,
-            const vector& existingConstraintMoment,
-            scalar deltaT,
-            vector& constraintPosition,
-            vector& constraintForceIncrement,
-            vector& constraintMomentIncrement
-        ) const;
+        //- Apply and accumulate translational constraints
+        virtual void constrainTranslation(pointConstraint&) const;
+
+        //- Apply and accumulate rotational constraints
+        virtual void constrainRotation(pointConstraint&) const;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& sDoFRBMCCoeff);
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C
similarity index 77%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C
index 22299c7f3c8cc3f2dc1b7fca94aa554a3c975665..6a497d4de5abe9bf8407e1c5c2634f98362ea1c1 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C
@@ -29,9 +29,8 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(sixDoFRigidBodyMotionConstraint, 0);
-
-defineRunTimeSelectionTable(sixDoFRigidBodyMotionConstraint, dictionary);
+    defineTypeNameAndDebug(sixDoFRigidBodyMotionConstraint, 0);
+    defineRunTimeSelectionTable(sixDoFRigidBodyMotionConstraint, dictionary);
 }
 
 
@@ -39,9 +38,11 @@ defineRunTimeSelectionTable(sixDoFRigidBodyMotionConstraint, dictionary);
 
 Foam::sixDoFRigidBodyMotionConstraint::sixDoFRigidBodyMotionConstraint
 (
+    const word& name,
     const dictionary& sDoFRBMCDict
 )
 :
+    name_(name),
     sDoFRBMCCoeffs_
     (
         sDoFRBMCDict.subDict
@@ -49,11 +50,6 @@ Foam::sixDoFRigidBodyMotionConstraint::sixDoFRigidBodyMotionConstraint
             word(sDoFRBMCDict.lookup("sixDoFRigidBodyMotionConstraint"))
           + "Coeffs"
         )
-    ),
-    tolerance_(readScalar(sDoFRBMCDict.lookup("tolerance"))),
-    relaxationFactor_
-    (
-        sDoFRBMCDict.lookupOrDefault<scalar>("relaxationFactor", 1)
     )
 {}
 
@@ -71,14 +67,6 @@ bool Foam::sixDoFRigidBodyMotionConstraint::read
     const dictionary& sDoFRBMCDict
 )
 {
-    tolerance_ = (readScalar(sDoFRBMCDict.lookup("tolerance")));
-
-    relaxationFactor_ = sDoFRBMCDict.lookupOrDefault<scalar>
-    (
-        "relaxationFactor",
-        1
-    );
-
     sDoFRBMCCoeffs_ = sDoFRBMCDict.subDict(type() + "Coeffs");
 
     return true;
@@ -86,12 +74,7 @@ bool Foam::sixDoFRigidBodyMotionConstraint::read
 
 
 void Foam::sixDoFRigidBodyMotionConstraint::write(Ostream& os) const
-{
-    os.writeKeyword("tolerance")
-        << tolerance_ << token::END_STATEMENT << nl;
+{}
 
-    os.writeKeyword("relaxationFactor")
-        << relaxationFactor_ << token::END_STATEMENT << nl;
-}
 
 // ************************************************************************* //
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H
similarity index 75%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H
index b1f5438948cf5ad3bd255f1e403db4a28d4f90ad..b4eed214f9db1ad4eb8c883862f0127cc35c283d 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H
@@ -46,7 +46,7 @@ SourceFiles
 #include "Time.H"
 #include "dictionary.H"
 #include "autoPtr.H"
-#include "vector.H"
+#include "pointConstraint.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -69,16 +69,12 @@ protected:
 
     // Protected data
 
+        //- Name of the constraint
+        word name_;
+
         //- Constraint model specific coefficient dictionary
         dictionary sDoFRBMCCoeffs_;
 
-        //- Solution tolerance.  Meaning depends on model, usually an
-        //  absolute distance or angle.
-        scalar tolerance_;
-
-        //- Relaxation factor for solution, default to one
-        scalar relaxationFactor_;
-
 
 public:
 
@@ -93,8 +89,8 @@ public:
             autoPtr,
             sixDoFRigidBodyMotionConstraint,
             dictionary,
-            (const dictionary& sDoFRBMCDict),
-            (sDoFRBMCDict)
+            (const word& name, const dictionary& sDoFRBMCDict),
+            (name, sDoFRBMCDict)
         );
 
 
@@ -103,6 +99,7 @@ public:
         //- Construct from the sDoFRBMCDict dictionary and Time
         sixDoFRigidBodyMotionConstraint
         (
+            const word& name,
             const dictionary& sDoFRBMCDict
         );
 
@@ -115,6 +112,7 @@ public:
         //- Select constructed from the sDoFRBMCDict dictionary and Time
         static autoPtr<sixDoFRigidBodyMotionConstraint> New
         (
+            const word& name,
             const dictionary& sDoFRBMCDict
         );
 
@@ -125,19 +123,17 @@ public:
 
     // Member Functions
 
-        //- Calculate the constraint position, force and moment.
-        //  Global reference frame vectors.  Returns boolean stating
-        //  whether the constraint been converged to tolerance.
-        virtual bool constrain
-        (
-            const sixDoFRigidBodyMotion& motion,
-            const vector& existingConstraintForce,
-            const vector& existingConstraintMoment,
-            scalar deltaT,
-            vector& constraintPosition,
-            vector& constraintForceIncrement,
-            vector& constraintMomentIncrement
-        ) const = 0;
+        //- Return the name
+        const word& name() const
+        {
+            return name_;
+        }
+
+        //- Apply and accumulate translational constraints
+        virtual void constrainTranslation(pointConstraint&) const = 0;
+
+        //- Apply and accumulate rotational constraints
+        virtual void constrainRotation(pointConstraint&) const = 0;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& sDoFRBMCDict);
@@ -150,18 +146,6 @@ public:
                 return sDoFRBMCCoeffs_;
             }
 
-            //- Return access to the tolerance
-            inline scalar tolerance() const
-            {
-                return tolerance_;
-            }
-
-            //- Return access to the relaxationFactor
-            inline scalar relaxationFactor() const
-            {
-                return relaxationFactor_;
-            }
-
         //- Write
         virtual void write(Ostream&) const;
 };
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraintNew.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraintNew.C
similarity index 89%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraintNew.C
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraintNew.C
index fa74999aa5d8d3ec9dd785422758ade35fba4de7..b6f3e8fb1f195cee52efe82a3e65f2e993a6c592 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraintNew.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraintNew.C
@@ -28,16 +28,17 @@ License
 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
 
 Foam::autoPtr<Foam::sixDoFRigidBodyMotionConstraint>
-Foam::sixDoFRigidBodyMotionConstraint::New(const dictionary& sDoFRBMCDict)
+Foam::sixDoFRigidBodyMotionConstraint::New
+(
+    const word& name,
+    const dictionary& sDoFRBMCDict
+)
 {
     const word constraintType
     (
         sDoFRBMCDict.lookup("sixDoFRigidBodyMotionConstraint")
     );
 
-    // Info<< "Selecting sixDoFRigidBodyMotionConstraint function "
-    //     << constraintType << endl;
-
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(constraintType);
 
@@ -56,7 +57,10 @@ Foam::sixDoFRigidBodyMotionConstraint::New(const dictionary& sDoFRBMCDict)
             << exit(FatalError);
     }
 
-    return autoPtr<sixDoFRigidBodyMotionConstraint>(cstrIter()(sDoFRBMCDict));
+    return autoPtr<sixDoFRigidBodyMotionConstraint>
+    (
+        cstrIter()(name, sDoFRBMCDict)
+    );
 }
 
 
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/linearAxialAngularSpring/linearAxialAngularSpring.C
similarity index 98%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/linearAxialAngularSpring/linearAxialAngularSpring.C
index 76805d60d47a8bf1772957292ca5a1a88ef9169a..3701dd8bf951461d88b77fbfce32fa001bfc0aea 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/linearAxialAngularSpring/linearAxialAngularSpring.C
@@ -51,10 +51,11 @@ namespace sixDoFRigidBodyMotionRestraints
 Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::
 linearAxialAngularSpring
 (
+    const word& name,
     const dictionary& sDoFRBMRDict
 )
 :
-    sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
+    sixDoFRigidBodyMotionRestraint(name, sDoFRBMRDict),
     refQ_(),
     axis_(),
     stiffness_(),
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/linearAxialAngularSpring/linearAxialAngularSpring.H
similarity index 99%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.H
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/linearAxialAngularSpring/linearAxialAngularSpring.H
index d7f8c5e2805bc572d86f902dade04204c4efd074..4a30e30c65c14d576a49e57a93118057630370fb 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/linearAxialAngularSpring/linearAxialAngularSpring.H
@@ -81,6 +81,7 @@ public:
         //- Construct from components
         linearAxialAngularSpring
         (
+            const word& name,
             const dictionary& sDoFRBMRDict
         );
 
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/linearSpring/linearSpring.C
similarity index 98%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.C
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/linearSpring/linearSpring.C
index c6b1c10ca063592136c08861faa11022d9006566..a9ed4ce7a3f413eaa936e88be8a6ceeb0c8970fa 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/linearSpring/linearSpring.C
@@ -49,10 +49,11 @@ namespace sixDoFRigidBodyMotionRestraints
 
 Foam::sixDoFRigidBodyMotionRestraints::linearSpring::linearSpring
 (
+    const word& name,
     const dictionary& sDoFRBMRDict
 )
 :
-    sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
+    sixDoFRigidBodyMotionRestraint(name, sDoFRBMRDict),
     anchor_(),
     refAttachmentPt_(),
     stiffness_(),
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/linearSpring/linearSpring.H
similarity index 99%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.H
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/linearSpring/linearSpring.H
index c0610d974b39d9843580e78c81c473206e6370cd..dbd12372c3a928556a1ccdc6760be3b0194b925b 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/linearSpring/linearSpring.H
@@ -84,6 +84,7 @@ public:
         //- Construct from components
         linearSpring
         (
+            const word& name,
             const dictionary& sDoFRBMRDict
         );
 
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C
similarity index 98%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C
index b06c5baf465f9262ec34ae7d69934c60c3584b38..31c16f85acdcc2d3e3870aa8378697038371d0f4 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C
@@ -39,9 +39,11 @@ defineRunTimeSelectionTable(sixDoFRigidBodyMotionRestraint, dictionary);
 
 Foam::sixDoFRigidBodyMotionRestraint::sixDoFRigidBodyMotionRestraint
 (
+    const word& name,
     const dictionary& sDoFRBMRDict
 )
 :
+    name_(name),
     sDoFRBMRCoeffs_
     (
         sDoFRBMRDict.subDict
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H
similarity index 92%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H
index d8bdf1469a3d60b404ce8ab0f25157002ec1723a..41216a0c2ba55a29eff5f340dd5d9155e647ff36 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H
@@ -69,6 +69,9 @@ protected:
 
     // Protected data
 
+        //- Name of the restraint
+        word name_;
+
         //- Restraint model specific coefficient dictionary
         dictionary sDoFRBMRCoeffs_;
 
@@ -86,8 +89,8 @@ public:
             autoPtr,
             sixDoFRigidBodyMotionRestraint,
             dictionary,
-            (const dictionary& sDoFRBMRDict),
-            (sDoFRBMRDict)
+            (const word& name, const dictionary& sDoFRBMRDict),
+            (name, sDoFRBMRDict)
         );
 
 
@@ -96,6 +99,7 @@ public:
         //- Construct from the sDoFRBMRDict dictionary and Time
         sixDoFRigidBodyMotionRestraint
         (
+            const word& name,
             const dictionary& sDoFRBMRDict
         );
 
@@ -108,6 +112,7 @@ public:
         //- Select constructed from the sDoFRBMRDict dictionary and Time
         static autoPtr<sixDoFRigidBodyMotionRestraint> New
         (
+            const word& name,
             const dictionary& sDoFRBMRDict
         );
 
@@ -118,6 +123,12 @@ public:
 
     // Member Functions
 
+        //- Return the name
+        const word& name() const
+        {
+            return name_;
+        }
+
         //- Calculate the restraint position, force and moment.
         //  Global reference frame vectors.
         virtual void restrain
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraintNew.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraintNew.C
similarity index 89%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraintNew.C
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraintNew.C
index 9361e0ab4df1c3b788435bb944d02d6ca50563da..1bb047a77d6d270f40a6fd121bbbee03c70099b5 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraintNew.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraintNew.C
@@ -28,16 +28,17 @@ License
 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
 
 Foam::autoPtr<Foam::sixDoFRigidBodyMotionRestraint>
-Foam::sixDoFRigidBodyMotionRestraint::New(const dictionary& sDoFRBMRDict)
+Foam::sixDoFRigidBodyMotionRestraint::New
+(
+    const word& name,
+    const dictionary& sDoFRBMRDict
+)
 {
     const word restraintType
     (
         sDoFRBMRDict.lookup("sixDoFRigidBodyMotionRestraint")
     );
 
-    // Info<< "Selecting sixDoFRigidBodyMotionRestraint function "
-    //     << restraintType << endl;
-
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(restraintType);
 
@@ -56,7 +57,10 @@ Foam::sixDoFRigidBodyMotionRestraint::New(const dictionary& sDoFRBMRDict)
             << exit(FatalError);
     }
 
-    return autoPtr<sixDoFRigidBodyMotionRestraint>(cstrIter()(sDoFRBMRDict));
+    return autoPtr<sixDoFRigidBodyMotionRestraint>
+    (
+        cstrIter()(name, sDoFRBMRDict)
+    );
 }
 
 
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sphericalAngularSpring/sphericalAngularSpring.C
similarity index 98%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.C
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sphericalAngularSpring/sphericalAngularSpring.C
index dbb480b93aa692727e6ee5ca9b45cdc572809e0c..f3dab425d7c7351ae7c67c61faa73298e9991d5f 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sphericalAngularSpring/sphericalAngularSpring.C
@@ -50,10 +50,11 @@ namespace sixDoFRigidBodyMotionRestraints
 Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::
 sphericalAngularSpring
 (
+    const word& name,
     const dictionary& sDoFRBMRDict
 )
 :
-    sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
+    sixDoFRigidBodyMotionRestraint(name, sDoFRBMRDict),
     refQ_(),
     stiffness_(),
     damping_()
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sphericalAngularSpring/sphericalAngularSpring.H
similarity index 99%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.H
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sphericalAngularSpring/sphericalAngularSpring.H
index a36ef5d4ab104f3e1f2683184cd8496b73cb33f6..f1bfa784d371d1e78d0eb77890fd5c7736919da4 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/sphericalAngularSpring/sphericalAngularSpring.H
@@ -78,6 +78,7 @@ public:
         //- Construct from components
         sphericalAngularSpring
         (
+            const word& name,
             const dictionary& sDoFRBMRDict
         );
 
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C
similarity index 98%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C
index 252aa471f7958fab956c2887b1aa6e209e4537f0..037639dde8560abb21d27b0d3a9f0625db0f1fc9 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C
@@ -52,10 +52,11 @@ namespace sixDoFRigidBodyMotionRestraints
 Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::
 tabulatedAxialAngularSpring
 (
+    const word& name,
     const dictionary& sDoFRBMRDict
 )
 :
-    sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
+    sixDoFRigidBodyMotionRestraint(name, sDoFRBMRDict),
     refQ_(),
     axis_(),
     moment_(),
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.H
similarity index 99%
rename from src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.H
rename to src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.H
index bb45f00d519361353833a4cb02574ccac7219757..bcd2c2c5b3f05362479d26028d2b698bcbe3aba9 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/restraints/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.H
@@ -88,6 +88,7 @@ public:
         //- Construct from components
         tabulatedAxialAngularSpring
         (
+            const word& name,
             const dictionary& sDoFRBMRDict
         );
 
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
index 6f002c5e2fec80cf2d5113c0de2b817a89db08b5..b3a013328c2e1c59a490453bfd6f39670fe26c58 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
@@ -41,7 +41,7 @@ void Foam::sixDoFRigidBodyMotion::applyRestraints()
         {
             if (report_)
             {
-                Info<< "Restraint " << restraintNames_[rI] << ": ";
+                Info<< "Restraint " << restraints_[rI].name() << ": ";
             }
 
             // restraint position
@@ -65,101 +65,6 @@ void Foam::sixDoFRigidBodyMotion::applyRestraints()
 }
 
 
-void Foam::sixDoFRigidBodyMotion::applyConstraints(scalar deltaT)
-{
-    if (constraints_.empty())
-    {
-        return;
-    }
-
-    if (Pstream::master())
-    {
-        label iteration = 0;
-
-        bool allConverged = true;
-
-        // constraint force accumulator
-        vector cFA = vector::zero;
-
-        // constraint moment accumulator
-        vector cMA = vector::zero;
-
-        do
-        {
-            allConverged = true;
-
-            forAll(constraints_, cI)
-            {
-                if (sixDoFRigidBodyMotionConstraint::debug)
-                {
-                    Info<< "Constraint " << constraintNames_[cI] << ": ";
-                }
-
-                // constraint position
-                point cP = vector::zero;
-
-                // constraint force
-                vector cF = vector::zero;
-
-                // constraint moment
-                vector cM = vector::zero;
-
-                bool constraintConverged = constraints_[cI].constrain
-                (
-                    *this,
-                    cFA,
-                    cMA,
-                    deltaT,
-                    cP,
-                    cF,
-                    cM
-                );
-
-                allConverged = allConverged && constraintConverged;
-
-                // Accumulate constraint force
-                cFA += cF;
-
-                // Accumulate constraint moment
-                cMA += cM + ((cP - centreOfMass()) ^ cF);
-            }
-
-        } while(++iteration < maxConstraintIterations_ && !allConverged);
-
-        if (iteration >= maxConstraintIterations_)
-        {
-            FatalErrorIn
-            (
-                "Foam::sixDoFRigidBodyMotion::applyConstraints(scalar)"
-            )
-                << nl << "Maximum number of sixDoFRigidBodyMotion constraint "
-                << "iterations ("
-                << maxConstraintIterations_
-                << ") exceeded." << nl
-                << exit(FatalError);
-        }
-
-        Info<< "sixDoFRigidBodyMotion constraints converged in "
-            << iteration << " iterations" << endl;
-
-        if (report_)
-        {
-            Info<< "Constraint force: " << cFA << nl
-                << "Constraint moment: " << cMA
-                << endl;
-        }
-
-        // Add the constraint forces and moments to the motion state variables
-        a() += cFA/mass_;
-
-        // The moment of constraint forces has already been added
-        // during accumulation.  Moments are returned in global axes,
-        // transforming to body local
-        tau() += Q().T() & cMA;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
@@ -167,10 +72,9 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
     motionState_(),
     motionState0_(),
     restraints_(),
-    restraintNames_(),
     constraints_(),
-    constraintNames_(),
-    maxConstraintIterations_(0),
+    tConstraints_(tensor::I),
+    rConstraints_(tensor::I),
     initialCentreOfMass_(vector::zero),
     initialQ_(I),
     momentOfInertia_(diagTensor::one*VSMALL),
@@ -209,10 +113,9 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     ),
     motionState0_(motionState_),
     restraints_(),
-    restraintNames_(),
     constraints_(),
-    constraintNames_(),
-    maxConstraintIterations_(0),
+    tConstraints_(tensor::I),
+    rConstraints_(tensor::I),
     initialCentreOfMass_(initialCentreOfMass),
     initialQ_(initialQ),
     momentOfInertia_(momentOfInertia),
@@ -232,10 +135,9 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     motionState_(stateDict),
     motionState0_(motionState_),
     restraints_(),
-    restraintNames_(),
     constraints_(),
-    constraintNames_(),
-    maxConstraintIterations_(0),
+    tConstraints_(tensor::I),
+    rConstraints_(tensor::I),
     initialCentreOfMass_
     (
         dict.lookupOrDefault
@@ -271,10 +173,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     motionState_(sDoFRBM.motionState_),
     motionState0_(sDoFRBM.motionState0_),
     restraints_(sDoFRBM.restraints_),
-    restraintNames_(sDoFRBM.restraintNames_),
     constraints_(sDoFRBM.constraints_),
-    constraintNames_(sDoFRBM.constraintNames_),
-    maxConstraintIterations_(sDoFRBM.maxConstraintIterations_),
     initialCentreOfMass_(sDoFRBM.initialCentreOfMass_),
     initialQ_(sDoFRBM.initialQ_),
     momentOfInertia_(sDoFRBM.momentOfInertia_),
@@ -306,29 +205,23 @@ void Foam::sixDoFRigidBodyMotion::addRestraints
 
         restraints_.setSize(restraintDict.size());
 
-        restraintNames_.setSize(restraintDict.size());
-
         forAllConstIter(IDLList<entry>, restraintDict, iter)
         {
             if (iter().isDict())
             {
-                // Info<< "Adding restraint: " << iter().keyword() << endl;
-
                 restraints_.set
                 (
-                    i,
-                    sixDoFRigidBodyMotionRestraint::New(iter().dict())
+                    i++,
+                    sixDoFRigidBodyMotionRestraint::New
+                    (
+                        iter().keyword(),
+                        iter().dict()
+                    )
                 );
-
-                restraintNames_[i] = iter().keyword();
-
-                i++;
             }
         }
 
         restraints_.setSize(i);
-
-        restraintNames_.setSize(i);
     }
 }
 
@@ -346,21 +239,25 @@ void Foam::sixDoFRigidBodyMotion::addConstraints
 
         constraints_.setSize(constraintDict.size());
 
-        constraintNames_.setSize(constraintDict.size());
+        pointConstraint pct;
+        pointConstraint pcr;
 
         forAllConstIter(IDLList<entry>, constraintDict, iter)
         {
             if (iter().isDict())
             {
-                // Info<< "Adding constraint: " << iter().keyword() << endl;
-
                 constraints_.set
                 (
                     i,
-                    sixDoFRigidBodyMotionConstraint::New(iter().dict())
+                    sixDoFRigidBodyMotionConstraint::New
+                    (
+                        iter().keyword(),
+                        iter().dict()
+                    )
                 );
 
-                constraintNames_[i] = iter().keyword();
+                constraints_[i].constrainTranslation(pct);
+                constraints_[i].constrainRotation(pcr);
 
                 i++;
             }
@@ -368,15 +265,11 @@ void Foam::sixDoFRigidBodyMotion::addConstraints
 
         constraints_.setSize(i);
 
-        constraintNames_.setSize(i);
+        tConstraints_ = pct.constraintTransformation();
+        rConstraints_ = pcr.constraintTransformation();
 
-        if (!constraints_.empty())
-        {
-            maxConstraintIterations_ = readLabel
-            (
-                constraintDict.lookup("maxIterations")
-            );
-        }
+        Info<< "Translational constraint tensor " << tConstraints_ << nl
+            << "Rotational constraint tensor " << rConstraints_ << endl;
     }
 }
 
@@ -392,8 +285,8 @@ void Foam::sixDoFRigidBodyMotion::updatePosition
 
     if (Pstream::master())
     {
-        v() = v0() + aDamp_*0.5*deltaT0*a();
-        pi() = pi0() + aDamp_*0.5*deltaT0*tau();
+        v() = tConstraints_ & aDamp_*(v0() + 0.5*deltaT0*a());
+        pi() = rConstraints_ & aDamp_*(pi0() + 0.5*deltaT0*tau());
 
         // Leapfrog move part
         centreOfMass() = centreOfMass0() + deltaT*v();
@@ -401,7 +294,7 @@ void Foam::sixDoFRigidBodyMotion::updatePosition
         // Leapfrog orientation adjustment
         Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
         Q() = Qpi.first();
-        pi() = Qpi.second();
+        pi() = rConstraints_ & Qpi.second();
     }
 
     Pstream::scatter(motionState_);
@@ -439,12 +332,9 @@ void Foam::sixDoFRigidBodyMotion::updateAcceleration
         }
         first = false;
 
-        // Apply constraints after relaxation
-        applyConstraints(deltaT);
-
         // Correct velocities
-        v() += aDamp_*0.5*deltaT*a();
-        pi() += aDamp_*0.5*deltaT*tau();
+        v() += tConstraints_ & aDamp_*0.5*deltaT*a();
+        pi() += rConstraints_ & aDamp_*0.5*deltaT*tau();
 
         if (report_)
         {
@@ -463,8 +353,8 @@ void Foam::sixDoFRigidBodyMotion::updateVelocity(scalar deltaT)
 
     if (Pstream::master())
     {
-        v() += aDamp_*0.5*deltaT*a();
-        pi() += aDamp_*0.5*deltaT*tau();
+        v() += tConstraints_ & aDamp_*0.5*deltaT*a();
+        pi() += rConstraints_ & aDamp_*0.5*deltaT*tau();
 
         if (report_)
         {
@@ -532,7 +422,10 @@ Foam::vector Foam::sixDoFRigidBodyMotion::predictedOrientation
     vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment));
     Tuple2<tensor, vector> QpiTemp = rotate(Q0(), piTemp, deltaT);
 
-    return (QpiTemp.first() & initialQ_.T() & vInitial);
+    vector o(QpiTemp.first() & initialQ_.T() & vInitial);
+    o /= mag(o);
+
+    return o;
 }
 
 
@@ -564,8 +457,6 @@ Foam::tmp<Foam::pointField> Foam::sixDoFRigidBodyMotion::scaledPosition
     const scalarField& scale
 ) const
 {
-    Info<< "initialCentreOfMass " << initialCentreOfMass() << endl;
-
     // Calculate the transformation septerion from the initial state
     septernion s
     (
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
index 670004ba693e0eba489e2d7204ee1d6e8cd37b0b..0e4ab39da268d4d82762f72be90fd1de5640065d 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
@@ -90,21 +90,17 @@ class sixDoFRigidBodyMotion
         //- Motion state data object for previous time-step
         sixDoFRigidBodyMotionState motionState0_;
 
-        //- Restraints on the motion
+        //- Motion restraints
         PtrList<sixDoFRigidBodyMotionRestraint> restraints_;
 
-        //- Names of the restraints
-        wordList restraintNames_;
-
-        //- Constaints on the motion
+        //- Motion constaints
         PtrList<sixDoFRigidBodyMotionConstraint> constraints_;
 
-        //- Names of the constraints
-        wordList constraintNames_;
+        //- Translational constraint tensor
+        tensor tConstraints_;
 
-        //- Maximum number of iterations allowed to attempt to obey
-        //  constraints
-        label maxConstraintIterations_;
+        //- Rotational constraint tensor
+        tensor rConstraints_;
 
         //- Centre of mass of initial state
         point initialCentreOfMass_;
@@ -155,9 +151,6 @@ class sixDoFRigidBodyMotion
         //- Apply the restraints to the object
         void applyRestraints();
 
-        //- Apply the constraints to the object
-        void applyConstraints(scalar deltaT);
-
 
         // Access functions retained as private because of the risk of
         // confusion over what is a body local frame vector and what is global
@@ -168,20 +161,10 @@ class sixDoFRigidBodyMotion
             inline const PtrList<sixDoFRigidBodyMotionRestraint>&
                 restraints() const;
 
-            //- Return access to the restraintNames
-            inline const wordList& restraintNames() const;
-
             //- Return access to the constraints
             inline const PtrList<sixDoFRigidBodyMotionConstraint>&
                 constraints() const;
 
-            //- Return access to the constraintNames
-            inline const wordList& constraintNames() const;
-
-            //- Return access to the maximum allowed number of
-            //  constraint iterations
-            inline label maxConstraintIterations() const;
-
             //- Return access to the initial centre of mass
             inline const point& initialCentreOfMass() const;
 
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
index f70ba68ced1cc79c20dcc26b009130d4734201b4..3186ca2fa1251f190104e49dcbfbb368e59df12d 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
@@ -111,12 +111,6 @@ Foam::sixDoFRigidBodyMotion::restraints() const
 }
 
 
-inline const Foam::wordList& Foam::sixDoFRigidBodyMotion::restraintNames() const
-{
-    return restraintNames_;
-}
-
-
 inline const Foam::PtrList<Foam::sixDoFRigidBodyMotionConstraint>&
 Foam::sixDoFRigidBodyMotion::constraints() const
 {
@@ -124,19 +118,6 @@ Foam::sixDoFRigidBodyMotion::constraints() const
 }
 
 
-inline const Foam::wordList&
-Foam::sixDoFRigidBodyMotion::constraintNames() const
-{
-    return constraintNames_;
-}
-
-
-inline Foam::label Foam::sixDoFRigidBodyMotion::maxConstraintIterations() const
-{
-    return maxConstraintIterations_;
-}
-
-
 inline const Foam::point&
 Foam::sixDoFRigidBodyMotion::initialCentreOfMass() const
 {
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
index 99bdba3ef022f370b9c096ab7482db75b2ede074..11670ba6231ceb64a269a506cbbaaac5d2397d35 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
@@ -54,7 +54,7 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
         {
             word restraintType = restraints_[rI].type();
 
-            os  << indent << restraintNames_[rI] << nl
+            os  << indent << restraints_[rI].name() << nl
                 << indent << token::BEGIN_BLOCK << incrIndent << endl;
 
             os.writeKeyword("sixDoFRigidBodyMotionRestraint")
@@ -79,14 +79,11 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
         os  << indent << "constraints" << nl
             << indent << token::BEGIN_BLOCK << incrIndent << nl;
 
-        os.writeKeyword("maxIterations")
-            << maxConstraintIterations_ << token::END_STATEMENT << nl;
-
         forAll(constraints_, rI)
         {
             word constraintType = constraints_[rI].type();
 
-            os  << indent << constraintNames_[rI] << nl
+            os  << indent << constraints_[rI].name() << nl
                 << indent << token::BEGIN_BLOCK << incrIndent << endl;
 
             os.writeKeyword("sixDoFRigidBodyMotionConstraint")
diff --git a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C
index 0327dd1da758b485c23cef211752c8b97c248f3e..4b51db9f17b1d1e0e57712a502e7ab078a962d4e 100644
--- a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C
@@ -263,9 +263,9 @@ void Foam::externalWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
                 forAll (thicknessLayers_, iLayer)
                 {
                     const scalar l = thicknessLayers_[iLayer];
-                    if (l > 0.0)
+                    if (kappaLayers_[iLayer] > 0.0)
                     {
-                        totalSolidRes += kappaLayers_[iLayer]/l;
+                        totalSolidRes += l/kappaLayers_[iLayer];
                     }
                 }
             }
diff --git a/tutorials/incompressible/pimpleFoam/pitzDaily/0/U b/tutorials/incompressible/pimpleFoam/pitzDaily/0/U
new file mode 100644
index 0000000000000000000000000000000000000000..885c5b9de439fa3983f723c56d7734369a6248ad
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/pitzDaily/0/U
@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform (10 0 0);
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    upperWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    lowerWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/pitzDaily/0/epsilon b/tutorials/incompressible/pimpleFoam/pitzDaily/0/epsilon
new file mode 100644
index 0000000000000000000000000000000000000000..d82c45e6290be093bee3d9a3d6af3f84072a543d
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/pitzDaily/0/epsilon
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 14.855;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 14.855;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    upperWall
+    {
+        type            epsilonWallFunction;
+        value           uniform 14.855;
+    }
+    lowerWall
+    {
+        type            epsilonWallFunction;
+        value           uniform 14.855;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/pitzDaily/0/k b/tutorials/incompressible/pimpleFoam/pitzDaily/0/k
new file mode 100644
index 0000000000000000000000000000000000000000..7de1adf4b5d6b4518dc306b6d4ab6d0f40ea2312
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/pitzDaily/0/k
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.375;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.375;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    upperWall
+    {
+        type            kqRWallFunction;
+        value           uniform 0.375;
+    }
+    lowerWall
+    {
+        type            kqRWallFunction;
+        value           uniform 0.375;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/pitzDaily/0/nuTilda b/tutorials/incompressible/pimpleFoam/pitzDaily/0/nuTilda
new file mode 100644
index 0000000000000000000000000000000000000000..cbe04420b91623c7b0736f7f818d8ac8c59634d8
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/pitzDaily/0/nuTilda
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      nuTilda;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    upperWall
+    {
+        type            zeroGradient;
+    }
+
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/pitzDaily/0/nut b/tutorials/incompressible/pimpleFoam/pitzDaily/0/nut
new file mode 100644
index 0000000000000000000000000000000000000000..585d8ec0f4c35a1a0d5c797222096d9483947cce
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/pitzDaily/0/nut
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    upperWall
+    {
+        type            nutkWallFunction;
+        value           uniform 0;
+    }
+    lowerWall
+    {
+        type            nutkWallFunction;
+        value           uniform 0;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/pitzDaily/0/p b/tutorials/incompressible/pimpleFoam/pitzDaily/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..0fdd33ec47f28523f2ae5b0964b19321c1d92f4d
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/pitzDaily/0/p
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    upperWall
+    {
+        type            zeroGradient;
+    }
+
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/pitzDaily/constant/RASProperties b/tutorials/incompressible/pimpleFoam/pitzDaily/constant/RASProperties
new file mode 100644
index 0000000000000000000000000000000000000000..a4937b503a46850b2626f0d301e4a07b9f691507
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/pitzDaily/constant/RASProperties
@@ -0,0 +1,25 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel        kEpsilon;
+
+turbulence      on;
+
+printCoeffs     on;
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/pitzDaily/constant/polyMesh/blockMeshDict b/tutorials/incompressible/pimpleFoam/pitzDaily/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..0e0fcf1c9855da788effb6cfe4f7773c99bf2594
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/pitzDaily/constant/polyMesh/blockMeshDict
@@ -0,0 +1,173 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.001;
+
+vertices
+(
+    (-20.6 0 -0.5)
+    (-20.6 3 -0.5)
+    (-20.6 12.7 -0.5)
+    (-20.6 25.4 -0.5)
+    (0 -25.4 -0.5)
+    (0 -5 -0.5)
+    (0 0 -0.5)
+    (0 3 -0.5)
+    (0 12.7 -0.5)
+    (0 25.4 -0.5)
+    (206 -25.4 -0.5)
+    (206 -8.5 -0.5)
+    (206 0 -0.5)
+    (206 6.5 -0.5)
+    (206 17 -0.5)
+    (206 25.4 -0.5)
+    (290 -16.6 -0.5)
+    (290 -6.3 -0.5)
+    (290 0 -0.5)
+    (290 4.5 -0.5)
+    (290 11 -0.5)
+    (290 16.6 -0.5)
+    (-20.6 0 0.5)
+    (-20.6 3 0.5)
+    (-20.6 12.7 0.5)
+    (-20.6 25.4 0.5)
+    (0 -25.4 0.5)
+    (0 -5 0.5)
+    (0 0 0.5)
+    (0 3 0.5)
+    (0 12.7 0.5)
+    (0 25.4 0.5)
+    (206 -25.4 0.5)
+    (206 -8.5 0.5)
+    (206 0 0.5)
+    (206 6.5 0.5)
+    (206 17 0.5)
+    (206 25.4 0.5)
+    (290 -16.6 0.5)
+    (290 -6.3 0.5)
+    (290 0 0.5)
+    (290 4.5 0.5)
+    (290 11 0.5)
+    (290 16.6 0.5)
+);
+
+blocks
+(
+    hex (0 6 7 1 22 28 29 23) (18 7 1) simpleGrading (0.5 1.8 1)
+    hex (1 7 8 2 23 29 30 24) (18 10 1) simpleGrading (0.5 4 1)
+    hex (2 8 9 3 24 30 31 25) (18 13 1) simpleGrading (0.5 0.25 1)
+    hex (4 10 11 5 26 32 33 27) (180 18 1) simpleGrading (4 1 1)
+    hex (5 11 12 6 27 33 34 28) (180 9 1) edgeGrading (4 4 4 4 0.5 1 1 0.5 1 1 1 1)
+    hex (6 12 13 7 28 34 35 29) (180 7 1) edgeGrading (4 4 4 4 1.8 1 1 1.8 1 1 1 1)
+    hex (7 13 14 8 29 35 36 30) (180 10 1) edgeGrading (4 4 4 4 4 1 1 4 1 1 1 1)
+    hex (8 14 15 9 30 36 37 31) (180 13 1) simpleGrading (4 0.25 1)
+    hex (10 16 17 11 32 38 39 33) (25 18 1) simpleGrading (2.5 1 1)
+    hex (11 17 18 12 33 39 40 34) (25 9 1) simpleGrading (2.5 1 1)
+    hex (12 18 19 13 34 40 41 35) (25 7 1) simpleGrading (2.5 1 1)
+    hex (13 19 20 14 35 41 42 36) (25 10 1) simpleGrading (2.5 1 1)
+    hex (14 20 21 15 36 42 43 37) (25 13 1) simpleGrading (2.5 0.25 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 22 23 1)
+            (1 23 24 2)
+            (2 24 25 3)
+        );
+    }
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (16 17 39 38)
+            (17 18 40 39)
+            (18 19 41 40)
+            (19 20 42 41)
+            (20 21 43 42)
+        );
+    }
+    upperWall
+    {
+        type wall;
+        faces
+        (
+            (3 25 31 9)
+            (9 31 37 15)
+            (15 37 43 21)
+        );
+    }
+    lowerWall
+    {
+        type wall;
+        faces
+        (
+            (0 6 28 22)
+            (6 5 27 28)
+            (5 4 26 27)
+            (4 10 32 26)
+            (10 16 38 32)
+        );
+    }
+    frontAndBack
+    {
+        type empty;
+        faces
+        (
+            (22 28 29 23)
+            (23 29 30 24)
+            (24 30 31 25)
+            (26 32 33 27)
+            (27 33 34 28)
+            (28 34 35 29)
+            (29 35 36 30)
+            (30 36 37 31)
+            (32 38 39 33)
+            (33 39 40 34)
+            (34 40 41 35)
+            (35 41 42 36)
+            (36 42 43 37)
+            (0 1 7 6)
+            (1 2 8 7)
+            (2 3 9 8)
+            (4 5 11 10)
+            (5 6 12 11)
+            (6 7 13 12)
+            (7 8 14 13)
+            (8 9 15 14)
+            (10 11 17 16)
+            (11 12 18 17)
+            (12 13 19 18)
+            (13 14 20 19)
+            (14 15 21 20)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/pitzDaily/constant/polyMesh/boundary b/tutorials/incompressible/pimpleFoam/pitzDaily/constant/polyMesh/boundary
new file mode 100644
index 0000000000000000000000000000000000000000..221823b308d5063df215d77e90a11eae26d6c1c6
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/pitzDaily/constant/polyMesh/boundary
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.2.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+5
+(
+    inlet
+    {
+        type            patch;
+        nFaces          30;
+        startFace       24170;
+    }
+    outlet
+    {
+        type            patch;
+        nFaces          57;
+        startFace       24200;
+    }
+    upperWall
+    {
+        type            wall;
+        nFaces          223;
+        startFace       24257;
+    }
+    lowerWall
+    {
+        type            wall;
+        nFaces          250;
+        startFace       24480;
+    }
+    frontAndBack
+    {
+        type            empty;
+        inGroups        1(empty);
+        nFaces          24450;
+        startFace       24730;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/pitzDaily/constant/transportProperties b/tutorials/incompressible/pimpleFoam/pitzDaily/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..b40b7d66cd884b7a54d4c7a61b50b1e39a466150
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/pitzDaily/constant/transportProperties
@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              nu [ 0 2 -1 0 0 0 0 ] 1e-05;
+
+CrossPowerLawCoeffs
+{
+    nu0             nu0 [ 0 2 -1 0 0 0 0 ] 1e-06;
+    nuInf           nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
+    m               m [ 0 0 1 0 0 0 0 ] 1;
+    n               n [ 0 0 0 0 0 0 0 ] 1;
+}
+
+BirdCarreauCoeffs
+{
+    nu0             nu0 [ 0 2 -1 0 0 0 0 ] 1e-06;
+    nuInf           nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
+    k               k [ 0 0 1 0 0 0 0 ] 0;
+    n               n [ 0 0 0 0 0 0 0 ] 1;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/pitzDaily/constant/turbulenceProperties b/tutorials/incompressible/pimpleFoam/pitzDaily/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..e7fa28c74a5fead3fbcdd79b5587ef684e8bacb4
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/pitzDaily/constant/turbulenceProperties
@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType RASModel;
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/pitzDaily/system/controlDict b/tutorials/incompressible/pimpleFoam/pitzDaily/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..a8ea900b3c47194aaf01f6803065fe4306bc6395
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/pitzDaily/system/controlDict
@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     pimpleFoam;
+
+startFrom       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1;
+
+deltaT          0.0001;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.001;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           5;
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/pitzDaily/system/fvSchemes b/tutorials/incompressible/pimpleFoam/pitzDaily/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..926d489aae55eb3d059e41849197d58c6129a2b3
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/pitzDaily/system/fvSchemes
@@ -0,0 +1,71 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+    grad(U)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      bounded Gauss linearUpwind grad(U);
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,R)      bounded Gauss upwind;
+    div(R)          Gauss linear;
+    div(phi,nuTilda) bounded Gauss upwind;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         none;
+    laplacian(nuEff,U) Gauss linear corrected;
+    laplacian(rAUf,p)  Gauss linear corrected;
+    laplacian(DkEff,k) Gauss linear corrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+    laplacian(DREff,R) Gauss linear corrected;
+    laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+    interpolate(U)  linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/pitzDaily/system/fvSolution b/tutorials/incompressible/pimpleFoam/pitzDaily/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..880d7bef7b27665cfb7675d2036a2b2faf848267
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/pitzDaily/system/fvSolution
@@ -0,0 +1,62 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver           GAMG;
+        tolerance        1e-7;
+        relTol           0.01;
+
+        smoother         GaussSeidel;
+
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator     faceAreaPair;
+        mergeLevels      1;
+    }
+
+    pFinal
+    {
+        $p;
+        relTol          0;
+    }
+
+    "(U|k|epsilon)"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    "(U|k|epsilon)Final"
+    {
+        $U;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    nCorrectors         2;
+}
+
+
+// ************************************************************************* //