Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/utilities/mesh/manipulation/createBaffles/createBafflesDict

@@ -19,13 +19,18 @@ FoamFile
 // - or converting boundary faces into a boundary face
 //   (internalFacesOnly=false)(though should use really createPatch
 //    to do this)
-// - can also create duplicate (overlapping) sets of baffles:
+// - specification in one of two modes:
+//      - patchPairs : create two patches of same type, same input
+//      - patches    : create patches separately, full control over what
+//                     to create on what side
+//                     (this mode can also create duplicate (overlapping)
+//                      sets of baffles:
 //                      - internalFacesOnly = false
 //                      - have 4 entries in patches:
 //                          - master
 //                          - slave
 //                          - additional master
-//          - additional slave
+//                          - additional slave)
 
 
 // Whether to convert internal faces only (so leave boundary faces intact).

doc/Doxygen/css/doxyTabs.css

@@ -12,7 +12,8 @@ div.navigation
 
 
 div.tabs,
-div.tabs2
+div.tabs2,
+div.tabs3
 {
     width           : 100%;
     padding-top     : 5px;
@@ -22,7 +23,8 @@ div.tabs2
 }
 
 div.tabs ul,
-div.tabs2 ul
+div.tabs2 ul,
+div.tabs3 ul
 {
     margin: 0px;
     padding-left: 10px;
@@ -33,7 +35,8 @@ div.tabs2 ul
 
 
 div.tabs li, div.tabs form,
-div.tabs2 li, div.tabs2 form
+div.tabs2 li, div.tabs2 form,
+div.tabs3 li, div.tabs3 form
 {
     display         : inline;
     margin          : 0px;
@@ -41,7 +44,8 @@ div.tabs2 li, div.tabs2 form
 }
 
 div.tabs ul li,
-div.tabs2 ul li
+div.tabs2 ul li,
+div.tabs3 ul li
 {
     padding-top: 10px;
     padding-bottom: 10px;
@@ -77,14 +81,19 @@ div.tabs a:hover,
 div.tabs2 a:link,
 div.tabs2 a:visited,
 div.tabs2 a:active,
-div.tabs2 a:hover
+div.tabs2 a:hover,
+div.tabs3 a:link,
+div.tabs3 a:visited,
+div.tabs3 a:active,
+div.tabs3 a:hover
 {
     color           : #000000;
     color:#555;
 }
 
 div.tabs span,
-div.tabs2 span
+div.tabs2 span,
+div.tabs3 span
 {
     display         : inline;
     padding         : 0px 9px;

doc/Doxygen/css/doxygen.css

@@ -122,11 +122,13 @@ a.elRef {
 }
 
 a.code {
-	color: #3030f0;
+        color: #3030f0;;
+               font-family: "Courier New", Courier, monospace, fixed;
 }
 
 a.codeRef {
-	color: #3030f0;
+        color: #3030f0;;
+               font-family: "Courier New", Courier, monospace, fixed;
 }
 
 /* @end */
@@ -136,8 +138,8 @@ dl.el {
 }
 
 .fragment {
-	font-family: monospace, fixed;
         font-size: 105%;
+        font-family: "Courier New", Courier, monospace, fixed;
 }
 
 pre.fragment {
@@ -228,47 +230,58 @@ img.footer {
 /* @group Code Colorization */
 
 span.keyword {
-	color: #008000
+        color: #008000;
+               font-family: "Courier New", Courier, monospace, fixed;
 }
 
 span.keywordtype {
-	color: #604020
+        color: #604020;
+        font-family: "Courier New", Courier, monospace, fixed;
 }
 
 span.keywordflow {
-	color: #e08000
+        color: #e08000;
+        font-family: "Courier New", Courier, monospace, fixed;
 }
 
 span.comment {
-	color: #800000
+        color: #800000;
+        font-family: "Courier New", Courier, monospace, fixed;
 }
 
 span.preprocessor {
-	color: #806020
+        color: #806020;
+               font-family: "Courier New", Courier, monospace, fixed;
 }
 
 span.stringliteral {
-	color: #002080
+        color: #002080;
+               font-family: "Courier New", Courier, monospace, fixed;
 }
 
 span.charliteral {
-	color: #008080
+        color: #008080;
+               font-family: "Courier New", Courier, monospace, fixed;
 }
 
 span.vhdldigit {
-	color: #ff00ff
+        color: #ff00ff;
+               font-family: "Courier New", Courier, monospace, fixed;
 }
 
 span.vhdlchar {
-	color: #000000
+        color: #000000;
+               font-family: "Courier New", Courier, monospace, fixed;
 }
 
 span.vhdlkeyword {
-	color: #700070
+        color: #700070;
+               font-family: "Courier New", Courier, monospace, fixed;
 }
 
 span.vhdllogic {
-	color: #ff0000
+        color: #ff0000;
+               font-family: "Courier New", Courier, monospace, fixed;
 }
 
 /* @end */

src/OpenFOAM/db/functionObjects/outputFilterOutputControl/outputFilterOutputControl.H

@@ -82,7 +82,7 @@ private:
         //  a value <= 1 means execute at every time step
         label outputInterval_;
 
-        //- Dumping counter for ocOutputTime
+        //- Dumping counter for ocOutputTime or index dump for ocAdjustableTime
         label outputTimeLastDump_;
 
         //- Dump each deltaT (adjust Ttime)
@@ -123,19 +123,19 @@ public:
         bool output();
 
         //- Return outputControl
-        outputControls outputControl()
+        outputControls outputControl() const
         {
             return outputControl_;
         }
 
         //- Return writeInterval
-        scalar writeInterval()
+        scalar writeInterval() const
         {
             return writeInterval_;
         }
 
         //- Return outputTimeLastDump
-        label outputTimeLastDump()
+        label outputTimeLastDump() const
         {
             return outputTimeLastDump_;
         }

src/OpenFOAM/meshes/polyMesh/polyMesh.C

@@ -1351,8 +1351,10 @@ bool Foam::polyMesh::pointInCell
 
         case FACECENTRETETS:
         {
-            const point& cc = cellCentres()[cellI];
+            // only test that point is on inside of plane defined by cell face
+            // triangles
             const cell& cFaces = cells()[cellI];
+
             forAll(cFaces, cFaceI)
             {
                 label faceI = cFaces[cFaceI];
@@ -1376,31 +1378,62 @@ bool Foam::polyMesh::pointInCell
                         nextPointI = f[fp];
                     }
 
-                    if
+                    triPointRef faceTri
                     (
-                        tetPointRef
-                        (
-                            points()[nextPointI],
                         points()[pointI],
-                            fc,
+                        points()[nextPointI],
-                            cc
+                        fc
-                        ).inside(p)
+                    );
-                    )
+
+                    vector proj = p - faceTri.centre();
+
+                    if ((faceTri.normal() & proj) > 0)
                     {
-                        return true;
-                    }
-                }
-            }
                         return false;
                     }
+                }
+            }
+            return true;
+        }
         break;
 
         case FACEDIAGTETS:
         {
-            label tetFaceI, tetPtI;
+            // only test that point is on inside of plane defined by cell face
-            findTetFacePt(cellI, p, tetFaceI, tetPtI);
+            // triangles
+            const cell& cFaces = cells()[cellI];
 
-            return tetFaceI != -1;
+            forAll(cFaces, cFaceI)
+            {
+                label faceI = cFaces[cFaceI];
+                const face& f = faces_[faceI];
+
+                for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
+                {
+                    // Get tetIndices of face triangle
+                    tetIndices faceTetIs
+                    (
+                        polyMeshTetDecomposition::triangleTetIndices
+                        (
+                            *this,
+                            faceI,
+                            cellI,
+                            tetPtI
+                        )
+                    );
+
+                    triPointRef faceTri = faceTetIs.faceTri(*this);
+
+                    vector proj = p - faceTri.centre();
+
+                    if ((faceTri.normal() & proj) > 0)
+                    {
+                        return false;
+                    }
+                }
+            }
+
+            return true;
         }
         break;
     }

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H

@@ -85,12 +85,12 @@ class tetIndices
     label faceBasePtI_;
 
     //- point on the face such that the right-hand circulation
-    //      {faceBasePtI_, facePtIA_, facePtBI_}
+    //      {faceBasePtI_, facePtAI_, facePtBI_}
     //  forms a triangle that points out of the tet
     label facePtAI_;
 
     //- point on the face such that the right-hand circulation
-    //      {faceBasePtI_, facePtIA_, facePtBI_}
+    //      {faceBasePtI_, facePtAI_, facePtBI_}
     //  forms a triangle that points out of the tet
     label facePtBI_;
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.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
@@ -66,8 +66,6 @@ bool Foam::primitiveMesh::pointInCell(const point& p, label celli) const
     const vectorField& cf = faceCentres();
     const vectorField& Sf = faceAreas();
 
-    bool inCell = true;
-
     forAll(f, facei)
     {
         label nFace = f[facei];
@@ -77,10 +75,14 @@ bool Foam::primitiveMesh::pointInCell(const point& p, label celli) const
         {
             normal = -normal;
         }
-        inCell = inCell && ((normal & proj) <= 0);
+
+        if ((normal & proj) > 0)
+        {
+            return false;
+        }
     }
 
-    return inCell;
+    return true;
 }
 
 

src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntryNew.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
@@ -34,7 +34,7 @@ Foam::autoPtr<Foam::DataEntry<Type> > Foam::DataEntry<Type>::New
     const dictionary& dict
 )
 {
-    Istream& is(dict.lookup(entryName));
+    Istream& is(dict.lookup(entryName, false));
 
     token firstToken(is);
 

src/dynamicMesh/Make/files

@@ -84,7 +84,9 @@ polyMeshAdder/polyMeshAdder.C
 fvMeshTools/fvMeshTools.C
 
 motionSmoother/motionSmoother.C
-motionSmoother/motionSmootherCheck.C
+motionSmoother/motionSmootherAlgo.C
+motionSmoother/motionSmootherAlgoCheck.C
+motionSmoother/motionSmootherData.C
 motionSmoother/polyMeshGeometry/polyMeshGeometry.C
 motionSmoother/badQualityToCell/badQualityToCell.C
 motionSmoother/badQualityToFace/badQualityToFace.C

src/dynamicMesh/motionSmoother/motionSmoother.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) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,246 +69,32 @@ Note
 
 SourceFiles
     motionSmoother.C
-    motionSmootherTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef motionSmoother_H
 #define motionSmoother_H
 
-#include "pointFields.H"
+#include "motionSmootherData.H"
-#include "HashSet.H"
+#include "motionSmootherAlgo.H"
-#include "PackedBoolList.H"
-#include "indirectPrimitivePatch.H"
-#include "className.H"
-#include "twoDPointCorrector.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-class polyMeshGeometry;
-class faceSet;
-
 /*---------------------------------------------------------------------------*\
                            Class motionSmoother Declaration
 \*---------------------------------------------------------------------------*/
 
 class motionSmoother
+:
+    public motionSmootherData,
+    public motionSmootherAlgo
+
 {
     // Private class
 
-        //- To synchronise displacements. We want max displacement since
-        //  this is what is specified on pp and internal mesh will have
-        //  zero displacement.
-        class maxMagEqOp
-        {
-
-        public:
-
-            void operator()(vector& x, const vector& y) const
-            {
-                for (direction i = 0; i < vector::nComponents; i++)
-                {
-                    scalar magX = mag(x[i]);
-                    scalar magY = mag(y[i]);
-
-                    if (magX < magY)
-                    {
-                        x[i] = y[i];
-                    }
-                    else if (magX == magY)
-                    {
-                        if (y[i] > x[i])
-                        {
-                            x[i] = y[i];
-                        }
-                    }
-                }
-            }
-        };
-
-
-    // Private data
-
-        //- Reference to polyMesh. Non-const since we move mesh.
-        polyMesh& mesh_;
-
-        //- Reference to pointMesh
-        pointMesh& pMesh_;
-
-        //- Reference to face subset of all adaptPatchIDs
-        indirectPrimitivePatch& pp_;
-
-        //- Indices of fixedValue patches that we're allowed to modify the
-        // displacement on.
-        const labelList adaptPatchIDs_;
-
-
-        // Smoothing and checking parameters
-        dictionary paramDict_;
-
-        // Internal data
-
-        //- Displacement field
-        pointVectorField displacement_;
-
-        //- Scale factor for displacement
-        pointScalarField scale_;
-
-        //- Starting mesh position
-        pointField oldPoints_;
-
-        //- Is mesh point on boundary or not
-        PackedBoolList isInternalPoint_;
-
-        //- Is edge master (always except if on coupled boundary and on
-        //  lower processor)
-        PackedBoolList isMasterEdge_;
-
-        //- 2-D motion corrector
-        twoDPointCorrector twoDCorrector_;
-
-        // Muli-patch constraints (from pointPatchInterpolation)
-
-            labelList patchPatchPointConstraintPoints_;
-            tensorField patchPatchPointConstraintTensors_;
-
-
-    // Private Member Functions
-
-        //- Average of connected points.
-        template<class Type>
-        tmp<GeometricField<Type, pointPatchField, pointMesh> > avg
-        (
-            const GeometricField<Type, pointPatchField, pointMesh>& fld,
-            const scalarField& edgeWeight
-        ) const;
-
-        //- Average postion of connected points.
-        template<class Type>
-        tmp<GeometricField<Type, pointPatchField, pointMesh> > avgPositions
-        (
-            const GeometricField<Type, pointPatchField, pointMesh>& fld,
-            const scalarField& edgeWeight
-        ) const;
-
-        //- Check constraints
-        template<class Type>
-        static void checkConstraints
-        (
-            GeometricField<Type, pointPatchField, pointMesh>&
-        );
-
-        //- Multi-patch constraints
-        template<class Type>
-        void applyCornerConstraints
-        (
-            GeometricField<Type, pointPatchField, pointMesh>&
-        ) const;
-
-        //- Test synchronisation of generic field (not positions!) on points
-        template<class Type, class CombineOp>
-        void testSyncField
-        (
-            const Field<Type>&,
-            const CombineOp& cop,
-            const Type& zero,
-            const scalar maxMag
-        ) const;
-
-        //- Test synchronisation of points
-        void testSyncPositions(const pointField&, const scalar maxMag) const;
-
-        //- Assemble tensors for multi-patch constraints
-        void makePatchPatchAddressing();
-
-        static void checkFld(const pointScalarField&);
-
-        //- Get points used by given faces
-        labelHashSet getPoints(const labelHashSet&) const;
-
-        //- Calculate per-edge weight
-        tmp<scalarField> calcEdgeWeights(const pointField&) const;
-
-        //- explicit smoothing and min on all affected internal points
-        void minSmooth
-        (
-            const scalarField& edgeWeights,
-            const PackedBoolList& isAffectedPoint,
-            const pointScalarField& fld,
-            pointScalarField& newFld
-        ) const;
-
-        //- same but only on selected points (usually patch points)
-        void minSmooth
-        (
-            const scalarField& edgeWeights,
-            const PackedBoolList& isAffectedPoint,
-            const labelList& meshPoints,
-            const pointScalarField& fld,
-            pointScalarField& newFld
-        ) const;
-
-        //- Scale certain (internal) points of a field
-        void scaleField
-        (
-            const labelHashSet& pointLabels,
-            const scalar scale,
-            pointScalarField&
-        ) const;
-
-        //- As above but points have to be in meshPoints as well
-        //  (usually to scale patch points)
-        void scaleField
-        (
-            const labelList& meshPoints,
-            const labelHashSet& pointLabels,
-            const scalar scale,
-            pointScalarField&
-        ) const;
-
-        //- Lower on internal points
-        void subtractField
-        (
-            const labelHashSet& pointLabels,
-            const scalar f,
-            pointScalarField&
-        ) const;
-
-        //- As above but points have to be in meshPoints as well
-        //  (usually to scale patch points)
-        void subtractField
-        (
-            const labelList& meshPoints,
-            const labelHashSet& pointLabels,
-            const scalar scale,
-            pointScalarField&
-        ) const;
-
-        //- Helper function. Is point internal?
-        bool isInternalPoint(const label pointI) const;
-
-        //- Given a set of faces that cause smoothing and a number of
-        //  iterations determine the maximum set of points who are affected
-        //  and the accordingly affected faces.
-        void getAffectedFacesAndPoints
-        (
-            const label nPointIter,
-            const faceSet& wrongFaces,
-
-            labelList& affectedFaces,
-            PackedBoolList& isAffectedPoint
-        ) const;
-
-        //- Disallow default bitwise copy construct
-        motionSmoother(const motionSmoother&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const motionSmoother&);
-
-
 public:
 
     ClassName("motionSmoother");
@@ -337,212 +123,15 @@ public:
             const dictionary& paramDict
         );
 
-
-    //- Destructor
-    ~motionSmoother();
-
-
-    // Member Functions
-
-        // Access
-
-            //- Reference to mesh
-            const polyMesh& mesh() const;
-
-            //- Reference to pointMesh
-            const pointMesh& pMesh() const;
-
-            //- Reference to patch
-            const indirectPrimitivePatch& patch() const;
-
-            //- Patch labels that are being adapted
-            const labelList& adaptPatchIDs() const;
-
-            const dictionary& paramDict() const;
-
-            //- Reference to displacement field
-            pointVectorField& displacement();
-
-            //- Reference to displacement field
-            const pointVectorField& displacement() const;
-
-            //- Reference to scale field
-            const pointScalarField& scale() const;
-
-            //- Starting mesh position
-            const pointField& oldPoints() const;
-
-            //- Return reference to 2D point motion correction
-            twoDPointCorrector& twoDCorrector()
-            {
-                return twoDCorrector_;
-            }
-
-
-
-        // Edit
-
-            //- Take over existing mesh position.
-            void correct();
-
-            //- Set patch fields on displacement to be consistent with
-            //  internal values.
-            void setDisplacementPatchFields();
-
-            //- Set displacement field from displacement on patch points.
-            //  Modify provided displacement to be consistent with actual
-            //  boundary conditions on displacement. Note: resets the
-            //  displacement to be 0 on coupled patches beforehand
-            //  to make sure shared points
-            //  partially on pp (on some processors) and partially not
-            //  (on other processors) get the value from pp.
-            void setDisplacement(pointField& patchDisp);
-
-            //- Special correctBoundaryConditions which evaluates fixedValue
-            //  patches first so they get overwritten with any constraint
-            //  bc's.
-            void correctBoundaryConditions(pointVectorField&) const;
-
-            //- Move mesh. Does 2D correction (modifies passed pointField) and
-            //  polyMesh::movePoints. Returns swept volumes.
-            tmp<scalarField> movePoints(pointField&);
-
-            //- Set the errorReduction (by how much to scale the displacement
-            //  at error locations) parameter. Returns the old value.
-            //  Set to 0 (so revert to old mesh) grows out one cell layer
-            //  from error faces.
-            scalar setErrorReduction(const scalar);
-
-            //- Move mesh with given scale. Return true if mesh ok or has
-            //  less than nAllow errors, false
-            //  otherwise and locally update scale. Smoothmesh=false means only
-            //  patch points get moved.
-            //  Parallel ok (as long as displacement field is consistent
-            //  across patches)
-            bool scaleMesh
-            (
-                labelList& checkFaces,
-                const bool smoothMesh = true,
-                const label nAllow = 0
-            );
-
-            //- Move mesh (with baffles) with given scale.
-            bool scaleMesh
-            (
-                labelList& checkFaces,
-                const List<labelPair>& baffles,
-                const bool smoothMesh = true,
-                const label nAllow = 0
-            );
-
-            //- Move mesh with externally provided mesh constraints
-            bool scaleMesh
-            (
-                labelList& checkFaces,
-                const List<labelPair>& baffles,
-                const dictionary& paramDict,
-                const dictionary& meshQualityDict,
-                const bool smoothMesh = true,
-                const label nAllow = 0
-            );
-
-            //- Update topology
-            void updateMesh();
-
-            //- Check mesh with mesh settings in dict. Collects incorrect faces
-            //  in set. Returns true if one or more faces in error.
-            //  Parallel ok.
-            static bool checkMesh
-            (
-                const bool report,
-                const polyMesh& mesh,
-                const dictionary& dict,
-                labelHashSet& wrongFaces
-            );
-
-            //- Check (subset of mesh) with mesh settings in dict.
-            //  Collects incorrect faces in set. Returns true if one
-            //  or more faces in error. Parallel ok.
-            static bool checkMesh
-            (
-                const bool report,
-                const polyMesh& mesh,
-                const dictionary& dict,
-                const labelList& checkFaces,
-                labelHashSet& wrongFaces
-            );
-
-            //- Check (subset of mesh including baffles) with mesh settings
-            //  in dict. Collects incorrect faces in set. Returns true if one
-            //  or more faces in error. Parallel ok.
-            static bool checkMesh
-            (
-                const bool report,
-                const polyMesh& mesh,
-                const dictionary& dict,
-                const labelList& checkFaces,
-                const List<labelPair>& baffles,
-                labelHashSet& wrongFaces
-            );
-
-            //- Check part of mesh with mesh settings in dict.
-            //  Collects incorrect faces in set. Returns true if one or
-            //  more faces in error. Parallel ok.
-            static bool checkMesh
-            (
-                const bool report,
-                const dictionary& dict,
-                const polyMeshGeometry&,
-                const labelList& checkFaces,
-                labelHashSet& wrongFaces
-            );
-
-            //- Check part of mesh including baffles with mesh settings in dict.
-            //  Collects incorrect faces in set. Returns true if one or
-            //  more faces in error. Parallel ok.
-            static bool checkMesh
-            (
-                const bool report,
-                const dictionary& dict,
-                const polyMeshGeometry&,
-                const labelList& checkFaces,
-                const List<labelPair>& baffles,
-                labelHashSet& wrongFaces
-            );
-
-            // Helper functions to manipulate displacement vector.
-
-                //- Fully explicit smoothing of fields (not positions)
-                //  of internal points with varying diffusivity.
-                template<class Type>
-                void smooth
-                (
-                    const GeometricField<Type, pointPatchField, pointMesh>& fld,
-                    const scalarField& edgeWeight,
-                    GeometricField<Type, pointPatchField, pointMesh>& newFld
-                ) const;
 };
 
 
