Merge branch 'master' of github.com:OpenFOAM/OpenFOAM-2.3.x

README.html

@@ -26,10 +26,10 @@
   .center {margin-left:auto; margin-right:auto; text-align:center;}
   p.verse { margin-left: 3% }
   pre {
-	border: 1pt solid #AEBDCC;
+        border: 1pt solid #AEBDCC;
-	background-color: #F3F5F7;
+        background-color: #F3F5F7;
-	padding: 5pt;
+        padding: 5pt;
-	font-family: courier, monospace;
+        font-family: courier, monospace;
         font-size: 90%;
         overflow:auto;
   }
@@ -137,7 +137,7 @@ for the JavaScript code in this tag.
 <h2 id="sec-1"><span class="section-number-2">1</span> About OpenFOAM</h2>
 <div class="outline-text-2" id="text-1">
 
-<p>  OpenFOAM is a free, open source computational fluid dynamcis (CFD) software
+<p>  OpenFOAM is a free, open source computational fluid dynamics (CFD) software
   package released by the OpenFOAM Foundation. It has a large user base across
   most areas of engineering and science, from both commercial and academic
   organisations. OpenFOAM has an extensive range of features to solve anything

README.org

@@ -8,7 +8,7 @@
 # Copyright (c) 2014 OpenFOAM Foundation.
 
 * About OpenFOAM
-  OpenFOAM is a free, open source computational fluid dynamcis (CFD) software
+  OpenFOAM is a free, open source computational fluid dynamics (CFD) software
   package released by the OpenFOAM Foundation. It has a large user base across
   most areas of engineering and science, from both commercial and academic
   organisations. OpenFOAM has an extensive range of features to solve anything

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/TomiyamaAspectRatio/TomiyamaAspectRatio.C

@@ -52,7 +52,7 @@ Foam::aspectRatioModels::TomiyamaAspectRatio::TomiyamaAspectRatio
     const orderedPhasePair& pair
 )
 :
-    aspectRatioModel(dict, pair),
+    VakhrushevEfremov(dict, pair),
     yWall_(pair.phase1().mesh().lookupObject<volScalarField>("yWall"))
 {}
 