-template<>
-void motionSmoother::applyCornerConstraints<scalar>
-(
-    GeometricField<scalar, pointPatchField, pointMesh>& pf
-) const;
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-#   include "motionSmootherTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/dynamicMesh/motionSmoother/motionSmootherAlgo.H

@@ -0,0 +1,534 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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::motionSmootherAlgo
+
+Description
+    Given a displacement moves the mesh by scaling the displacement back
+    until there are no more mesh errors.
+
+    Holds displacement field (read upon construction since need boundary
+    conditions) and scaling factor and optional patch number on which to
+    scale back displacement.
+
+    E.g.
+    \verbatim
+        // Construct iterative mesh mover.
+        motionSmoother meshMover(mesh, labelList(1, patchI));
+
+        // Set desired displacement:
+        meshMover.displacement() = ..
+
+        for (label iter = 0; iter < maxIter; iter++)
+        {
+            if (meshMover.scaleMesh(true))
+            {
+                Info<< "Successfully moved mesh" << endl;
+                return true;
+            }
+        }
+    \endverbatim
+
+Note
+    - Shared points (parallel): a processor can have points which are part of
+    pp on another processor but have no pp itself (i.e. it has points
+    and/or edges but no faces of pp). Hence we have to be careful when e.g.
+    synchronising displacements that the value from the processor which has
+    faces of pp get priority. This is currently handled in setDisplacement
+    by resetting the internal displacement to zero before doing anything
+    else. The combine operator used will give preference to non-zero
+    values.
+
+    - Various routines take baffles. These are sets of boundary faces that
+    are treated as a single internal face. This is a hack used to apply
+    movement to internal faces.
+
+    - Mesh constraints are looked up from the supplied dictionary. (uses
+    recursive lookup)
+
+SourceFiles
+    motionSmootherAlgo.C
+    motionSmootherAlgoTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef motionSmootherAlgo_H
+#define motionSmootherAlgo_H
+
+#include "pointFields.H"
+#include "HashSet.H"
+#include "PackedBoolList.H"
+#include "indirectPrimitivePatch.H"
+#include "className.H"
+#include "twoDPointCorrector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMeshGeometry;
+class faceSet;
+
+/*---------------------------------------------------------------------------*\
+                           Class motionSmootherAlgo Declaration
+\*---------------------------------------------------------------------------*/
+
+class motionSmootherAlgo
+{
+    // Private class
+
+        //- To synchronise displacements. We want max displacement since
+        //  this is what is specified on pp and internal mesh will have
+        //  zero displacement.
+        class maxMagEqOp
+        {
+
+        public:
+
+            void operator()(vector& x, const vector& y) const
+            {
+                for (direction i = 0; i < vector::nComponents; i++)
+                {
+                    scalar magX = mag(x[i]);
+                    scalar magY = mag(y[i]);
+
+                    if (magX < magY)
+                    {
+                        x[i] = y[i];
+                    }
+                    else if (magX == magY)
+                    {
+                        if (y[i] > x[i])
+                        {
+                            x[i] = y[i];
+                        }
+                    }
+                }
+            }
+        };
+
+
+    // Private data
+
+        //- Reference to polyMesh. Non-const since we move mesh.
+        polyMesh& mesh_;
+
+        //- Reference to pointMesh
+        pointMesh& pMesh_;
+
+        //- Reference to face subset of all adaptPatchIDs
+        indirectPrimitivePatch& pp_;
+
+        //- Displacement field
+        pointVectorField& displacement_;
+
+        //- Scale factor for displacement
+        pointScalarField& scale_;
+
+        //- Starting mesh position
+        pointField& oldPoints_;
+
+
+        // Internal data
+
+            //- Indices of fixedValue patches that we're allowed to modify the
+            // displacement on.
+            const labelList adaptPatchIDs_;
+
+            // Smoothing and checking parameters
+            dictionary paramDict_;
+
+            //- Is mesh point on boundary or not
+            PackedBoolList isInternalPoint_;
+
+            //- Is edge master (always except if on coupled boundary and on
+            //  lower processor)
+            PackedBoolList isMasterEdge_;
+
+            //- 2-D motion corrector
+            twoDPointCorrector twoDCorrector_;
+
+            // Muli-patch constraints (from pointPatchInterpolation)
+
+                labelList patchPatchPointConstraintPoints_;
+                tensorField patchPatchPointConstraintTensors_;
+
+
+    // Private Member Functions
+
+        //- Average of connected points.
+        template<class Type>
+        tmp<GeometricField<Type, pointPatchField, pointMesh> > avg
+        (
+            const GeometricField<Type, pointPatchField, pointMesh>& fld,
+            const scalarField& edgeWeight
+        ) const;
+
+        //- Average postion of connected points.
+        template<class Type>
+        tmp<GeometricField<Type, pointPatchField, pointMesh> > avgPositions
+        (
+            const GeometricField<Type, pointPatchField, pointMesh>& fld,
+            const scalarField& edgeWeight
+        ) const;
+
+        //- Check constraints
+        template<class Type>
+        static void checkConstraints
+        (
+            GeometricField<Type, pointPatchField, pointMesh>&
+        );
+
+        //- Multi-patch constraints
+        template<class Type>
+        void applyCornerConstraints
+        (
+            GeometricField<Type, pointPatchField, pointMesh>&
+        ) const;
+
+        //- Test synchronisation of generic field (not positions!) on points
+        template<class Type, class CombineOp>
+        void testSyncField
+        (
+            const Field<Type>&,
+            const CombineOp& cop,
+            const Type& zero,
+            const scalar maxMag
+        ) const;
+
+        //- Test synchronisation of points
+        void testSyncPositions(const pointField&, const scalar maxMag) const;
+
+        //- Assemble tensors for multi-patch constraints
+        void makePatchPatchAddressing();
+
+        static void checkFld(const pointScalarField&);
+
+        //- Get points used by given faces
+        labelHashSet getPoints(const labelHashSet&) const;
+
+        //- Calculate per-edge weight
+        tmp<scalarField> calcEdgeWeights(const pointField&) const;
+
+        //- explicit smoothing and min on all affected internal points
+        void minSmooth
+        (
+            const scalarField& edgeWeights,
+            const PackedBoolList& isAffectedPoint,
+            const pointScalarField& fld,
+            pointScalarField& newFld
+        ) const;
+
+        //- same but only on selected points (usually patch points)
+        void minSmooth
+        (
+            const scalarField& edgeWeights,
+            const PackedBoolList& isAffectedPoint,
+            const labelList& meshPoints,
+            const pointScalarField& fld,
+            pointScalarField& newFld
+        ) const;
+
+        //- Scale certain (internal) points of a field
+        void scaleField
+        (
+            const labelHashSet& pointLabels,
+            const scalar scale,
+            pointScalarField&
+        ) const;
+
+        //- As above but points have to be in meshPoints as well
+        //  (usually to scale patch points)
+        void scaleField
+        (
+            const labelList& meshPoints,
+            const labelHashSet& pointLabels,
+            const scalar scale,
+            pointScalarField&
+        ) const;
+
+        //- Lower on internal points
+        void subtractField
+        (
+            const labelHashSet& pointLabels,
+            const scalar f,
+            pointScalarField&
+        ) const;
+
+        //- As above but points have to be in meshPoints as well
+        //  (usually to scale patch points)
+        void subtractField
+        (
+            const labelList& meshPoints,
+            const labelHashSet& pointLabels,
+            const scalar scale,
+            pointScalarField&
+        ) const;
+
+        //- Helper function. Is point internal?
+        bool isInternalPoint(const label pointI) const;
+
+        //- Given a set of faces that cause smoothing and a number of
+        //  iterations determine the maximum set of points who are affected
+        //  and the accordingly affected faces.
+        void getAffectedFacesAndPoints
+        (
+            const label nPointIter,
+            const faceSet& wrongFaces,
+
+            labelList& affectedFaces,
+            PackedBoolList& isAffectedPoint
+        ) const;
+
+        //- Disallow default bitwise copy construct
+        motionSmootherAlgo(const motionSmootherAlgo&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const motionSmootherAlgo&);
+
+
+public:
+
+    ClassName("motionSmootherAlgo");
+
+    // Constructors
+
+        //- Construct from mesh, patches to work on and smoothing parameters.
+        motionSmootherAlgo
+        (
+            polyMesh&,
+            pointMesh&,
+            indirectPrimitivePatch& pp,         // 'outside' points
+            pointVectorField& displacement,
+            pointScalarField& scale,
+            pointField& oldPoints,
+            const labelList& adaptPatchIDs,     // patches forming 'outside'
+            const dictionary& paramDict
+        );
+
+
+    //- Destructor
+    ~motionSmootherAlgo();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Reference to mesh
+            const polyMesh& mesh() const;
+
+            //- Reference to pointMesh
+            const pointMesh& pMesh() const;
+
+            //- Reference to patch
+            const indirectPrimitivePatch& patch() const;
+
+            //- Patch labels that are being adapted
+            const labelList& adaptPatchIDs() const;
+
+            const dictionary& paramDict() const;
+
+            //- Return reference to 2D point motion correction
+            twoDPointCorrector& twoDCorrector()
+            {
+                return twoDCorrector_;
+            }
+
+
+
+        // Edit
+
+            //- Take over existing mesh position.
+            void correct();
+
+            //- Set patch fields on displacement to be consistent with
+            //  internal values.
+            void setDisplacementPatchFields();
+
+            //- Set displacement field from displacement on patch points.
+            //  Modify provided displacement to be consistent with actual
+            //  boundary conditions on displacement. Note: resets the
+            //  displacement to be 0 on coupled patches beforehand
+            //  to make sure shared points
+            //  partially on pp (on some processors) and partially not
+            //  (on other processors) get the value from pp.
+            void setDisplacement(pointField& patchDisp);
+
+            //- Special correctBoundaryConditions which evaluates fixedValue
+            //  patches first so they get overwritten with any constraint
+            //  bc's.
+            void correctBoundaryConditions(pointVectorField&) const;
+
+            //- Apply optional point constraint (2d correction)
+            void modifyMotionPoints(pointField& newPoints) const;
+
+            //- Get the current points (oldPoints+scale*displacement)
+            tmp<pointField> curPoints() const;
+
+            //- Set the errorReduction (by how much to scale the displacement
+            //  at error locations) parameter. Returns the old value.
+            //  Set to 0 (so revert to old mesh) grows out one cell layer
+            //  from error faces.
+            scalar setErrorReduction(const scalar);
+
+            //- Move mesh with given scale. Return true if mesh ok or has
+            //  less than nAllow errors, false
+            //  otherwise and locally update scale. Smoothmesh=false means only
+            //  patch points get moved.
+            //  Parallel ok (as long as displacement field is consistent
+            //  across patches)
+            bool scaleMesh
+            (
+                labelList& checkFaces,
+                const bool smoothMesh = true,
+                const label nAllow = 0
+            );
+
+            //- Move mesh (with baffles) with given scale.
+            bool scaleMesh
+            (
+                labelList& checkFaces,
+                const List<labelPair>& baffles,
+                const bool smoothMesh = true,
+                const label nAllow = 0
+            );
+
+            //- Move mesh with externally provided mesh constraints
+            bool scaleMesh
+            (
+                labelList& checkFaces,
+                const List<labelPair>& baffles,
+                const dictionary& paramDict,
+                const dictionary& meshQualityDict,
+                const bool smoothMesh = true,
+                const label nAllow = 0
+            );
+
+
+            //- Update for new mesh geometry
+            void movePoints();
+
+            //- Update for new mesh topology
+            void updateMesh();
+
+
+            //- Check mesh with mesh settings in dict. Collects incorrect faces
+            //  in set. Returns true if one or more faces in error.
+            //  Parallel ok.
+            static bool checkMesh
+            (
+                const bool report,
+                const polyMesh& mesh,
+                const dictionary& dict,
+                labelHashSet& wrongFaces
+            );
+
+            //- Check (subset of mesh) with mesh settings in dict.
+            //  Collects incorrect faces in set. Returns true if one
+            //  or more faces in error. Parallel ok.
+            static bool checkMesh
+            (
+                const bool report,
+                const polyMesh& mesh,
+                const dictionary& dict,
+                const labelList& checkFaces,
+                labelHashSet& wrongFaces
+            );
+
+            //- Check (subset of mesh including baffles) with mesh settings
+            //  in dict. Collects incorrect faces in set. Returns true if one
+            //  or more faces in error. Parallel ok.
+            static bool checkMesh
+            (
+                const bool report,
+                const polyMesh& mesh,
+                const dictionary& dict,
+                const labelList& checkFaces,
+                const List<labelPair>& baffles,
+                labelHashSet& wrongFaces
+            );
+
+            //- Check part of mesh with mesh settings in dict.
+            //  Collects incorrect faces in set. Returns true if one or
+            //  more faces in error. Parallel ok.
+            static bool checkMesh
+            (
+                const bool report,
+                const dictionary& dict,
+                const polyMeshGeometry&,
+                const labelList& checkFaces,
+                labelHashSet& wrongFaces
+            );
+
+            //- Check part of mesh including baffles with mesh settings in dict.
+            //  Collects incorrect faces in set. Returns true if one or
+            //  more faces in error. Parallel ok.
+            static bool checkMesh
+            (
+                const bool report,
+                const dictionary& dict,
+                const polyMeshGeometry&,
+                const labelList& checkFaces,
+                const List<labelPair>& baffles,
+                labelHashSet& wrongFaces
+            );
+
+            // Helper functions to manipulate displacement vector.
+
+                //- Fully explicit smoothing of fields (not positions)
+                //  of internal points with varying diffusivity.
+                template<class Type>
+                void smooth
+                (
+                    const GeometricField<Type, pointPatchField, pointMesh>& fld,
+                    const scalarField& edgeWeight,
+                    GeometricField<Type, pointPatchField, pointMesh>& newFld
+                ) const;
+};
+
+
+template<>
+void motionSmootherAlgo::applyCornerConstraints<scalar>
+(
+    GeometricField<scalar, pointPatchField, pointMesh>& pf
+) const;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "motionSmootherAlgoTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C

@@ -23,13 +23,13 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "motionSmoother.H"
+#include "motionSmootherAlgo.H"
 #include "polyMeshGeometry.H"
 #include "IOmanip.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::motionSmoother::checkMesh
+bool Foam::motionSmootherAlgo::checkMesh
 (
     const bool report,
     const polyMesh& mesh,
@@ -50,7 +50,7 @@ bool Foam::motionSmoother::checkMesh
     );
 }
 
-bool Foam::motionSmoother::checkMesh
+bool Foam::motionSmootherAlgo::checkMesh
 (
     const bool report,
     const polyMesh& mesh,
@@ -411,7 +411,7 @@ bool Foam::motionSmoother::checkMesh
 }
 
 
-bool Foam::motionSmoother::checkMesh
+bool Foam::motionSmootherAlgo::checkMesh
 (
     const bool report,
     const polyMesh& mesh,
@@ -429,7 +429,7 @@ bool Foam::motionSmoother::checkMesh
     );
 }
 
-bool Foam::motionSmoother::checkMesh
+bool Foam::motionSmootherAlgo::checkMesh
 (
     const bool report,
     const dictionary& dict,
@@ -452,7 +452,7 @@ bool Foam::motionSmoother::checkMesh
 }
 
 