@@ -69,7 +69,7 @@ Foam::tmp<Foam::volScalarField>
 Foam::aspectRatioModels::TomiyamaAspectRatio::E() const
 {
     return
-        pair_.Eo()
+        VakhrushevEfremov::E()
        *max
        (
            scalar(1) - 0.35*yWall_/pair_.dispersed().d(),

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/TomiyamaAspectRatio/TomiyamaAspectRatio.H

@@ -44,7 +44,7 @@ SourceFiles
 #ifndef TomiyamaAspectRatio_H
 #define TomiyamaAspectRatio_H
 
-#include "aspectRatioModel.H"
+#include "VakhrushevEfremov.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,7 +59,7 @@ namespace aspectRatioModels
 
 class TomiyamaAspectRatio
 :
-    public aspectRatioModel
+    public VakhrushevEfremov
 {
 private:
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/TomiyamaAnalytic/TomiyamaAnalytic.C

@@ -69,17 +69,17 @@ Foam::dragModels::TomiyamaAnalytic::CdRe() const
     volScalarField Eo(max(pair_.Eo(), residualEo_));
     volScalarField E(max(pair_.E(), residualE_));
 
-    volScalarField OmEsq(max(scalar(1) - sqr(E), residualE_));
+    volScalarField OmEsq(max(scalar(1) - sqr(E), sqr(residualE_)));
     volScalarField rtOmEsq(sqrt(OmEsq));
 
-    volScalarField F((asin(rtOmEsq) - E*rtOmEsq)/OmEsq);
+    volScalarField F(max(asin(rtOmEsq) - E*rtOmEsq, residualE_)/OmEsq);
 
     return
         (8.0/3.0)
        *Eo
        /(
             Eo*pow(E, 2.0/3.0)/OmEsq
-          + 16*pow(E, 0.75)
+          + 16*pow(E, 4.0/3.0)
         )
        /sqr(F)
        *max(pair_.Re(), residualRe_);

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/Moraga/Moraga.C

@@ -62,15 +62,39 @@ Foam::liftModels::Moraga::~Moraga()
 
 Foam::tmp<Foam::volScalarField> Foam::liftModels::Moraga::Cl() const
 {
-    volScalarField ReSqrSr
+    volScalarField Re(pair_.Re());
+
+    volScalarField sqrSr
     (
-        pair_.Re()
+        sqr(pair_.dispersed().d())
-       *sqr(pair_.dispersed().d())
        /pair_.continuous().nu()
        *mag(fvc::grad(pair_.continuous().U()))
     );
 
-    return 0.2*exp(- ReSqrSr/3.6e5 - 0.12)*exp(ReSqrSr/3.0e7);
+    if
+    (
+        min(Re).value() < 1200.0
+     || max(Re).value() > 18800.0
+     || min(sqrSr).value() < 0.0016
+     || max(sqrSr).value() > 0.04
+    )
+    {
+        WarningIn
+        (
+            "Foam::tmp<Foam::volScalarField> "
+            "Foam::liftModels::Moraga::Cl() const"
+        )   << "Re and/or Sr are out of the range of applicability of the "
+            << "Moraga model. Clamping to range bounds"
+            << endl;
+    }
+
+    Re.min(1200.0);
+    Re.max(18800.0);
+
+    sqrSr.min(0.0016);
+    sqrSr.max(0.04);
+
+    return 0.2*exp(- Re*sqrSr/3.6e5 - 0.12)*exp(Re*sqrSr/3.0e7);
 }
 
 

applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C

@@ -19,7 +19,7 @@ Foam::label Foam::findOppositeWedge
 {
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-    scalar wppCosAngle = wpp.centreNormal()&wpp.patchNormal();
+    scalar wppCosAngle = wpp.cosAngle();
 
     forAll(patches, patchI)
     {
@@ -30,13 +30,11 @@ Foam::label Foam::findOppositeWedge
          && isA<wedgePolyPatch>(patches[patchI])
         )
         {
-            const wedgePolyPatch& pp = refCast<const wedgePolyPatch>
+            const wedgePolyPatch& pp =
-            (
+                refCast<const wedgePolyPatch>(patches[patchI]);
-                patches[patchI]
-            );
 
             // Calculate (cos of) angle to wpp (not pp!) centre normal
-            scalar ppCosAngle = wpp.centreNormal()&pp.patchNormal();
+            scalar ppCosAngle = wpp.centreNormal() & pp.n();
 
             if
             (
@@ -73,12 +71,10 @@ bool Foam::checkWedges
     {
         if (patches[patchI].size() && isA<wedgePolyPatch>(patches[patchI]))
         {
-            const wedgePolyPatch& pp = refCast<const wedgePolyPatch>
+            const wedgePolyPatch& pp =
-            (
+                refCast<const wedgePolyPatch>(patches[patchI]);
-                patches[patchI]
-            );
 
-            scalar wedgeAngle = acos(pp.centreNormal()&pp.patchNormal());
+            scalar wedgeAngle = acos(pp.cosAngle());
 
             if (report)
             {
@@ -100,10 +96,8 @@ bool Foam::checkWedges
                 return true;
             }
 
-            const wedgePolyPatch& opp = refCast<const wedgePolyPatch>
+            const wedgePolyPatch& opp =
-            (
+                refCast<const wedgePolyPatch>(patches[oppositePatchI]);
-                patches[oppositePatchI]
-            );
 
 
             if (mag(opp.axis() & pp.axis()) < (1-1e-3))
@@ -140,7 +134,7 @@ bool Foam::checkWedges
             forAll(pp.meshPoints(), i)
             {
                 const point& pt = p[pp.meshPoints()[i]];
-                scalar d = mag((pt-p0) & pp.patchNormal());
+                scalar d = mag((pt - p0) & pp.n());
 
                 if (d > sqrt(SMALL))
                 {
@@ -385,10 +379,8 @@ bool Foam::checkCoupledPoints
     {
         if (patches[patchI].coupled())
         {
-            const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
+            const coupledPolyPatch& cpp =
-            (
+                refCast<const coupledPolyPatch>(patches[patchI]);
-                patches[patchI]
-            );
 
             if (cpp.owner())
             {

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,9 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    Wedge front and back plane patch
-
 \*---------------------------------------------------------------------------*/
 
 #include "wedgePointPatch.H"
@@ -46,6 +43,19 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::wedgePointPatch::wedgePointPatch
+(
+    const polyPatch& patch,
+    const pointBoundaryMesh& bm
+)
+:
+    facePointPatch(patch, bm),
+    wedgePolyPatch_(refCast<const wedgePolyPatch>(patch))
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::wedgePointPatch::applyConstraint
@@ -54,7 +64,7 @@ void Foam::wedgePointPatch::applyConstraint
     pointConstraint& pc
 ) const
 {
-    pc.applyConstraint(pointNormals()[pointi]);
+    pc.applyConstraint(wedgePolyPatch_.n());
 }
 
 

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,6 +51,11 @@ class wedgePointPatch
 :
     public facePointPatch
 {
+    // Private data
+
+        //- Local reference cast into the symmetryPlane patch
+        const wedgePolyPatch& wedgePolyPatch_;
+
 
 public:
 
@@ -65,10 +70,7 @@ public:
         (
             const polyPatch& patch,
             const pointBoundaryMesh& bm
-        )
+        );
-        :
-            facePointPatch(patch, bm)
-        {}
 
 
     // Member Functions
@@ -85,6 +87,12 @@ public:
             const label pointi,
             pointConstraint&
         ) const;
+
+        //- Return symmetry plane normal
+        const vector& n() const
+        {
+            return wedgePolyPatch_.n();
+        }
 };
 
 

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/faceMapper/faceMapper.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,9 +28,6 @@ License
 #include "polyMesh.H"
 #include "mapPolyMesh.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::faceMapper::calcAddressing() const
@@ -168,7 +165,7 @@ void Foam::faceMapper::calcAddressing() const
         }
 
 
-        // Grab inserted points (for them the size of addressing is still zero)
+        // Grab inserted faces (for them the size of addressing is still zero)
 
         insertedFaceLabelsPtr_ = new labelList(mesh_.nFaces());
         labelList& insertedFaces = *insertedFaceLabelsPtr_;
@@ -413,13 +410,4 @@ const Foam::labelList& Foam::faceMapper::oldPatchSizes() const
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyMesh.C

@@ -921,7 +921,6 @@ Foam::polyMesh::cellTree() const
 }
 
 
-// Add boundary patches. Constructor helper
 void Foam::polyMesh::addPatches
 (
     const List<polyPatch*>& p,
@@ -968,7 +967,6 @@ void Foam::polyMesh::addPatches
 }
 
 
-// Add mesh zones. Constructor helper
 void Foam::polyMesh::addZones
 (
     const List<pointZone*>& pz,
@@ -1084,7 +1082,6 @@ const Foam::labelList& Foam::polyMesh::faceNeighbour() const
 }
 
 
-// Return old mesh motion points
 const Foam::pointField& Foam::polyMesh::oldPoints() const
 {
     if (oldPointsPtr_.empty())
@@ -1129,11 +1126,14 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
 
     points_ = newPoints;
 
+    bool moveError = false;
     if (debug)
     {
         // Check mesh motion
         if (checkMeshMotion(points_, true))
         {
+            moveError = true;
+
             Info<< "tmp<scalarField> polyMesh::movePoints"
                 << "(const pointField&) : "
                 << "Moving the mesh with given points will "
@@ -1176,6 +1176,12 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
 
     const_cast<Time&>(time()).functionObjects().movePoints(*this);
 
+
+    if (debug && moveError)
+    {
+        write();
+    }
+
     return sweptVols;
 }
 
@@ -1219,7 +1225,6 @@ Foam::label& Foam::polyMesh::comm()
 }
 
 
-// Remove all files and some subdirs (eg, sets)
 void Foam::polyMesh::removeFiles(const fileName& instanceDir) const
 {
     fileName meshFilesPath = thisDb().time().path()/instanceDir/meshDir();

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C

@@ -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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,56 +40,90 @@ namespace Foam
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-void Foam::wedgePolyPatch::initTransforms()
+void Foam::wedgePolyPatch::calcGeometry(PstreamBuffers&)
 {
-    if (size() > 0)
+    if (axis_ != vector::rootMax)
     {
-        const pointField& points = this->points();
+        return;
+    }
 
-        patchNormal_ = operator[](0).normal(points);
+    if (returnReduce(size(), sumOp<label>()))
-        patchNormal_ /= mag(patchNormal_);
+    {
+        const vectorField& nf(faceNormals());
+        n_ = gAverage(nf);
+
+        if (debug)
+        {
+            Info<< "Patch " << name() << " calculated average normal "
+                << n_ << endl;
+        }
+
+
+        // Check the wedge is planar
+        forAll(nf, faceI)
+        {
+            if (magSqr(n_ - nf[faceI]) > SMALL)
+            {
+                // only issue warning instead of error so that the case can
+                // still be read for post-processing
+                WarningIn
+                (
+                    "wedgePolyPatch::calcGeometry(PstreamBuffers&)"
+                )
+                    << "Wedge patch '" << name() << "' is not planar." << nl
+                    << "At local face at "
+                    << primitivePatch::faceCentres()[faceI]
+                    << " the normal " << nf[faceI]
+                    << " differs from the average normal " << n_
+                    << " by " << magSqr(n_ - nf[faceI]) << nl
+                    << "Either correct the patch or split it into planar parts"
+                    << endl;
+            }
+        }
 
         centreNormal_ =
             vector
             (
-                sign(patchNormal_.x())*(max(mag(patchNormal_.x()), 0.5) - 0.5),
+                sign(n_.x())*(max(mag(n_.x()), 0.5) - 0.5),
-                sign(patchNormal_.y())*(max(mag(patchNormal_.y()), 0.5) - 0.5),
+                sign(n_.y())*(max(mag(n_.y()), 0.5) - 0.5),
-                sign(patchNormal_.z())*(max(mag(patchNormal_.z()), 0.5) - 0.5)
+                sign(n_.z())*(max(mag(n_.z()), 0.5) - 0.5)
             );
         centreNormal_ /= mag(centreNormal_);
 
+        cosAngle_ = centreNormal_ & n_;
+
         const scalar cnCmptSum =
             centreNormal_.x() + centreNormal_.y() + centreNormal_.z();
 
         if (mag(cnCmptSum) < (1 - SMALL))
         {
-            FatalErrorIn("wedgePolyPatch::initTransforms()")
+            FatalErrorIn("wedgePolyPatch::calcGeometry(PstreamBuffers&)")
                 << "wedge " << name()
                 << " centre plane does not align with a coordinate plane by "
                 << 1 - mag(cnCmptSum)
                 << exit(FatalError);
         }
 
-        axis_ = centreNormal_ ^ patchNormal_;
+        axis_ = centreNormal_ ^ n_;
         scalar magAxis = mag(axis_);
 
         if (magAxis < SMALL)
         {
-            FatalErrorIn("wedgePolyPatch::initTransforms()")
+            FatalErrorIn("wedgePolyPatch::calcGeometry(PstreamBuffers&)")
                 << "wedge " << name()
                 << " plane aligns with a coordinate plane." << nl
                 << "    The wedge plane should make a small angle (~2.5deg)"
                    " with the coordinate plane" << nl
                 << "    and the the pair of wedge planes should be symmetric"
                 << " about the coordinate plane." << nl
-                << "    Normal of face " << 0 << " is " << patchNormal_
+                << "    Normal of wedge plane is " << n_
                 << " , implied coordinate plane direction is " << centreNormal_
                 << exit(FatalError);
         }
 
         axis_ /= magAxis;
 
-        faceT_ = rotationTensor(centreNormal_, patchNormal_);
+        faceT_ = rotationTensor(centreNormal_, n_);
         cellT_ = faceT_ & faceT_;
     }
 }
@@ -107,10 +141,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
     const word& patchType
 )
 :
-    polyPatch(name, size, start, index, bm, patchType)
+    polyPatch(name, size, start, index, bm, patchType),
-{
+    axis_(vector::rootMax),
-    initTransforms();
+    centreNormal_(vector::rootMax),
-}
+    n_(vector::rootMax),
+    cosAngle_(0.0),
+    faceT_(tensor::zero),
+    cellT_(tensor::zero)
+{}
 
 
 Foam::wedgePolyPatch::wedgePolyPatch
@@ -122,10 +160,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
     const word& patchType
 )
 :
-    polyPatch(name, dict, index, bm, patchType)
+    polyPatch(name, dict, index, bm, patchType),
-{
+    axis_(vector::rootMax),
-    initTransforms();
+    centreNormal_(vector::rootMax),
-}
+    n_(vector::rootMax),
+    cosAngle_(0.0),
+    faceT_(tensor::zero),
+    cellT_(tensor::zero)
+{}
 
 
 Foam::wedgePolyPatch::wedgePolyPatch
@@ -134,10 +176,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
     const polyBoundaryMesh& bm
 )
 :
-    polyPatch(pp, bm)
+    polyPatch(pp, bm),
-{
+    axis_(pp.axis_),
-    initTransforms();
+    centreNormal_(pp.centreNormal_),
-}
+    n_(pp.n_),
+    cosAngle_(pp.cosAngle_),
+    faceT_(pp.faceT_),
+    cellT_(pp.cellT_)
+{}
 
 
 Foam::wedgePolyPatch::wedgePolyPatch
@@ -149,10 +195,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
     const label newStart
 )
 :
-    polyPatch(pp, bm, index, newSize, newStart)
+    polyPatch(pp, bm, index, newSize, newStart),
-{
+    axis_(pp.axis_),
-    initTransforms();
+    centreNormal_(pp.centreNormal_),
-}
+    n_(pp.n_),
+    cosAngle_(pp.cosAngle_),
+    faceT_(pp.faceT_),
+    cellT_(pp.cellT_)
+{}
 
 
 Foam::wedgePolyPatch::wedgePolyPatch
@@ -164,10 +214,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
     const label newStart
 )
 :
-    polyPatch(pp, bm, index, mapAddressing, newStart)
+    polyPatch(pp, bm, index, mapAddressing, newStart),
-{
+    axis_(pp.axis_),
-    initTransforms();
+    centreNormal_(pp.centreNormal_),
-}
+    n_(pp.n_),
+    cosAngle_(pp.cosAngle_),
+    faceT_(pp.faceT_),
+    cellT_(pp.cellT_)
+{}
 
 
 // ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,7 +59,10 @@ class wedgePolyPatch
         vector centreNormal_;
 
         //- Normal to the patch
-        vector patchNormal_;
+        vector n_;
+
+        //- Cosine of the wedge angle
+        scalar cosAngle_;
 
         //- Face transformation tensor
         tensor faceT_;
@@ -67,8 +70,13 @@ class wedgePolyPatch
         //- Neighbour-cell transformation tensor
         tensor cellT_;
 
-        //- Calculate the above tensors
+
-        void initTransforms();
+protected:
+
+    // Protected Member Functions
+
+        //- Calculate the patch geometry
+        virtual void calcGeometry(PstreamBuffers&);
 
 
 public:
@@ -180,9 +188,15 @@ public:
             }
 
             //- Return the normal to the patch
-            const vector& patchNormal() const
+            const vector& n() const
             {
-                return patchNormal_;
+                return n_;
+            }
+
+            //- Return the cosine of the wedge angle
+            scalar cosAngle() const
+            {
+                return cosAngle_;
             }
 
             //- Return face transformation tensor

src/dynamicMesh/layerAdditionRemoval/addCellLayer.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,6 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-//- Calculate the offset to the next layer
 Foam::tmp<Foam::vectorField> Foam::layerAdditionRemoval::extrusionDir() const
 {
     const polyMesh& mesh = topoChanger().mesh();
@@ -79,6 +78,7 @@ Foam::tmp<Foam::vectorField> Foam::layerAdditionRemoval::extrusionDir() const
 
         extrusionDir = minLayerThickness_*masterFaceLayer.pointNormals();
     }
+
     return textrusionDir;
 }
 
@@ -120,7 +120,7 @@ void Foam::layerAdditionRemoval::addCellLayer
 
     // Get the extrusion direction for the added points
 
-    tmp<vectorField> tpointOffsets = extrusionDir();
+    const vectorField pointOffsets(extrusionDir());
 
     // Add the new points
     labelList addedPoints(mp.size());
@@ -135,7 +135,7 @@ void Foam::layerAdditionRemoval::addCellLayer
                 polyAddPoint
                 (
                     points[mp[pointI]]                  // point
-                  + addDelta_*tpointOffsets()[pointI],
+                  + addDelta_*pointOffsets[pointI],
                     mp[pointI],                         // master point
                     -1,                                 // zone for point
                     true                                // supports a cell
@@ -143,14 +143,18 @@ void Foam::layerAdditionRemoval::addCellLayer
             );
     }
 
-    // Pout<< "mp: " << mp << " addedPoints: " << addedPoints << endl;
+    if (debug > 1)
+    {
+        Pout<< "mp: " << mp << " addedPoints: " << addedPoints << endl;
+    }
+
     // Create the cells
 
     const labelList& mc =
         mesh.faceZones()[faceZoneID_.index()].masterCells();
     const labelList& sc =
         mesh.faceZones()[faceZoneID_.index()].slaveCells();
-    // Pout<< "mc: " << mc << " sc: " << sc << endl;
+
     const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
     const boolList& mfFlip = mesh.faceZones()[faceZoneID_.index()].flipMap();
 
@@ -231,10 +235,13 @@ void Foam::layerAdditionRemoval::addCellLayer
             )
         );
 
-        // Pout<< "adding face: " << newFace
+        if (debug > 1)
-        //     << " own: " << mc[faceI]
+        {
-        //     << " nei: " << addedCells[faceI]
+            Pout<< "adding face: " << newFace
-        //     << endl;
+                << " own: " << mc[faceI]
+                << " nei: " << addedCells[faceI]
+                << endl;
+        }
     }
 
     // Modify the faces from the master zone for the new neighbour
@@ -267,10 +274,14 @@ void Foam::layerAdditionRemoval::addCellLayer
                 )
             );
 
-            // Pout<< "Modifying a boundary face. Face: " << curfaceID
+            if (debug > 1)
-            //     << " flip: " << mfFlip[faceI]
+            {
-            //     << endl;
+                Pout<< "Modifying a boundary face. Face: " << curfaceID
+                    << " flip: " << mfFlip[faceI]
+                    << endl;
+            }
         }
+
         // If slave cell is owner, the face remains the same (but with
         // a new neighbour - the newly created cell).  Otherwise, the
         // face is flipped.
@@ -293,10 +304,13 @@ void Foam::layerAdditionRemoval::addCellLayer
                 )
             );
 
-            // Pout<< "modify face, no flip " << curfaceID
+            if (debug > 1)
-            //     << " own: " << own[curfaceID]
+            {
-            //     << " nei: " << addedCells[faceI]
+                Pout<< "modify face, no flip " << curfaceID
-            //     << endl;
+                    << " own: " << own[curfaceID]
+                    << " nei: " << addedCells[faceI]
+                    << endl;
+            }
         }
         else
         {
@@ -317,10 +331,13 @@ void Foam::layerAdditionRemoval::addCellLayer
                 )
             );
 
-            // Pout<< "modify face, with flip " << curfaceID
+            if (debug > 1)
-            //     << " own: " << own[curfaceID]
+            {
-            //     << " nei: " << addedCells[faceI]
+                Pout<< "modify face, with flip " << curfaceID
-            //     << endl;
+                    << " own: " << own[curfaceID]
+                    << " nei: " << addedCells[faceI]
+                    << endl;
+            }
         }
     }
 
@@ -362,10 +379,13 @@ void Foam::layerAdditionRemoval::addCellLayer
             )
         );
 
-        // Pout<< "Add internal face off edge: " << newFace
+        if (debug > 1)
-        //     << " own: " << addedCells[edgeFaces[curEdgeID][0]]
+        {
-        //     << " nei: " << addedCells[edgeFaces[curEdgeID][1]]
+            Pout<< "Add internal face off edge: " << newFace
-        //     << endl;
+                << " own: " << addedCells[edgeFaces[curEdgeID][0]]
+                << " nei: " << addedCells[edgeFaces[curEdgeID][1]]
+                << endl;
+        }
     }
 
     // Prepare creation of faces from boundary edges.
@@ -448,10 +468,13 @@ void Foam::layerAdditionRemoval::addCellLayer
             )
         );
 
-        // Pout<< "add boundary face: " << newFace
+        if (debug > 1)
-        //     << " into patch " << patchID
+        {
-        //     << " own: " << addedCells[edgeFaces[curEdgeID][0]]
+            Pout<< "add boundary face: " << newFace
-        //     << endl;
+                << " into patch " << patchID
+                << " own: " << addedCells[edgeFaces[curEdgeID][0]]
+                << endl;
+        }
     }
 
     // Modify the remaining faces of the master cells to reconnect to the new
@@ -562,12 +585,15 @@ void Foam::layerAdditionRemoval::addCellLayer
                     )
                 );
 
-                // Pout<< "modifying stick-out face. Internal Old face: "
+                if (debug > 1)
-                //     << oldFace
+                {
-                //     << " new face: " << newFace
+                    Pout<< "modifying stick-out face. Internal Old face: "
-                //     << " own: " << own[curFaceID]
+                        << oldFace
-                //     << " nei: " << nei[curFaceID]
+                        << " new face: " << newFace
-                //     << endl;
+                        << " own: " << own[curFaceID]
+                        << " nei: " << nei[curFaceID]
+                        << endl;
+                }
             }
             else
             {
@@ -588,22 +614,24 @@ void Foam::layerAdditionRemoval::addCellLayer
                     )
                 );
 