-bool Foam::motionSmoother::checkMesh
+bool Foam::motionSmootherAlgo::checkMesh
 (
     const bool report,
     const dictionary& dict,

src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C

@@ -23,7 +23,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "motionSmoother.H"
+#include "motionSmootherAlgo.H"
 #include "meshTools.H"
 #include "processorPointPatchFields.H"
 #include "pointConstraint.H"
@@ -32,7 +32,7 @@ License
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 template<class Type>
-void Foam::motionSmoother::checkConstraints
+void Foam::motionSmootherAlgo::checkConstraints
 (
     GeometricField<Type, pointPatchField, pointMesh>& pf
 )
@@ -119,7 +119,7 @@ void Foam::motionSmoother::checkConstraints
                 {
                     FatalErrorIn
                     (
-                        "motionSmoother::checkConstraints"
+                        "motionSmootherAlgo::checkConstraints"
                         "(GeometricField<Type, pointPatchField, pointMesh>&)"
                     )   << "Patch fields are not consistent on mesh point "
                         << ppp << " coordinate " << mesh.points()[ppp]
@@ -136,7 +136,7 @@ void Foam::motionSmoother::checkConstraints
 
 
 template<class Type>
-void Foam::motionSmoother::applyCornerConstraints
+void Foam::motionSmootherAlgo::applyCornerConstraints
 (
     GeometricField<Type, pointPatchField, pointMesh>& pf
 ) const
@@ -155,7 +155,7 @@ void Foam::motionSmoother::applyCornerConstraints
 // Average of connected points.
 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::pointPatchField, Foam::pointMesh> >
-Foam::motionSmoother::avg
+Foam::motionSmootherAlgo::avg
 (
     const GeometricField<Type, pointPatchField, pointMesh>& fld,
     const scalarField& edgeWeight
@@ -253,7 +253,7 @@ Foam::motionSmoother::avg
 
 // smooth field (point-jacobi)
 template<class Type>
-void Foam::motionSmoother::smooth
+void Foam::motionSmootherAlgo::smooth
 (
     const GeometricField<Type, pointPatchField, pointMesh>& fld,
     const scalarField& edgeWeight,
@@ -278,7 +278,7 @@ void Foam::motionSmoother::smooth
 
 //- Test synchronisation of generic field (not positions!) on points
 template<class Type, class CombineOp>
-void Foam::motionSmoother::testSyncField
+void Foam::motionSmootherAlgo::testSyncField
 (
     const Field<Type>& fld,
     const CombineOp& cop,
@@ -308,7 +308,7 @@ void Foam::motionSmoother::testSyncField
         {
             FatalErrorIn
             (
-                "motionSmoother::testSyncField"
+                "motionSmootherAlgo::testSyncField"
                 "(const Field<Type>&, const CombineOp&"
                 ", const Type&, const bool)"
             )   << "On element " << i << " value:" << fld[i]

src/dynamicMesh/motionSmoother/motionSmootherData.C

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "motionSmootherData.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::motionSmootherData::motionSmootherData
+(
+    const pointMesh& pMesh
+)
+:
+    displacement_
+    (
+        IOobject
+        (
+            "displacement",
+            pMesh.time().timeName(),
+            pMesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pMesh
+    ),
+    scale_
+    (
+        IOobject
+        (
+            "scale",
+            pMesh.time().timeName(),
+            pMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pMesh,
+        dimensionedScalar("scale", dimless, 1.0)
+    ),
+    oldPoints_(pMesh().points())
+{}
+
+
+Foam::motionSmootherData::motionSmootherData
+(
+    const pointVectorField& displacement
+)
+:
+    displacement_
+    (
+        IOobject
+        (
+            "displacement",
+            displacement.time().timeName(),
+            displacement.mesh()(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        displacement
+    ),
+    scale_
+    (
+        IOobject
+        (
+            "scale",
+            displacement.time().timeName(),
+            displacement.mesh()(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        displacement.mesh(),
+        dimensionedScalar("scale", dimless, 1.0)
+    ),
+    oldPoints_(displacement.mesh()().points())
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::pointVectorField& Foam::motionSmootherData::displacement()
+{
+    return displacement_;
+}
+
+
+const Foam::pointVectorField& Foam::motionSmootherData::displacement() const
+{
+    return displacement_;
+}
+
+
+const Foam::pointScalarField& Foam::motionSmootherData::scale() const
+{
+    return scale_;
+}
+
+
+const Foam::pointField& Foam::motionSmootherData::oldPoints() const
+{
+    return oldPoints_;
+}
+
+
+// ************************************************************************* //

src/dynamicMesh/motionSmoother/motionSmootherData.H

@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::motionSmootherData
+
+Description
+
+SourceFiles
+    motionSmootherData.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef motionSmootherData_H
+#define motionSmootherData_H
+
+#include "pointFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class motionSmootherData Declaration
+\*---------------------------------------------------------------------------*/
+
+class motionSmootherData
+{
+protected:
+
+    // Private data
+
+        //- Displacement field
+        pointVectorField displacement_;
+
+        //- Scale factor for displacement
+        pointScalarField scale_;
+
+        //- Starting mesh position
+        pointField oldPoints_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct read
+        motionSmootherData
+        (
+            const pointMesh&
+        );
+
+        //- Construct from pointDisplacement
+        motionSmootherData
+        (
+            const pointVectorField&
+        );
+
+
+    // Member Functions
+
+        //- Reference to displacement field
+        pointVectorField& displacement();
+
+        //- Reference to displacement field
+        const pointVectorField& displacement() const;
+
+        //- Reference to scale field
+        const pointScalarField& scale() const;
+
+        //- Starting mesh position
+        const pointField& oldPoints() const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C

@@ -3185,7 +3185,8 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
 
     forAll(refineCell, cellI)
     {
-        if (refineCell.get(cellI))
+//        if (refineCell.get(cellI))
+        if (refineCell[cellI])
         {
             nRefined++;
         }
@@ -3196,7 +3197,8 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
 
     forAll(refineCell, cellI)
     {
-        if (refineCell.get(cellI))
+//        if (refineCell.get(cellI))
+        if (refineCell[cellI])
         {
             newCellsToRefine[nRefined++] = cellI;
         }

src/finiteVolume/Make/files

@@ -189,9 +189,11 @@ $(derivedFvPatchFields)/uniformJumpAMI/uniformJumpAMIFvPatchFields.C
 $(derivedFvPatchFields)/uniformTotalPressure/uniformTotalPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/variableHeightFlowRate/variableHeightFlowRateFvPatchField.C
 $(derivedFvPatchFields)/variableHeightFlowRateInletVelocity/variableHeightFlowRateInletVelocityFvPatchVectorField.C
+$(derivedFvPatchFields)/outletPhaseMeanVelocity/outletPhaseMeanVelocityFvPatchVectorField.C
 $(derivedFvPatchFields)/waveTransmissive/waveTransmissiveFvPatchFields.C
 $(derivedFvPatchFields)/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
 $(derivedFvPatchFields)/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.C
+$(derivedFvPatchFields)/prghPressure/prghPressureFvPatchScalarField.C
 
 fvsPatchFields = fields/fvsPatchFields
 $(fvsPatchFields)/fvsPatchField/fvsPatchFields.C

src/finiteVolume/cfdTools/general/include/alphaControls.H

@@ -10,3 +10,11 @@ bool alphaOuterCorrectors
 (
     alphaControls.lookupOrDefault<Switch>("alphaOuterCorrectors", false)
 );
+
+// Apply the compression correction from the previous iteration
+// Improves efficiency for steady-simulations but can only be applied
+// once the alpha field is reasonably steady, i.e. fully developed
+bool alphaApplyPrevCorr
+(
+    alphaControls.lookupOrDefault<Switch>("alphaApplyPrevCorr", false)
+);

src/finiteVolume/fields/fvPatchFields/derived/outletPhaseMeanVelocity/outletPhaseMeanVelocityFvPatchVectorField.C

@@ -0,0 +1,188 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "outletPhaseMeanVelocityFvPatchVectorField.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::outletPhaseMeanVelocityFvPatchVectorField
+::outletPhaseMeanVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    mixedFvPatchField<vector>(p, iF),
+    Umean_(0),
+    alphaName_("none")
+{
+    refValue() = vector::zero;
+    refGrad() = vector::zero;
+    valueFraction() = 0.0;
+}
+
+
+Foam::outletPhaseMeanVelocityFvPatchVectorField
+::outletPhaseMeanVelocityFvPatchVectorField
+(
+    const outletPhaseMeanVelocityFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchField<vector>(ptf, p, iF, mapper),
+    Umean_(ptf.Umean_),
+    alphaName_(ptf.alphaName_)
+{}
+
+
+Foam::outletPhaseMeanVelocityFvPatchVectorField
+::outletPhaseMeanVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchField<vector>(p, iF),
+    Umean_(readScalar(dict.lookup("Umean"))),
+    alphaName_(dict.lookup("alpha"))
+{
+    refValue() = vector::zero;
+    refGrad() = vector::zero;
+    valueFraction() = 0.0;
+
+    if (dict.found("value"))
+    {
+        fvPatchVectorField::operator=
+        (
+            vectorField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        fvPatchVectorField::operator=(patchInternalField());
+    }
+}
+
+
+Foam::outletPhaseMeanVelocityFvPatchVectorField
+::outletPhaseMeanVelocityFvPatchVectorField
+(
+    const outletPhaseMeanVelocityFvPatchVectorField& ptf
+)
+:
+    mixedFvPatchField<vector>(ptf),
+    Umean_(ptf.Umean_),
+    alphaName_(ptf.alphaName_)
+{}
+
+
+Foam::outletPhaseMeanVelocityFvPatchVectorField
+::outletPhaseMeanVelocityFvPatchVectorField
+(
+    const outletPhaseMeanVelocityFvPatchVectorField& ptf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    mixedFvPatchField<vector>(ptf, iF),
+    Umean_(ptf.Umean_),
+    alphaName_(ptf.alphaName_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::outletPhaseMeanVelocityFvPatchVectorField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    scalarField alphap =
+        patch().lookupPatchField<volScalarField, scalar>(alphaName_);
+
+    alphap = max(alphap, scalar(0));
+    alphap = min(alphap, scalar(1));
+
+    // Get the patchInternalField (zero-gradient field)
+    vectorField Uzg(patchInternalField());
+
+    // Calculate the phase mean zero-gradient velocity
+    scalar Uzgmean =
+        gSum(alphap*(patch().Sf() & Uzg))
+       /gSum(alphap*patch().magSf());
+
+    // Set the refValue and valueFraction to adjust the boundary field
+    // such that the phase mean is Umean_
+    if (Uzgmean >= Umean_)
+    {
+        refValue() = vector::zero;
+        valueFraction() = 1.0 - Umean_/Uzgmean;
+    }
+    else
+    {
+        refValue() = (Umean_ + Uzgmean)*patch().nf();
+        valueFraction() = 1.0 - Uzgmean/Umean_;
+    }
+
+    mixedFvPatchField<vector>::updateCoeffs();
+}
+
+
+void Foam::outletPhaseMeanVelocityFvPatchVectorField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchField<vector>::write(os);
+
+    os.writeKeyword("Umean") << Umean_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("alpha") << alphaName_
+        << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+   makePatchTypeField
+   (
+       fvPatchVectorField,
+       outletPhaseMeanVelocityFvPatchVectorField
+   );
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/outletPhaseMeanVelocity/outletPhaseMeanVelocityFvPatchVectorField.H

@@ -0,0 +1,197 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::outletPhaseMeanVelocityFvPatchVectorField
+
+Group
+    grpOutletBoundaryConditions
+
+Description
+    This boundary condition adjusts the velocity for the given phase to achieve
+    the specified mean thus causing the phase-fraction to adjust according to
+    the mass flow rate.
+
+    Typical usage is as the outlet condition for a towing-tank ship simulation
+    to maintain the outlet water level at the level as the inlet.
+
+    \heading Patch usage
+    \table
+        Property     | Description             | Required    | Default value
+        Umean        | mean velocity normal to the boundary [m/s] | yes |
+        alpha        | phase-fraction field    | yes |
+    \endtable
+
+    Example of the boundary condition specification:
+    \verbatim
+    myPatch
+    {
+        type            outletPhaseMeanVelocity;
+        Umean           1.2;
+        alpha           alpha.water;
+        value           uniform (1.2 0 0);
+    }
+    \endverbatim
+
+SeeAlso
+    Foam::mixedFvPatchField
+    Foam::variableHeightFlowRateInletVelocityFvPatchVectorField
+
+SourceFiles
+    outletPhaseMeanVelocityFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef outletPhaseMeanVelocityFvPatchVectorField_H
+#define outletPhaseMeanVelocityFvPatchVectorField_H
+
+#include "mixedFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+    Class outletPhaseMeanVelocityFvPatchVectorField Declaration
+\*---------------------------------------------------------------------------*/
+
+class outletPhaseMeanVelocityFvPatchVectorField
+:
+    public mixedFvPatchVectorField
+{
+    // Private data
+
+        //- Inlet integral flow rate
+        scalar Umean_;
+
+        //- Name of the phase-fraction field
+        word alphaName_;
+
+
+public:
+
+   //- Runtime type information
+   TypeName("outletPhaseMeanVelocity");
+
+
+   // Constructors
+
+        //- Construct from patch and internal field
+        outletPhaseMeanVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        outletPhaseMeanVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  outletPhaseMeanVelocityFvPatchVectorField
+        //  onto a new patch
+        outletPhaseMeanVelocityFvPatchVectorField
+        (
+            const outletPhaseMeanVelocityFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        outletPhaseMeanVelocityFvPatchVectorField
+        (
+            const outletPhaseMeanVelocityFvPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new outletPhaseMeanVelocityFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        outletPhaseMeanVelocityFvPatchVectorField
+        (
+            const outletPhaseMeanVelocityFvPatchVectorField&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchVectorField> clone
+        (
+            const DimensionedField<vector, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new outletPhaseMeanVelocityFvPatchVectorField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+
+    // Member functions
+
+        // Access
+
+            //- Return the flux
+            scalar Umean() const
+            {
+                return Umean_;
+            }
+
+            //- Return reference to the flux to allow adjustment
+            scalar& Umean()
+            {
+                return Umean_;
+            }
+
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/prghPressure/prghPressureFvPatchScalarField.C

@@ -0,0 +1,185 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "prghPressureFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "uniformDimensionedFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::prghPressureFvPatchScalarField::
+prghPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    rhoName_("rho"),
+    p_(p.size(), 0.0)
+{}
+
+
+Foam::prghPressureFvPatchScalarField::
+prghPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    rhoName_(dict.lookupOrDefault<word>("rhoName", "rho")),
+    p_("p", dict, p.size())
+{
+    if (dict.found("value"))
+    {
+        fvPatchScalarField::operator=
+        (
+            scalarField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        fvPatchField<scalar>::operator=(p_);
+    }
+}
+
+
+Foam::prghPressureFvPatchScalarField::
+prghPressureFvPatchScalarField
+(
+    const prghPressureFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    rhoName_(ptf.rhoName_),
+    p_(ptf.p_, mapper)
+{}
+
+
+Foam::prghPressureFvPatchScalarField::
+prghPressureFvPatchScalarField
+(
+    const prghPressureFvPatchScalarField& ptf
+)
+:
+    fixedValueFvPatchScalarField(ptf),
+    rhoName_(ptf.rhoName_),
+    p_(ptf.p_)
+{}
+
+
+Foam::prghPressureFvPatchScalarField::
+prghPressureFvPatchScalarField
+(
+    const prghPressureFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(ptf, iF),
+    rhoName_(ptf.rhoName_),
+    p_(ptf.p_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::prghPressureFvPatchScalarField::autoMap
+(
+    const fvPatchFieldMapper& m
+)
+{
+    fixedValueFvPatchScalarField::autoMap(m);
+    p_.autoMap(m);
+}
+
+
+void Foam::prghPressureFvPatchScalarField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList& addr
+)
+{
+    fixedValueFvPatchScalarField::rmap(ptf, addr);
+
+    const prghPressureFvPatchScalarField& tiptf =
+        refCast<const prghPressureFvPatchScalarField>(ptf);
+
+    p_.rmap(tiptf.p_, addr);
+}
+
+
+void Foam::prghPressureFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const scalarField& rhop = patch().lookupPatchField<volScalarField, scalar>
+    (
+        rhoName_
+    );
+
+    const uniformDimensionedVectorField& g =
+        db().lookupObject<uniformDimensionedVectorField>("g");
+
+    operator==(p_ - rhop*((g.value() & patch().Cf())));
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::prghPressureFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    if (rhoName_ != "rho")
+    {
+        os.writeKeyword("rhoName")
+            << rhoName_ << token::END_STATEMENT << nl;
+    }
+    p_.writeEntry("p", os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        prghPressureFvPatchScalarField
+    );
+}
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/prghPressure/prghPressureFvPatchScalarField.H

@@ -0,0 +1,236 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::prghPressureFvPatchScalarField
+
+Group
+    grpGenericBoundaryConditions
+
+Description
+    This boundary condition provides static pressure condition for p_rgh,
+    calculated as:
+
+        \f[
+            p_rgh = p - \rho g h
+        \f]
+
+    where
+    \vartable
+        p_rgh   | Pseudo hydrostatic pressure [Pa]
+        p       | Static pressure [Pa]
+        h       | Height in the opposite direction to gravity
+        \rho    | density
+        g       | acceleration due to gravity [m/s2]
+    \endtable
+
+    \heading Patch usage
+
+    \table
+        Property     | Description             | Required    | Default value
+        rhoName      | rho field name          | no          | rho
+        p            | static pressure         | yes         |
+    \endtable
+
+    Example of the boundary condition specification:
+    \verbatim
+    myPatch
+    {
+        type            prghPressure;
+        rhoName         rho;
+        p               uniform 0;
+        value           uniform 0; // optional initial value
+    }
+    \endverbatim
+
+SeeAlso
+    Foam::fixedValueFvPatchScalarField
+
+SourceFiles
+    prghPressureFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef prghPressureFvPatchScalarField_H
+#define prghPressureFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+         Class prghPressureFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class prghPressureFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+
+protected:
+
+    // Protected data
+
+        //- Name of phase-fraction field
+        word rhoName_;
+
+        //- Static pressure
+        scalarField p_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("prghPressure");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        prghPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        prghPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  prghPressureFvPatchScalarField onto a new patch
+        prghPressureFvPatchScalarField
+        (
+            const prghPressureFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        prghPressureFvPatchScalarField
+        (
+            const prghPressureFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField >
+            (
+                new prghPressureFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        prghPressureFvPatchScalarField
+        (
+            const prghPressureFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new prghPressureFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Access
+
+            //- Return the rhoName
+            const word& rhoName() const
+            {
+                return rhoName_;
+            }
+
+            //- Return reference to the rhoName to allow adjustment
+            word& rhoName()
+            {
+                return rhoName_;
+            }
+
+            //- Return the static pressure
+            const scalarField& p() const
+            {
+                return p_;
+            }
+
+            //- Return reference to the static pressure to allow adjustment
+            scalarField& p()
+            {
+                return p_;
+            }
+
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
+
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C

@@ -120,7 +120,11 @@ timeVaryingMappedFixedValueFvPatchField
     }
     else
     {
-        updateCoeffs();
+        // Note: we use evaluate() here to trigger updateCoeffs followed
+        //       by re-setting of fvatchfield::updated_ flag. This is
+        //       so if first use is in the next time step it retriggers
+        //       a new update.
+        this->evaluate(Pstream::blocking);
     }
 }
 

src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.H

@@ -33,47 +33,47 @@ Description
 
     1. incompressible subsonic:
         \f[
-            p_T = p_0 + 0.5 |U|^2
+            p_p = p_0 - 0.5 |U|^2
         \f]
         where
         \vartable
-            p_T     | incompressible total pressure [m2/s2]
+            p_p     | incompressible pressure at patch [m2/s2]
-            p_0     | incompressible reference pressure [m2/s2]
+            p_0     | incompressible total pressure [m2/s2]
             U       | velocity
         \endvartable
 
     2. compressible subsonic:
         \f[
-            p_T = p_0 + 0.5 \rho |U|^2
+            p_p = p_0 - 0.5 \rho |U|^2
         \f]
         where
         \vartable
-            p_T     | total pressure [Pa]
+            p_p     | pressure at patch [Pa]
-            p_0     | reference pressure [Pa]
+            p_0     | total pressure [Pa]
             \rho    | density [kg/m3]
             U       | velocity
         \endvartable
 
     3. compressible transonic (\gamma <= 1):
         \f[
-            p_T = \frac{p_0}{1 + 0.5 \psi |U|^2}
+            p_p = \frac{p_0}{1 + 0.5 \psi |U|^2}
         \f]
         where
         \vartable
-            p_T     | total pressure [Pa]
+            p_p     | pressure at patch [Pa]
-            p_0     | reference pressure [Pa]
+            p_0     | total pressure [Pa]
             G       | coefficient given by \f$\frac{\gamma}{1-\gamma}\f$
         \endvartable
 
     4. compressible supersonic (\gamma > 1):
         \f[
-            p_T = \frac{p_0}{(1 + 0.5 \psi G |U|^2)^{\frac{1}{G}}}
+            p_p = \frac{p_0}{(1 + 0.5 \psi G |U|^2)^{\frac{1}{G}}}
         \f]
         where
         \vartable
+            p_p     | pressure at patch [Pa]
+            p_0     | total pressure [Pa]
             \gamma  | ratio of specific heats (Cp/Cv)
-            p_T     | total pressure [Pa]
-            p_0     | reference pressure [Pa]
             \psi    | compressibility [m2/s2]
             G       | coefficient given by \f$\frac{\gamma}{1-\gamma}\f$
         \endvartable
@@ -98,7 +98,7 @@ Description
         rho          | density field name      | no          | none
         psi          | compressibility field name | no       | none
         gamma        | ratio of specific heats (Cp/Cv) | yes |
-        p0           | static pressure reference | yes       |
+        p0           | total pressure          | yes       |
     \endtable
 
     Example of the boundary condition specification:

src/finiteVolume/finiteVolume/fvc/fvcSimpleReconstruct.C

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "fvcReconstruct.H"
+#include "fvMesh.H"
+#include "zeroGradientFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fvc
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+tmp
+<
+    GeometricField
+    <
+        typename outerProduct<vector,Type>::type, fvPatchField, volMesh
+    >
+>
+reconstruct
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
+)
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+
+    const fvMesh& mesh = ssf.mesh();
+
+    const labelUList& owner = mesh.owner();
+    const labelUList& neighbour = mesh.neighbour();
+
+    const volVectorField& C = mesh.C();
+    const surfaceVectorField& Cf = mesh.Cf();
+
+    tmp<GeometricField<GradType, fvPatchField, volMesh> > treconField
+    (
+        new GeometricField<GradType, fvPatchField, volMesh>
+        (
+            IOobject
+            (
+                "reconstruct("+ssf.name()+')',
+                ssf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensioned<GradType>
+            (
+                "0",
+                ssf.dimensions()/dimArea,
+                pTraits<GradType>::zero
+            ),
+            zeroGradientFvPatchField<GradType>::typeName
+        )
+    );
+
+    Field<GradType>& rf = treconField();
+
+    forAll(owner, facei)
+    {
+        label own = owner[facei];
+        label nei = neighbour[facei];
+
+        rf[own] += (Cf[facei] - C[own])*ssf[facei];
+        rf[nei] -= (Cf[facei] - C[nei])*ssf[facei];
+    }
+
+    const typename GeometricField<Type, fvsPatchField, surfaceMesh>::
+    GeometricBoundaryField& bsf = ssf.boundaryField();
+
+    forAll(bsf, patchi)
+    {
+        const fvsPatchField<Type>& psf = bsf[patchi];
+
+        const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
+        const vectorField& pCf = Cf.boundaryField()[patchi];
+
+        forAll(pOwner, pFacei)
+        {
+            label own = pOwner[pFacei];
+            rf[own] += (pCf[pFacei] - C[own])*psf[pFacei];
+        }
+    }
+
+    rf /= mesh.V();
+
+    treconField().correctBoundaryConditions();
+
+    return treconField;
+}
+
+
+template<class Type>
+tmp
+<
+    GeometricField
+    <
+        typename outerProduct<vector, Type>::type, fvPatchField, volMesh
+    >
+>
+reconstruct
+(
+    const tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >& tssf
+)
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+    tmp<GeometricField<GradType, fvPatchField, volMesh> > tvf
+    (
+        fvc::reconstruct(tssf())
+    );
+    tssf.clear();
+    return tvf;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fvc
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/invDistLeastSquaresVectors.C

@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "leastSquaresVectors.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(leastSquaresVectors, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::leastSquaresVectors::leastSquaresVectors(const fvMesh& mesh)
+:
+    MeshObject<fvMesh, Foam::MoveableMeshObject, leastSquaresVectors>(mesh),
+    pVectors_
+    (
+        IOobject
+        (
+            "LeastSquaresP",
+            mesh_.pointsInstance(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh_,
+        dimensionedVector("zero", dimless/dimLength, vector::zero)
+    ),
+    nVectors_
+    (
+        IOobject
+        (
+            "LeastSquaresN",
+            mesh_.pointsInstance(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh_,
+        dimensionedVector("zero", dimless/dimLength, vector::zero)
+    )
+{
+    calcLeastSquaresVectors();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::leastSquaresVectors::~leastSquaresVectors()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::leastSquaresVectors::calcLeastSquaresVectors()
+{
+    if (debug)
+    {
+        Info<< "leastSquaresVectors::calcLeastSquaresVectors() :"
+            << "Calculating least square gradient vectors"
+            << endl;
+    }
+
+    const fvMesh& mesh = mesh_;
+
+    // Set local references to mesh data
+    const labelUList& owner = mesh_.owner();
+    const labelUList& neighbour = mesh_.neighbour();
+
+    const volVectorField& C = mesh.C();
+
+    // Set up temporary storage for the dd tensor (before inversion)
+    symmTensorField dd(mesh_.nCells(), symmTensor::zero);
+
+    forAll(owner, facei)
+    {
+        label own = owner[facei];
+        label nei = neighbour[facei];
+
+        vector d = C[nei] - C[own];
+        symmTensor wdd = sqr(d)/magSqr(d);
+        dd[own] += wdd;
+        dd[nei] += wdd;
+    }
+
+
+    surfaceVectorField::GeometricBoundaryField& blsP =
+        pVectors_.boundaryField();
+
+    forAll(blsP, patchi)
+    {
+        const fvsPatchVectorField& patchLsP = blsP[patchi];
+
+        const fvPatch& p = patchLsP.patch();
+        const labelUList& faceCells = p.patch().faceCells();
+
+        // Build the d-vectors
+        vectorField pd(p.delta());
+
+        forAll(pd, patchFacei)
+        {
+            const vector& d = pd[patchFacei];
+
+            dd[faceCells[patchFacei]] += sqr(d)/magSqr(d);
+        }
+    }
+
+
+    // Invert the dd tensor
+    const symmTensorField invDd(inv(dd));
+
+
+    // Revisit all faces and calculate the pVectors_ and nVectors_ vectors
+    forAll(owner, facei)
+    {
+        label own = owner[facei];
+        label nei = neighbour[facei];
+
+        vector d = C[nei] - C[own];
+
+        pVectors_[facei] = (invDd[own] & d)/magSqr(d);
+        nVectors_[facei] = -(invDd[nei] & d)/magSqr(d);
+    }
+
+    forAll(blsP, patchi)
+    {
+        fvsPatchVectorField& patchLsP = blsP[patchi];
+
+        const fvPatch& p = patchLsP.patch();
+        const labelUList& faceCells = p.faceCells();
+
+        // Build the d-vectors
+        vectorField pd(p.delta());
+
+        forAll(pd, patchFacei)
+        {
+            const vector& d = pd[patchFacei];
+
+            patchLsP[patchFacei] = (invDd[faceCells[patchFacei]] & d)/magSqr(d);
+        }
+    }
+
+    if (debug)
+    {
+        Info<< "leastSquaresVectors::calcLeastSquaresVectors() :"
+            << "Finished calculating least square gradient vectors"
+            << endl;
+    }
+}
+
+
+bool Foam::leastSquaresVectors::movePoints()
+{
+    calcLeastSquaresVectors();
+    return true;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/unweightedLeastSquaresVectors.C

@@ -0,0 +1,183 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "leastSquaresVectors.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(leastSquaresVectors, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::leastSquaresVectors::leastSquaresVectors(const fvMesh& mesh)
+:
+    MeshObject<fvMesh, Foam::MoveableMeshObject, leastSquaresVectors>(mesh),
+    pVectors_
+    (
+        IOobject
+        (
+            "LeastSquaresP",
+            mesh_.pointsInstance(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh_,
+        dimensionedVector("zero", dimless/dimLength, vector::zero)
+    ),
+    nVectors_
+    (
+        IOobject
+        (
+            "LeastSquaresN",
+            mesh_.pointsInstance(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh_,
+        dimensionedVector("zero", dimless/dimLength, vector::zero)
+    )
+{
+    calcLeastSquaresVectors();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::leastSquaresVectors::~leastSquaresVectors()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::leastSquaresVectors::calcLeastSquaresVectors()
+{
+    if (debug)
+    {
+        Info<< "leastSquaresVectors::calcLeastSquaresVectors() :"
+            << "Calculating least square gradient vectors"
+            << endl;
+    }
+
+    const fvMesh& mesh = mesh_;
+
+    // Set local references to mesh data
+    const labelUList& owner = mesh_.owner();
+    const labelUList& neighbour = mesh_.neighbour();
+
+    const volVectorField& C = mesh.C();
+
+    // Set up temporary storage for the dd tensor (before inversion)
+    symmTensorField dd(mesh_.nCells(), symmTensor::zero);
+
+    forAll(owner, facei)
+    {
+        label own = owner[facei];
+        label nei = neighbour[facei];
+
+        symmTensor wdd = sqr(C[nei] - C[own]);
+        dd[own] += wdd;
+        dd[nei] += wdd;
+    }
+
+
+    surfaceVectorField::GeometricBoundaryField& blsP =
+        pVectors_.boundaryField();
+
+    forAll(blsP, patchi)
+    {
+        const fvsPatchVectorField& patchLsP = blsP[patchi];
+
+        const fvPatch& p = patchLsP.patch();
+        const labelUList& faceCells = p.patch().faceCells();
+
+        // Build the d-vectors
+        vectorField pd(p.delta());
+
+        forAll(pd, patchFacei)
+        {
+            dd[faceCells[patchFacei]] += sqr(pd[patchFacei]);
+        }
+    }
+
+
+    // Invert the dd tensor
+    const symmTensorField invDd(inv(dd));
+
+
+    // Revisit all faces and calculate the pVectors_ and nVectors_ vectors
+    forAll(owner, facei)
+    {
+        label own = owner[facei];
+        label nei = neighbour[facei];
+
+        vector d = C[nei] - C[own];
+
+        pVectors_[facei] = (invDd[own] & d);
+        nVectors_[facei] = -(invDd[nei] & d);
+    }
+
+    forAll(blsP, patchi)
+    {
+        fvsPatchVectorField& patchLsP = blsP[patchi];
+
+        const fvPatch& p = patchLsP.patch();
+        const labelUList& faceCells = p.faceCells();
+
+        // Build the d-vectors
+        vectorField pd(p.delta());
+
+        forAll(pd, patchFacei)
+        {
+            patchLsP[patchFacei] =
+                (invDd[faceCells[patchFacei]] & pd[patchFacei]);
+        }
+    }
+
+    if (debug)
+    {
+        Info<< "leastSquaresVectors::calcLeastSquaresVectors() :"
+            << "Finished calculating least square gradient vectors"
+            << endl;
+    }
+}
+
+
+bool Foam::leastSquaresVectors::movePoints()
+{
+    calcLeastSquaresVectors();
+    return true;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C

@@ -30,6 +30,7 @@ License
 void Foam::MULES::correct
 (
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiPsiCorr,
     const scalar psiMax,
     const scalar psiMin
@@ -39,6 +40,7 @@ void Foam::MULES::correct
     (
         geometricOneField(),
         psi,
+        phi,
         phiPsiCorr,
         zeroField(), zeroField(),
         psiMax, psiMin
@@ -49,6 +51,7 @@ void Foam::MULES::correct
 void Foam::MULES::LTScorrect
 (
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiPsiCorr,
     const scalar psiMax,
     const scalar psiMin
@@ -58,6 +61,7 @@ void Foam::MULES::LTScorrect
     (
         geometricOneField(),
         psi,
+        phi,
         phiPsiCorr,
         zeroField(), zeroField(),
         psiMax, psiMin

src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H

@@ -66,6 +66,7 @@ void correct
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su
@@ -76,6 +77,7 @@ void correct
 (
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -86,6 +88,7 @@ void correct
 void correct
 (
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const scalar psiMax,
     const scalar psiMin
@@ -96,6 +99,7 @@ void LTScorrect
 (
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -106,6 +110,7 @@ void LTScorrect
 void LTScorrect
 (
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const scalar psiMax,
     const scalar psiMin
@@ -119,6 +124,7 @@ void limiterCorr
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
+    const surfaceScalarField& phi,
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -133,6 +139,7 @@ void limitCorr
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,

src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C

@@ -37,6 +37,7 @@ void Foam::MULES::correct
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su
@@ -77,6 +78,7 @@ void Foam::MULES::correct
 (
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -94,8 +96,18 @@ void Foam::MULES::correct
         readLabel(MULEScontrols.lookup("nLimiterIter"))
     );
 
-    limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
+    limitCorr
-    correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
+    (
+        rDeltaT,
+        rho,
+        psi,
+        phi,
+        phiCorr,
+        Sp, Su,
+        psiMax, psiMin,
+        nLimiterIter
+    );
+    correct(rDeltaT, rho, psi, phi, phiCorr, Sp, Su);
 }
 
 
@@ -104,6 +116,7 @@ void Foam::MULES::LTScorrect
 (
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -123,8 +136,18 @@ void Foam::MULES::LTScorrect
         readLabel(MULEScontrols.lookup("nLimiterIter"))
     );
 
-    limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
+    limitCorr
-    correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
+    (
+        rDeltaT,
+        rho,
+        psi,
+        phi,
+        phiCorr,
+        Sp, Su,
+        psiMax, psiMin,
+        nLimiterIter
+    );
+    correct(rDeltaT, rho, psi, phi, phiCorr, Sp, Su);
 }
 
 
@@ -135,6 +158,7 @@ void Foam::MULES::limiterCorr
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
+    const surfaceScalarField& phi,
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -153,6 +177,9 @@ void Foam::MULES::limiterCorr
     tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
     const scalarField& V = tVsc();
 
+    const surfaceScalarField::GeometricBoundaryField& phiBf =
+        phi.boundaryField();
+
     const scalarField& phiCorrIf = phiCorr;
     const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
         phiCorr.boundaryField();
@@ -408,12 +435,12 @@ void Foam::MULES::limiterCorr
             {
                 const labelList& pFaceCells =
                     mesh.boundary()[patchi].faceCells();
-                const scalarField& phiCorrPf = phiCorrBf[patchi];
+                const scalarField& phiPf = phiBf[patchi];
 
                 forAll(lambdaPf, pFacei)
                 {
                     // Limit outlet faces only
-                    if (phiCorrPf[pFacei] > SMALL*SMALL)
+                    if (phiPf[pFacei] > SMALL*SMALL)
                     {
                         label pfCelli = pFaceCells[pFacei];
 
@@ -443,6 +470,7 @@ void Foam::MULES::limitCorr
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -478,6 +506,7 @@ void Foam::MULES::limitCorr
         rDeltaT,
         rho,
         psi,
+        phi,
         phiCorr,
         Sp,
         Su,

src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C

@@ -103,7 +103,22 @@ timeVaryingMappedFixedValuePointPatchField
     endAverage_(pTraits<Type>::zero)
 {
     dict.readIfPresent("fieldTableName", fieldTableName_);
-    updateCoeffs();
+
+    if (dict.found("value"))
+    {
+        fixedValuePointPatchField<Type>::operator==
+        (
+            Field<Type>("value", dict, p.size())
+        );
+    }
+    else
+    {
+        // Note: use evaluate to do updateCoeffs followed by a reset
+        //       of the pointPatchField::updated_ flag. This is
+        //       so if first use is in the next time step it retriggers
+        //       a new update.
+        pointPatchField<Type>::evaluate(Pstream::blocking);
+    }
 }
 
 

src/lagrangian/basic/particle/particleTemplates.C

@@ -424,7 +424,7 @@ Foam::scalar Foam::particle::trackToFace
                 (
                     position_,
                     endPosition,
-                    triI,
+                    tI,
                     tetAreas[tI],
                     tetPlaneBasePtIs[tI],
                     cellI_,

src/mesh/autoMesh/Make/files

@@ -26,5 +26,11 @@ $(autoHexMesh)/shellSurfaces/shellSurfaces.C
 $(autoHexMesh)/trackedParticle/trackedParticle.C
 $(autoHexMesh)/trackedParticle/trackedParticleCloud.C
 
+meshMover = $(autoHexMesh)/externalDisplacementMeshMover
+$(meshMover)/displacementMeshMoverMotionSolver.C
+$(meshMover)/externalDisplacementMeshMover.C
+$(meshMover)/medialAxisMeshMover.C
+
+
 LIB = $(FOAM_LIBBIN)/libautoMesh
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -53,6 +53,7 @@ Description
 #include "calculatedPointPatchFields.H"
 #include "cyclicSlipPointPatchFields.H"
 #include "fixedValueFvPatchFields.H"
+#include "localPointRegion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -833,7 +834,7 @@ Foam::autoLayerDriver::makeLayerDisplacementField
 
     forAll(numLayers, patchI)
     {
-        //  0 layers: do not allow lslip so fixedValue 0
+        //  0 layers: do not allow slip so fixedValue 0
         // >0 layers: fixedValue which gets adapted
         if (numLayers[patchI] >= 0)
         {
@@ -2809,14 +2810,7 @@ void Foam::autoLayerDriver::mergePatchFacesUndo
 
     const fvMesh& mesh = meshRefiner_.mesh();
 
-    List<labelPair> couples
+    List<labelPair> couples(localPointRegion::findDuplicateFacePairs(mesh));
-    (
-        meshRefiner_.getDuplicateFaces   // get all baffles
-        (
-            identity(mesh.nFaces()-mesh.nInternalFaces())
-          + mesh.nInternalFaces()
-        )
-    );
 
     labelList duplicateFace(mesh.nFaces(), -1);
     forAll(couples, i)
@@ -3522,7 +3516,9 @@ void Foam::autoLayerDriver::addLayers
         }
 
         // Reset mesh points and start again
-        meshMover().movePoints(oldPoints);
+        mesh.movePoints(oldPoints);
+        // Update meshmover for change of mesh geometry
+        meshMover().movePoints();
         meshMover().correct();
 
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H

@@ -58,7 +58,9 @@ class layerParameters;
 
 class autoLayerDriver
 {
-    // Static data members
+public:
+
+    // Public data types
 
         //- Extrusion controls
         enum extrudeMode
@@ -69,6 +71,7 @@ class autoLayerDriver
                             /*!< faces locally */
         };
 
+private:
 
     // Private classes
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C

@@ -1756,19 +1756,10 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
         Info<< "Writing wanted-displacement mesh (possibly illegal) to "
             << meshRefiner_.timeName() << endl;
         pointField oldPoints(mesh.points());
-        vectorField totalDisp
+
-        (
+        meshRefiner_.mesh().movePoints(meshMover.curPoints());
-            meshMover.scale().internalField()
+        // Warn meshMover for changed geometry
-          * displacement.internalField()
+        meshMover.movePoints();
-        );
-        syncTools::syncPointList
-        (
-            mesh,
-            totalDisp,
-            minMagSqrEqOp<point>(),
-            vector(GREAT, GREAT, GREAT)
-        );
-        meshMover.movePoints((mesh.points()+totalDisp)());
 
         // Above move will have changed the instance only on the points (which
         // is correct).
@@ -1791,7 +1782,11 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
         dispVec.write();
         medialDist.write();
         medialRatio.write();
-        meshMover.movePoints(oldPoints);
+
+        // Move mesh back
+        meshRefiner_.mesh().movePoints(oldPoints);
+        // Warn meshMover for changed geometry
+        meshMover.movePoints();
     }
 
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C

@@ -36,6 +36,7 @@ License
 #include "mapDistributePolyMesh.H"
 #include "unitConversion.H"
 #include "snapParameters.H"
+#include "localPointRegion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1036,19 +1037,11 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
 
 
     // Merge all baffles that are still remaining after duplicating points.
-    List<labelPair> couples
+    List<labelPair> couples(localPointRegion::findDuplicateFacePairs(mesh));
-    (
-        meshRefiner_.getDuplicateFaces   // get all baffles
-        (
-            identity(mesh.nFaces()-mesh.nInternalFaces())
-          + mesh.nInternalFaces()
-        )
-    );
 
     label nCouples = returnReduce(couples.size(), sumOp<label>());
 
-    Info<< "Detected unsplittable baffles : "
+    Info<< "Detected unsplittable baffles : " << nCouples << endl;
-        << nCouples << endl;
 
     if (nCouples > 0)
     {

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMeshMoverMotionSolver.C

@@ -0,0 +1,136 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "displacementMeshMoverMotionSolver.H"
+#include "addToRunTimeSelectionTable.H"
+#include "localPointRegion.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(displacementMeshMoverMotionSolver, 0);
+
+    addToRunTimeSelectionTable
+    (
+        motionSolver,
+        displacementMeshMoverMotionSolver,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::displacementMeshMoverMotionSolver::displacementMeshMoverMotionSolver
+(
+    const polyMesh& mesh,
+    const IOdictionary& dict
+)
+:
+    displacementMotionSolver(mesh, dict, typeName) // read pointDisplacement
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::displacementMeshMoverMotionSolver::
+~displacementMeshMoverMotionSolver()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+Foam::externalDisplacementMeshMover&
+Foam::displacementMeshMoverMotionSolver::meshMover() const
+{
+    if (!meshMoverPtr_.valid())
+    {
+        const word moverType(coeffDict().lookup("meshMover"));
+
+        meshMoverPtr_ = externalDisplacementMeshMover::New
+        (
+            moverType,
+            coeffDict().subDict(moverType + "Coeffs"),
+            localPointRegion::findDuplicateFacePairs(mesh()),
+            pointDisplacement_
+        );
+    }
+    return meshMoverPtr_();
+}
+
+
+Foam::tmp<Foam::pointField>
+Foam::displacementMeshMoverMotionSolver::curPoints() const
+{
+    // Return actual points. Cannot do a reference since complains about
+    // assignment to self in polyMesh::movePoints
+    return tmp<pointField>(new pointField(mesh().points()));
+}
+
+
+void Foam::displacementMeshMoverMotionSolver::solve()
+{
+    // The points have moved so before calculation update
+    // the mesh and motionSolver accordingly
+    movePoints(mesh().points());
+
+    // Update any point motion bcs (e.g. timevarying)
+    pointDisplacement().boundaryField().updateCoeffs();
+
+    label nAllowableErrors = 0;
+    meshMover().move(nAllowableErrors);
+}
+
+
+void Foam::displacementMeshMoverMotionSolver::movePoints(const pointField& p)
+{
+    displacementMotionSolver::movePoints(p);
+
+    // Update meshMover for new geometry
+    if (meshMoverPtr_.valid())
+    {
+        meshMover().movePoints(p);
+    }
+}
+
+
+void Foam::displacementMeshMoverMotionSolver::updateMesh
+(
+    const mapPolyMesh& map
+)
+{
+    displacementMotionSolver::updateMesh(map);
+
+    // Update meshMover for new topology
+    meshMoverPtr_.clear();
+}
+
+
+// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMeshMoverMotionSolver.H

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::displacementMeshMoverMotionSolver
+
+Description
+    Mesh motion solver for an fvMesh.  Based on solving the cell-centre
+    Laplacian for the motion displacement.
+
+SourceFiles
+    displacementMeshMoverMotionSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef displacementMeshMoverMotionSolver_H
+#define displacementMeshMoverMotionSolver_H
+
+#include "displacementMotionSolver.H"
+#include "externalDisplacementMeshMover.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+             Class displacementMeshMoverMotionSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class displacementMeshMoverMotionSolver
+:
+    public displacementMotionSolver
+{
+    // Private data
+
+        mutable autoPtr<externalDisplacementMeshMover> meshMoverPtr_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        displacementMeshMoverMotionSolver
+        (
+            const displacementMeshMoverMotionSolver&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const displacementMeshMoverMotionSolver&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("displacementMeshMover");
+
+
+    // Constructors
+
+        //- Construct from polyMesh and IOdictionary
+        displacementMeshMoverMotionSolver(const polyMesh&, const IOdictionary&);
+
+
+    //- Destructor
+    ~displacementMeshMoverMotionSolver();
+
+
+    // Member Functions
+
+        externalDisplacementMeshMover& meshMover() const;
+
+        //- Return point location obtained from the current motion field
+        virtual tmp<pointField> curPoints() const;
+
+        //- Solve for motion
+        virtual void solve();
+
+        //- Update local data for geometry changes
+        virtual void movePoints(const pointField&);
+
+        //- Update topology
+        virtual void updateMesh(const mapPolyMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "externalDisplacementMeshMover.H"
+#include "mapPolyMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(externalDisplacementMeshMover, 0);
+    defineRunTimeSelectionTable(externalDisplacementMeshMover, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::externalDisplacementMeshMover::externalDisplacementMeshMover
+(
+    const dictionary& dict,
+    const List<labelPair>& baffles,
+    pointVectorField& pointDisplacement
+)
+:
+    dict_(dict),
+    baffles_(baffles),
+    pointDisplacement_(pointDisplacement)
+{}
+
+
+
+// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::externalDisplacementMeshMover>
+Foam::externalDisplacementMeshMover::New
+(
+    const word& type,
+    const dictionary& dict,
+    const List<labelPair>& baffles,
+    pointVectorField& pointDisplacement
+)
+{
+    Info<< "Selecting externalDisplacementMeshMover " << type << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(type);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "externalDisplacementMeshMover::New(const word&"
+            ", pointVectorField&, const List<labelPair>&"
+            ", const dictionary&)"
+        )   << "Unknown externalDisplacementMeshMover type "
+            << type << nl << nl
+            << "Valid externalDisplacementMeshMover types:" << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<externalDisplacementMeshMover>
+    (
+        cstrIter()(dict, baffles, pointDisplacement)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::externalDisplacementMeshMover::~externalDisplacementMeshMover()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::externalDisplacementMeshMover::movePoints(const pointField&)
+{
+    // No local data to update
+}
+
+
+void Foam::externalDisplacementMeshMover::updateMesh(const mapPolyMesh& mpm)
+{
+    // Renumber baffles
+    DynamicList<labelPair> newBaffles(baffles_.size());
+    forAll(baffles_, i)
+    {
+        label f0 = mpm.reverseFaceMap()[baffles_[i].first()];
+        label f1 = mpm.reverseFaceMap()[baffles_[i].second()];
+
+        if (f0 >= 0 && f1 >= 0)
+        {
+            newBaffles.append(labelPair(f0, f1));
+        }
+    }
+    newBaffles.shrink();
+    baffles_.transfer(newBaffles);
+}
+
+
+// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.H

@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::externalDisplacementMeshMover
+
+Description
+    Virtual base class for mesh movers with externally provided displacement
+    field giving the boundary conditions. Move the mesh from the current
+    location to a new location (so modify the mesh; v.s. motionSolver that
+    only returns the new location)
+
+SourceFiles
+    externalDisplacementMeshMover.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef externalDisplacementMeshMover_H
+#define externalDisplacementMeshMover_H
+
+#include "pointFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class mapPolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                   Class externalDisplacementMeshMover Declaration
+\*---------------------------------------------------------------------------*/
+
+class externalDisplacementMeshMover
+{
+protected:
+
+    // Protected data
+
+        //- Settings
+        dictionary dict_;
+
+        //- Baffles in the mesh
+        List<labelPair> baffles_;
+
+        //- Reference to point motion field
+        pointVectorField& pointDisplacement_;
+
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        externalDisplacementMeshMover
+        (
+            const externalDisplacementMeshMover&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const externalDisplacementMeshMover&);
+
+public:
+
+    //- Runtime type information
+    TypeName("externalDisplacementMeshMover");
+
+
+    // Declare run-time New selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            externalDisplacementMeshMover,
+            dictionary,
+            (
+                const dictionary& dict,
+                const List<labelPair>& baffles,
+                pointVectorField& pointDisplacement
+            ),
+            (dict, baffles, pointDisplacement)
+        );
+
+
+    // Constructors
+
+        //- Construct from dictionary and displacement field
+        externalDisplacementMeshMover
+        (
+            const dictionary& dict,
+            const List<labelPair>& baffles,
+            pointVectorField& pointDisplacement
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected meshMover model
+        static autoPtr<externalDisplacementMeshMover> New
+        (
+            const word& type,
+            const dictionary& dict,
+            const List<labelPair>& baffles,
+            pointVectorField& pointDisplacement
+        );
+
+
+    //- Destructor
+    virtual ~externalDisplacementMeshMover();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return reference to the point motion displacement field
+            pointVectorField& pointDisplacement()
+            {
+                return pointDisplacement_;
+            }
+
+            //- Return const reference to the point motion displacement field
+            const pointVectorField& pointDisplacement() const
+            {
+                return pointDisplacement_;
+            }
+
+            const pointMesh& pMesh() const
+            {
+                return pointDisplacement_.mesh();
+            }
+
+            const polyMesh& mesh() const
+            {
+                return pMesh()();
+            }
+
+
+        // Mesh mover
+
+            //- Move mesh. Return true if succesful
+            virtual bool move(const label nAllowableErrors) = 0;
+
+            //- Update local data for geometry changes
+            virtual void movePoints(const pointField&);
+
+            //-  Update local data for topology changes
+            virtual void updateMesh(const mapPolyMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H

@@ -0,0 +1,355 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::medialAxisMeshMover
+
+Description
+    Mesh motion solver that uses a medial axis algorithm to work
+    out a fraction between the (nearest point on a) moving surface
+    and the (nearest point on a) fixed surface.
+    This fraction is then used to scale the motion.
+
+SourceFiles
+    medialAxisMeshMover.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef medialAxisMeshMover_H
+#define medialAxisMeshMover_H
+
+#include "externalDisplacementMeshMover.H"
+#include "motionSmootherAlgo.H"
+#include "autoLayerDriver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class pointData;
+
+/*---------------------------------------------------------------------------*\
+             Class medialAxisMeshMover Declaration
+\*---------------------------------------------------------------------------*/
+
+class medialAxisMeshMover
+:
+    public externalDisplacementMeshMover
+{
+    // Private data
+
+        const labelList adaptPatchIDs_;
+
+        autoPtr<indirectPrimitivePatch> adaptPatchPtr_;
+
+        //- Scale factor for displacement
+        pointScalarField scale_;
+
+        //- Starting mesh position
+        pointField oldPoints_;
+
+        //- Mesh mover algorithm
+        motionSmootherAlgo meshMover_;
+
+
+    // Settings
+
+        //- Smoothing wanted patch thickness
+        label nSmoothPatchThickness_;
+
+        //- (lambda-mu) smoothing of internal displacement
+        label nSmoothDisplacement_;
+
+        //- Layer thickness too big
+        scalar maxThicknessToMedialRatio_;
+
+        //- Feature angle when to stop adding layers
+        scalar featureAngle_;
+
+        //- Stop layer growth where mesh wraps around sharp edge
+        scalar minCosLayerTermination_;
+
+        //- Smooth internal normals
+        label nSmoothNormals_;
+
+        //- Smooth surface normals
+        label nSmoothSurfaceNormals_;
+
+        //- When is medial axis
+        scalar minMedianAxisAngleCos_;
+
+        //- When to slip along wall
+        scalar slipFeatureAngle_;
+
+        //- Number of attempts shrinking the mesh
+        label nSnap_;
+
+
+    // Pre-calculated medial axis information
+
+        //- normal of nearest wall. Where this normal changes direction
+        //  defines the medial axis
+        pointVectorField dispVec_;
+
+        //- ratio of medial distance to wall distance.
+        //  (1 at wall, 0 at medial axis)
+        pointScalarField medialRatio_;
+
+        //- Distance to nearest medial axis point
+        pointScalarField medialDist_;
+
+        //- Location on nearest medial axis point
+        pointVectorField medialVec_;
+
+
+    // Private Member Functions
+
+        const dictionary& coeffDict() const;
+
+        const dictionary& meshQualityDict() const;
+
+        //- Read model coefficients
+        void read();
+
+        //- Extract fixed-value patchfields
+        static labelList getFixedValueBCs(const pointVectorField&);
+
+        //- Extract bc types. Replace fixedValue derivatives with fixedValue
+        wordList getPatchFieldTypes(const pointVectorField& fld);
+
+        //- Construct patch on selected patches
+        static autoPtr<indirectPrimitivePatch> getPatch
+        (
+            const polyMesh&,
+            const labelList&
+        );
+
+        //- Weighted sum (over all subset of mesh points) by
+        //  summing contribution from (master) edges
+        template<class Type>
+        static void weightedSum
+        (
+            const polyMesh& mesh,
+            const PackedBoolList& isMasterEdge,
+            const labelList& meshEdges,
+            const labelList& meshPoints,
+            const edgeList& edges,
+            const scalarField& edgeWeights,
+            const Field<Type>& data,
+            Field<Type>& sum
+        );
+
+        void calculateEdgeWeights
+        (
+            const PackedBoolList& isMasterEdge,
+            const labelList& meshEdges,
+            const labelList& meshPoints,
+            const edgeList& edges,
+            scalarField& edgeWeights,
+            scalarField& invSumWeight
+        ) const;
+
+
+        // Calculation of medial axis information
+
+            //- Smooth normals on patch
+            void smoothPatchNormals
+            (
+                const label nSmoothDisp,
+                const PackedBoolList& isMasterPoint,
+                const PackedBoolList& isMasterEdge,
+                const labelList& meshEdges,
+                pointField& normals
+            ) const;
+
+            //- Smooth normals on interior
+            void smoothNormals
+            (
+                const label nSmoothDisp,
+                const PackedBoolList& isMasterPoint,
+                const PackedBoolList& isMasterEdge,
+                const labelList& fixedPoints,
+                pointVectorField& normals
+            ) const;
+
+            //- Is mesh edge on a cusp of displacement
+            bool isMaxEdge
+            (
+                const List<pointData>& pointWallDist,
+                const label edgeI,
+                const scalar minCos
+            ) const;
+
+            //- Initialise medial axis
+            void update();
+
+
+        // Calculation of mesh movement
+
+            //- Unmark extrusion at given point
+            static bool unmarkExtrusion
+            (
+                const label patchPointI,
+                pointField& patchDisp,
+                List<autoLayerDriver::extrudeMode>& extrudeStatus
+            );
+
+            //- Synchronise extrusion
+            void syncPatchDisplacement
+            (
+                const scalarField& minThickness,
+                pointField& patchDisp,
+                List<autoLayerDriver::extrudeMode>& extrudeStatus
+            ) const;
+
+            void smoothLambdaMuDisplacement
+            (
+                const label nSmoothDisp,
+                const PackedBoolList& isMasterPoint,
+                const PackedBoolList& isMasterEdge,
+                vectorField& displacement
+            ) const;
+
+            void minSmoothField
+            (
+                const label nSmoothDisp,
+                const PackedBoolList& isMasterPoint,
+                const PackedBoolList& isMasterEdge,
+                const labelList& meshEdges,
+                const scalarField& fieldMin,
+                scalarField& field
+            ) const;
+
+            //- Stop layer growth at feature edges
+            void handleFeatureAngleLayerTerminations
+            (
+                const scalar minCos,
+                const PackedBoolList& isMasterPoint,
+                const labelList& meshEdges,
+                List<autoLayerDriver::extrudeMode>& extrudeStatus,
+                pointField& patchDisp,
+                label& nPointCounter
+            ) const;
+
+            //- Find isolated islands (points, edges and faces and layer
+            //  terminations) in the layer mesh and stop any layer growth
+            //  at these points
+            void findIsolatedRegions
+            (
+                const scalar minCosLayerTermination,
+                const PackedBoolList& isMasterPoint,
+                const PackedBoolList& isMasterEdge,
+                const labelList& meshEdges,
+                const scalarField& minThickness,
+                List<autoLayerDriver::extrudeMode>& extrudeStatus,
+                pointField& patchDisp
+            ) const;
+
+
+        //- Calculate desired displacement. Modifies per-patch displacement
+        //  and calculates displacement for whole mesh
+        //  (in pointDisplacement)
+        void calculateDisplacement
+        (
+            const scalarField& minThickness,
+            List<autoLayerDriver::extrudeMode>& extrudeStatus,
+            pointField& patchDisp
+        );
+
+        //- Move mesh according to calculated displacement
+        bool shrinkMesh
+        (
+            const dictionary& meshQualityDict,
+            const label nAllowableErrors
+        );
+
+
+        //- Disallow default bitwise copy construct
+        medialAxisMeshMover
+        (
+            const medialAxisMeshMover&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const medialAxisMeshMover&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("displacementMedialAxis");
+
+
+    // Constructors
+
+        //- Construct from dictionary and displacement field
+        medialAxisMeshMover
+        (
+            const dictionary& dict,
+            const List<labelPair>& baffles,
+            pointVectorField& pointDisplacement
+        );
+
+
+    // Destructor
+
+        virtual ~medialAxisMeshMover();
+
+
+    // Member Functions
+
+        //- Move mesh. Return true if succesful
+        virtual bool move(const label nAllowableErrors);
+
+        //- Update local data for geometry changes
+        virtual void movePoints(const pointField&);
+
+        //-  Update local data for topology changes
+        virtual void updateMesh(const mapPolyMesh&)
+        {
+            notImplemented
+            (
+                "medialAxisMeshMover::updateMesh"
+                "(const mapPolyMesh&)"
+            );
+        }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "medialAxisMeshMoverTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMoverTemplates.C

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "medialAxisMeshMover.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::medialAxisMeshMover::weightedSum
+(
+    const polyMesh& mesh,
+    const PackedBoolList& isMasterEdge,
+    const labelList& meshEdges,
+    const labelList& meshPoints,
+    const edgeList& edges,
+    const scalarField& edgeWeights,
+    const Field<Type>& pointData,
+    Field<Type>& sum
+)
+{
+    if
+    (
+        mesh.nEdges() != isMasterEdge.size()
+     || edges.size() != meshEdges.size()
+     || edges.size() != edgeWeights.size()
+     || meshPoints.size() != pointData.size()
+    )
+    {
+        FatalErrorIn("medialAxisMeshMover::weightedSum(..)")
+            << "Inconsistent sizes for edge or point data:"
+            << " meshEdges:" << meshEdges.size()
+            << " isMasterEdge:" << isMasterEdge.size()
+            << " edgeWeights:" << edgeWeights.size()
+            << " edges:" << edges.size()
+            << " pointData:" << pointData.size()
+            << " meshPoints:" << meshPoints.size()
+            << abort(FatalError);
+    }
+
+    sum.setSize(meshPoints.size());
+    sum = pTraits<Type>::zero;
+
+    forAll(edges, edgeI)
+    {
+        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
+        {
+            const edge& e = edges[edgeI];
+
+            scalar eWeight = edgeWeights[edgeI];
+
+            label v0 = e[0];
+            label v1 = e[1];
+
+            sum[v0] += eWeight*pointData[v1];
+            sum[v1] += eWeight*pointData[v0];
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh,
+        meshPoints,
+        sum,
+        plusEqOp<Type>(),
+        pTraits<Type>::zero     // null value
+    );
+}
+
+
+// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -1519,10 +1519,9 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                 labelList coupledFace(mesh_.nFaces(), -1);
                 {
                     // Get boundary baffles that need to stay together
-                    List<labelPair> allCouples = getDuplicateFaces
+                    List<labelPair> allCouples
                     (
-                        identity(nBnd)
+                        localPointRegion::findDuplicateFacePairs(mesh_)
-                       +mesh_.nInternalFaces()
                     );
 
                     // Merge with any couples from
@@ -2229,6 +2228,36 @@ void Foam::meshRefinement::selectSeparatedCoupledFaces(boolList& selected) const
 }
 
 
+Foam::label Foam::meshRefinement::findRegion
+(
+    const polyMesh& mesh,
+    const labelList& cellToRegion,
+    const vector& perturbVec,
+    const point& p
+)
+{
+    label regionI = -1;
+    label cellI = mesh.findCell(p);
+    if (cellI != -1)
+    {
+        regionI = cellToRegion[cellI];
+    }
+    reduce(regionI, maxOp<label>());
+
+    if (regionI == -1)
+    {
+        // See if we can perturb a bit
+        cellI = mesh.findCell(p+perturbVec);
+        if (cellI != -1)
+        {
+            regionI = cellToRegion[cellI];
+        }
+        reduce(regionI, maxOp<label>());
+    }
+    return regionI;
+}
+
+
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions
 (
     const labelList& globalToMasterPatch,
@@ -2247,16 +2276,13 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions
 
     regionSplit cellRegion(mesh_, blockedFace);
 
-    label regionI = -1;
+    label regionI = findRegion
-
+    (
-    label cellI = mesh_.findCell(keepPoint);
+        mesh_,
-
+        cellRegion,
-    if (cellI != -1)
+        mergeDistance_*vector(1,1,1), // note:1,1,1 should really be normalised
-    {
+        keepPoint
-        regionI = cellRegion[cellI];
+    );
-    }
-
-    reduce(regionI, maxOp<label>());
 
     if (regionI == -1)
     {

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -853,10 +853,6 @@ public:
                 List<labelPair>&
             );
 
-            //- Return a list of coupled face pairs, i.e. faces that
-            //  use the same vertices.
-            List<labelPair> getDuplicateFaces(const labelList& testFaces) const;
-
             //- Merge baffles. Gets pairs of faces.
             autoPtr<mapPolyMesh> mergeBaffles(const List<labelPair>&);
 
@@ -894,6 +890,15 @@ public:
             //- Select coupled faces that are not collocated
             void selectSeparatedCoupledFaces(boolList&) const;
 
+            //- Find region cell is in. Uses optional perturbation to re-test.
+            static label findRegion
+            (
+                const polyMesh&,
+                const labelList& cellRegion,
+                const vector& perturbVec,
+                const point& p
+            );
+
             //- Split mesh. Keep part containing point.
             autoPtr<mapPolyMesh> splitMeshRegions
             (

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -621,89 +621,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
 }
 
 
-// Return a list of coupled face pairs, i.e. faces that use the same vertices.
-// (this information is recalculated instead of maintained since would be too
-// hard across splitMeshRegions).
-Foam::List<Foam::labelPair> Foam::meshRefinement::getDuplicateFaces
-(
-    const labelList& testFaces
-) const
-{
-    labelList duplicateFace
-    (
-        localPointRegion::findDuplicateFaces
-        (
-            mesh_,
-            testFaces
-        )
-    );
-
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-    // Convert into list of coupled face pairs (mesh face labels).
-    List<labelPair> duplicateFaces(testFaces.size());
-    label dupI = 0;
-
-    forAll(duplicateFace, i)
-    {
-        label otherFaceI = duplicateFace[i];
-
-        if (otherFaceI != -1 && i < otherFaceI)
-        {
-            label meshFace0 = testFaces[i];
-            label patch0 = patches.whichPatch(meshFace0);
-            label meshFace1 = testFaces[otherFaceI];
-            label patch1 = patches.whichPatch(meshFace1);
-
-            if
-            (
-                (patch0 != -1 && isA<processorPolyPatch>(patches[patch0]))
-             || (patch1 != -1 && isA<processorPolyPatch>(patches[patch1]))
-            )
-            {
-                FatalErrorIn
-                (
-                    "meshRefinement::getDuplicateFaces"
-                    "(const bool, const labelList&)"
-                )   << "One of two duplicate faces is on"
-                    << " processorPolyPatch."
-                    << "This is not allowed." << nl
-                    << "Face:" << meshFace0
-                    << " is on patch:" << patches[patch0].name()
-                    << nl
-                    << "Face:" << meshFace1
-                    << " is on patch:" << patches[patch1].name()
-                    << abort(FatalError);
-            }
-
-            duplicateFaces[dupI++] = labelPair(meshFace0, meshFace1);
-        }
-    }
-    duplicateFaces.setSize(dupI);
-
-    Info<< "getDuplicateFaces : found " << returnReduce(dupI, sumOp<label>())
-        << " pairs of duplicate faces." << nl << endl;
-
-
-    if (debug&MESH)
-    {
-        faceSet duplicateFaceSet(mesh_, "duplicateFaces", 2*dupI);
-
-        forAll(duplicateFaces, i)
-        {
-            duplicateFaceSet.insert(duplicateFaces[i][0]);
-            duplicateFaceSet.insert(duplicateFaces[i][1]);
-        }
-        Pout<< "Writing duplicate faces (baffles) to faceSet "
-            << duplicateFaceSet.name() << nl << endl;
-        duplicateFaceSet.instance() = timeName();
-        duplicateFaceSet.write();
-    }
-
-    return duplicateFaces;
-}
-
-
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
 (
     const labelList& globalToMasterPatch,
@@ -1464,18 +1381,16 @@ void Foam::meshRefinement::findCellZoneInsideWalk
             << endl;
 
         // Find the region containing the insidePoint
-        label keepRegionI = -1;
+        label keepRegionI = findRegion
-
+        (
-        label cellI = mesh_.findCell(insidePoint);
+            mesh_,
-
+            cellRegion,
-        if (cellI != -1)
+            mergeDistance_*vector(1,1,1),
-        {
+            insidePoint
-            keepRegionI = cellRegion[cellI];
+        );
-        }
-        reduce(keepRegionI, maxOp<label>());
 
         Info<< "For surface " << surfaces_.names()[surfI]
-            << " found point " << insidePoint << " in cell " << cellI
+            << " found point " << insidePoint
             << " in global region " << keepRegionI
             << " out of " << cellRegion.nRegions() << " regions." << endl;
 
@@ -1631,19 +1546,16 @@ void Foam::meshRefinement::findCellZoneTopo
     }
 
 
-
     // Find the region containing the keepPoint
-    label keepRegionI = -1;
+    label keepRegionI = findRegion
+    (
+        mesh_,
+        cellRegion,
+        mergeDistance_*vector(1,1,1),
+        keepPoint
+    );
 
-    label cellI = mesh_.findCell(keepPoint);
+    Info<< "Found point " << keepPoint
-
-    if (cellI != -1)
-    {
-        keepRegionI = cellRegion[cellI];
-    }
-    reduce(keepRegionI, maxOp<label>());
-
-    Info<< "Found point " << keepPoint << " in cell " << cellI
         << " in global region " << keepRegionI
         << " out of " << cellRegion.nRegions() << " regions." << endl;
 
@@ -2621,11 +2533,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
         (
             freeStandingBaffles    // filter out freestanding baffles
             (
-                getDuplicateFaces   // get all baffles
+                localPointRegion::findDuplicateFacePairs(mesh_),
-                (
-                    identity(mesh_.nFaces()-mesh_.nInternalFaces())
-                   +mesh_.nInternalFaces()
-                ),
                 planarAngle
             )
         );
@@ -2713,17 +2621,15 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
     blockedFace.clear();
 
     // Find the region containing the keepPoint
-    label keepRegionI = -1;
+    label keepRegionI = findRegion
+    (
+        mesh_,
+        cellRegion,
+        mergeDistance_*vector(1,1,1),
+        keepPoint
+    );
 
-    label cellI = mesh_.findCell(keepPoint);
+    Info<< "Found point " << keepPoint
-
-    if (cellI != -1)
-    {
-        keepRegionI = cellRegion[cellI];
-    }
-    reduce(keepRegionI, maxOp<label>());
-
-    Info<< "Found point " << keepPoint << " in cell " << cellI
         << " in global region " << keepRegionI
         << " out of " << cellRegion.nRegions() << " regions." << endl;
 

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C

@@ -1204,6 +1204,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
 
 
     // Bit of statistics
+    if (debug)
     {
         label nSet = 0;
         label nNormals = 9;

src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C

@@ -90,16 +90,16 @@ void Foam::cyclicACMIPolyPatch::resetAMI
             AMIPatchToPatchInterpolation::imPartialFaceAreaWeight
         );
 
-        const scalarField& srcWeightSum = AMI().srcWeightsSum();
+        srcMask_ =
+            min(1.0 - tolerance_, max(tolerance_, AMI().srcWeightsSum()));
+
+        tgtMask_ =
+            min(1.0 - tolerance_, max(tolerance_, AMI().tgtWeightsSum()));
 
-        // set patch face areas based on sum of AMI weights per face
         forAll(Sf, faceI)
         {
-            scalar w = srcWeightSum[faceI];
+            Sf[faceI] *= srcMask_[faceI];
-            w = min(1.0 - tolerance_, max(tolerance_, w));
+            noSf[faceI] *= 1.0 - srcMask_[faceI];
-
-            Sf[faceI] *= w;
-            noSf[faceI] *= 1.0 - w;
         }
 
         setNeighbourFaceAreas();
@@ -116,8 +116,6 @@ void Foam::cyclicACMIPolyPatch::setNeighbourFaceAreas() const
         refCast<const cyclicACMIPolyPatch>(this->neighbPatch());
     const polyPatch& pp = cp.nonOverlapPatch();
 
-    const scalarField& tgtWeightSum = AMI().tgtWeightsSum();
-
     const vectorField& faceAreas0 = cp.faceAreas0();
 
     vectorField::subField Sf = cp.faceAreas();
@@ -125,11 +123,8 @@ void Foam::cyclicACMIPolyPatch::setNeighbourFaceAreas() const
 
     forAll(Sf, faceI)
     {
-        scalar w = tgtWeightSum[faceI];
+        Sf[faceI] = tgtMask_[faceI]*faceAreas0[faceI];
-        w = min(1.0 - tolerance_, max(tolerance_, w));
+        noSf[faceI] = (1.0 - tgtMask_[faceI])*faceAreas0[faceI];
-
-        Sf[faceI] = w*faceAreas0[faceI];
-        noSf[faceI] = (1.0 - w)*faceAreas0[faceI];
     }
 }
 
@@ -184,6 +179,18 @@ void Foam::cyclicACMIPolyPatch::clearGeom()
 }
 
 
+const Foam::scalarField& Foam::cyclicACMIPolyPatch::srcMask() const
+{
+    return srcMask_;
+}
+
+
+const Foam::scalarField& Foam::cyclicACMIPolyPatch::tgtMask() const
+{
+    return tgtMask_;
+}
+
+
 // * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * //
 
 Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
@@ -201,6 +208,8 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     faceAreas0_(),
     nonOverlapPatchName_(word::null),
     nonOverlapPatchID_(-1),
+    srcMask_(),
+    tgtMask_(),
     updated_(false)
 {
     // Non-overlapping patch might not be valid yet so cannot determine
@@ -221,6 +230,8 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     faceAreas0_(),
     nonOverlapPatchName_(dict.lookup("nonOverlapPatch")),
     nonOverlapPatchID_(-1),
+    srcMask_(),
+    tgtMask_(),
     updated_(false)
 {
     if (nonOverlapPatchName_ == name)
@@ -256,6 +267,8 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     faceAreas0_(),
     nonOverlapPatchName_(pp.nonOverlapPatchName_),
     nonOverlapPatchID_(-1),
+    srcMask_(),
+    tgtMask_(),
     updated_(false)
 {
     // Non-overlapping patch might not be valid yet so cannot determine
@@ -278,6 +291,8 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     faceAreas0_(),
     nonOverlapPatchName_(nonOverlapPatchName),
     nonOverlapPatchID_(-1),
+    srcMask_(),
+    tgtMask_(),
     updated_(false)
 {
     if (nonOverlapPatchName_ == name())
@@ -314,6 +329,8 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     faceAreas0_(),
     nonOverlapPatchName_(pp.nonOverlapPatchName_),
     nonOverlapPatchID_(-1),
+    srcMask_(),
+    tgtMask_(),
     updated_(false)
 {}
 
@@ -326,6 +343,13 @@ Foam::cyclicACMIPolyPatch::~cyclicACMIPolyPatch()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+const Foam::cyclicACMIPolyPatch& Foam::cyclicACMIPolyPatch::neighbPatch() const
+{
+    const polyPatch& pp = this->boundaryMesh()[neighbPatchID()];
+    return refCast<const cyclicACMIPolyPatch>(pp);
+}
+
+
 Foam::label Foam::cyclicACMIPolyPatch::nonOverlapPatchID() const
 {
     if (nonOverlapPatchID_ == -1)

src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.H

@@ -66,6 +66,12 @@ private:
         //- Index of non-overlapping patch
         mutable label nonOverlapPatchID_;
 
+        //- Mask/weighting for source patch
+        mutable scalarField srcMask_;
+
+        //- Mask/weighting for target patch
+        mutable scalarField tgtMask_;
+
         //- Flag to indicate that AMI has been updated
         mutable bool updated_;
 
@@ -111,6 +117,12 @@ protected:
         //- Clear geometry
         virtual void clearGeom();
 
+        //- Return the mask/weighting for the source patch
+        virtual const scalarField& srcMask() const;
+
+        //- Return the mask/weighting for the target patch
+        virtual const scalarField& tgtMask() const;
+
 
 public:
 
@@ -245,6 +257,9 @@ public:
             //- Return access to the original patch face areas
             inline const vectorField& faceAreas0() const;
 
+            //- Return a reference to the neighbour patch
+            virtual const cyclicACMIPolyPatch& neighbPatch() const;
+
             //- Non-overlapping patch name
             inline const word& nonOverlapPatchName() const;
 

src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchI.H

@@ -71,11 +71,11 @@ inline const Foam::scalarField& Foam::cyclicACMIPolyPatch::mask() const
 {
     if (owner())
     {
-        return AMI().srcWeightsSum();
+        return srcMask_;
     }
     else
     {
-        return neighbPatch().AMI().tgtWeightsSum();
+        return neighbPatch().tgtMask();
     }
 }
 

src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchTemplates.C

@@ -34,7 +34,7 @@ Foam::tmp<Foam::Field<Type> > Foam::cyclicACMIPolyPatch::interpolate
 {
     if (owner())
     {
-        const scalarField& w = AMI().srcWeightsSum();
+        const scalarField& w = srcMask_;
 
         return
             w*AMI().interpolateToSource(fldCouple)
@@ -42,7 +42,7 @@ Foam::tmp<Foam::Field<Type> > Foam::cyclicACMIPolyPatch::interpolate
     }
     else
     {
-        const scalarField& w = neighbPatch().AMI().tgtWeightsSum();
+        const scalarField& w = neighbPatch().tgtMask();
 
         return
             w*neighbPatch().AMI().interpolateToTarget(fldCouple)
@@ -73,14 +73,14 @@ void Foam::cyclicACMIPolyPatch::interpolate
 {
     if (owner())
     {
-        const scalarField& w = AMI().srcWeightsSum();
+        const scalarField& w = srcMask_;
 
         AMI().interpolateToSource(fldCouple, cop, result);
         result = w*result + (1.0 - w)*fldNonOverlap;
     }
     else
     {
-        const scalarField& w = neighbPatch().AMI().tgtWeightsSum();
+        const scalarField& w = neighbPatch().tgtMask();
 
         neighbPatch().AMI().interpolateToTarget(fldCouple, cop, result);
         result = w*result + (1.0 - w)*fldNonOverlap;

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C

@@ -642,6 +642,13 @@ bool Foam::cyclicAMIPolyPatch::owner() const
 }
 
 
+const Foam::cyclicAMIPolyPatch& Foam::cyclicAMIPolyPatch::neighbPatch() const
+{
+    const polyPatch& pp = this->boundaryMesh()[neighbPatchID()];
+    return refCast<const cyclicAMIPolyPatch>(pp);
+}
+
+
 const Foam::autoPtr<Foam::searchableSurface>&
 Foam::cyclicAMIPolyPatch::surfPtr() const
 {

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H

@@ -284,7 +284,7 @@ public:
             virtual bool owner() const;
 
             //- Return a reference to the neighbour patch
-            inline const cyclicAMIPolyPatch& neighbPatch() const;
+            virtual const cyclicAMIPolyPatch& neighbPatch() const;
 
             //- Return a reference to the projection surface
             const autoPtr<searchableSurface>& surfPtr() const;

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchI.H

@@ -38,14 +38,6 @@ inline const Foam::word& Foam::cyclicAMIPolyPatch::neighbPatchName() const
 }
 
 
-inline const Foam::cyclicAMIPolyPatch&
-Foam::cyclicAMIPolyPatch::neighbPatch() const
-{
-    const polyPatch& pp = this->boundaryMesh()[neighbPatchID()];
-    return refCast<const cyclicAMIPolyPatch>(pp);
-}
-
-
 inline const Foam::vector& Foam::cyclicAMIPolyPatch::rotationAxis() const
 {
     return rotationAxis_;

src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C

@@ -203,7 +203,7 @@ void Foam::mappedPatchBase::findSamples
     {
         case NEARESTCELL:
         {
-            if (samplePatch().size() && samplePatch() != "none")
+            if (samplePatch_.size() && samplePatch_ != "none")
             {
                 FatalErrorIn
                 (
@@ -495,7 +495,6 @@ void Foam::mappedPatchBase::findSamples
 void Foam::mappedPatchBase::calcMapping() const
 {
     static bool hasWarned = false;
-
     if (mapPtr_.valid())
     {
         FatalErrorIn("mappedPatchBase::calcMapping() const")
@@ -509,10 +508,8 @@ void Foam::mappedPatchBase::calcMapping() const
     // Get offsetted points
     const pointField offsettedPoints(samplePoints(patchPoints()));
 
-
+    // Do a sanity check - am I sampling my own patch?
-    // Do a sanity check
+    // This only makes sense for a non-zero offset.
-    // Am I sampling my own patch? This only makes sense for a non-zero
-    // offset.
     bool sampleMyself =
     (
         mode_ == NEARESTPATCHFACE
@@ -550,7 +547,6 @@ void Foam::mappedPatchBase::calcMapping() const
             << "offsetMode_:" << offsetModeNames_[offsetMode_] << endl;
     }
 
-
     // Get global list of all samples and the processor and face they come from.
     pointField samples;
     labelList patchFaceProcs;
@@ -565,7 +561,6 @@ void Foam::mappedPatchBase::calcMapping() const
         patchFc
     );
 
-
     // Find processor and cell/face samples are in and actual location.
     labelList sampleProcs;
     labelList sampleIndices;
@@ -641,7 +636,6 @@ void Foam::mappedPatchBase::calcMapping() const
         }
     }
 
-
     // Now we have all the data we need:
     // - where sample originates from (so destination when mapping):
     //   patchFaces, patchFaceProcs.
@@ -688,7 +682,6 @@ void Foam::mappedPatchBase::calcMapping() const
         }
     }
 
-
     // Determine schedule.
     mapPtr_.reset(new mapDistribute(sampleProcs, patchFaceProcs));
 
@@ -806,17 +799,15 @@ void Foam::mappedPatchBase::calcAMI() const
     }
 
     AMIPtr_.clear();
-/*
-    const polyPatch& nbr = samplePolyPatch();
-
-//    pointField nbrPoints(offsettedPoints());
-    pointField nbrPoints(nbr.localPoints());
 
     if (debug)
     {
+        const polyPatch& nbr = samplePolyPatch();
+
+        pointField nbrPoints(nbr.localPoints());
+
         OFstream os(patch_.name() + "_neighbourPatch-org.obj");
         meshTools::writeOBJ(os, samplePolyPatch().localFaces(), nbrPoints);
-    }
 
         // transform neighbour patch to local system
         primitivePatch nbrPatch0
@@ -829,15 +820,13 @@ void Foam::mappedPatchBase::calcAMI() const
             nbrPoints
         );
 
-    if (debug)
-    {
         OFstream osN(patch_.name() + "_neighbourPatch-trans.obj");
         meshTools::writeOBJ(osN, nbrPatch0, nbrPoints);
 
         OFstream osO(patch_.name() + "_ownerPatch.obj");
         meshTools::writeOBJ(osO, patch_.localFaces(), patch_.localPoints());
     }
-*/
+
     // Construct/apply AMI interpolation to determine addressing and weights
     AMIPtr_.reset
     (

src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H

@@ -392,6 +392,7 @@ public:
             //- Get the patch on the region
             const polyPatch& samplePolyPatch() const;
 
+
         // Helpers
 
             //- Get the sample points

src/meshTools/regionSplit/localPointRegion.C

@@ -585,6 +585,65 @@ Foam::labelList Foam::localPointRegion::findDuplicateFaces
 }
 
 
+Foam::List<Foam::labelPair> Foam::localPointRegion::findDuplicateFacePairs
+(
+    const polyMesh& mesh
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Faces to test: all boundary faces
+    labelList testFaces
+    (
+        identity(mesh.nFaces()-mesh.nInternalFaces())
+      + mesh.nInternalFaces()
+    );
+
+    // Find correspondencing baffle face (or -1)
+    const labelList duplicateFace(findDuplicateFaces(mesh, testFaces));
+
+    // Convert into list of coupled face pairs (mesh face labels).
+    DynamicList<labelPair> baffles(testFaces.size());
+
+    forAll(duplicateFace, i)
+    {
+        label otherFaceI = duplicateFace[i];
+
+        if (otherFaceI != -1 && i < otherFaceI)
+        {
+            label meshFace0 = testFaces[i];
+            label patch0 = patches.whichPatch(meshFace0);
+            label meshFace1 = testFaces[otherFaceI];
+            label patch1 = patches.whichPatch(meshFace1);
+
+            // Check for illegal topology. Should normally not happen!
+            if
+            (
+                (patch0 != -1 && isA<processorPolyPatch>(patches[patch0]))
+             || (patch1 != -1 && isA<processorPolyPatch>(patches[patch1]))
+            )
+            {
+                FatalErrorIn
+                (
+                    "localPointRegion::findDuplicateFacePairs(const polyMesh&)"
+                )   << "One of two duplicate faces is on"
+                    << " processorPolyPatch."
+                    << "This is not allowed." << nl
+                    << "Face:" << meshFace0
+                    << " is on patch:" << patches[patch0].name()
+                    << nl
+                    << "Face:" << meshFace1
+                    << " is on patch:" << patches[patch1].name()
+                    << abort(FatalError);
+            }
+
+            baffles.append(labelPair(meshFace0, meshFace1));
+        }
+    }
+    return baffles.shrink();
+}
+
+
 void Foam::localPointRegion::updateMesh(const mapPolyMesh& map)
 {
     {

src/meshTools/regionSplit/localPointRegion.H

@@ -48,6 +48,7 @@ SourceFiles
 #include "HashSet.H"
 #include "faceList.H"
 #include "boolList.H"
+#include "labelPair.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -149,6 +150,9 @@ public:
                 const labelList&
             );
 
+            //- Helper routine to find all baffles (two boundary faces
+            //  using the same points but in reverse order)
+            static labelPairList findDuplicateFacePairs(const polyMesh&);
 
         // Access
 

src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C

@@ -209,6 +209,38 @@ void Foam::pointToPointPlanarInterpolation::calcWeights
         nearestVertex_,
         nearestVertexWeight_
     );
+
+    //if (debug)
+    //{
+    //    forAll(sourcePoints, i)
+    //    {
+    //        Pout<< "source:" << i << " at:" << sourcePoints[i]
+    //            << " 2d:" << localVertices[i]
+    //            << endl;
+    //    }
+    //
+    //
+    //    forAll(destPoints, i)
+    //    {
+    //        label v0 = nearestVertex_[i][0];
+    //        label v1 = nearestVertex_[i][1];
+    //        label v2 = nearestVertex_[i][2];
+    //
+    //        Pout<< "For location " << destPoints[i]
+    //            << " 2d:" << localFaceCentres[i]
+    //            << " sampling vertices" << nl
+    //            << "    " << v0
+    //            << " at:" << sourcePoints[v0]
+    //            << " weight:" << nearestVertexWeight_[i][0] << nl
+    //            << "    " << v1
+    //            << " at:" << sourcePoints[v1]
+    //            << " weight:" << nearestVertexWeight_[i][1] << nl
+    //            << "    " << v2
+    //            << " at:" << sourcePoints[v2]
+    //            << " weight:" << nearestVertexWeight_[i][2] << nl
+    //            << endl;
+    //    }
+    //}
 }
 
 

src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -2745,6 +2745,12 @@ void Foam::triSurfaceTools::calcInterpolationWeights
 
                     calcInterpolationWeights(tri, nearest.rawPoint(), weights);
 
+                    //Pout<< "calcScalingFactors : samplePt:" << samplePt
+                    //    << " distance:" << nearest.distance()
+                    //    << " to verts:" << verts
+                    //    << " weights:" << weights
+                    //    << endl;
+
                     break;
                 }
             }

src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C

@@ -1103,75 +1103,12 @@ void Foam::decompositionMethod::setConstraints
             << "Keeping owner of faces in baffles "
             << " on same processor." << endl;
 
-        // Faces to test: all boundary faces
+        explicitConnections = localPointRegion::findDuplicateFacePairs(mesh);
-        labelList testFaces
+        forAll(explicitConnections, i)
-        (
-            identity(mesh.nFaces()-mesh.nInternalFaces())
-          + mesh.nInternalFaces()
-        );
-
-        // Find correspondencing baffle face (or -1)
-        labelList duplicateFace
-        (
-            localPointRegion::findDuplicateFaces
-            (
-                mesh,
-                testFaces
-            )
-        );
-
-        const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-        // Convert into list of coupled face pairs (mesh face labels).
-        explicitConnections.setSize(testFaces.size());
-        label dupI = 0;
-
-        forAll(duplicateFace, i)
         {
-            label otherFaceI = duplicateFace[i];
+            blockedFace[explicitConnections[i].first()] = false;
-
+            blockedFace[explicitConnections[i].second()] = false;
-            if (otherFaceI != -1 && i < otherFaceI)
-            {
-                label meshFace0 = testFaces[i];
-                label patch0 = patches.whichPatch(meshFace0);
-                label meshFace1 = testFaces[otherFaceI];
-                label patch1 = patches.whichPatch(meshFace1);
-
-                // Check for illegal topology. Should normally not happen!
-                if
-                (
-                    (patch0 != -1 && isA<processorPolyPatch>(patches[patch0]))
-                 || (patch1 != -1 && isA<processorPolyPatch>(patches[patch1]))
-                )
-                {
-                    FatalErrorIn
-                    (
-                        "decompositionMethod::decompose(const polyMesh&)"
-                    )   << "One of two duplicate faces is on"
-                        << " processorPolyPatch."
-                        << "This is not allowed." << nl
-                        << "Face:" << meshFace0
-                        << " is on patch:" << patches[patch0].name()
-                        << nl
-                        << "Face:" << meshFace1
-                        << " is on patch:" << patches[patch1].name()
-                        << abort(FatalError);
         }
-
-                explicitConnections[dupI++] = labelPair(meshFace0, meshFace1);
-                if (blockedFace[meshFace0])
-                {
-                    blockedFace[meshFace0] = false;
-                    //nUnblocked++;
-                }
-                if (blockedFace[meshFace1])
-                {
-                    blockedFace[meshFace1] = false;
-                    //nUnblocked++;
-                }
-            }
-        }
-        explicitConnections.setSize(dupI);
     }
 
     if
@@ -1256,161 +1193,6 @@ Foam::labelList Foam::decompositionMethod::decompose
     const scalarField& cellWeights
 )
 {
-    //labelHashSet sameProcFaces;
-    //
-    //if (decompositionDict_.found("preservePatches"))
-    //{
-    //    wordList pNames(decompositionDict_.lookup("preservePatches"));
-    //
-    //    Info<< nl
-    //        << "Keeping owner of faces in patches " << pNames
-    //        << " on same processor. This only makes sense for cyclics."
-    //        << endl;
-    //
-    //    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-    //
-    //    forAll(pNames, i)
-    //    {
-    //        const label patchI = patches.findPatchID(pNames[i]);
-    //
-    //        if (patchI == -1)
-    //        {
-    //            FatalErrorIn
-    //            (
-    //                "decompositionMethod::decompose(const polyMesh&)")
-    //                << "Unknown preservePatch " << pNames[i]
-    //                << endl << "Valid patches are " << patches.names()
-    //                << exit(FatalError);
-    //        }
-    //
-    //        const polyPatch& pp = patches[patchI];
-    //
-    //        forAll(pp, i)
-    //        {
-    //            sameProcFaces.insert(pp.start() + i);
-    //        }
-    //    }
-    //}
-    //if (decompositionDict_.found("preserveFaceZones"))
-    //{
-    //    wordList zNames(decompositionDict_.lookup("preserveFaceZones"));
-    //
-    //    Info<< nl
-    //        << "Keeping owner and neighbour of faces in zones " << zNames
-    //        << " on same processor" << endl;
-    //
-    //    const faceZoneMesh& fZones = mesh.faceZones();
-    //
-    //    forAll(zNames, i)
-    //    {
-    //        label zoneI = fZones.findZoneID(zNames[i]);
-    //
-    //        if (zoneI == -1)
-    //        {
-    //            FatalErrorIn
-    //            ("decompositionMethod::decompose(const polyMesh&)")
-    //                << "Unknown preserveFaceZone " << zNames[i]
-    //                << endl << "Valid faceZones are " << fZones.names()
-    //                << exit(FatalError);
-    //        }
-    //
-    //        const faceZone& fz = fZones[zoneI];
-    //
-    //        forAll(fz, i)
-    //        {
-    //            sameProcFaces.insert(fz[i]);
-    //        }
-    //    }
-    //}
-    //
-    //
-    //// Specified processor for group of cells connected to faces
-    //
-    ////- Sets of faces to move together
-    //PtrList<labelList> specifiedProcessorFaces;
-    ////- Destination processor
-    //labelList specifiedProcessor;
-    //
-    //label nProcSets = 0;
-    //if (decompositionDict_.found("singleProcessorFaceSets"))
-    //{
-    //    List<Tuple2<word, label> > zNameAndProcs
-    //    (
-    //        decompositionDict_.lookup("singleProcessorFaceSets")
-    //    );
-    //
-    //    specifiedProcessorFaces.setSize(zNameAndProcs.size());
-    //    specifiedProcessor.setSize(zNameAndProcs.size());
-    //
-    //    forAll(zNameAndProcs, setI)
-    //    {
-    //        Info<< "Keeping all cells connected to faceSet "
-    //            << zNameAndProcs[setI].first()
-    //            << " on processor " << zNameAndProcs[setI].second() << endl;
-    //
-    //        // Read faceSet
-    //        faceSet fz(mesh, zNameAndProcs[setI].first());
-    //
-    //        specifiedProcessorFaces.set(setI, new labelList(fz.sortedToc()));
-    //        specifiedProcessor[setI] = zNameAndProcs[setI].second();
-    //        nProcSets += fz.size();
-    //    }
-    //    reduce(nProcSets, sumOp<label>());
-    //}
-    //
-    //
-    //label nUnblocked = returnReduce(sameProcFaces.size(), sumOp<label>());
-    //
-    //if (nProcSets+nUnblocked > 0)
-    //{
-    //    Info<< "Constrained decomposition:" << endl
-    //        << "    faces with same owner and neighbour processor : "
-    //        << nUnblocked << endl
-    //        << "    faces all on same processor                   : "
-    //        << nProcSets << endl << endl;
-    //}
-    //
-    //
-    //// Faces where owner and neighbour are not 'connected' (= all except
-    //// sameProcFaces)
-    //boolList blockedFace(mesh.nFaces(), true);
-    //{
-    //    forAllConstIter(labelHashSet, sameProcFaces, iter)
-    //    {
-    //        blockedFace[iter.key()] = false;
-    //    }
-    //    syncTools::syncFaceList(mesh, blockedFace, andEqOp<bool>());
-    //
-    //    // Add all point connected faces
-    //    boolList procFacePoint(mesh.nPoints(), false);
-    //    forAll(specifiedProcessorFaces, setI)
-    //    {
-    //        const labelList& set = specifiedProcessorFaces[setI];
-    //        forAll(set, fI)
-    //        {
-    //            const face& f = mesh.faces()[set[fI]];
-    //            forAll(f, fp)
-    //            {
-    //                procFacePoint[f[fp]] = true;
-    //            }
-    //        }
-    //    }
-    //    syncTools::syncPointList(mesh, procFacePoint, orEqOp<bool>(), false);
-    //
-    //    forAll(procFacePoint, pointI)
-    //    {
-    //        if (procFacePoint[pointI])
-    //        {
-    //            const labelList& pFaces = mesh.pointFaces()[pointI];
-    //            forAll(pFaces, i)
-    //            {
-    //                blockedFace[pFaces[i]] = false;
-    //            }
-    //        }
-    //    }
-    //    syncTools::syncFaceList(mesh, blockedFace, andEqOp<bool>());
-    //}
-
     boolList blockedFace;
     PtrList<labelList> specifiedProcessorFaces;
     labelList specifiedProcessor;

src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H

@@ -47,7 +47,7 @@ Description
     - alpha_liquidCore : alpha with outside liquid core set to 0
     - alpha_background : alpha with outside background set to 0.
 
-    Histogram:
+    %Histogram:
     - determine histogram of diameter (given minDiameter, maxDiameter, nBins)
     - write graph of number of droplets per bin
     - write graph of sum, average and deviation of droplet volume per bin

src/postProcessing/functionObjects/forces/forceCoeffs/forceCoeffs.H

@@ -46,10 +46,14 @@ Description
         pitchAxis   (0 0 1);
         magUInf     100;
         lRef        3.5;
-        ARef        2.2;
+        Aref        2.2;
+
+        binData
+        {
             nBin        20;
-        binDir      (1 0 0);
+            direction   (1 0 0);
-        binFormat   gnuplot;
+            cumulative  yes;
+        }
     }
     \endverbatim
 
@@ -64,10 +68,15 @@ Description
         pitchAxis    | picth axis              | yes         |
         magUInf      | free stream velocity magnitude | yes  |
         lRef         | reference length scale for moment calculations | yes |
-        ARef         | reference area          | yes |
+        Aref         | reference area          | yes |
-        nBin         | number of data bins     | no          |
+    \endtable
-        binDir       | direction along which bins are defined | no |
+
-        binFormat    | output format for bin data | no       |
+    Bin data is optional, but if the dictionary is present, the entries must
+    be defined according o
+    \table
+        nBin         | number of data bins     | yes         |
+        direction    | direction along which bins are defined | yes |
+        cumulative   | bin data accumulated with incresing distance | yes |
     \endtable
 
 SeeAlso

src/postProcessing/functionObjects/forces/forces/forces.H

@@ -44,8 +44,13 @@ Description
         ...
         log         yes;
         patches     (walls);
+
+        binData
+        {
             nBin        20;
-        binDir      (1 0 0);
+            direction   (1 0 0);
+            cumulative  yes;
+        }
     }
     \endverbatim
 
@@ -55,8 +60,6 @@ Description
         type         | type name: forces       | yes         |
         log          | write force data to standard output | no | no
         patches      | patches included in the forces calculation | yes |
-        nBin         | number of data bins     | no          |
-        binDir       | direction along which bins are defined | no |
         pName        | pressure field name     | no          | p
         UName        | velocity field name     | no          | U
         rhoName      | density field name (see below) | no   | rho
@@ -65,6 +68,14 @@ Description
         fDName       | name of force density field (see below) | no | fD
     \endtable
 
+    Bin data is optional, but if the dictionary is present, the entries must
+    be defined according o
+    \table
+        nBin         | number of data bins     | yes         |
+        direction    | direction along which bins are defined | yes |
+        cumulative   | bin data accumulated with incresing distance | yes |
+    \endtable
+
 Note
   - For incompressible cases, set \c rhoName to \c rhoInf.  You will then be
     required to provide a \c rhoInf value corresponding to the free-stream

src/postProcessing/functionObjects/utilities/Make/files

@@ -36,4 +36,6 @@ yPlusLES/yPlusLESFunctionObject.C
 yPlusRAS/yPlusRAS.C
 yPlusRAS/yPlusRASFunctionObject.C
 
+setTimeStep/setTimeStepFunctionObject.C
+
 LIB = $(FOAM_LIBBIN)/libutilityFunctionObjects

src/postProcessing/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.C

@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "setTimeStepFunctionObject.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(setTimeStepFunctionObject, 0);
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        setTimeStepFunctionObject,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::setTimeStepFunctionObject::setTimeStepFunctionObject
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    functionObject(name),
+    time_(runTime),
+    enabled_(true)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::setTimeStepFunctionObject::on()
+{
+    enabled_ = true;
+}
+
+
+void Foam::setTimeStepFunctionObject::off()
+{
+    enabled_ = false;
+}
+
+
+bool Foam::setTimeStepFunctionObject::start()
+{
+    return true;
+}
+
+
+bool Foam::setTimeStepFunctionObject::execute(const bool forceWrite)
+{
+    return true;
+}
+
+
+bool Foam::setTimeStepFunctionObject::end()
+{
+    return true;
+}
+
+
+bool Foam::setTimeStepFunctionObject::timeSet()
+{
+    return true;
+}
+
+
+bool Foam::setTimeStepFunctionObject::adjustTimeStep()
+{
+    if (enabled())
+    {
+        // Wanted timestep
+        scalar newDeltaT = timeStepPtr_().value(time_.timeOutputValue());
+
+        const_cast<Time&>(time()).setDeltaT(newDeltaT, false);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+bool Foam::setTimeStepFunctionObject::read(const dictionary& dict)
+{
+    enabled_ = dict.lookupOrDefault("enabled", true);
+
+    if (enabled_)
+    {
+        timeStepPtr_ = DataEntry<scalar>::New("deltaT", dict);
+
+        // Check that time has adjustTimeStep
+        const dictionary& controlDict = time_.controlDict();
+
+        Switch adjust;
+        if
+        (
+            !controlDict.readIfPresent<Switch>("adjustTimeStep", adjust)
+         || !adjust
+        )
+        {
+            FatalIOErrorIn("setTimeStep::read(const dictionary&)", dict)
+                << "Need to have 'adjustTimeStep' true to enable external"
+                << " timestep control" << exit(FatalIOError);
+        }
+    }
+    return true;
+}
+
+
+void Foam::setTimeStepFunctionObject::updateMesh(const mapPolyMesh& mpm)
+{}
+
+
+void Foam::setTimeStepFunctionObject::movePoints(const polyMesh& mesh)
+{}
+
+
+// ************************************************************************* //

src/postProcessing/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.H

@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::setTimeStepFunctionObject
+
+Group
+    grpUtilitiesFunctionObjects
+
+Description
+    Overrides the timeStep. Can only be used with
+    solvers with adjustTimeStep control (e.g. pimpleFoam). Makes no attempt
+    to cooperate with other timeStep 'controllers' (maxCo, other
+    functionObjects). Supports 'enabled' flag but none of othe other ones
+    'timeStart', 'timeEnd', 'outputControl' etc.
+
+SourceFiles
+    setTimeStepFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef setTimeStepFunctionObject_H
+#define setTimeStepFunctionObject_H
+
+#include "functionObject.H"
+#include "dictionary.H"
+#include "DataEntry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class setTimeStepFunctionObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class setTimeStepFunctionObject
+:
+    public functionObject
+{
+    // Private data
+
+        //- Reference to the time database
+        const Time& time_;
+
+
+        // Optional user inputs
+
+            //- Switch for the execution - defaults to 'yes/on'
+            bool enabled_;
+
+            //- Time step
+            autoPtr<DataEntry<scalar> > timeStepPtr_;
+
+
+        //- Disallow default bitwise copy construct
+        setTimeStepFunctionObject(const setTimeStepFunctionObject&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const setTimeStepFunctionObject&);
+
+
+public:
+    //- Runtime type information
+    TypeName("setTimeStep");
+
+    // Constructors
+
+        //- Construct from components
+        setTimeStepFunctionObject
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return time database
+            virtual const Time& time() const
+            {
+                return time_;
+            }
+
+            //- Return the enabled flag
+            virtual bool enabled() const
+            {
+                return enabled_;
+            }
+
+
+        // Function object control
+
+            //- Switch the function object on
+            virtual void on();
+
+            //- Switch the function object off
+            virtual void off();
+
+
+            //- Called at the start of the time-loop
+            virtual bool start();
+
+            //- Called at each ++ or += of the time-loop
+            virtual bool execute(const bool forceWrite);
+
+            //- Called when Time::run() determines that the time-loop exits
+            virtual bool end();
+
+            //- Called when time was set at the end of the Time::operator++
+            virtual bool timeSet();
+
+            //- Called at the end of Time::adjustDeltaT() if adjustTime is true
+            virtual bool adjustTimeStep();
+
+            //- Read and set the function object if its data have changed
+            virtual bool read(const dictionary&);
+
+            //- Update for changes of mesh
+            virtual void updateMesh(const mapPolyMesh& mpm);
+
+            //- Update for changes of mesh
+            virtual void movePoints(const polyMesh& mesh);
+
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/sampledSurface/sampledPatch/sampledPatch.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
@@ -248,6 +248,7 @@ Foam::tmp<Foam::vectorField> Foam::sampledPatch::sample
     return sampleField(vField);
 }
 
+
 Foam::tmp<Foam::sphericalTensorField> Foam::sampledPatch::sample
 (
     const volSphericalTensorField& vField
@@ -292,6 +293,7 @@ Foam::tmp<Foam::vectorField> Foam::sampledPatch::sample
     return sampleField(sField);
 }
 
+
 Foam::tmp<Foam::sphericalTensorField> Foam::sampledPatch::sample
 (
     const surfaceSphericalTensorField& sField
@@ -318,6 +320,7 @@ Foam::tmp<Foam::tensorField> Foam::sampledPatch::sample
     return sampleField(sField);
 }
 
+
 Foam::tmp<Foam::scalarField> Foam::sampledPatch::interpolate
 (
     const interpolation<scalar>& interpolator
@@ -335,6 +338,7 @@ Foam::tmp<Foam::vectorField> Foam::sampledPatch::interpolate
     return interpolateField(interpolator);
 }
 
+
 Foam::tmp<Foam::sphericalTensorField> Foam::sampledPatch::interpolate
 (
     const interpolation<sphericalTensor>& interpolator

src/sampling/sampledSurface/sampledPatch/sampledPatch.H

@@ -55,6 +55,7 @@ class sampledPatch
     //- Private typedefs for convenience
     typedef MeshedSurface<face> MeshStorage;
 
+
     // Private data
 
         //- Name of patches
@@ -78,16 +79,17 @@ class sampledPatch
         //- Start indices (in patchFaceLabels_) of patches
         labelList patchStart_;
 
+
     // Private Member Functions
 
-        //- sample field on faces
+        //- Sample field on faces
         template<class Type>
         tmp<Field<Type> > sampleField
         (
             const GeometricField<Type, fvPatchField, volMesh>& vField
         ) const;
 
-        //- sample surface field on faces
+        //- Sample surface field on faces
         template<class Type>
         tmp<Field<Type> > sampleField
         (
@@ -95,10 +97,9 @@ class sampledPatch
         ) const;
 
         template<class Type>
-        tmp<Field<Type> >
+        tmp<Field<Type> > interpolateField(const interpolation<Type>&) const;
-        interpolateField(const interpolation<Type>&) const;
 
-        //- remap action on triangulation or cleanup
+        //- Re-map action on triangulation or cleanup
         virtual void remapFaces(const labelUList& faceMap);
 
 
@@ -180,31 +181,32 @@ public:
         }
 
 
-        //- sample field on surface
+        // Sample
+            //- Sample field on surface
             virtual tmp<scalarField> sample
             (
                 const volScalarField&
             ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
             virtual tmp<vectorField> sample
             (
                 const volVectorField&
             ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
             virtual tmp<sphericalTensorField> sample
             (
                 const volSphericalTensorField&
             ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
             virtual tmp<symmTensorField> sample
             (
                 const volSymmTensorField&
             ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
             virtual tmp<tensorField> sample
             (
                 const volTensorField&
@@ -240,32 +242,34 @@ public:
                 const surfaceTensorField&
             ) const;
 
-        //- interpolate field on surface
+
+        // Interpolate
+
+            //- Interpolate field on surface
             virtual tmp<scalarField> interpolate
             (
                 const interpolation<scalar>&
             ) const;
 
-
+            //- Interpolate field on surface
-        //- interpolate field on surface
             virtual tmp<vectorField> interpolate
             (
                 const interpolation<vector>&
             ) const;
 
-        //- interpolate field on surface
+            //- Interpolate field on surface
             virtual tmp<sphericalTensorField> interpolate
             (
                 const interpolation<sphericalTensor>&
             ) const;
 
-        //- interpolate field on surface
+            //- Interpolate field on surface
             virtual tmp<symmTensorField> interpolate
             (
                 const interpolation<symmTensor>&
             ) const;
 
-        //- interpolate field on surface
+            //- Interpolate field on surface
             virtual tmp<tensorField> interpolate
             (
                 const interpolation<tensor>&

src/sampling/sampledSurface/sampledPatchInternalField/sampledPatchInternalField.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
@@ -204,6 +204,7 @@ Foam::tmp<Foam::vectorField> Foam::sampledPatchInternalField::interpolate
     return interpolateField(interpolator);
 }
 
+
 Foam::tmp<Foam::sphericalTensorField>
 Foam::sampledPatchInternalField::interpolate
 (

src/sampling/sampledSurface/sampledPatchInternalField/sampledPatchInternalField.H

@@ -74,6 +74,7 @@ class sampledPatchInternalField
         template<class Type>
         tmp<Field<Type> > interpolateField(const interpolation<Type>&) const;
 
+
 public:
 
     //- Runtime type information
@@ -97,68 +98,72 @@ public:
 
     // Member Functions
 
-        //- sample field on surface
+        // Sample
+
+            //- Sample field on surface
             virtual tmp<scalarField> sample
             (
                 const volScalarField&
             ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
             virtual tmp<vectorField> sample
             (
                 const volVectorField&
             ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
             virtual tmp<sphericalTensorField> sample
             (
                 const volSphericalTensorField&
             ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
             virtual tmp<symmTensorField> sample
             (
                 const volSymmTensorField&
             ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
             virtual tmp<tensorField> sample
             (
                 const volTensorField&
             ) const;
 
 
-        //- interpolate field on surface
+        // Interpolate
+
+            //- Interpolate field on surface
             virtual tmp<scalarField> interpolate
             (
                 const interpolation<scalar>&
             ) const;
 
-
+            //- Interpolate field on surface
-        //- interpolate field on surface
             virtual tmp<vectorField> interpolate
             (
                 const interpolation<vector>&
             ) const;
 
-        //- interpolate field on surface
+            //- Interpolate field on surface
             virtual tmp<sphericalTensorField> interpolate
             (
                 const interpolation<sphericalTensor>&
             ) const;
 
-        //- interpolate field on surface
+            //- Interpolate field on surface
             virtual tmp<symmTensorField> interpolate
             (
                 const interpolation<symmTensor>&
             ) const;
 
-        //- interpolate field on surface
+            //- Interpolate field on surface
             virtual tmp<tensorField> interpolate
             (
                 const interpolation<tensor>&
             ) const;
 
+
         //- Write
         virtual void print(Ostream&) const;
 };

src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C

@@ -113,7 +113,7 @@ void Foam::sampledSurfaces::sampleAndWrite
 )
 {
     // interpolator for this field
-    autoPtr< interpolation<Type> > interpolator;
+    autoPtr<interpolation<Type> > interpolatorPtr;
 
     const word& fieldName = vField.name();
     const fileName outputDir = outputPath_/vField.time().timeName();
@@ -126,16 +126,16 @@ void Foam::sampledSurfaces::sampleAndWrite
 
         if (s.interpolate())
         {
-            if (interpolator.empty())
+            if (interpolatorPtr.empty())
             {
-                interpolator = interpolation<Type>::New
+                interpolatorPtr = interpolation<Type>::New
                 (
                     interpolationScheme_,
                     vField
                 );
             }
 
-            values = s.interpolate(interpolator());
+            values = s.interpolate(interpolatorPtr());
         }
         else
         {

src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C

@@ -281,9 +281,8 @@ void Foam::chemkinReader::addPressureDependentReaction
                     thirdBodyEfficiencies(speciesTable_, efficiencies)
                 )
             );
-        }
             break;
-
+        }
         case Troe:
         {
             scalarList TroeCoeffs
@@ -337,9 +336,8 @@ void Foam::chemkinReader::addPressureDependentReaction
                     thirdBodyEfficiencies(speciesTable_, efficiencies)
                 )
             );
-        }
             break;
-
+        }
         case SRI:
         {
             scalarList SRICoeffs
@@ -395,12 +393,9 @@ void Foam::chemkinReader::addPressureDependentReaction
                     thirdBodyEfficiencies(speciesTable_, efficiencies)
                 )
             );
-        }
             break;
-
+        }
         default:
-        {
-            if (fofType < 4)
         {
             FatalErrorIn("chemkinReader::addPressureDependentReaction")
                 << "Fall-off function type "
@@ -409,14 +404,6 @@ void Foam::chemkinReader::addPressureDependentReaction
                 << " not implemented"
                 << exit(FatalError);
         }
-            else
-            {
-                FatalErrorIn("chemkinReader::addPressureDependentReaction")
-                    << "Unknown fall-off function type " << fofType
-                    << " on line " << lineNo_-1
-                    << exit(FatalError);
-            }
-        }
     }
 }
 
@@ -537,9 +524,8 @@ void Foam::chemkinReader::addReaction
                     )
                 );
             }
-        }
             break;
-
+        }
         case thirdBodyArrhenius:
         {
             if (rType == nonEquilibriumReversible)
@@ -597,9 +583,8 @@ void Foam::chemkinReader::addReaction
                     )
                 );
             }
-        }
             break;
-
+        }
         case unimolecularFallOff:
         {
             addPressureDependentReaction<FallOffReactionRate>
@@ -616,9 +601,8 @@ void Foam::chemkinReader::addReaction
                 Afactor,
                 RR
             );
-        }
             break;
-
+        }
         case chemicallyActivatedBimolecular:
         {
             addPressureDependentReaction<ChemicallyActivatedReactionRate>
@@ -635,9 +619,8 @@ void Foam::chemkinReader::addReaction
                 Afactor/concFactor,
                 RR
             );
-        }
             break;
-
+        }
         case LandauTeller:
         {
             const scalarList& LandauTellerCoeffs =
@@ -705,9 +688,8 @@ void Foam::chemkinReader::addReaction
                     )
                 );
             }
-        }
             break;
-
+        }
         case Janev:
         {
             const scalarList& JanevCoeffs =
@@ -727,9 +709,8 @@ void Foam::chemkinReader::addReaction
                     FixedList<scalar, 9>(JanevCoeffs)
                 )
             );
-        }
             break;
-
+        }
         case powerSeries:
         {
             const scalarList& powerSeriesCoeffs =
@@ -749,21 +730,17 @@ void Foam::chemkinReader::addReaction
                     FixedList<scalar, 4>(powerSeriesCoeffs)
                 )
             );
-        }
             break;
-
+        }
         case unknownReactionRateType:
         {
             FatalErrorIn("chemkinReader::addReaction")
                 << "Internal error on line " << lineNo_-1
                 << ": reaction rate type has not been set"
                 << exit(FatalError);
-        }
             break;
-
+        }
         default:
-        {
-            if (rrType < 9)
         {
             FatalErrorIn("chemkinReader::addReaction")
                 << "Reaction rate type " << reactionRateTypeNames[rrType]
@@ -771,14 +748,6 @@ void Foam::chemkinReader::addReaction
                 << " not implemented"
                 << exit(FatalError);
         }
-            else
-            {
-                FatalErrorIn("chemkinReader::addReaction")
-                    << "Unknown reaction rate type " << rrType
-                    << " on line " << lineNo_-1
-                    << exit(FatalError);
-            }
-        }
     }
 
 

src/transportModels/interfaceProperties/interfaceProperties.C

@@ -112,7 +112,7 @@ void Foam::interfaceProperties::calculateK()
     const surfaceVectorField& Sf = mesh.Sf();
 
     // Cell gradient of alpha
-    const volVectorField gradAlpha(fvc::grad(alpha1_));
+    const volVectorField gradAlpha(fvc::grad(alpha1_, "nHat"));
 
     // Interpolated face-gradient of alpha
     surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
@@ -123,6 +123,11 @@ void Foam::interfaceProperties::calculateK()
 
     // Face unit interface normal
     surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
+    // surfaceVectorField nHatfv
+    // (
+    //     (gradAlphaf + deltaN_*vector(0, 0, 1)
+    //    *sign(gradAlphaf.component(vector::Z)))/(mag(gradAlphaf) + deltaN_)
+    // );
     correctContactAngle(nHatfv.boundaryField(), gradAlphaf.boundaryField());
 
     // Face unit interface normal flux

src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    externalWallHeatFluxTemperatureFvPatchScalarField
+    Foam::externalWallHeatFluxTemperatureFvPatchScalarField
 
 Group
     grpThermoBoundaryConditions grpWallBoundaryConditions
@@ -34,38 +34,53 @@ Description
     fixed heat transfer coefficient mode.
 
     The condition can operate in two modes:
-        \li fixed heat transfer coefficient: supply h and Ta
+    - fixed heat transfer coefficient: supply h and Ta
-        \li fixed heat flux: supply q
+    - fixed heat flux: supply q
 
-    where
+    where:
-        \var h  = heat transfer coefficient [W/m2/K]
+    \vartable
-        \var Ta = ambient temperature [K]
+        h  | heat transfer coefficient [W/m2/K]
-        \var q  = heat flux [W/m2]
+        Ta | ambient temperature [K]
+        q  | heat flux [W/m2]
+    \endvartable
 
     The thermal conductivity, \c kappa, can either be retrieved from the
     mesh database using the \c lookup option, or from a \c solidThermo
     thermophysical package.
 
+
+    \heading Patch usage
+
+    \table
+        Property     | Description             | Required    | Default value
+        kappa        | thermal conductivity option | yes     |
+        q            | heat flux [W/m2]            | yes*    |
+        Ta           | ambient temperature [K]     | yes*    |
+        h            | heat transfer coefficient [W/m/K] | yes*|
+        thicknessLayers | list of thicknesses per layer [m] | yes |
+        kappaLayers  | list of thermal conductivites per layer [W/m/K] | yes |
+        kappaName    | name of thermal conductivity field | yes  |
+    \endtable
+
     Example of the boundary condition specification:
     \verbatim
     myPatch
     {
         type            externalWallHeatFluxTemperature;
-           kappa           fluidThermo;       // fluidThermo, solidThermo or
+        kappa           fluidThermo;
-                                               // lookup
+        q               uniform 1000;
-            q               uniform 1000;      // heat flux / [W/m2]
+        Ta              uniform 300.0;
-            Ta              uniform 300.0;     // ambient temperature /[K]
+        h               uniform 10.0;
-            h               uniform 10.0;      // heat transfer coeff /[W/Km2]
+        thicknessLayers (0.1 0.2 0.3 0.4);
-            thicknessLayers (0.1 0.2 0.3 0.4); // thickness of layer [m]
+        kappaLayers     (1 2 3 4)
-            kappaLayers     (1 2 3 4)          // thermal conductivity of
+        value           uniform 300.0;
-                                               // layer [W/m/K]
-            value           uniform 300.0;     // initial temperature / [K]
         kappaName       none;
     }
     \endverbatim
 
 Note
-    Only supply \c h and \c Ta, or \c q in the dictionary (see above)
+  - Only supply \c h and \c Ta, or \c q in the dictionary (see above)
+  - kappa entries can be: fluidThermo, solidThermo or lookup
 
 SourceFiles
     externalWallHeatFluxTemperatureFvPatchScalarField.C

src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C

@@ -231,15 +231,15 @@ baffleThickness() const
     {
         if (thickness_.size() != patch().size())
         {
-            FatalErrorIn
+            FatalIOErrorIn
             (
                 " template<class solidType>"
-                " tmp<scalarField> thermalBaffle1DFvPatchScalarField<solidType>
+                " tmp<scalarField> thermalBaffle1DFvPatchScalarField<solidType>"
-                " baffleThickness() const"
+                " baffleThickness() const",
+                solidDict_
             )<< " Field thickness has not been specified "
             << " for patch " << this->patch().name()
-            << " in dictionary " <<  solidDict_
+            << exit(FatalIOError);
-            << abort(FatalError);
         }
 
         return thickness_;
@@ -302,10 +302,15 @@ void thermalBaffle1DFvPatchScalarField<solidType>::autoMap
 )
 {
     mixedFvPatchScalarField::autoMap(m);
+
+    if (this->owner())
+    {
         thickness_.autoMap(m);
         Qs_.autoMap(m);
+    }
 }
 
+
 template<class solidType>
 void thermalBaffle1DFvPatchScalarField<solidType>::rmap
 (

tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict

@@ -159,7 +159,7 @@ castellatedMeshControls
     // section reachable from the locationInMesh is kept.
     // NOTE: This point should never be on a face, always inside a cell, even
     // after refinement.
-    locationInMesh (3 3 0.43);
+    locationInMesh (3.0001 3.0001 0.43);
 
 
     // Whether any faceZones (as specified in the refinementSurfaces)
@@ -311,12 +311,11 @@ meshQualityControls
 
 // Advanced
 
-// Flags for optional output
+// Write flags
-// 0 : only write final meshes
+writeFlags
-// 1 : write intermediate meshes
+(
-// 2 : write volScalarField with cellLevel for postprocessing
+    layerFields     // write volScalarField for layer coverage
-// 4 : write current intersections as .obj files
+);
-debug 0;
 
 
 // Merge tolerance. Is fraction of overall bounding box of initial mesh.

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/0/T

@@ -36,7 +36,6 @@ boundaryField
         type            inletOutlet;
         inletValue      uniform 293;
         value           uniform 293;
-
     }
     frontAndBack
     {

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/0/U

@@ -33,7 +33,8 @@ boundaryField
     }
     outlet
     {
-        type            zeroGradient;
+        type            pressureInletOutletVelocity;
+        value           $internalField;
     }
     frontAndBack
     {

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/0/alphat

@@ -1,45 +0,0 @@
-/*--------------------------------*- 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      alphat;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [1 -1 -1 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    fuel
-    {
-        type            fixedValue;
-        value           uniform 0;
-    }
-    air
-    {
-        type            fixedValue;
-        value           uniform 0;
-    }
-    outlet
-    {
-        type            zeroGradient;
-    }
-    frontAndBack
-    {
-        type            empty;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/0/epsilon

@@ -1,45 +0,0 @@
-/*--------------------------------*- 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 50;
-
-boundaryField
-{
-    fuel
-    {
-        type            fixedValue;
-        value           uniform 1e-05;
-    }
-    air
-    {
-        type            fixedValue;
-        value           uniform 1e-05;
-    }
-    outlet
-    {
-        type            zeroGradient;
-    }
-    frontAndBack
-    {
-        type            empty;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/0/p

@@ -29,14 +29,16 @@ boundaryField
     {
         type            zeroGradient;
     }
-    inlet2
-    {
-        type            zeroGradient;
-    }
     outlet
     {
-        type            fixedValue;
+        type            totalPressure;
-        value           uniform 1e5;
+        p0              $internalField;
+        U               U;
+        phi             phi;
+        rho             none;
+        psi             none;
+        gamma           1;
+        value           $internalField;
     }
     frontAndBack
     {

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/constant/chemistryProperties

@@ -17,14 +17,25 @@ FoamFile
 
 chemistryType
 {
-    chemistrySolver   noChemistrySolver;
+    chemistrySolver   EulerImplicit;
     chemistryThermo   rho;
 }
 
-chemistry           on; // off;
+chemistry           on;
 
-chemCalcFreq        1;
+initialChemicalTimeStep 1e-07;
 
-initialChemicalTimeStep  1e-8; // NOT USED
+EulerImplicitCoeffs
+{
+    cTauChem        1;
+    equilibriumRateLimiter off;
+}
+
+odeCoeffs
+{
+    solver          Rosenbrock43;
+    absTol          1e-12;
+    relTol          0.01;
+}
 
 // ************************************************************************* //

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/constant/combustionProperties

@@ -15,15 +15,12 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  PaSR<rhoChemistryCombustion>;
+combustionModel  laminar<rhoChemistryCombustion>;
 
 active  true;
 
-PaSRCoeffs
+laminarCoeffs
 {
-    Cmix                0.1;
-    turbulentReaction   off;
-    useReactionRate     true;
 }
 
 

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/constant/polyMesh/boundary

@@ -9,37 +9,39 @@ FoamFile
 {
     version     2.0;
     format      ascii;
-    class       volScalarField;
+    class       polyBoundaryMesh;
-    location    "0";
+    location    "constant/polyMesh";
-    object      k;
+    object      boundary;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 2 -2 0 0 0 0];
+4
-
+(
-internalField   uniform 3.75e-05;
-
-boundaryField
-{
     fuel
     {
-        type            fixedValue;
+        type            patch;
-        value           uniform 3.75e-05;
+        nFaces          40;
+        startFace       7860;
     }
     air
     {
-        type            fixedValue;
+        type            patch;
-        value           uniform 3.75e-05;
+        nFaces          40;
+        startFace       7900;
     }
     outlet
     {
-        type            zeroGradient;
+        type            patch;
+        nFaces          200;
+        startFace       7940;
     }
     frontAndBack
     {
         type            empty;
+        inGroups        1(empty);
+        nFaces          8000;
+        startFace       8140;
     }
-}
+)
-
 
 // ************************************************************************* //

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/constant/turbulenceProperties

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-simulationType  RASModel;
+simulationType  laminar;
 
 
 // ************************************************************************* //

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/controlDict

@@ -23,7 +23,7 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         10000;
+endTime         2000;
 
 deltaT          1;
 

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvSolution

@@ -17,79 +17,69 @@ FoamFile
 
 solvers
 {
-    rho
+    "rho.*"
     {
-        solver          PCG;
+        solver          diagonal;
-        preconditioner  DIC;
-        tolerance       0;
-        relTol          0.1;
-    }
-
-    rhoFinal
-    {
-        solver          PCG;
-        preconditioner  DIC;
-        tolerance       1e-4;
-        relTol          0;
     }
 
     p
     {
-        solver           GAMG;
+        solver           PCG;
-
+        preconditioner   DIC;
-        tolerance        0;
+        tolerance        1e-6;
-        relTol           0.05;
+        relTol           0.1;
-
+    }
-        smoother         DICGaussSeidel;
-        nPreSweeps       0;
-        nPostSweeps      2;
-
-        cacheAgglomeration true;
-
-        nCellsInCoarsestLevel 10;
-        agglomerator     faceAreaPair;
-        mergeLevels      1;
-
-        maxIter          50;
-    };
 
     pFinal
     {
         $p;
-
-        tolerance        1e-4;
         relTol           0;
-    };
+    }
 
-    "(U|Yi|h|k|epsilon)"
+    "(U|h|k|epsilon)"
     {
-        solver          smoothSolver;
+        solver          PBiCG;
-        smoother        GaussSeidel;
+        preconditioner  DILU;
-        tolerance       0;
+        tolerance       1e-6;
         relTol          0.1;
     }
 
-    "(U|Yi|h|k|epsilon)Final"
+    "(U|h|k|epsilon)Final"
     {
-        solver          smoothSolver;
+        $U;
-        smoother        GaussSeidel;
+        relTol          0.1;
-        tolerance       1e-4;
+    }
-        relTol          0;
+
+    Yi
+    {
+        $U;
+        relTol          0.1;
     }
 }
 
 PIMPLE
 {
-    nCorrectors     2;
+    momentumPredictor no;
+    nOuterCorrectors  1;
+    nCorrectors     1;
     nNonOrthogonalCorrectors 0;
-    momentumPredictor yes;
-    rhoMin          0.1;
-    rhoMax          1.5;
 
-    maxCo           0.1;
+    maxDeltaT       1e-2;
-    rDeltaTSmoothingCoeff 0.1;
+    maxCo           1;
-    maxDeltaT       2e-4;
+    alphaTemp       0.05;
-    alphaTemp       0.005;
+    rDeltaTSmoothingCoeff 1;
+    rDeltaTDampingCoeff 1;
+}
+
+relaxationFactors
+{
+    fields
+    {
+    }
+    equations
+    {
+        ".*" 1;
+    }
 }
 
 

tutorials/lagrangian/reactingParcelFoam/filter/constant/chemistryProperties

@@ -17,7 +17,7 @@ FoamFile
 
 chemistryType
 {
-    chemistrySolver   ode;
+    chemistrySolver   noChemistrySolver;
     chemistryThermo   rho;
 }
 

tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/chemistryProperties

@@ -17,25 +17,12 @@ FoamFile
 
 chemistryType
 {
-    chemistrySolver   ode;
+    chemistrySolver   noChemistrySolver;
     chemistryThermo   rho;
 }
 
 chemistry       off;
 
-initialChemicalTimeStep 1e-07;
+initialChemicalTimeStep 1e-7;
-
-EulerImplicitCoeffs
-{
-    cTauChem        0.05;
-    equilibriumRateLimiter off;
-}
-
-odeCoeffs
-{
-    solver          RK;
-    eps             0.05;
-}
-
 
 // ************************************************************************* //

tutorials/multiphase/LTSInterFoam/wigleyHull/system/fvSchemes

@@ -27,7 +27,7 @@ gradSchemes
 
 divSchemes
 {
-    div(rho*phi,U)  Gauss linearUpwind grad(U);
+    div(rhoPhi,U)   Gauss linearUpwind grad(U);
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
     div(phi,k)      Gauss upwind;