-                // Pout<< "modifying stick-out face. Boundary Old face: "
+                if (debug > 1)
-                //     << oldFace
+                {
-                //     << " new face: " << newFace
+                    Pout<< "modifying stick-out face. Boundary Old face: "
-                //     << " own: " << own[curFaceID]
+                        << oldFace
-                //     << " patch: "
+                        << " new face: " << newFace
-                //     << mesh.boundaryMesh().whichPatch(curFaceID)
+                        << " own: " << own[curFaceID]
-                //     << endl;
+                        << " patch: "
+                        << mesh.boundaryMesh().whichPatch(curFaceID)
+                        << endl;
+                }
             }
         }
     }
 
     if (debug)
     {
-        Pout<< "void layerAdditionRemoval::addCellLayer("
+        Pout<< "void layerAdditionRemoval::addCellLayer(polyTopoChange&) const "
-            << "polyTopoChange& ref) const "
+            << " for object " << name() << ": "
-            << " for object " << name() << " : "
             << "Finished adding cell layer" << endl;
     }
 }

src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,9 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    Remove a layer of cells and prepare addressing data
-
 \*---------------------------------------------------------------------------*/
 
 #include "layerAdditionRemoval.H"
@@ -87,6 +84,7 @@ bool Foam::layerAdditionRemoval::validCollapse() const
     }
 }
 
+
 void Foam::layerAdditionRemoval::removeCellLayer
 (
     polyTopoChange& ref
@@ -115,7 +113,7 @@ void Foam::layerAdditionRemoval::removeCellLayer
 
     // Remove all the cells from the master layer
     const labelList& mc =
-        topoChanger().mesh().faceZones()[faceZoneID_.index()].masterCells();
+        mesh.faceZones()[faceZoneID_.index()].masterCells();
 
     forAll(mc, faceI)
     {
@@ -205,7 +203,10 @@ void Foam::layerAdditionRemoval::removeCellLayer
 
     labelList ftm = facesToModify.toc();
 
-//Pout<< "faces to modify: " << ftm << endl;
+    if (debug > 1)
+    {
+        Pout<< "faces to modify: " << ftm << endl;
+    }
 
     forAll(ftm, faceI)
     {
@@ -228,9 +229,12 @@ void Foam::layerAdditionRemoval::removeCellLayer
             }
         }
 
-//Pout<< "face label: " << curFaceID
+        if (debug > 1)
-//    << " old face: " << faces[curFaceID]
+        {
-//    << " new face: " << newFace << endl;
+            Pout<< "face label: " << curFaceID
+                << " old face: " << faces[curFaceID]
+                << " new face: " << newFace << endl;
+        }
 
         // Get face zone and its flip
         label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
@@ -238,22 +242,15 @@ void Foam::layerAdditionRemoval::removeCellLayer
 
         if (modifiedFaceZone >= 0)
         {
-            modifiedFaceZoneFlip =
+            const faceZone& fz = mesh.faceZones()[modifiedFaceZone];
-                mesh.faceZones()[modifiedFaceZone].flipMap()
+            modifiedFaceZoneFlip = fz.flipMap()[fz.whichFace(curFaceID)];
-                [
-                    mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
-                ];
         }
 
-        label newNei;
+        label newNeighbour = -1;
 
         if (curFaceID < mesh.nInternalFaces())
         {
-            newNei = nei[curFaceID];
+            newNeighbour = nei[curFaceID];
-        }
-        else
-        {
-            newNei = -1;
         }
 
         // Modify the face
@@ -264,7 +261,7 @@ void Foam::layerAdditionRemoval::removeCellLayer
                 newFace,                // modified face
                 curFaceID,              // label of face being modified
                 own[curFaceID],         // owner
-                newNei,                 // neighbour
+                newNeighbour,           // neighbour
                 false,                  // face flip
                 mesh.boundaryMesh().whichPatch(curFaceID),// patch for face
                 false,                  // remove from zone
@@ -285,20 +282,34 @@ void Foam::layerAdditionRemoval::removeCellLayer
         // of the cell to be removed
         label masterSideCell = own[mf[faceI]];
 
-        if (mesh.isInternalFace(mf[faceI]) && masterSideCell == mc[faceI])
+        if (masterSideCell == mc[faceI])
         {
-            // Owner cell of the face is being removed.
+            if (mesh.isInternalFace(mf[faceI]))
-            // Grab the neighbour instead
+            {
-            masterSideCell = nei[mf[faceI]];
+                // Owner cell of the face is being removed.
+                // Grab the neighbour instead
+                masterSideCell = nei[mf[faceI]];
+            }
+            else
+            {
+                masterSideCell = -1;
+            }
         }
 
         label slaveSideCell = own[ftc[faceI]];
 
-        if (mesh.isInternalFace(ftc[faceI]) && slaveSideCell == mc[faceI])
+        if (slaveSideCell == mc[faceI])
         {
-            // Owner cell of the face is being removed.
+            if (mesh.isInternalFace(ftc[faceI]))
-            // Grab the neighbour instead
+            {
-            slaveSideCell = nei[ftc[faceI]];
+                // Owner cell of the face is being removed.
+                // Grab the neighbour instead
+                slaveSideCell = nei[ftc[faceI]];
+            }
+            else
+            {
+                slaveSideCell = -1;
+            }
         }
 
         // Find out if the face needs to be flipped
@@ -374,15 +385,18 @@ void Foam::layerAdditionRemoval::removeCellLayer
             zoneFlip = !zoneFlip;
         }
 
-// Pout<< "Modifying face " << mf[faceI]
+        if (debug > 1)
-//     << " newFace: " << newFace << nl
+        {
-//     << " newOwner: " << newOwner
+            Pout<< "Modifying face " << mf[faceI]
-//     << " newNeighbour: " << newNeighbour
+                << " newFace: " << newFace << nl
-//     << " flipFace: " << flipFace
+                << " newOwner: " << newOwner
-//     << " newPatchID: " << newPatchID
+                << " newNeighbour: " << newNeighbour
-//     << " newZoneID: " << newZoneID << nl
+                << " flipFace: " << flipFace
-//     << " oldOwn: " << own[mf[faceI]]
+                << " newPatchID: " << newPatchID
-//     << " oldNei: " << nei[mf[faceI]] << endl;
+                << " newZoneID: " << newZoneID << nl
+                << " oldOwn: " << own[mf[faceI]]
+                << " oldNei: " << nei[mf[faceI]] << endl;
+        }
 
         ref.setAction
         (

src/dynamicMesh/layerAdditionRemoval/setLayerPairing.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -88,10 +88,14 @@ bool Foam::layerAdditionRemoval::setLayerPairing() const
 
     facesPairingPtr_ = new labelList(mf.size(), -1);
     labelList& ftc = *facesPairingPtr_;
-    // Pout<< "meshPoints: " << meshPoints << nl
+
-    //      << "localPoints: "
+    if (debug > 1)
-    //     << mesh.faceZones()[faceZoneID_.index()]().localPoints()
+    {
-    //     << endl;
+        Pout<< "meshPoints: " << meshPoints << nl
+            << "localPoints: "
+            << mesh.faceZones()[faceZoneID_.index()]().localPoints()
+            << endl;
+    }
 
     // For all faces, create the mapping
     label nPointErrors = 0;
@@ -119,12 +123,15 @@ bool Foam::layerAdditionRemoval::setLayerPairing() const
             continue;
         }
 
-// Pout<< "curMasterFace: " << faces[mf[faceI]] << nl
+        if (debug > 1)
-//     << "cell shape: " << mesh.cellShapes()[mc[faceI]] << nl
+        {
-//     << "curLocalFace: " << curLocalFace << nl
+            Pout<< "curMasterFace: " << faces[mf[faceI]] << nl
-//     << "lidFace: " << lidFace
+                << "cell shape: " << mesh.cellShapes()[mc[faceI]] << nl
-//     << " master index: " << lidFace.masterIndex()
+                << "curLocalFace: " << curLocalFace << nl
-//     << " oppositeIndex: " << lidFace.oppositeIndex() << endl;
+                << "lidFace: " << lidFace
+                << " master index: " << lidFace.masterIndex()
+                << " oppositeIndex: " << lidFace.oppositeIndex() << endl;
+        }
 
         // Grab the opposite face for face collapse addressing
         ftc[faceI] = lidFace.oppositeIndex();
@@ -145,17 +152,20 @@ bool Foam::layerAdditionRemoval::setLayerPairing() const
                 if (ptc[clp] != lidFace[pointI])
                 {
                     nPointErrors++;
-//                     Pout<< "Topological error in cell layer pairing.  "
+
-//                         << "This mesh is either topologically incorrect "
+                    if (debug > 1)
-//                         << "or the master afce layer is not defined "
+                    {
-//                         << "consistently.  Please check the "
+                        Pout<< "Topological error in cell layer pairing.  "
-//                         << "face zone flip map." << nl
+                            << "This mesh is either topologically incorrect "
-//                         << "First index: " << ptc[clp]
+                            << "or the master face layer is not defined "
-//                         << " new index: " << lidFace[pointI] << endl;
+                            << "consistently.  Please check the "
+                            << "face zone flip map." << nl
+                            << "First index: " << ptc[clp]
+                            << " new index: " << lidFace[pointI] << endl;
+                    }
                 }
             }
         }
-//         Pout<< "ptc: " << ptc << endl;
     }
 
     reduce(nPointErrors, sumOp<label>());

src/dynamicMesh/motionSolver/displacement/displacementMotionSolver.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,76 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * Protected Data Members * * * * * * * * * * * * * //
+
+Foam::IOobject Foam::displacementMotionSolver::points0IO
+(
+    const polyMesh& mesh
+) const
+{
+    const word instance =
+        time().findInstance
+        (
+            mesh.meshDir(),
+            "points0",
+            IOobject::READ_IF_PRESENT
+        );
+
+    if (instance != time().constant())
+    {
+        // points0 written to a time folder
+
+        return
+            IOobject
+            (
+                "points0",
+                instance,
+                polyMesh::meshSubDir,
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            );
+    }
+    else
+    {
+        // check that points0 are actually in constant directory
+
+        IOobject io
+        (
+            "points0",
+            instance,
+            polyMesh::meshSubDir,
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        );
+
+        if (io.headerOk())
+        {
+            return io;
+        }
+        else
+        {
+            // copy of original mesh points
+
+            return
+                IOobject
+                (
+                    "points",
+                    instance,
+                    polyMesh::meshSubDir,
+                    mesh,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE,
+                    false
+                );
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::displacementMotionSolver::displacementMotionSolver
@@ -49,29 +119,14 @@ Foam::displacementMotionSolver::displacementMotionSolver
         IOobject
         (
             "pointDisplacement",
-            mesh.time().timeName(),
+            time().timeName(),
             mesh,
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
         pointMesh::New(mesh)
     ),
-    points0_
+    points0_(pointIOField(points0IO(mesh)))
-    (
-        pointIOField
-        (
-            IOobject
-            (
-                "points",
-                mesh.time().constant(),
-                polyMesh::meshSubDir,
-                mesh,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE,
-                false
-            )
-        )
-    )
 {
     if (points0_.size() != mesh.nPoints())
     {
@@ -81,7 +136,8 @@ Foam::displacementMotionSolver::displacementMotionSolver
             "displacementMotionSolver\n"
             "(\n"
             "    const polyMesh&,\n"
-            "    const IOdictionary&\n"
+            "    const IOdictionary&,\n"
+            "    const word&\n"
             ")"
         )   << "Number of points in mesh " << mesh.nPoints()
             << " differs from number of points " << points0_.size()
@@ -90,7 +146,7 @@ Foam::displacementMotionSolver::displacementMotionSolver
                 IOobject
                 (
                     "points",
-                    mesh.time().constant(),
+                    time().constant(),
                     polyMesh::meshSubDir,
                     mesh,
                     IOobject::MUST_READ,
@@ -118,7 +174,7 @@ void Foam::displacementMotionSolver::movePoints(const pointField&)
 
 void Foam::displacementMotionSolver::updateMesh(const mapPolyMesh& mpm)
 {
-    // pointMesh already updates pointFields.
+    // pointMesh already updates pointFields
 
     motionSolver::updateMesh(mpm);
 
@@ -136,7 +192,7 @@ void Foam::displacementMotionSolver::updateMesh(const mapPolyMesh& mpm)
 
     // Note: boundBox does reduce
     const vector span0 = boundBox(points0_).span();
-    const vector span  = boundBox(points).span();
+    const vector span = boundBox(points).span();
 
     vector scaleFactors(cmptDivide(span0, span));
 
@@ -156,11 +212,11 @@ void Foam::displacementMotionSolver::updateMesh(const mapPolyMesh& mpm)
             }
             else
             {
-                // New point. Assume motion is scaling.
+                // New point - assume motion is scaling
                 newPoints0[pointI] = points0_[oldPointI] + cmptMultiply
                 (
                     scaleFactors,
-                    points[pointI]-points[masterPointI]
+                    points[pointI] - points[masterPointI]
                 );
             }
         }
@@ -170,12 +226,21 @@ void Foam::displacementMotionSolver::updateMesh(const mapPolyMesh& mpm)
             (
                 "displacementMotionSolver::updateMesh"
                 "(const mapPolyMesh&)"
-            )   << "Cannot work out coordinates of introduced vertices."
+            )   << "Cannot determine co-ordinates of introduced vertices."
-                << " New vertex " << pointI << " at coordinate "
+                << " New vertex " << pointI << " at co-ordinate "
                 << points[pointI] << exit(FatalError);
         }
     }
+
+    twoDCorrectPoints(newPoints0);
+
     points0_.transfer(newPoints0);
+
+    // points0 changed - set to write and check-in to database
+    points0_.rename("points0");
+    points0_.writeOpt() = IOobject::AUTO_WRITE;
+    points0_.instance() = time().timeName();
+    points0_.checkIn();
 }
 
 

src/dynamicMesh/motionSolver/displacement/displacementMotionSolver.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,6 +40,7 @@ SourceFiles
 
 #include "motionSolver.H"
 #include "pointFields.H"
+#include "pointIOField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -63,21 +64,22 @@ protected:
         //- Point motion field
         mutable pointVectorField pointDisplacement_;
 
-private:
-
-    // Private data
-
         //- Starting points
-        pointField points0_;
+        pointIOField points0_;
+
+
+    // Protected Member Functions
+
+        //- Return IO object for points0
+        IOobject points0IO(const polyMesh& mesh) const;
+
+private:
 
 
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
-        displacementMotionSolver
+        displacementMotionSolver(const displacementMotionSolver&);
-        (
-            const displacementMotionSolver&
-        );
 
         //- Disallow default bitwise assignment
         void operator=(const displacementMotionSolver&);

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,6 +30,8 @@ License
 #include "PointEdgeWave.H"
 #include "syncTools.H"
 #include "interpolationTable.H"
+#include "mapPolyMesh.H"
+#include "pointConstraints.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -111,11 +113,13 @@ void Foam::displacementLayeredMotionMotionSolver::calcZoneMask
             0
         );
 
-
+        if (debug)
-        Info<< "On cellZone " << cz.name()
+        {
-            << " marked " << returnReduce(nPoints, sumOp<label>())
+            Info<< "On cellZone " << cz.name()
-            << " points and " << returnReduce(nEdges, sumOp<label>())
+                << " marked " << returnReduce(nPoints, sumOp<label>())
-            << " edges." << endl;
+                << " points and " << returnReduce(nEdges, sumOp<label>())
+                << " edges." << endl;
+        }
     }
 }
 
@@ -239,14 +243,21 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
         {
             FatalIOErrorIn
             (
-                "displacementLayeredMotionMotionSolver::faceZoneEvaluate(..)",
+                "displacementLayeredMotionMotionSolver::faceZoneEvaluate"
+                "("
+                    "const faceZone&, "
+                    "const labelList&, "
+                    "const dictionary&, "
+                    "const PtrList<pointVectorField>&, "
+                    "const label"
+                ") const",
                 *this
-            )   << "slip can only be used on second faceZonePatch of pair."
+            )   << "slip can only be used on second faceZone patch of pair.  "
                 << "FaceZone:" << fz.name()
                 << exit(FatalIOError);
         }
         // Use field set by previous bc
-        fld = vectorField(patchDisp[patchI-1], meshPoints);
+        fld = vectorField(patchDisp[patchI - 1], meshPoints);
     }
     else if (type == "follow")
     {
@@ -269,7 +280,14 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
     {
         FatalIOErrorIn
         (
-            "displacementLayeredMotionMotionSolver::faceZoneEvaluate(..)",
+            "displacementLayeredMotionMotionSolver::faceZoneEvaluate"
+            "("
+                "const faceZone&, "
+                "const labelList&, "
+                "const dictionary&, "
+                "const PtrList<pointVectorField>&, "
+                "const label"
+            ") const",
             *this
         )   << "Unknown faceZonePatch type " << type << " for faceZone "
             << fz.name() << exit(FatalIOError);
@@ -295,9 +313,9 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
         FatalIOErrorIn
         (
             "displacementLayeredMotionMotionSolver::"
-            "correctBoundaryConditions(..)",
+            "cellZoneSolve(const label, const dictionary&)",
             *this
-        )   << "Can only handle 2 faceZones (= patches) per cellZone. "
+        )   << "Two faceZones (patches) must be specifed per cellZone. "
             << " cellZone:" << cellZoneI
             << " patches:" << patchesDict.toc()
             << exit(FatalIOError);
@@ -317,7 +335,7 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
             FatalIOErrorIn
             (
                 "displacementLayeredMotionMotionSolver::"
-                "correctBoundaryConditions(..)",
+                "cellZoneSolve(const label, const dictionary&)",
                 *this
             )   << "Cannot find faceZone " << faceZoneName
                 << endl << "Valid zones are " << mesh().faceZones().names()
@@ -386,14 +404,17 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
             patchI
         );
 
-
+        if (debug)
-Info<< "For cellZone:" << cellZoneI
+        {
-    << " for faceZone:" << fz.name() << " nPoints:" << tseed().size()
+            Info<< "For cellZone:" << cellZoneI
-    << " have patchField:"
+                << " for faceZone:" << fz.name()
-    << " max:" << gMax(tseed())
+                << " nPoints:" << tseed().size()
-    << " min:" << gMin(tseed())
+                << " have patchField:"
-    << " avg:" << gAverage(tseed())
+                << " max:" << gMax(tseed())
-    << endl;
+                << " min:" << gMin(tseed())
+                << " avg:" << gAverage(tseed())
+                << endl;
+        }
 
         // Set distance and transported value
         walkStructured
@@ -417,16 +438,15 @@ Info<< "For cellZone:" << cellZoneI
 
     // Solve
     // ~~~~~
-    // solving the interior is just interpolating
 
     if (debug)
     {
-        // Get normalised distance
+        // Normalised distance
         pointScalarField distance
         (
             IOobject
             (
-                "distance",
+                mesh().cellZones()[cellZoneI].name() + ":distance",
                 mesh().time().timeName(),
                 mesh(),
                 IOobject::NO_READ,
@@ -434,40 +454,66 @@ Info<< "For cellZone:" << cellZoneI
                 false
             ),
             pointMesh::New(mesh()),
-            dimensionedScalar("distance", dimLength, 0.0)
+            dimensionedScalar("zero", dimLength, 0.0)
         );
+
         forAll(distance, pointI)
+        {
+            scalar d1 = patchDist[0][pointI];
+            scalar d2 = patchDist[1][pointI];
+            if (d1 + d2 > SMALL)
+            {
+                scalar s = d1/(d1 + d2);
+                distance[pointI] = s;
+            }
+        }
+
+        Info<< "Writing " << pointScalarField::typeName
+            << distance.name() << " to "
+            << mesh().time().timeName() << endl;
+        distance.write();
+    }
+
+
+    const word interpolationScheme = zoneDict.lookup("interpolationScheme");
+
+    if (interpolationScheme == "oneSided")
+    {
+        forAll(pointDisplacement_, pointI)
+        {
+            if (isZonePoint[pointI])
+            {
+                pointDisplacement_[pointI] = patchDisp[0][pointI];
+            }
+        }
+    }
+    else if (interpolationScheme == "linear")
+    {
+        forAll(pointDisplacement_, pointI)
         {
             if (isZonePoint[pointI])
             {
                 scalar d1 = patchDist[0][pointI];
                 scalar d2 = patchDist[1][pointI];
-                if (d1+d2 > SMALL)
+                scalar s = d1/(d1 + d2 + VSMALL);
-                {
+
-                    scalar s = d1/(d1+d2);
+                const vector& pd1 = patchDisp[0][pointI];
-                    distance[pointI] = s;
+                const vector& pd2 = patchDisp[1][pointI];
-                }
+
+                pointDisplacement_[pointI] = (1 - s)*pd1 + s*pd2;
             }
         }
-        Info<< "Writing distance pointScalarField to "
-            << mesh().time().timeName() << endl;
-        distance.write();
     }
-
+    else
-    // Average
-    forAll(pointDisplacement_, pointI)
     {
-        if (isZonePoint[pointI])
+        FatalErrorIn
-        {
+        (
-            scalar d1 = patchDist[0][pointI];
+            "displacementLayeredMotionMotionSolver::"
-            scalar d2 = patchDist[1][pointI];
+            "cellZoneSolve(const label, const dictionary&)"
-
+        )
-            scalar s = d1/(d1+d2+VSMALL);
+            << "Invalid interpolationScheme: " << interpolationScheme
-
+            << ". Valid schemes are 'oneSided' and 'linear'"
-            pointDisplacement_[pointI] =
+            << exit(FatalError);
-                (1-s)*patchDisp[0][pointI]
-              + s*patchDisp[1][pointI];
-        }
     }
 }
 
@@ -502,23 +548,19 @@ Foam::displacementLayeredMotionMotionSolver::curPoints() const
         points0() + pointDisplacement_.internalField()
     );
 
-    twoDCorrectPoints(tcurPoints());
-
     return tcurPoints;
 }
 
 
 void Foam::displacementLayeredMotionMotionSolver::solve()
 {
-    // The points have moved so before interpolation update
+    // The points have moved so before interpolation update the motionSolver
-    // the motionSolver accordingly
     movePoints(mesh().points());
 
     // Apply boundary conditions
     pointDisplacement_.boundaryField().updateCoeffs();
 
-    // Apply all regions (=cellZones)
+    // Solve motion on all regions (=cellZones)
-
     const dictionary& regionDicts = coeffDict().subDict("regions");
     forAllConstIter(dictionary, regionDicts, regionIter)
     {
@@ -527,14 +569,13 @@ void Foam::displacementLayeredMotionMotionSolver::solve()
 
         label zoneI = mesh().cellZones().findZoneID(cellZoneName);
 
-        Info<< "solve : zone:" << cellZoneName << " index:" << zoneI
+        Info<< "solving for zone: " << cellZoneName << endl;
-            << endl;
 
         if (zoneI == -1)
         {
             FatalIOErrorIn
             (
-                "displacementLayeredMotionMotionSolver::solve(..)",
+                "displacementLayeredMotionMotionSolver::solve()",
                 *this
             )   << "Cannot find cellZone " << cellZoneName
                 << endl << "Valid zones are " << mesh().cellZones().names()
@@ -545,7 +586,9 @@ void Foam::displacementLayeredMotionMotionSolver::solve()
     }
 
     // Update pointDisplacement for solved values
-    pointDisplacement_.correctBoundaryConditions();
+    const pointConstraints& pcs =
+        pointConstraints::New(pointDisplacement_.mesh());
+    pcs.constrainDisplacement(pointDisplacement_, false);
 }
 
 
@@ -555,6 +598,27 @@ void Foam::displacementLayeredMotionMotionSolver::updateMesh
 )
 {
     displacementMotionSolver::updateMesh(mpm);
+
+    const vectorField displacement(this->newPoints() - points0_);
+
+    forAll(points0_, pointI)
+    {
+        label oldPointI = mpm.pointMap()[pointI];
+
+        if (oldPointI >= 0)
+        {
+            label masterPointI = mpm.reversePointMap()[oldPointI];
+
+            if ((masterPointI != pointI))
+            {
+                // newly inserted point in this cellZone
+
+                // need to set point0 so that it represents the position that
+                // it would have had if it had existed for all time
+                points0_[pointI] -= displacement[pointI];
+            }
+        }
+    }
 }
 
 

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -149,7 +149,7 @@ public:
                     (-rAB + rA - rB)
                    *(-rAB - rA + rB)
                    *(-rAB + rA + rB)
-                   *( rAB + rA + rA)
+                   *( rAB + rA + rB)
                 );
         }
 

src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,7 +29,7 @@ License
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
 template<class Type>
-Foam::autoPtr<labelList> Foam::AveragingMethods::Dual<Type>::size
+Foam::autoPtr<Foam::labelList> Foam::AveragingMethods::Dual<Type>::size
 (
     const fvMesh& mesh
 )

src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/makeAveragingMethods.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,7 +23,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "fvCFD.H"
+#include "Field.H"
+#include "fvcGrad.H"
 #include "polyMeshTetDecomposition.H"
 
 #include "Basic.H"

src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/CorrectionLimitingMethod/CorrectionLimitingMethod.H

@@ -35,7 +35,7 @@ SourceFiles
 #ifndef CorrectionLimitingMethod_H
 #define CorrectionLimitingMethod_H
 
-#include "fvCFD.H"
+#include "volFieldsFwd.H"
 #include "dictionary.H"
 #include "runTimeSelectionTables.H"
 

src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,10 @@ License
 
 #include "Implicit.H"
 #include "fixedValueFvsPatchField.H"
+#include "fvmDdt.H"
+#include "fvmDiv.H"
+#include "fvmLaplacian.H"
+#include "fvcReconstruct.H"
 #include "volPointInterpolation.H"                                              
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/ParticleStressModel/ParticleStressModel.H

@@ -35,7 +35,8 @@ SourceFiles
 #ifndef ParticleStressModel_H
 #define ParticleStressModel_H
 
-#include "fvCFD.H"
+#include "Field.H"
+#include "FieldField.H"
 #include "dictionary.H"
 #include "runTimeSelectionTables.H"
 

src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/TimeScaleModel/TimeScaleModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,9 @@ SourceFiles
 #ifndef TimeScaleModel_H
 #define TimeScaleModel_H
 
-#include "fvCFD.H"
+#include "constants.H"
+#include "Field.H"
+#include "FieldField.H"
 #include "dictionary.H"
 #include "runTimeSelectionTables.H"
 

src/mesh/blockMesh/Make/files

@@ -7,6 +7,7 @@ curvedEdges/curvedEdge.C
 curvedEdges/lineEdge.C
 curvedEdges/polyLineEdge.C
 curvedEdges/lineDivide.C
+curvedEdges/BSplineEdge.C
 curvedEdges/splineEdge.C
 
 blockDescriptor/blockDescriptor.C

src/mesh/blockMesh/curvedEdges/BSplineEdge.C

@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 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 "BSplineEdge.H"
+#include "addToRunTimeSelectionTable.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(BSplineEdge, 0);
+
+    addToRunTimeSelectionTable
+    (
+        curvedEdge,
+        BSplineEdge,
+        Istream
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::BSplineEdge::BSplineEdge
+(
+    const pointField& points,
+    const label start,
+    const label end,
+    const pointField& internalPoints
+)
+:
+    curvedEdge(points, start, end),
+    BSpline(appendEndPoints(points, start, end, internalPoints))
+{}
+
+
+Foam::BSplineEdge::BSplineEdge(const pointField& points, Istream& is)
+:
+    curvedEdge(points, is),
+    BSpline(appendEndPoints(points, start_, end_, pointField(is)))
+{
+    token t(is);
+    is.putBack(t);
+
+    // discard unused start/end tangents
+    if (t == token::BEGIN_LIST)
+    {
+        vector tangent0Ignored(is);
+        vector tangent1Ignored(is);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::BSplineEdge::~BSplineEdge()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::point Foam::BSplineEdge::position(const scalar mu) const
+{
+    return BSpline::position(mu);
+}
+
+
+Foam::scalar Foam::BSplineEdge::length() const
+{
+    return BSpline::length();
+}
+
+
+// ************************************************************************* //

src/mesh/blockMesh/curvedEdges/BSplineEdge.H

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 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::BSplineEdge
+
+Description
+    A curvedEdge interface for B-splines.
+
+SourceFiles
+    BSplineEdge.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BSplineEdge_H
+#define BSplineEdge_H
+
+#include "curvedEdge.H"
+#include "BSpline.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class BSplineEdge Declaration
+\*---------------------------------------------------------------------------*/
+
+class BSplineEdge
+:
+    public curvedEdge,
+    public BSpline
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        BSplineEdge(const BSplineEdge&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const BSplineEdge&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("BSpline");
+
+
+    // Constructors
+
+        //- Construct from components
+        BSplineEdge
+        (
+            const pointField&,
+            const label start,
+            const label end,
+            const pointField& internalPoints
+        );
+
+        //- Construct from Istream, setting pointsList
+        BSplineEdge(const pointField&, Istream&);
+
+
+    //- Destructor
+    virtual ~BSplineEdge();
+
+
+    // Member Functions
+
+        //- Return the point position corresponding to the curve parameter
+        //  0 <= lambda <= 1
+        virtual point position(const scalar) const;
+
+        //- Return the length of the spline curve (not implemented)
+        virtual scalar length() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/twoDPointCorrector/twoDPointCorrector.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,26 +28,20 @@ License
 #include "wedgePolyPatch.H"
 #include "emptyPolyPatch.H"
 #include "SubField.H"
-
+#include "meshTools.H"
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-const scalar twoDPointCorrector::edgeOrthogonalityTol = 1.0 - 1e-4;
+const Foam::scalar Foam::twoDPointCorrector::edgeOrthogonalityTol = 1.0 - 1e-4;
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void twoDPointCorrector::calcAddressing() const
+void Foam::twoDPointCorrector::calcAddressing() const
 {
     // Find geometry normal
     planeNormalPtr_ = new vector(0, 0, 0);
     vector& pn = *planeNormalPtr_;
 
-    bool isWedge = false;
-
     // Algorithm:
     // Attempt to find wedge patch and work out the normal from it.
     // If not found, find an empty patch with faces in it and use the
@@ -61,9 +55,15 @@ void twoDPointCorrector::calcAddressing() const
     {
         if (isA<wedgePolyPatch>(patches[patchI]))
         {
-            isWedge = true;
+            isWedge_ = true;
 
-            pn = refCast<const wedgePolyPatch>(patches[patchI]).centreNormal();
+            const wedgePolyPatch& wp =
+                refCast<const wedgePolyPatch>(patches[patchI]);
+
+            pn = wp.centreNormal();
+
+            wedgeAxis_ = wp.axis();
+            wedgeAngle_ = mag(acos(wp.cosAngle()));
 
             if (polyMesh::debug)
             {
@@ -75,7 +75,7 @@ void twoDPointCorrector::calcAddressing() const
     }
 
     // Try to find an empty patch with faces
-    if (!isWedge)
+    if (!isWedge_)
     {
         forAll(patches, patchI)
         {
@@ -96,11 +96,8 @@ void twoDPointCorrector::calcAddressing() const
 
     if (mag(pn) < VSMALL)
     {
-        FatalErrorIn
+        FatalErrorIn("twoDPointCorrector::calcAddressing()")
-        (
+            << "Cannot determine normal vector from patches."
-            "twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh, "
-            "const vector& n)"
-        )   << "Cannot determine normal vector from patches."
             << abort(FatalError);
     }
     else
@@ -125,9 +122,9 @@ void twoDPointCorrector::calcAddressing() const
 
     forAll(meshEdges, edgeI)
     {
-        vector edgeVector =
+        const edge& e = meshEdges[edgeI];
-            meshEdges[edgeI].vec(meshPoints)/
+
-            (meshEdges[edgeI].mag(meshPoints) + VSMALL);
+        vector edgeVector = e.vec(meshPoints)/(e.mag(meshPoints) + VSMALL);
 
         if (mag(edgeVector & pn) > edgeOrthogonalityTol)
         {
@@ -141,15 +138,12 @@ void twoDPointCorrector::calcAddressing() const
     // Construction check: number of points in a read 2-D or wedge geometry
     // should be odd and the number of edges normal to the plane should be
     // exactly half the number of points
-    if (!isWedge)
+    if (!isWedge_)
     {
         if (meshPoints.size() % 2 != 0)
         {
-            WarningIn
+            WarningIn("twoDPointCorrector::calcAddressing()")
-            (
+                << "the number of vertices in the geometry "
-                "twoDPointCorrector::twoDPointCorrector("
-                "const polyMesh& mesh, const vector& n)"
-            )   << "the number of vertices in the geometry "
                 << "is odd - this should not be the case for a 2-D case. "
                 << "Please check the geometry."
                 << endl;
@@ -157,11 +151,8 @@ void twoDPointCorrector::calcAddressing() const
 
         if (2*nNormalEdges != meshPoints.size())
         {
-            WarningIn
+            WarningIn("twoDPointCorrector::calcAddressing()")
-            (
+                << "The number of points in the mesh is "
-                "twoDPointCorrector::twoDPointCorrector("
-                "const polyMesh& mesh, const vector& n)"
-            )   << "The number of points in the mesh is "
                 << "not equal to twice the number of edges normal to the plane "
                 << "- this may be OK only for wedge geometries.\n"
                 << "    Please check the geometry or adjust "
@@ -174,21 +165,38 @@ void twoDPointCorrector::calcAddressing() const
 }
 
 
-void twoDPointCorrector::clearAddressing() const
+void Foam::twoDPointCorrector::clearAddressing() const
 {
     deleteDemandDrivenData(planeNormalPtr_);
     deleteDemandDrivenData(normalEdgeIndicesPtr_);
 }
 
 
+void Foam::twoDPointCorrector::snapToWedge
+(
+    const vector& n,
+    const point& A,
+    point& p
+) const
+{
+    scalar ADash = mag(A - wedgeAxis_*(wedgeAxis_ & A));
+    vector pDash = ADash*tan(wedgeAngle_)*planeNormal();
+
+    p = A + sign(n & p)*pDash;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh)
+Foam::twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh)
 :
     MeshObject<polyMesh, Foam::UpdateableMeshObject, twoDPointCorrector>(mesh),
     required_(mesh_.nGeometricD() == 2),
     planeNormalPtr_(NULL),
-    normalEdgeIndicesPtr_(NULL)
+    normalEdgeIndicesPtr_(NULL),
+    isWedge_(false),
+    wedgeAxis_(vector::zero),
+    wedgeAngle_(0.0)
 {}
 
 
@@ -203,7 +211,7 @@ Foam::twoDPointCorrector::~twoDPointCorrector()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-direction twoDPointCorrector::normalDir() const
+Foam::direction Foam::twoDPointCorrector::normalDir() const
 {
     const vector& pn = planeNormal();
 
@@ -231,8 +239,7 @@ direction twoDPointCorrector::normalDir() const
 }
 
 
-// Return plane normal
+const Foam::vector& Foam::twoDPointCorrector::planeNormal() const
-const vector& twoDPointCorrector::planeNormal() const
 {
     if (!planeNormalPtr_)
     {
@@ -243,8 +250,7 @@ const vector& twoDPointCorrector::planeNormal() const
 }
 
 
-// Return indices of normal edges.
+const Foam::labelList& Foam::twoDPointCorrector::normalEdgeIndices() const
-const labelList& twoDPointCorrector::normalEdgeIndices() const
 {
     if (!normalEdgeIndicesPtr_)
     {
@@ -255,7 +261,7 @@ const labelList& twoDPointCorrector::normalEdgeIndices() const
 }
 
 
-void twoDPointCorrector::correctPoints(pointField& p) const
+void Foam::twoDPointCorrector::correctPoints(pointField& p) const
 {
     if (!required_) return;
 
@@ -277,16 +283,25 @@ void twoDPointCorrector::correctPoints(pointField& p) const
         point& pEnd = p[meshEdges[neIndices[edgeI]].end()];
 
         // calculate average point position
-        const point A = 0.5*(pStart + pEnd);
+        point A = 0.5*(pStart + pEnd);
+        meshTools::constrainToMeshCentre(mesh_, A);
 
-        // correct point locations
+        if (isWedge_)
-        pStart = A + pn*(pn & (pStart - A));
+        {
-        pEnd = A + pn*(pn & (pEnd - A));
+            snapToWedge(pn, A, pStart);
+            snapToWedge(pn, A, pEnd);
+        }
+        else
+        {
+            // correct point locations
+            pStart = A + pn*(pn & (pStart - A));
+            pEnd = A + pn*(pn & (pEnd - A));
+        }
     }
 }
 
 
-void twoDPointCorrector::correctDisplacement
+void Foam::twoDPointCorrector::correctDisplacement
 (
     const pointField& p,
     vectorField& disp
@@ -316,29 +331,36 @@ void twoDPointCorrector::correctDisplacement
         point pEnd = p[endPointI] + disp[endPointI];
 
         // calculate average point position
-        const point A = 0.5*(pStart + pEnd);
+        point A = 0.5*(pStart + pEnd);
+        meshTools::constrainToMeshCentre(mesh_, A);
 
-        // correct point locations
+        if (isWedge_)
-        disp[startPointI] = (A + pn*(pn & (pStart - A))) - p[startPointI];
+        {
-        disp[endPointI] = (A + pn*(pn & (pEnd - A))) - p[endPointI];
+            snapToWedge(pn, A, pStart);
+            snapToWedge(pn, A, pEnd);
+            disp[startPointI] = pStart - p[startPointI];
+            disp[endPointI] = pEnd - p[endPointI];
+        }
+        else
+        {
+            // correct point locations
+            disp[startPointI] = (A + pn*(pn & (pStart - A))) - p[startPointI];
+            disp[endPointI] = (A + pn*(pn & (pEnd - A))) - p[endPointI];
+        }
     }
 }
 
 
-void twoDPointCorrector::updateMesh(const mapPolyMesh&)
+void Foam::twoDPointCorrector::updateMesh(const mapPolyMesh&)
 {
     clearAddressing();
 }
 
 
-bool twoDPointCorrector::movePoints()
+bool Foam::twoDPointCorrector::movePoints()
 {
     return true;
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/meshTools/twoDPointCorrector/twoDPointCorrector.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,7 @@ namespace Foam
 class polyMesh;
 
 /*---------------------------------------------------------------------------*\
-                        Class twoDPointCorrector Declaration
+                     Class twoDPointCorrector Declaration
 \*---------------------------------------------------------------------------*/
 
 class twoDPointCorrector
@@ -74,6 +74,15 @@ class twoDPointCorrector
         //- Indices of edges normal to plane
         mutable labelList* normalEdgeIndicesPtr_;
 
+        //- Flag to indicate a wedge geometry
+        mutable bool isWedge_;
+
+        //- Wedge axis (if wedge geometry)
+        mutable vector wedgeAxis_;
+
+        //- Wedge angle (if wedge geometry)
+        mutable scalar wedgeAngle_;
+
 
     // Private Member Functions
 
@@ -90,6 +99,9 @@ class twoDPointCorrector
         //- Clear addressing
         void clearAddressing() const;
 
+        //- Snap a point to the wedge patch(es)
+        void snapToWedge(const vector& n, const point& A, point& p) const;
+
 
     // Static data members
 

src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C

@@ -30,8 +30,8 @@ License
 
 bool Foam::sixDoFRigidBodyMotion::read(const dictionary& dict)
 {
-    dict.lookup("momentOfInertia") >> momentOfInertia_;
     dict.lookup("mass") >> mass_;
+    dict.lookup("momentOfInertia") >> momentOfInertia_;
     aRelax_ = dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0);
     aDamp_ = dict.lookupOrDefault<scalar>("accelerationDamping", 1.0);
     report_ = dict.lookupOrDefault<Switch>("report", false);
@@ -50,16 +50,18 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
 {
     motionState_.write(os);
 
-    os.writeKeyword("initialCentreOfRotation")
+    os.writeKeyword("centreOfMass")
-        << initialCentreOfRotation_ << token::END_STATEMENT << nl;
+        << initialCentreOfMass_ << token::END_STATEMENT << nl;
-    os.writeKeyword("initialOrientation")
+    os.writeKeyword("orientation")
         << initialQ_ << token::END_STATEMENT << nl;
-    os.writeKeyword("momentOfInertia")
-        << momentOfInertia_ << token::END_STATEMENT << nl;
     os.writeKeyword("mass")
         << mass_ << token::END_STATEMENT << nl;
+    os.writeKeyword("momentOfInertia")
+        << momentOfInertia_ << token::END_STATEMENT << nl;
     os.writeKeyword("accelerationRelaxation")
         << aRelax_ << token::END_STATEMENT << nl;
+    os.writeKeyword("accelerationDamping")
+        << aDamp_ << token::END_STATEMENT << nl;
     os.writeKeyword("report")
         << report_ << token::END_STATEMENT << nl;
 

src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C

@@ -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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -113,7 +113,7 @@ void Foam::Reaction<ReactionThermo>::setThermo
 {
     if (rhs_.size() > 0)
     {
-        ReactionThermo::operator=
+        ReactionThermo::thermoType::operator=
         (
             rhs_[0].stoichCoeff*(*thermoDatabase[species_[rhs_[0].index]])
         );
@@ -149,7 +149,7 @@ Foam::Reaction<ReactionThermo>::Reaction
     const HashPtrTable<ReactionThermo>& thermoDatabase
 )
 :
-    ReactionThermo(*thermoDatabase[species[0]]),
+    ReactionThermo::thermoType(*thermoDatabase[species[0]]),
     name_("un-named-reaction-" + Foam::name(getNewReactionID())),
     species_(species),
     lhs_(lhs),
@@ -166,7 +166,7 @@ Foam::Reaction<ReactionThermo>::Reaction
     const speciesTable& species
 )
 :
-    ReactionThermo(r),
+    ReactionThermo::thermoType(r),
     name_(r.name() + "Copy"),
     species_(species),
     lhs_(r.lhs_),
@@ -320,7 +320,7 @@ Foam::Reaction<ReactionThermo>::Reaction
     Istream& is
 )
 :
-    ReactionThermo(*thermoDatabase[species[0]]),
+    ReactionThermo::thermoType(*thermoDatabase[species[0]]),
     name_("un-named-reaction" + Foam::name(getNewReactionID())),
     species_(species)
 {
@@ -337,7 +337,7 @@ Foam::Reaction<ReactionThermo>::Reaction
     const dictionary& dict
 )
 :
-    ReactionThermo(*thermoDatabase[species[0]]),
+    ReactionThermo::thermoType(*thermoDatabase[species[0]]),
     name_(dict.dictName()),
     species_(species)
 {

src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -64,7 +64,7 @@ inline Ostream& operator<<(Ostream&, const Reaction<ReactionThermo>&);
 template<class ReactionThermo>
 class Reaction
 :
-    public ReactionThermo
+    public ReactionThermo::thermoType
 {
 protected:
 

src/thermophysicalModels/specie/thermo/thermo/thermo.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -122,6 +122,10 @@ class thermo
 
 public:
 
+    //- The thermodynamics of the individual species'
+    typedef thermo<Thermo, Type> thermoType;
+
+
     // Constructors
 
         //- construct from components