Merge branch 'master' into molecularDynamics

applications/solvers/compressible/rhoCentralFoam/Allwmake

@@ -1,5 +1,8 @@
 #!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
 
 wmake libso BCs 
 wmake
 
+# ----------------------------------------------------------------- end-of-file

applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C

@@ -77,13 +77,13 @@ smoluchowskiJumpTFvPatchScalarField::smoluchowskiJumpTFvPatchScalarField
 :
     mixedFvPatchScalarField(p, iF),
     accommodationCoeff_(readScalar(dict.lookup("accommodationCoeff"))),
-    Twall_("Twall", dict, p.size())
+    Twall_("Twall", dict, p.size()),
+    gamma_(dict.lookupOrDefault<scalar>("gamma", 1.4))
 {
     if
     (
         mag(accommodationCoeff_) < SMALL
-        ||
+     || mag(accommodationCoeff_) > 2.0
-        mag(accommodationCoeff_) > 2.0
     )
     {
         FatalIOErrorIn
@@ -96,8 +96,8 @@ smoluchowskiJumpTFvPatchScalarField::smoluchowskiJumpTFvPatchScalarField
             "    const dictionary&"
             ")",
             dict
-        )   << "unphysical accommodationCoeff_ specified"
+        )   << "unphysical accommodationCoeff specified"
-            << "(0 < accommodationCoeff_ <= 1)" << endl
+            << "(0 < accommodationCoeff <= 1)" << endl
             << exit(FatalError);
     }
 
@@ -113,15 +113,6 @@ smoluchowskiJumpTFvPatchScalarField::smoluchowskiJumpTFvPatchScalarField
         fvPatchField<scalar>::operator=(patchInternalField());
     }
 
-    if (dict.found("gamma"))
-    {
-        gamma_ = readScalar(dict.lookup("gamma"));
-    }
-    else
-    {
-        gamma_ = 1.4;
-    }
-
     refValue() = *this;
     refGrad() = 0.0;
     valueFraction() = 0.0;

applications/solvers/compressible/rhopSonicFoam/Allwmake

@@ -1,5 +1,8 @@
 #!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake libso BCs 
 wmake
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/incompressible/boundaryFoam/boundaryFoam.C

@@ -82,11 +82,14 @@ int main(int argc, char *argv[])
         U += (Ubar - UbarStar);
         gradP += (Ubar - UbarStar)/(1.0/UEqn.A())().weightedAverage(mesh.V());
 
+        label id = y.size() - 1;
+
         scalar wallShearStress =
-            flowDirection & turbulence->R()()[0] & wallNormal;
+            flowDirection & turbulence->R()()[id] & wallNormal;
 
         scalar yplusWall
-            = ::sqrt(mag(wallShearStress))*y[0]/laminarTransport.nu()()[0];
+//            = ::sqrt(mag(wallShearStress))*y[id]/laminarTransport.nu()()[id];
+            = ::sqrt(mag(wallShearStress))*y[id]/turbulence->nuEff()()[id];
 
         Info<< "Uncorrected Ubar = " << (flowDirection & UbarStar.value())<< tab
             << "pressure gradient = " << (flowDirection & gradP.value()) << tab

applications/solvers/incompressible/simpleFoam/Make/options

@@ -7,5 +7,4 @@ EXE_LIBS = \
     -lincompressibleRASModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -lmeshTools \
+    -lmeshTools
-    /* $(LIB_WM_OPTIONS_DIR)/libfbsdmalloc.o */

applications/solvers/multiphase/lesInterFoam/lesInterFoam.C

@@ -85,6 +85,8 @@ int main(int argc, char *argv[])
 
         #include "continuityErrs.H"
 
+        p = pd + rho*gh;
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake

@@ -1,7 +1,10 @@
 #!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake libso phaseModel
 wmake libso interfacialModels
 wmake libso kineticTheoryModels
 wmake
+
+# ----------------------------------------------------------------- end-of-file

applications/utilities/mesh/advanced/refineWallLayer/refineWallLayer.C

@@ -184,10 +184,6 @@ int main(int argc, char *argv[])
         << endl;
 
     // Transfer DynamicLists to straight ones.
-    labelList cutEdges;
-    cutEdges.transfer(allCutEdges);
-    allCutEdges.clear();
-
     scalarField cutEdgeWeights;
     cutEdgeWeights.transfer(allCutEdgeWeights);
     allCutEdgeWeights.clear();
@@ -199,7 +195,7 @@ int main(int argc, char *argv[])
         mesh,
         cutCells.toc(),     // cells candidate for cutting
         labelList(0),       // cut vertices
-        cutEdges,           // cut edges
+        allCutEdges,        // cut edges
         cutEdgeWeights      // weight on cut edges
     );
 

applications/utilities/mesh/conversion/Optional/Allwmake

@@ -2,7 +2,6 @@
 #
 # Build optional components (eg, may depend on third-party libraries)
 # -----------------------------------------------------------------------------
-
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
@@ -18,4 +17,4 @@ then
    wmake ccm26ToFoam
 fi
 
-# end
+# ----------------------------------------------------------------- end-of-file

applications/utilities/mesh/generation/extrudeMesh/Allwmake

@@ -1,6 +1,8 @@
 #!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake libso extrudeModel
 wmake
 
+# ----------------------------------------------------------------- end-of-file

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C

@@ -433,7 +433,7 @@ int main(int argc, char *argv[])
     Info<< "Finished meshing in = "
         << runTime.elapsedCpuTime() << " s." << endl;
 
-    Pout<< "End\n" << endl;
+    Info<< "End\n" << endl;
 
     return(0);
 }

applications/utilities/mesh/manipulation/setSet/Allwmake

@@ -1,4 +1,5 @@
 #!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
 
 READLINE=0
 if [ -f /usr/include/readline/readline.h ]; then
@@ -9,3 +10,5 @@ if [ -f /usr/include/readline/readline.h ]; then
 fi
 export READLINE
 wmake
+
+# ----------------------------------------------------------------- end-of-file

applications/utilities/parallelProcessing/decomposePar/decomposeMesh.C

@@ -35,11 +35,8 @@ Description
 #include "IOstreams.H"
 #include "SLPtrList.H"
 #include "boolList.H"
-#include "cellList.H"
 #include "primitiveMesh.H"
-#include "cyclicFvPatch.H"
+#include "cyclicPolyPatch.H"
-#include "fvPatchList.H"
-#include "DynamicList.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -33,7 +33,7 @@ License
 #include "OSspecific.H"
 #include "Map.H"
 #include "globalMeshData.H"
-
+#include "DynamicList.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 

applications/utilities/postProcessing/dataConversion/foamToFieldview9/foamToFieldview9.C

@@ -45,7 +45,6 @@ Description
 #include "IOobjectList.H"
 #include "boolList.H"
 #include "stringList.H"
-#include "DynamicList.H"
 #include "cellModeller.H"
 
 #include "floatScalar.H"

applications/utilities/postProcessing/dataConversion/foamToVTK/internalWriter.C

@@ -27,9 +27,6 @@ License
 #include "internalWriter.H"
 #include "writeFuns.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components

applications/utilities/postProcessing/dataConversion/foamToVTK/lagrangianWriter.C

@@ -29,9 +29,6 @@ License
 #include "Cloud.H"
 #include "passiveParticle.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components

applications/utilities/postProcessing/dataConversion/foamToVTK/writePatchGeom.C

@@ -78,6 +78,7 @@ void writePatchGeom
     writeFuns::write(pStream, binary, vertLabels);
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

applications/utilities/postProcessing/dataConversion/foamToVTK/writePointSet.C

@@ -103,6 +103,7 @@ void writePointSet
     writeFuns::write(pStream, binary, pointIDs);
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

applications/utilities/postProcessing/graphics/PV3FoamReader/Allwmake

@@ -17,3 +17,5 @@ then
       ;;
    esac
 fi
+
+# ----------------------------------------------------------------- end-of-file

applications/utilities/postProcessing/graphics/PVFoamReader/Allwmake

@@ -1,4 +1,5 @@
 #!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 if [ -d "$ParaView_DIR" -a -r "$ParaView_DIR" ]
@@ -16,3 +17,5 @@ then
       ;;
    esac
 fi
+
+# ----------------------------------------------------------------- end-of-file

applications/utilities/postProcessing/graphics/ensightFoamReader/Allwmake

@@ -1,2 +1,6 @@
 #!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
 wmake libso
+
+# ----------------------------------------------------------------- end-of-file

applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C

@@ -75,7 +75,7 @@ int main(int argc, char *argv[])
             }
 
             // Give patch area
-            Info<< "    Patch area = " << sum(mesh.Sf()) << endl;
+            Info<< "    Patch area = " << sum(mesh.Sf().boundaryField()[patchi]) << endl;
 
             if (fieldHeader.headerClassName() == "volScalarField")
             {

applications/utilities/preProcessing/applyWallFunctionBounaryConditions/Make/files

@@ -0,0 +1,3 @@
+applyWallFunctionBounaryConditions.C
+
+EXE = $(FOAM_APPBIN)/applyWallFunctionBounaryConditions

applications/utilities/preProcessing/applyWallFunctionBounaryConditions/Make/options

@@ -0,0 +1,6 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume
+

applications/utilities/preProcessing/applyWallFunctionBounaryConditions/applyWallFunctionBounaryConditions.C

@@ -0,0 +1,279 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    applyWallFunctionBounaryConditions
+
+Description
+    Updates OpenFOAM RAS cases to use the new wall function framework
+    Attempts to determine whether case is compressible or incompressible, or
+    can be supplied with -compressible command line argument
+
+\*---------------------------------------------------------------------------*/
+
+
+#include "argList.H"
+#include "fvMesh.H"
+#include "Time.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+#include "wallPolyPatch.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+bool caseIsCompressible(const fvMesh& mesh)
+{
+    // Attempt flux field
+    IOobject phiHeader
+    (
+        "phi",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (phiHeader.headerOk())
+    {
+        surfaceScalarField phi(phiHeader, mesh);
+        if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
+        {
+            return true;
+        }
+    }
+
+    // Attempt density field
+    IOobject rhoHeader
+    (
+        "rho",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (rhoHeader.headerOk())
+    {
+        volScalarField rho(rhoHeader, mesh);
+        if (rho.dimensions() == dimDensity)
+        {
+            return true;
+        }
+    }
+
+    // Attempt pressure field
+    IOobject pHeader
+    (
+        "p",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (pHeader.headerOk())
+    {
+        volScalarField p(pHeader, mesh);
+        if (p.dimensions() == dimMass/sqr(dimTime)/dimLength)
+        {
+            return true;
+        }
+    }
+
+    // Attempt hydrostatic pressure field
+    IOobject pdHeader
+    (
+        "pd",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (pdHeader.headerOk())
+    {
+        volScalarField pd(pdHeader, mesh);
+        if (pd.dimensions() == dimMass/sqr(dimTime)/dimLength)
+        {
+            return true;
+        }
+    }
+
+    // If none of the above are true, assume that the case is incompressible
+    return false;
+}
+
+
+void createVolScalarField
+(
+    const fvMesh& mesh,
+    const word& fieldName,
+    const dimensionSet& dims
+)
+{
+    IOobject fieldHeader
+    (
+        fieldName,
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (!fieldHeader.headerOk())
+    {
+        Info<< "Creating field " << fieldName << nl << endl;
+
+        volScalarField field
+        (
+            IOobject
+            (
+                fieldName,
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("zero", dims, 0.0)
+        );
+
+        field.write();
+    }
+}
+
+
+void replaceBoundaryType
+(
+    const fvMesh& mesh,
+    const word& fieldName,
+    const word& boundaryType,
+    const string& boundaryValue
+)
+{
+    IOobject header
+    (
+        fieldName,
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (!header.headerOk())
+    {
+        return;
+    }
+
+    Info<< "Updating boundary types for field " << header.name() << endl;
+
+    const word oldTypeName = IOdictionary::typeName;
+    const_cast<word&>(IOdictionary::typeName) = word::null;
+
+    IOdictionary dict(header);
+
+    const_cast<word&>(IOdictionary::typeName) = oldTypeName;
+    const_cast<word&>(dict.type()) = dict.headerClassName();
+
+    // Make a backup of the old field
+    word backupName(dict.name() + ".old");
+    Info<< "    copying " << dict.name() << " to "
+        << backupName << endl;
+    IOdictionary dictOld = dict;
+    dictOld.rename(backupName);
+    dictOld.regIOobject::write();
+
+    // Loop through boundary patches and update
+    const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
+    dictionary& boundaryDict = dict.subDict("boundaryField");
+    forAll(bMesh, patchI)
+    {
+        if (isType<wallPolyPatch>(bMesh[patchI]))
+        {
+            word patchName = bMesh[patchI].name();
+            dictionary& oldPatch = boundaryDict.subDict(patchName);
+
+            dictionary newPatch(dictionary::null);
+            newPatch.add("type", boundaryType);
+            newPatch.add("value", ("uniform " + boundaryValue).c_str());
+
+            oldPatch = newPatch;
+        }
+    }
+
+    Info<< "    writing updated " << dict.name() << nl << endl;
+    dict.regIOobject::write();
+}
+
+
+int main(int argc, char *argv[])
+{
+
+#   include "addTimeOptions.H"
+    argList::validOptions.insert("compressible", "");
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    bool compressible = args.options().found("compressible");
+
+    Info<< "Updating turbulence fields to operate using new run time "
+        << "selectable" << nl << "wall functions"
+        << nl << endl;
+
+    if (compressible || caseIsCompressible(mesh))
+    {
+        Info<< "Case treated as compressible" << nl << endl;
+        createVolScalarField
+        (
+            mesh,
+            "mut",
+            dimArea/dimTime*dimDensity
+        );
+        replaceBoundaryType(mesh, "mut", "mutWallFunction", "0");
+    }
+    else
+    {
+        Info<< "Case treated as incompressible" << nl << endl;
+        createVolScalarField(mesh, "nut", dimArea/dimTime);
+        replaceBoundaryType(mesh, "nut", "nutWallFunction", "0");
+    }
+
+    replaceBoundaryType(mesh, "epsilon", "epsilonWallFunction", "0");
+    replaceBoundaryType(mesh, "omega", "omegaWallFunction", "0");
+    replaceBoundaryType(mesh, "k", "kQRWallFunction", "0");
+    replaceBoundaryType(mesh, "q", "kQRWallFunction", "0");
+    replaceBoundaryType(mesh, "R", "kQRWallFunction", "(0 0 0 0 0 0)");
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

doc/Doxygen/Allwmake

@@ -1,4 +1,5 @@
 #!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 umask 22

etc/controlDict

@@ -808,7 +808,7 @@ OptimisationSwitches
 {
     fileModificationSkew 10;
     commsType       nonBlocking; //scheduled; //blocking;
-    floatTransfer   1;
+    floatTransfer   0;
     nProcsSimpleSum 0;
 }
 

etc/settings.csh

@@ -82,7 +82,10 @@ switch ("$WM_COMPILER_INST")
 case OpenFOAM:
     switch ("$WM_COMPILER")
     case Gcc:
-        setenv WM_COMPILER_DIR $WM_THIRD_PARTY_DIR/gcc-4.3.1/platforms/$WM_ARCH$WM_COMPILER_ARCH
+        #setenv WM_COMPILER_DIR $WM_THIRD_PARTY_DIR/gcc-4.3.1/platforms/$WM_ARCH$WM_COMPILER_ARCH
+        setenv WM_COMPILER_DIR $WM_THIRD_PARTY_DIR/gcc-4.3.2/platforms/$WM_ARCH$WM_COMPILER_ARCH
+        _foamAddLib $WM_THIRD_PARTY_DIR/mpfr-2.3.2/platforms/$WM_ARCH$WM_COMPILER_ARCH/lib
+        _foamAddLib $WM_THIRD_PARTY_DIR/gmp-4.2.3/platforms/$WM_ARCH$WM_COMPILER_ARCH/lib
     breaksw
     case Gcc42:
         setenv WM_COMPILER_DIR $WM_THIRD_PARTY_DIR/gcc-4.2.4/platforms/$WM_ARCH$WM_COMPILER_ARCH
@@ -105,13 +108,8 @@ case OpenFOAM:
 endsw
 
 if ($?WM_COMPILER_BIN) then
-    set path=($WM_COMPILER_BIN $path)
+    _foamAddPath $WM_COMPILER_BIN
-
+    _foamAddLib $WM_COMPILER_LIB
-    if ($?LD_LIBRARY_PATH) then
-        setenv LD_LIBRARY_PATH ${WM_COMPILER_LIB}:${LD_LIBRARY_PATH}
-    else
-        setenv LD_LIBRARY_PATH ${WM_COMPILER_LIB}
-    endif
 endif
 
 unset WM_COMPILER_BIN

etc/settings.sh

@@ -99,7 +99,10 @@ case "$WM_COMPILER_INST" in
 OpenFOAM)
     case "$WM_COMPILER" in
     Gcc)
-        export WM_COMPILER_DIR=$WM_THIRD_PARTY_DIR/gcc-4.3.1/platforms/$WM_ARCH$WM_COMPILER_ARCH
+        #export WM_COMPILER_DIR=$WM_THIRD_PARTY_DIR/gcc-4.3.1/platforms/$WM_ARCH$WM_COMPILER_ARCH
+        export WM_COMPILER_DIR=$WM_THIRD_PARTY_DIR/gcc-4.3.2/platforms/$WM_ARCH$WM_COMPILER_ARCH
+        _foamAddLib $WM_THIRD_PARTY_DIR/mpfr-2.3.2/platforms/$WM_ARCH$WM_COMPILER_ARCH/lib
+        _foamAddLib $WM_THIRD_PARTY_DIR/gmp-4.2.3/platforms/$WM_ARCH$WM_COMPILER_ARCH/lib
         ;;
     Gcc42)
         export WM_COMPILER_DIR=$WM_THIRD_PARTY_DIR/gcc-4.2.4/platforms/$WM_ARCH$WM_COMPILER_ARCH
@@ -123,8 +126,8 @@ OpenFOAM)
 esac
 
 if [ -d "$WM_COMPILER_BIN" ]; then
-    export PATH=$WM_COMPILER_BIN:$PATH
+    _foamAddPath $WM_COMPILER_BIN
-    export LD_LIBRARY_PATH=$WM_COMPILER_LIB:$LD_LIBRARY_PATH
+    _foamAddLib $WM_COMPILER_LIB
 fi
 
 unset WM_COMPILER_BIN WM_COMPILER_LIB

src/OpenFOAM/containers/Lists/DynamicList/DynamicList.C

@@ -33,7 +33,7 @@ template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList(Istream& is)
 :
     List<T>(is),
-    nextFree_(List<T>::size())
+    allocSize_(List<T>::size())
 {}
 
 
@@ -44,9 +44,6 @@ Foam::Ostream& Foam::operator<<
     const Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>& DL
 )
 {
-    const_cast<DynamicList<T, SizeInc, SizeMult, SizeDiv>&>(DL)
-        .setSize(DL.nextFree_);
-
     os << static_cast<const List<T>&>(DL);
     return os;
 }
@@ -60,7 +57,7 @@ Foam::Istream& Foam::operator>>
 )
 {
     is >> static_cast<List<T>&>(DL);
-    DL.nextFree_ = DL.List<T>::size();
+    DL.allocSize_ = DL.List<T>::size();
 
     return is;
 }

src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H

@@ -26,7 +26,7 @@ Class
     Foam::DynamicList
 
 Description
-    A 1D vector of objects of type \<T\> which resizes itself as necessary to
+    A 1D vector of objects of type \<T\> that resizes itself as necessary to
     accept the new objects.
 
     Internal storage is a compact array and the list can be shrunk to compact
@@ -81,24 +81,29 @@ class DynamicList
 {
     // Private data
 
-        //- Number of next free element
+        //- Allocated size for underlying List.
-        label nextFree_;
+        label allocSize_;
 
 
 public:
 
+    // Related types
+
+        //- Declare friendship with the List class
+        friend class List<T>;
+
     // Constructors
 
         //- Construct null
         inline DynamicList();
 
-        //- Construct given size
+        //- Construct given size.
         explicit inline DynamicList(const label);
 
-        //- Construct from UList. nextFree_ set to size().
+        //- Construct from UList. Size set to UList size.
         explicit inline DynamicList(const UList<T>&);
 
-        //- Construct from Istream. nextFree_ set to size().
+        //- Construct from Istream. Size set to size of read list.
         explicit DynamicList(Istream&);
 
 
@@ -106,22 +111,24 @@ public:
 
         // Access
 
-            //- Size of the active part of the list.
+            //- Size of the underlying storage.
-            //  Direct over-ride of list size member function
+            inline label allocSize() const;
-            inline label size() const;
 
 
         // Edit
 
             //- Reset size of List.
-            void setSize(const label);
+            inline void setSize(const label);
 
             //- Reset size of List and value for new elements.
-            void setSize(const label, const T&);
+            inline void setSize(const label, const T&);
 
-            //- Clear the list, i.e. set next free to zero.
+            //- Clear the list, i.e. set the size to zero.
             //  Allocated size does not change
-            void clear();
+            inline void clear();
+
+            //- Clear the list and delete storage.
+            inline void clearStorage();
 
             //- Shrink the List<T> to the number of elements used
             inline DynamicList<T, SizeInc, SizeMult, SizeDiv>& shrink();
@@ -129,12 +136,12 @@ public:
             //- Transfer the contents of the argument List into this List
             //  and annull the argument list. Is same as List::transfer except
             //  checks that you're not changing the underlying list to something
-            //  smaller than nextFree_.
+            //  smaller than allocSize_.
-            void transfer(List<T>&);
+            inline void transfer(List<T>&);
 
             //- Transfer the contents of the argument DynamicList into this
             //  DynamicList and annull the argument list.
-            void transfer(DynamicList<T, SizeInc, SizeMult, SizeDiv>&);
+            inline void transfer(DynamicList<T, SizeInc, SizeMult, SizeDiv>&);
 
 
     // Member Operators
@@ -142,7 +149,7 @@ public:
         //- Append an element at the end of the list
         inline void append(const T& e);
 
-        //- Return and remove the top element
+        //- Remove and return the top element
         inline T remove();
 
         //- Return non-const access to an element,
@@ -152,7 +159,7 @@ public:
         //- Assignment of all entries to the given value
         inline void operator=(const T&);
 
-        //- Assignment to List<T>
+        //- Assignment from List<T>. Also handles assignment from DynamicList.
         inline void operator=(const List<T>&);
 
 

src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H

@@ -31,8 +31,10 @@ template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList()
 :
     List<T>(SizeInc),
-    nextFree_(0)
+    allocSize_(SizeInc)
-{}
+{
+    List<T>::size() = 0;
+}
 
 
 //- Construct given size
@@ -43,8 +45,10 @@ inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList
 )
 :
     List<T>(s),
-    nextFree_(0)
+    allocSize_(s)
-{}
+{
+    List<T>::size() = 0;
+}
 
 
 //- Construct given size
@@ -55,17 +59,17 @@ inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList
 )
 :
     List<T>(s),
-    nextFree_(s.size())
+    allocSize_(s.size())
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
-inline Foam::label Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::size()
+inline Foam::label Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::allocSize()
 const
 {
-    return nextFree_;
+    return allocSize_;
 }
 
 
@@ -75,13 +79,16 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::setSize
     const label s
 )
 {
-    if (s < nextFree_)
+    if (s < List<T>::size())
     {
-        nextFree_ = s;
+        List<T>::size() = s;
     }
     else
     {
-        List<T>::setSize(s);
+        label nextFree = List<T>::size();
+        allocSize_ = s;
+        List<T>::setSize(allocSize_);
+        List<T>::size() = nextFree;
     }
 }
 
@@ -93,14 +100,16 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::setSize
     const T& t
 )
 {
-    if (s < nextFree_)
+    if (s < List<T>::size())
     {
-        nextFree_ = s;
+        List<T>::size() = s;
     }
     else
     {
-        List<T>::setSize(s, t);
+        label nextFree = List<T>::size();
-        nextFree_ = s;
+        allocSize_ = s;
+        List<T>::setSize(allocSize_, t);
+        List<T>::size() = nextFree;
     }
 }
 
@@ -108,7 +117,15 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::setSize
 template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::clear()
 {
-    nextFree_ = 0;
+    List<T>::size() = 0;
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::clearStorage()
+{
+    List<T>::clear();
+    allocSize_ = 0;
 }
 
 
@@ -116,7 +133,8 @@ template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>&
 Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::shrink()
 {
-    List<T>::setSize(nextFree_);
+    allocSize_ = List<T>::size();
+    List<T>::setSize(allocSize_);
     return *this;
 }
 
@@ -125,20 +143,20 @@ template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline void
 Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::transfer(List<T>& l)
 {
-    if (l.size() < nextFree_)
+    if (l.size() < List<T>::size())
     {
         FatalErrorIn
         (
             "void DynamicList<T, SizeInc, SizeMult"
             ", SizeDiv>::transfer(List<T>&)"
         )   << "Cannot replace the underlying storage of this DynamicList"
-            << " of which " << nextFree_ << " elements are used" << nl
+            << " of which " << List<T>::size() << " elements are used" << nl
             << "with a List of size " << l.size() << abort(FatalError);
     }
     else
     {
+        allocSize_ = l.size();
         List<T>::transfer(l);   // take over storage
-        l.clear();              // set nextFree of l to 0
     }
 }
 
@@ -150,48 +168,53 @@ Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::transfer
     DynamicList<T, SizeInc, SizeMult, SizeDiv>& l
 )
 {
-    List<T>::transfer(l);       // take over storage
+    allocSize_ = l.allocSize();
-    nextFree_ = l.size();       // take over used size
+    List<T>::transfer(l);       // take over storage. Null l.
-    l.clear();                  // set nextFree of l to 0
 }
 
 
 template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
-inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::append
+inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::append(const T& e)
-(
-    const T& e
-)
 {
-    nextFree_++;
+    // Work on copy free index since gets overwritten by setSize
+    label nextFree = List<T>::size();
 
-    if (nextFree_ > List<T>::size())
+    nextFree++;
+
+    if (nextFree > allocSize_)
     {
-        List<T>::setSize
+        allocSize_ = max
         (
-            max
+            nextFree,
-            (
+            label(SizeMult*allocSize_/SizeDiv + SizeInc)
-                nextFree_,
-                label(SizeMult*List<T>::size()/SizeDiv + SizeInc)
-            )
         );
+        List<T>::setSize(allocSize_);
     }
 
-    this->operator[](nextFree_ - 1) = e;
+    List<T>::size() = nextFree;
+
+    this->operator[](nextFree - 1) = e;
 }
 
 
 template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline T Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::remove()
 {
-    if (nextFree_ == 0)
+    if (List<T>::size() == 0)
     {
         FatalErrorIn
         (
-            "void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::remove()"
+            "Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::remove()"
         )   << "List is empty" << abort(FatalError);
     }
 
-    return List<T>::operator[](--nextFree_);
+    label nextFree = List<T>::size()-1;
+
+    const T& val = List<T>::operator[](nextFree);
+
+    List<T>::size() = nextFree;
+
+    return val;
 }
 
 
@@ -203,20 +226,22 @@ inline T& Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::operator()
     const label i
 )
 {
-    nextFree_ = max(nextFree_, i + 1);
+    label nextFree = List<T>::size();
 
-    if (nextFree_ > List<T>::size())
+    nextFree = max(nextFree, i + 1);
+
+    if (nextFree > allocSize_)
     {
-        List<T>::setSize
+        allocSize_ = max
         (
-            max
+            nextFree,
-            (
+            label(SizeMult*allocSize_/SizeDiv + SizeInc)
-                nextFree_,
-                label(SizeMult*List<T>::size()/SizeDiv + SizeInc)
-            )
         );
+        List<T>::setSize(allocSize_);
     }
 
+    List<T>::size() = nextFree;
+
     return this->operator[](i);
 }
 
@@ -228,7 +253,7 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::operator=
 )
 {
     List<T>::operator=(t);
-    nextFree_ = List<T>::size();
+    allocSize_ = List<T>::size();
 }
 
 
@@ -239,7 +264,7 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::operator=
 )
 {
     List<T>::operator=(l);
-    nextFree_ = l.size();
+    allocSize_ = List<T>::size();
 }
 
 

src/OpenFOAM/containers/Lists/List/List.C

@@ -38,14 +38,11 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * //
 
 // Construct with length specified
 template<class T>
-List<T>::List(const label s)
+Foam::List<T>::List(const label s)
 :
     UList<T>(NULL, s)
 {
@@ -69,7 +66,7 @@ List<T>::List(const label s)
 
 // Construct with length and single value specified
 template<class T>
-List<T>::List(const label s, const T& a)
+Foam::List<T>::List(const label s, const T& a)
 :
     UList<T>(NULL, s)
 {
@@ -98,7 +95,7 @@ List<T>::List(const label s, const T& a)
 
 // Construct as copy
 template<class T>
-List<T>::List(const List<T>& a)
+Foam::List<T>::List(const List<T>& a)
 :
     UList<T>(NULL, a.size_)
 {
@@ -130,7 +127,7 @@ List<T>::List(const List<T>& a)
 
 // Construct as copy or re-use as specified.
 template<class T>
-List<T>::List(List<T>& a, bool reUse)
+Foam::List<T>::List(List<T>& a, bool reUse)
 :
     UList<T>(NULL, a.size_)
 {
@@ -166,10 +163,10 @@ List<T>::List(List<T>& a, bool reUse)
 }
 
 
-// Construct given size and start and end iterators.
+// Construct given start and end iterators.
 template<class T>
 template<class InputIterator>
-List<T>::List(InputIterator first, InputIterator last)
+Foam::List<T>::List(InputIterator first, InputIterator last)
 {
     label s = 0;
     for
@@ -200,8 +197,8 @@ List<T>::List(InputIterator first, InputIterator last)
 
 // Construct as copy of FixedList<T, Size>
 template<class T>
-template<label Size>
+template<Foam::label Size>
-List<T>::List(const FixedList<T, Size>& fl)
+Foam::List<T>::List(const FixedList<T, Size>& fl)
 :
     UList<T>(NULL, Size)
 {
@@ -223,7 +220,7 @@ List<T>::List(const FixedList<T, Size>& fl)
 
 // Construct as copy of PtrList<T>
 template<class T>
-List<T>::List(const PtrList<T>& sptrl)
+Foam::List<T>::List(const PtrList<T>& sptrl)
 :
     UList<T>(NULL, sptrl.size())
 {
@@ -245,7 +242,7 @@ List<T>::List(const PtrList<T>& sptrl)
 
 // Construct as copy of SLList<T>
 template<class T>
-List<T>::List(const SLList<T>& sll)
+Foam::List<T>::List(const SLList<T>& sll)
 :
     UList<T>(NULL, sll.size())
 {
@@ -273,7 +270,7 @@ List<T>::List(const SLList<T>& sll)
 
 // Construct as copy of IndirectList<T>
 template<class T>
-List<T>::List(const IndirectList<T>& idl)
+Foam::List<T>::List(const IndirectList<T>& idl)
 :
     UList<T>(NULL, idl.size())
 {
@@ -295,7 +292,7 @@ List<T>::List(const IndirectList<T>& idl)
 
 // Construct as copy of BiIndirectList<T>
 template<class T>
-List<T>::List(const BiIndirectList<T>& idl)
+Foam::List<T>::List(const BiIndirectList<T>& idl)
 :
     UList<T>(NULL, idl.size())
 {
@@ -319,7 +316,7 @@ List<T>::List(const BiIndirectList<T>& idl)
 
 // Destroy list elements
 template<class T>
-List<T>::~List()
+Foam::List<T>::~List()
 {
     if (this->v_) delete[] this->v_;
 }
@@ -328,7 +325,7 @@ List<T>::~List()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class T>
-const List<T>& List<T>::null()
+const Foam::List<T>& Foam::List<T>::null()
 {
     List<T>* nullPtr = reinterpret_cast<List<T>*>(NULL);
     return *nullPtr;
@@ -336,7 +333,7 @@ const List<T>& List<T>::null()
 
 
 template<class T>
-void List<T>::setSize(const label newSize)
+void Foam::List<T>::setSize(const label newSize)
 {
     if (newSize < 0)
     {
@@ -382,7 +379,7 @@ void List<T>::setSize(const label newSize)
 
 
 template<class T>
-void List<T>::setSize(const label newSize, const T& a)
+void Foam::List<T>::setSize(const label newSize, const T& a)
 {
     label oldSize = this->size_;
     this->setSize(newSize);
@@ -397,7 +394,7 @@ void List<T>::setSize(const label newSize, const T& a)
 
 
 template<class T>
-void List<T>::clear()
+void Foam::List<T>::clear()
 {
     if (this->v_) delete[] this->v_;
     this->size_ = 0;
@@ -408,10 +405,9 @@ void List<T>::clear()
 // Transfer the contents of the argument List into this List
 // and anull the argument list
 template<class T>
-void List<T>::transfer(List<T>& a)
+void Foam::List<T>::transfer(List<T>& a)
 {
     if (this->v_) delete[] this->v_;
-
     this->size_ = a.size_;
     this->v_ = a.v_;
 
@@ -420,29 +416,45 @@ void List<T>::transfer(List<T>& a)
 }
 
 
+// Transfer the contents of the argument DynamicList into this List
+// and anull the argument list
 template<class T>
-void sort(List<T>& a)
+template<unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+void Foam::List<T>::transfer(DynamicList<T, SizeInc, SizeMult, SizeDiv>& a)
+{
+    if (this->v_) delete[] this->v_;
+    this->size_ = a.size_;
+    this->v_ = a.v_;
+
+    a.size_ = 0;
+    a.v_ = 0;
+    a.allocSize_ = 0;
+}
+
+
+template<class T>
+void Foam::sort(List<T>& a)
 {
     std::sort(a.begin(), a.end());
 }
 
 
 template<class T, class Cmp>
-void sort(List<T>& a, const Cmp& cmp)
+void Foam::sort(List<T>& a, const Cmp& cmp)
 {
     std::sort(a.begin(), a.end(), cmp);
 }
 
 
 template<class T>
-void stableSort(List<T>& a)
+void Foam::stableSort(List<T>& a)
 {
     std::stable_sort(a.begin(), a.end());
 }
 
 
 template<class T, class Cmp>
-void stableSort(List<T>& a, const Cmp& cmp)
+void Foam::stableSort(List<T>& a, const Cmp& cmp)
 {
     std::stable_sort(a.begin(), a.end(), cmp);
 }
@@ -452,7 +464,7 @@ void stableSort(List<T>& a, const Cmp& cmp)
 
 // Assignment to UList operator. Takes linear time.
 template<class T>
-void List<T>::operator=(const UList<T>& a)
+void Foam::List<T>::operator=(const UList<T>& a)
 {
     if (a.size_ != this->size_)
     {
@@ -484,7 +496,7 @@ void List<T>::operator=(const UList<T>& a)
 
 // Assignment operator. Takes linear time.
 template<class T>
-void List<T>::operator=(const List<T>& a)
+void Foam::List<T>::operator=(const List<T>& a)
 {
     if (this == &a)
     {
@@ -499,7 +511,7 @@ void List<T>::operator=(const List<T>& a)
 
 // Assignment operator. Takes linear time.
 template<class T>
-void List<T>::operator=(const SLList<T>& sll)
+void Foam::List<T>::operator=(const SLList<T>& sll)
 {
     if (sll.size() != this->size_)
     {
@@ -527,7 +539,7 @@ void List<T>::operator=(const SLList<T>& sll)
 
 // Assignment operator. Takes linear time.
 template<class T>
-void List<T>::operator=(const IndirectList<T>& idl)
+void Foam::List<T>::operator=(const IndirectList<T>& idl)
 {
     if (idl.size() != this->size_)
     {
@@ -549,7 +561,7 @@ void List<T>::operator=(const IndirectList<T>& idl)
 
 // Assignment operator. Takes linear time.
 template<class T>
-void List<T>::operator=(const BiIndirectList<T>& idl)
+void Foam::List<T>::operator=(const BiIndirectList<T>& idl)
 {
     if (idl.size() != this->size_)
     {
@@ -568,11 +580,6 @@ void List<T>::operator=(const BiIndirectList<T>& idl)
     }
 }
 
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 
 #include "ListIO.C"

src/OpenFOAM/containers/Lists/List/List.H

@@ -61,6 +61,8 @@ template<class T> Istream& operator>>(Istream&, List<T>&);
 template<class T, label Size> class FixedList;
 template<class T> class PtrList;
 template<class T> class SLList;
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+ class DynamicList;
 template<class T> class IndirectList;
 template<class T> class BiIndirectList;
 
@@ -94,7 +96,7 @@ public:
         //- Construct as copy or re-use as specified.
         List(List<T>&, bool reUse);
 
-        //- Construct given size and start and end iterators.
+        //- Construct given start and end iterators.
         template<class InputIterator>
         List(InputIterator first, InputIterator last);
 
@@ -140,6 +142,9 @@ public:
         //- Return the number of elements in the UList.
         inline label size() const;
 
+        //- Override size to be inconsistent with allocated storage.
+        //  Use with care.
+        inline label& size();
 
         // Edit
 
@@ -156,14 +161,14 @@ public:
             //  and annull the argument list.
             void transfer(List<T>&);
 
+            //- Transfer the contents of the argument List into this List
+            //  and annull the argument list.
+            template<unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+            void transfer(DynamicList<T, SizeInc, SizeMult, SizeDiv>&);
+
             //- Return subscript-checked element of UList.
             inline T& newElmt(const label);
 
-            //- Override size to be inconsistent with allocated storage.
-            //  Use with care.
-            inline label& size();
-
-
     // Member operators
 
         //- Assignment from UList operator. Takes linear time.

src/OpenFOAM/containers/Lists/List/ListIO.C

@@ -30,16 +30,11 @@ License
 #include "SLList.H"
 #include "contiguous.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 // Construct from Istream
 template<class T>
-List<T>::List(Istream& is)
+Foam::List<T>::List(Istream& is)
 :
     UList<T>(NULL, 0)
 {
@@ -48,7 +43,7 @@ List<T>::List(Istream& is)
 
 
 template<class T>
-Istream& operator>>(Istream& is, List<T>& L)
+Foam::Istream& Foam::operator>>(Istream& is, List<T>& L)
 {
     // Anull list
     L.setSize(0);
@@ -161,9 +156,4 @@ Istream& operator>>(Istream& is, List<T>& L)
     return is;
 }
 
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/OpenFOAM/containers/Lists/PackedList/PackedList.C

@@ -40,10 +40,7 @@ PackedList<nBits>::PackedList(const label size, const unsigned int val)
     List<unsigned int>(intSize(size)),
     size_(size)
 {
-    for (label i = 0; i < size; i++)
+    operator=(val);
-    {
-        set(i, val);
-    }
 }
 
 

src/OpenFOAM/containers/Lists/PackedList/PackedList.H

@@ -160,13 +160,15 @@ public:
             //- Number of packed elements
             inline label size() const;
 
-
             //- Get value at index I
             inline unsigned int get(const label i) const;
 
             //- Set value at index I. Return true if value changed.
             inline bool set(const label i, const unsigned int val);
 
+            //- Underlying storage
+            inline List<unsigned int>& storage();
+
 
     // Member operators
 
@@ -179,7 +181,8 @@ public:
             //- Assignment operator. Takes linear time.
             void operator=(const PackedList<nBits>&);
 
-            //- Assignment of all entries to the given value
+            //- Assignment of all entries to the given value. Does set on all
+            //  elements.
             inline void operator=(const unsigned int val);
 
             //- Return as labelList

src/OpenFOAM/containers/Lists/PackedList/PackedListI.H

@@ -83,7 +83,7 @@ inline void PackedList<nBits>::checkValue(const unsigned int val) const
 {
     if (val>=(1u << nBits))
     {
-        FatalErrorIn("PackedList<T>::set(const unsigned int)")
+        FatalErrorIn("PackedList<T>::checkValue(const unsigned int)")
             << "value " << label(val) << " out of range 0 ... "
             << label((1u << nBits)-1)
             << " representable by " << nBits << " bits"
@@ -183,6 +183,13 @@ inline bool PackedList<nBits>::set(const label i, const unsigned int val)
 }
 
 
+template<int nBits>
+inline List<unsigned int>& PackedList<nBits>::storage()
+{
+    return static_cast<List<unsigned int>&>(*this);
+}
+
+
 template<int nBits>
 inline ::Foam::reference PackedList<nBits>::operator[](const label i)
 {
@@ -213,9 +220,21 @@ inline void PackedList<nBits>::operator=(const unsigned int val)
 #   ifdef DEBUGList
     checkValue(val);
 #   endif
-    List<unsigned int>::operator=(val);
+
+    if (val == 0)
+    {
+        List<unsigned int>::operator=(val);
+    }
+    else
+    {
+        for (label i = 0; i < size_; i++)
+        {
+            set(i, val);
+        }
+    }
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/containers/Lists/PackedList/PackedListIO.C

@@ -1,152 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "PackedList.H"
-#include "Ostream.H"
-#include "token.H"
-#include "contiguous.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * Ostream Operator *  * * * * * * * * * * * * //
-
-template<class T>
-void PackedList<T>::writeEntry(Ostream& os) const
-{
-    if
-    (
-        size()
-     && token::compound::isCompound
-        (
-            "List<" + word(pTraits<T>::typeName) + '>'
-        )
-    )
-    {
-        os  << word("List<" + word(pTraits<T>::typeName) + '>') << " ";
-    }
-    
-    os << *this;
-}
-
-
-template<class T>
-void PackedList<T>::writeEntry(const word& keyword, Ostream& os) const
-{
-    os.writeKeyword(keyword);
-    writeEntry(os);
-    os << token::END_STATEMENT << endl;
-}
-
-
-template<class T>
-Ostream& operator<<(Ostream& os, const PackedList<T>& L)
-{
-    // Write list contents depending on data format
-    if (os.format() == IOstream::ASCII || !contiguous<T>())
-    {
-        bool uniform = false;
-
-        if (L.size() > 1 && contiguous<T>())
-        {
-            uniform = true;
-
-            unsigned int L0 = L.get(0);
-
-            forAll(L, i)
-            {
-                if (L.get(i) != L0)
-                {
-                    uniform = false;
-                    break;
-                }
-            }
-        }
-
-        if (uniform)
-        {
-            // Write size of list and start contents delimiter
-            os << L.size() << token::BEGIN_BLOCK;
-
-            // Write list contents
-            os << L.get(0);
-
-            // Write end of contents delimiter
-            os << token::END_BLOCK;
-        }
-        else if (L.size() < 11 && contiguous<T>())
-        {
-            // Write size of list and start contents delimiter
-            os << L.size() << token::BEGIN_LIST;
-
-            // Write list contents
-            forAll(L, i)
-            {
-                if (i > 0) os << token::SPACE;
-                os << L.get(i);
-            }
-
-            // Write end of contents delimiter
-            os << token::END_LIST;
-        }
-        else
-        {
-            // Write size of list and start contents delimiter
-            os << nl << L.size() << nl << token::BEGIN_LIST;
-
-            // Write list contents
-            forAll(L, i)
-            {
-                os << nl << L.get(i);
-            }
-
-            // Write end of contents delimiter
-            os << nl << token::END_LIST << nl;
-        }
-    }
-    else
-    {
-        os << nl << L.size() << nl;
-        if (L.size())
-        {
-            os.write((const char*)L.v_, L.size()*sizeof(T));
-        }
-    }
-
-    // Check state of IOstream
-    os.check("Ostream& operator<<(Ostream&, const PackedList&)");
-
-    return os;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/OpenFOAM/meshes/boundBox/boundBox.C

@@ -79,10 +79,9 @@ boundBox::boundBox(const pointField& points, const bool doReduce)
 
 
 boundBox::boundBox(Istream& is)
-:
+{
-    min_(is),
+    operator>>(is, *this);
-    max_(is)
+}
-{}
 
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
@@ -113,7 +112,7 @@ Istream& operator>>(Istream& is, boundBox& bb)
 {
     if (is.format() == IOstream::ASCII)
     {
-        return is >> bb.min_ >>  bb.max_;
+        return is >> bb.min_ >> bb.max_;
     }
     else
     {

src/OpenFOAM/meshes/meshShapes/edge/edge.H

@@ -141,6 +141,9 @@ inline label Hash<edge>::operator()(const edge& e) const
     return e[0]*e[1] + e[0]+e[1];
 }
 
+template<>
+inline bool contiguous<edge>()  {return true;}
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/meshShapes/triFace/triFace.H

@@ -166,6 +166,9 @@ inline label Hash<triFace>::operator()(const triFace& t) const
     return (t[0]*t[1]*t[2] + t[0]+t[1]+t[2]);
 }
 
+template<>
+inline bool contiguous<triFace>()  {return true;}
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.C

@@ -26,12 +26,6 @@ License
 
 #include "globalIndex.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::globalIndex::globalIndex(const label localSize)
@@ -68,15 +62,6 @@ Foam::globalIndex::globalIndex(Istream& is)
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
 Foam::Istream& Foam::operator>>(Istream& is, globalIndex& gi)

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C

@@ -397,13 +397,11 @@ void Foam::globalMeshData::calcSharedEdges() const
     sharedEdgeLabelsPtr_ = new labelList();
     labelList& sharedEdgeLabels = *sharedEdgeLabelsPtr_;
     sharedEdgeLabels.transfer(dynSharedEdgeLabels);
-    dynSharedEdgeLabels.clear();
 
     dynSharedEdgeAddr.shrink();
     sharedEdgeAddrPtr_ = new labelList();
     labelList& sharedEdgeAddr = *sharedEdgeAddrPtr_;
     sharedEdgeAddr.transfer(dynSharedEdgeAddr);
-    dynSharedEdgeAddr.clear();
 
     if (debug)
     {

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C

@@ -84,16 +84,32 @@ void Foam::mapDistribute::distribute
         // Receive sub field from neighbour
         for (label domain = 0; domain < Pstream::nProcs(); domain++)
         {
-            if
+            const labelList& map = constructMap[domain];
-            (
+
-                domain != Pstream::myProcNo()
+            if (domain != Pstream::myProcNo() && map.size() > 0)
-             && constructMap[domain].size() > 0
-            )
             {
                 IPstream fromNbr(Pstream::blocking, domain);
                 List<T> subField(fromNbr);
 
-                const labelList& map = constructMap[domain];
+                if (subField.size() != map.size())
+                {
+                    FatalErrorIn
+                    (
+                        "template<class T>\n"
+                        "void mapDistribute::distribute\n"
+                        "(\n"
+                        "    const Pstream::commsTypes commsType,\n"
+                        "    const List<labelPair>& schedule,\n"
+                        "    const label constructSize,\n"
+                        "    const labelListList& subMap,\n"
+                        "    const labelListList& constructMap,\n"
+                        "    List<T>& field\n"
+                        ")\n"
+                    )   << "Expected from processor " << domain
+                        << " " << map.size() << " but received "
+                        << subField.size() << " elements."
+                        << abort(FatalError);
+                }
 
                 forAll(map, i)
                 {
@@ -126,6 +142,7 @@ void Foam::mapDistribute::distribute
             newField[map[i]] = subField[i];
         }
 
+        // Schedule will already have pruned 0-sized comms
         forAll(schedule, i)
         {
             const labelPair& twoProcs = schedule[i];
@@ -154,6 +171,26 @@ void Foam::mapDistribute::distribute
 
                 const labelList& map = constructMap[sendProc];
 
+                if (subField.size() != map.size())
+                {
+                    FatalErrorIn
+                    (
+                        "template<class T>\n"
+                        "void mapDistribute::distribute\n"
+                        "(\n"
+                        "    const Pstream::commsTypes commsType,\n"
+                        "    const List<labelPair>& schedule,\n"
+                        "    const label constructSize,\n"
+                        "    const labelListList& subMap,\n"
+                        "    const labelListList& constructMap,\n"
+                        "    List<T>& field\n"
+                        ")\n"
+                    )   << "Expected from processor " << sendProc
+                        << " " << map.size() << " but received "
+                        << subField.size() << " elements."
+                        << abort(FatalError);
+                }
+
                 forAll(map, i)
                 {
                     newField[map[i]] = subField[i];
@@ -164,7 +201,74 @@ void Foam::mapDistribute::distribute
     }
     else if (commsType == Pstream::nonBlocking)
     {
-        List<T> newField(constructSize);
+        if (!contiguous<T>())
+        {
+            FatalErrorIn
+            (
+                "template<class T>\n"
+                "void mapDistribute::distribute\n"
+                "(\n"
+                "    const Pstream::commsTypes commsType,\n"
+                "    const List<labelPair>& schedule,\n"
+                "    const label constructSize,\n"
+                "    const labelListList& subMap,\n"
+                "    const labelListList& constructMap,\n"
+                "    List<T>& field\n"
+                ")\n"
+            )   << "Non-blocking only supported for contiguous data."
+                << exit(FatalError);
+        }
+
+        // Set up sends to neighbours
+
+        List<List<T > > sendFields(Pstream::nProcs());
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = subMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
+            {
+                List<T>& subField = sendFields[domain];
+                subField.setSize(map.size());
+                forAll(map, i)
+                {
+                    subField[i] = field[map[i]];
+                }
+
+                OPstream::write
+                (
+                    Pstream::nonBlocking,
+                    domain,
+                    reinterpret_cast<const char*>(subField.begin()),
+                    subField.size()*sizeof(T)
+                );
+            }
+        }
+
+        // Set up receives from neighbours
+
+        List<List<T > > recvFields(Pstream::nProcs());
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = constructMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
+            {
+                recvFields[domain].setSize(map.size());
+                IPstream::read
+                (
+                    Pstream::nonBlocking,
+                    domain,
+                    reinterpret_cast<char*>(recvFields[domain].begin()),
+                    recvFields[domain].size()*sizeof(T)
+                );
+            }
+        }
+
+
+        // Combine bits. Note that can reuse field storage
 
         // Subset myself
         const labelList& mySubMap = subMap[Pstream::myProcNo()];
@@ -175,54 +279,56 @@ void Foam::mapDistribute::distribute
             subField[i] = field[mySubMap[i]];
         }
 
+        field.setSize(constructSize);
+
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
 
         forAll(map, i)
         {
-            newField[map[i]] = subField[i];
+            field[map[i]] = subField[i];
-        }
-
-        // Send sub field to neighbour
-        for (label domain = 0; domain < Pstream::nProcs(); domain++)
-        {
-            const labelList& map = subMap[domain];
-
-            if (domain != Pstream::myProcNo() && map.size() > 0)
-            {
-
-                List<T> subField(map.size());
-                forAll(map, i)
-                {
-                    subField[i] = field[map[i]];
-                }
-
-                OPstream toNbr(Pstream::nonBlocking, domain);
-                toNbr << subField;
-            }
         }
 
 
-        // Receive sub field from neighbour
+        // Wait for all to finish
+
+        OPstream::waitRequests();
+        IPstream::waitRequests();
+        
+        // Collect neighbour fields
+
         for (label domain = 0; domain < Pstream::nProcs(); domain++)
         {
             const labelList& map = constructMap[domain];
 
             if (domain != Pstream::myProcNo() && map.size() > 0)
             {
-                IPstream fromNbr(Pstream::nonBlocking, domain);
+                if (recvFields[domain].size() != map.size())
-                List<T> subField(fromNbr);
+                {
+                    FatalErrorIn
+                    (
+                        "template<class T>\n"
+                        "void mapDistribute::distribute\n"
+                        "(\n"
+                        "    const Pstream::commsTypes commsType,\n"
+                        "    const List<labelPair>& schedule,\n"
+                        "    const label constructSize,\n"
+                        "    const labelListList& subMap,\n"
+                        "    const labelListList& constructMap,\n"
+                        "    List<T>& field\n"
+                        ")\n"
+                    )   << "Expected from processor " << domain
+                        << " " << map.size() << " but received "
+                        << recvFields[domain].size() << " elements."
+                        << abort(FatalError);
+                }
 
                 forAll(map, i)
                 {
-                    newField[map[i]] = subField[i];
+                    field[map[i]] = recvFields[domain][i];
                 }
             }
         }
-        OPstream::waitRequests();
-        IPstream::waitRequests();
-
-        field.transfer(newField);
     }
     else
     {

src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H

@@ -221,6 +221,41 @@ public:
                 const bool applySeparation
             );
 
+        // PackedList versions
+
+            template <int nBits, class CombineOp>
+            static void syncFaceList
+            (
+                const polyMesh& mesh,
+                PackedList<nBits>& faceValues,
+                const CombineOp& cop
+            );
+
+            template <int nBits>
+            static void swapFaceList
+            (
+                const polyMesh& mesh,
+                PackedList<nBits>& faceValues
+            );
+
+            template <int nBits, class CombineOp>
+            static void syncPointList
+            (
+                const polyMesh& mesh,
+                PackedList<nBits>& pointValues,
+                const CombineOp& cop,
+                const unsigned int nullValue
+            );
+
+            template <int nBits, class CombineOp>
+            static void syncEdgeList
+            (
+                const polyMesh& mesh,
+                PackedList<nBits>& edgeValues,
+                const CombineOp& cop,
+                const unsigned int nullValue
+            );
+
         // Other
 
             //- Get per point whether is it master (of a coupled set of points)

src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C

@@ -940,17 +940,7 @@ void Foam::syncTools::syncPointList
                     fromNbr >> nbrPatchInfo;
                 }
                 // Null any value which is not on neighbouring processor
-                label nbrSize = nbrPatchInfo.size();
+                nbrPatchInfo.setSize(procPatch.nPoints(), nullValue);
-                nbrPatchInfo.setSize(procPatch.nPoints());
-                for
-                (
-                    label i = nbrSize;
-                    i < procPatch.nPoints();
-                    i++
-                )
-                {
-                    nbrPatchInfo[i] = nullValue;
-                }
 
                 if (!procPatch.parallel())
                 {
@@ -1225,18 +1215,7 @@ void Foam::syncTools::syncEdgeList
                     fromNeighb >> nbrPatchInfo;
                 }
                 // Null any value which is not on neighbouring processor
-                label nbrSize = nbrPatchInfo.size();
+                nbrPatchInfo.setSize(procPatch.nEdges(), nullValue);
-                nbrPatchInfo.setSize(procPatch.nEdges());
-                for
-                (
-                    label i = nbrSize;
-                    i < procPatch.nEdges();
-                    i++
-                )
-                {
-                    nbrPatchInfo[i] = nullValue;
-                }
-
 
                 if (!procPatch.parallel())
                 {
@@ -1591,4 +1570,469 @@ void Foam::syncTools::swapFaceList
 }
 
 
+template <int nBits, class CombineOp>
+void Foam::syncTools::syncFaceList
+(
+    const polyMesh& mesh,
+    PackedList<nBits>& faceValues,
+    const CombineOp& cop
+)
+{
+    if (faceValues.size() != mesh.nFaces())
+    {
+        FatalErrorIn
+        (
+            "syncTools<int nBits, class CombineOp>::syncFaceList"
+            "(const polyMesh&, PackedList<nBits>&, const CombineOp&)"
+        )   << "Number of values " << faceValues.size()
+            << " is not equal to the number of faces in the mesh "
+            << mesh.nFaces() << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    if (!hasCouples(patches))
+    {
+        return;
+    }
+
+    // Patch data (proc patches only).
+    List<List<unsigned int> > patchValues(patches.size());
+
+    if (Pstream::parRun())
+    {
+        // Send
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].size() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                patchValues[patchI].setSize(procPatch.size());
+                forAll(procPatch, i)
+                {
+                    patchValues[patchI][i] =
+                        faceValues.get(procPatch.start()+i);
+                }
+
+                OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
+                toNbr << patchValues[patchI];
+            }
+        }
+
+
+        // Receive and combine.
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].size() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                {
+                    IPstream fromNbr
+                    (
+                        Pstream::blocking,
+                        procPatch.neighbProcNo()
+                    );
+                    fromNbr >> patchValues[patchI];
+                }
+
+                // Combine (bitwise)
+                forAll(procPatch, i)
+                {
+                    unsigned int patchVal = patchValues[patchI][i];
+                    label meshFaceI = procPatch.start()+i;
+                    unsigned int faceVal = faceValues.get(meshFaceI);
+                    cop(faceVal, patchVal);
+                    faceValues.set(meshFaceI, faceVal);
+                }
+            }
+        }
+    }
+
+    // Do the cyclics.
+    forAll(patches, patchI)
+    {
+        if (isA<cyclicPolyPatch>(patches[patchI]))
+        {
+            const cyclicPolyPatch& cycPatch =
+                refCast<const cyclicPolyPatch>(patches[patchI]);
+
+            label half = cycPatch.size()/2;
+
+            for (label i = 0; i < half; i++)
+            {
+                label meshFace0 = cycPatch.start()+i;
+                unsigned int val0 = faceValues.get(meshFace0);
+                label meshFace1 = meshFace0 + half;
+                unsigned int val1 = faceValues.get(meshFace1);
+
+                unsigned int t = val0;
+                cop(t, val1);
+                faceValues.set(meshFace0, t);
+
+                cop(val1, val0);                
+                faceValues.set(meshFace1, val1);
+            }
+        }
+    }
+}
+
+
+template <int nBits>
+void Foam::syncTools::swapFaceList
+(
+    const polyMesh& mesh,
+    PackedList<nBits>& faceValues
+)
+{
+    syncFaceList(mesh, faceValues, eqOp<unsigned int>());
+}
+
+
+template <int nBits, class CombineOp>
+void Foam::syncTools::syncPointList
+(
+    const polyMesh& mesh,
+    PackedList<nBits>& pointValues,
+    const CombineOp& cop,
+    const unsigned int nullValue
+)
+{
+    if (pointValues.size() != mesh.nPoints())
+    {
+        FatalErrorIn
+        (
+            "syncTools<int nBits, class CombineOp>::syncPointList"
+            "(const polyMesh&, PackedList<nBits>&, const CombineOp&"
+            ", const unsigned int&)"
+        )   << "Number of values " << pointValues.size()
+            << " is not equal to the number of points in the mesh "
+            << mesh.nPoints() << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    if (!hasCouples(patches))
+    {
+        return;
+    }
+
+    // Patch data (proc patches only).
+    List<List<unsigned int> > patchValues(patches.size());
+
+    if (Pstream::parRun())
+    {
+        // Send
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].nPoints() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                patchValues[patchI].setSize(procPatch.nPoints());
+                patchValues[patchI] = nullValue;
+
+                const labelList& meshPts = procPatch.meshPoints();
+                const labelList& nbrPts = procPatch.neighbPoints();
+
+                forAll(nbrPts, pointI)
+                {
+                    label nbrPointI = nbrPts[pointI];
+                    if (nbrPointI >= 0 && nbrPointI < procPatch.nPoints())
+                    {
+                        patchValues[patchI][nbrPointI] =
+                            pointValues.get(meshPts[pointI]);
+                    }
+                }
+
+                OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
+                toNbr << patchValues[patchI];
+            }
+        }
+
+
+        // Receive and combine.
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].nPoints() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                {
+                    // We do not know the number of points on the other side
+                    // so cannot use Pstream::read.
+                    IPstream fromNbr
+                    (
+                        Pstream::blocking,
+                        procPatch.neighbProcNo()
+                    );
+                    fromNbr >> patchValues[patchI];
+                }
+
+                // Null any value which is not on neighbouring processor
+                patchValues[patchI].setSize(procPatch.nPoints(), nullValue);
+
+                const labelList& meshPts = procPatch.meshPoints();
+
+                forAll(meshPts, pointI)
+                {
+                    label meshPointI = meshPts[pointI];
+                    unsigned int pointVal = pointValues.get(meshPointI);
+                    cop(pointVal, patchValues[patchI][pointI]);
+                    pointValues.set(meshPointI, pointVal);
+                }
+            }
+        }
+    }
+
+    // Do the cyclics.
+    forAll(patches, patchI)
+    {
+        if (isA<cyclicPolyPatch>(patches[patchI]))
+        {
+            const cyclicPolyPatch& cycPatch =
+                refCast<const cyclicPolyPatch>(patches[patchI]);
+
+            const edgeList& coupledPoints = cycPatch.coupledPoints();
+            const labelList& meshPts = cycPatch.meshPoints();
+
+            forAll(coupledPoints, i)
+            {
+                const edge& e = coupledPoints[i];
+
+                label point0 = meshPts[e[0]];
+                label point1 = meshPts[e[1]];
+
+                unsigned int val0 = pointValues.get(point0);
+                unsigned int t = val0;
+                unsigned int val1 = pointValues.get(point1);
+
+                cop(t, val1);
+                pointValues.set(point0, t);
+                cop(val1, val0);
+                pointValues.set(point1, val1);
+            }
+        }
+    }
+
+    // Synchronize multiple shared points.
+    const globalMeshData& pd = mesh.globalData();
+
+    if (pd.nGlobalPoints() > 0)
+    {
+        // Values on shared points. Use unpacked storage for ease!
+        List<unsigned int> sharedPts(pd.nGlobalPoints(), nullValue);
+
+        forAll(pd.sharedPointLabels(), i)
+        {
+            label meshPointI = pd.sharedPointLabels()[i];
+            // Fill my entries in the shared points
+            sharedPts[pd.sharedPointAddr()[i]] = pointValues.get(meshPointI);
+        }
+
+        // Combine on master.
+        Pstream::listCombineGather(sharedPts, cop);
+        Pstream::listCombineScatter(sharedPts);
+
+        // Now we will all have the same information. Merge it back with
+        // my local information.
+        forAll(pd.sharedPointLabels(), i)
+        {
+            label meshPointI = pd.sharedPointLabels()[i];
+            pointValues.set(meshPointI, sharedPts[pd.sharedPointAddr()[i]]);
+        }
+    }
+}
+
+
+template <int nBits, class CombineOp>
+void Foam::syncTools::syncEdgeList
+(
+    const polyMesh& mesh,
+    PackedList<nBits>& edgeValues,
+    const CombineOp& cop,
+    const unsigned int nullValue
+)
+{
+    if (edgeValues.size() != mesh.nEdges())
+    {
+        FatalErrorIn
+        (
+            "syncTools<int nBits, class CombineOp>::syncEdgeList"
+            "(const polyMesh&, PackedList<nBits>&, const CombineOp&"
+            ", const unsigned int&)"
+        )   << "Number of values " << edgeValues.size()
+            << " is not equal to the number of edges in the mesh "
+            << mesh.nEdges() << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    if (!hasCouples(patches))
+    {
+        return;
+    }
+
+    // Patch data (proc patches only).
+    List<List<unsigned int> > patchValues(patches.size());
+
+    if (Pstream::parRun())
+    {
+        // Send
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].nEdges() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                patchValues[patchI].setSize(procPatch.nEdges(), nullValue);
+
+                const labelList& meshEdges = procPatch.meshEdges();
+                const labelList& neighbEdges = procPatch.neighbEdges();
+
+                forAll(neighbEdges, edgeI)
+                {
+                    label nbrEdgeI = neighbEdges[edgeI];
+                    if (nbrEdgeI >= 0 && nbrEdgeI < procPatch.nEdges())
+                    {
+                        patchValues[patchI][nbrEdgeI] =
+                            edgeValues.get(meshEdges[edgeI]);
+                    }
+                }
+
+                OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
+                toNbr << patchValues[patchI];
+            }
+        }
+
+
+        // Receive and combine.
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].nEdges() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                {
+                    IPstream fromNeighb
+                    (
+                        Pstream::blocking,
+                        procPatch.neighbProcNo()
+                    );
+                    fromNeighb >> patchValues[patchI];
+                }
+
+                patchValues[patchI].setSize(procPatch.nEdges(), nullValue);
+
+                const labelList& meshEdges = procPatch.meshEdges();
+
+                forAll(meshEdges, edgeI)
+                {
+                    unsigned int patchVal = patchValues[patchI][edgeI];
+                    label meshEdgeI = meshEdges[edgeI];
+                    unsigned int edgeVal = edgeValues.get(meshEdgeI);
+                    cop(edgeVal, patchVal);
+                    edgeValues.set(meshEdgeI, edgeVal);
+                }
+            }
+        }
+    }
+
+    // Do the cyclics.
+    forAll(patches, patchI)
+    {
+        if (isA<cyclicPolyPatch>(patches[patchI]))
+        {
+            const cyclicPolyPatch& cycPatch =
+                refCast<const cyclicPolyPatch>(patches[patchI]);
+
+            const edgeList& coupledEdges = cycPatch.coupledEdges();
+            const labelList& meshEdges = cycPatch.meshEdges();
+
+            forAll(coupledEdges, i)
+            {
+                const edge& e = coupledEdges[i];
+
+                label edge0 = meshEdges[e[0]];
+                label edge1 = meshEdges[e[1]];
+
+                unsigned int val0 = edgeValues.get(edge0);
+                unsigned int t = val0;
+                unsigned int val1 = edgeValues.get(edge1);
+
+                cop(t, val1);
+                edgeValues.set(edge0, t);
+                cop(val1, val0);
+                edgeValues.set(edge1, val1);
+            }
+        }
+    }
+
+    // Synchronize multiple shared edges.
+    const globalMeshData& pd = mesh.globalData();
+
+    if (pd.nGlobalEdges() > 0)
+    {
+        // Values on shared edges. Use unpacked storage for ease!
+        List<unsigned int> sharedPts(pd.nGlobalEdges(), nullValue);
+
+        forAll(pd.sharedEdgeLabels(), i)
+        {
+            label meshEdgeI = pd.sharedEdgeLabels()[i];
+            // Fill my entries in the shared edges
+            sharedPts[pd.sharedEdgeAddr()[i]] = edgeValues.get(meshEdgeI);
+        }
+
+        // Combine on master.
+        Pstream::listCombineGather(sharedPts, cop);
+        Pstream::listCombineScatter(sharedPts);
+
+        // Now we will all have the same information. Merge it back with
+        // my local information.
+        forAll(pd.sharedEdgeLabels(), i)
+        {
+            label meshEdgeI = pd.sharedEdgeLabels()[i];
+            edgeValues.set(meshEdgeI, sharedPts[pd.sharedEdgeAddr()[i]]);
+        }
+    }
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/primitiveMesh/PrimitivePatchExtra/PrimitivePatchExtraSearch.C

@@ -62,7 +62,7 @@ void Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::markZone
 
             const labelList& fEdges = faceEs[faceI];
 
-            forAllfEdges, i)
+            forAll(fEdges, i)
             {
                 label edgeI = fEdges[i];
 
@@ -102,8 +102,8 @@ void Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::markZone
             break;
         }
 
-        changedFaces.transfer(newChangedFaces.shrink());
+        // New dynamicList: can leave dynamicList unshrunk
-        newChangedFaces.clear();
+        changedFaces.transfer(newChangedFaces);
     }
 }
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C

@@ -66,7 +66,6 @@ primitiveMesh::primitiveMesh()
     ppPtr_(NULL),
     cpPtr_(NULL),
 
-    allocSize_(0),
     labels_(0),
 
     cellCentresPtr_(NULL),
@@ -109,7 +108,6 @@ primitiveMesh::primitiveMesh
     ppPtr_(NULL),
     cpPtr_(NULL),
 
-    allocSize_(0),
     labels_(0),
 
     cellCentresPtr_(NULL),

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H

@@ -54,6 +54,7 @@ SourceFiles
 #ifndef primitiveMesh_H
 #define primitiveMesh_H
 
+#include "DynamicList.H"
 #include "edgeList.H"
 #include "pointField.H"
 #include "SubField.H"
@@ -157,10 +158,8 @@ class primitiveMesh
 
         // On-the-fly edge addresing storage
 
-            //- Temporary storage for addressing. allocSize is the real size
+            //- Temporary storage for addressing.
-            //  of the labelList.
+            mutable DynamicList<label> labels_;
-            mutable label allocSize_;
-            mutable labelList labels_;
 
             //- Temporary storage for addressing
             mutable labelHashSet labelSet_;
@@ -705,31 +704,80 @@ public:
 
             // On-the-fly addressing calculation. These functions return either
             // a reference to the full addressing (if already calculated) or
-            // a reference to member data labels_ so be careful when not storing
+            // a reference to the supplied storage. The one-argument ones
+            // use member DynamicList labels_ so be careful when not storing
             // result.
 
-            //- cellCells using cells
+            //- cellCells using cells.
+            const labelList& cellCells
+            (
+                const label cellI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& cellCells(const label cellI) const;
 
             //- cellPoints using cells
+            const labelList& cellPoints
+            (
+                const label cellI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& cellPoints(const label cellI) const;
 
             //- pointCells using pointFaces
+            const labelList& pointCells
+            (
+                const label pointI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& pointCells(const label pointI) const;
 
             //- pointPoints using edges, pointEdges
+            const labelList& pointPoints
+            (
+                const label pointI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& pointPoints(const label pointI) const;
 
             //- faceEdges using pointFaces, edges, pointEdges
+            const labelList& faceEdges
+            (
+                const label faceI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& faceEdges(const label faceI) const;
 
             //- edgeFaces using pointFaces, edges, pointEdges
+            const labelList& edgeFaces
+            (
+                const label edgeI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& edgeFaces(const label edgeI) const;
 
             //- edgeCells using pointFaces, edges, pointEdges
+            const labelList& edgeCells
+            (
+                const label edgeI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& edgeCells(const label edgeI) const;
 
             //- cellEdges using cells, pointFaces, edges, pointEdges
+            const labelList& cellEdges
+            (
+                const label cellI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& cellEdges(const label cellI) const;
 
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCells.C

@@ -41,6 +41,14 @@ void primitiveMesh::calcCellCells() const
     {
         Pout<< "primitiveMesh::calcCellCells() : calculating cellCells"
             << endl;
+
+        if (debug == -1)
+        {
+            // For checking calls:abort so we can quickly hunt down
+            // origin of call
+            FatalErrorIn("primitiveMesh::calcCellCells()")
+                << abort(FatalError);
+        }
     }
 
     // It is an error to attempt to recalculate cellCells
@@ -105,7 +113,11 @@ const labelListList& primitiveMesh::cellCells() const
 }
 
 
-const labelList& primitiveMesh::cellCells(const label cellI) const
+const labelList& primitiveMesh::cellCells
+(
+    const label cellI,
+    DynamicList<label>& storage
+) const
 {
     if (hasCellCells())
     {
@@ -117,16 +129,7 @@ const labelList& primitiveMesh::cellCells(const label cellI) const
         const labelList& nei = faceNeighbour();
         const cell& cFaces = cells()[cellI];
 
-        labels_.size() = allocSize_;
+        storage.clear();
-
-        if (cFaces.size() > allocSize_)
-        {
-            labels_.clear();
-            allocSize_ = cFaces.size();
-            labels_.setSize(allocSize_);
-        }
-
-        label n = 0;
 
         forAll(cFaces, i)
         {
@@ -136,22 +139,26 @@ const labelList& primitiveMesh::cellCells(const label cellI) const
             {
                 if (own[faceI] == cellI)
                 {
-                    labels_[n++] = nei[faceI];
+                    storage.append(nei[faceI]);
                 }
                 else
                 {
-                    labels_[n++] = own[faceI];
+                    storage.append(own[faceI]);
                 }
             }
         }
 
-        labels_.size() = n;
+        return storage;
-
-        return labels_;
     }
 }
 
 
+const labelList& primitiveMesh::cellCells(const label cellI) const
+{
+    return cellCells(cellI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellEdges.C

@@ -40,6 +40,14 @@ void Foam::primitiveMesh::calcCellEdges() const
         Pout<< "primitiveMesh::calcCellEdges() : "
             << "calculating cellEdges"
             << endl;
+
+        if (debug == -1)
+        {
+            // For checking calls:abort so we can quickly hunt down
+            // origin of call
+            FatalErrorIn("primitiveMesh::calcCellEdges()")
+                << abort(FatalError);
+        }
     }
 
     // It is an error to attempt to recalculate cellEdges

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellPoints.C

@@ -42,6 +42,14 @@ const labelListList& primitiveMesh::cellPoints() const
         {
             Pout<< "primitiveMesh::cellPoints() : "
                 << "calculating cellPoints" << endl;
+
+            if (debug == -1)
+            {
+                // For checking calls:abort so we can quickly hunt down
+                // origin of call
+                FatalErrorIn("primitiveMesh::cellPoints()")
+                    << abort(FatalError);
+            }
         }
 
         // Invert pointCells
@@ -53,7 +61,11 @@ const labelListList& primitiveMesh::cellPoints() const
 }
 
 
-const labelList& primitiveMesh::cellPoints(const label cellI) const
+const labelList& primitiveMesh::cellPoints
+(
+    const label cellI,
+    DynamicList<label>& storage
+) const
 {
     if (hasCellPoints())
     {
@@ -76,29 +88,28 @@ const labelList& primitiveMesh::cellPoints(const label cellI) const
             }
         }
 
-        labels_.size() = allocSize_;
+        storage.clear();
-
+        if (labelSet_.size() > storage.allocSize())
-        if (labelSet_.size() > allocSize_)
         {
-            labels_.clear();
+            storage.setSize(labelSet_.size());
-            allocSize_ = labelSet_.size();
-            labels_.setSize(allocSize_);
         }
 
-        label n = 0;
-
         forAllConstIter(labelHashSet, labelSet_, iter)
         {
-            labels_[n++] = iter.key();
+            storage.append(iter.key());
         }
 
-        labels_.size() = n;
+        return storage;
-
-        return labels_;
     }
 }
 
 
+const labelList& primitiveMesh::cellPoints(const label cellI) const
+{
+    return cellPoints(cellI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeCells.C

@@ -41,6 +41,14 @@ const labelListList& primitiveMesh::edgeCells() const
         if (debug)
         {
             Pout<< "primitiveMesh::edgeCells() : calculating edgeCells" << endl;
+
+            if (debug == -1)
+            {
+                // For checking calls:abort so we can quickly hunt down
+                // origin of call
+                FatalErrorIn("primitiveMesh::edgeCells()")
+                    << abort(FatalError);
+            }
         }
         // Invert cellEdges
         ecPtr_ = new labelListList(nEdges());
@@ -51,7 +59,11 @@ const labelListList& primitiveMesh::edgeCells() const
 }
 
 
-const labelList& primitiveMesh::edgeCells(const label edgeI) const
+const labelList& primitiveMesh::edgeCells
+(
+    const label edgeI,
+    DynamicList<label>& storage
+) const
 {
     if (hasEdgeCells())
     {
@@ -62,24 +74,11 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
         const labelList& own = faceOwner();
         const labelList& nei = faceNeighbour();
 
-        // edge faces can either return labels_ or reference in edgeLabels.
+        // Construct edgeFaces
-        labelList labelsCopy;
+        DynamicList<label> eFacesStorage;
-        if (!hasEdgeFaces())
+        const labelList& eFaces = edgeFaces(edgeI, eFacesStorage);
-        {
-            labelsCopy = edgeFaces(edgeI);
-        }
 
-        const labelList& eFaces =
+        storage.clear();
-        (
-            hasEdgeFaces()
-          ? edgeFaces()[edgeI]
-          : labelsCopy
-        );
-
-        labels_.size() = allocSize_;
-
-        // labels_ should certainly be big enough for edge cells.
-        label n = 0;
 
         // Do quadratic insertion.
         forAll(eFaces, i)
@@ -89,10 +88,10 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
             {
                 label ownCellI = own[faceI];
 
-                // Check if not already in labels_
+                // Check if not already in storage
-                for (label j = 0; j < n; j++)
+                forAll(storage, j)
                 {
-                    if (labels_[j] == ownCellI)
+                    if (storage[j] == ownCellI)
                     {
                         ownCellI = -1;
                         break;
@@ -101,7 +100,7 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
 
                 if (ownCellI != -1)
                 {
-                    labels_[n++] = ownCellI;
+                    storage.append(ownCellI);
                 }
             }
 
@@ -109,9 +108,9 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
             {
                 label neiCellI = nei[faceI];
 
-                for (label j = 0; j < n; j++)
+                forAll(storage, j)
                 {
-                    if (labels_[j] == neiCellI)
+                    if (storage[j] == neiCellI)
                     {
                         neiCellI = -1;
                         break;
@@ -120,18 +119,22 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
 
                 if (neiCellI != -1)
                 {
-                    labels_[n++] = neiCellI;
+                    storage.append(neiCellI);
                 }
             }
         }
 
-        labels_.size() = n;
+        return storage;
-
-        return labels_;
     }
 }
 
 
+const labelList& primitiveMesh::edgeCells(const label edgeI) const
+{
+    return edgeCells(edgeI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeFaces.C

@@ -41,6 +41,14 @@ const labelListList& primitiveMesh::edgeFaces() const
         if (debug)
         {
             Pout<< "primitiveMesh::edgeFaces() : calculating edgeFaces" << endl;
+
+            if (debug == -1)
+            {
+                // For checking calls:abort so we can quickly hunt down
+                // origin of call
+                FatalErrorIn("primitiveMesh::edgeFaces()")
+                    << abort(FatalError);
+            }
         }
 
         // Invert faceEdges
@@ -52,7 +60,11 @@ const labelListList& primitiveMesh::edgeFaces() const
 }
 
 
-const labelList& primitiveMesh::edgeFaces(const label edgeI) const
+const labelList& primitiveMesh::edgeFaces
+(
+    const label edgeI,
+    DynamicList<label>& storage
+) const
 {
     if (hasEdgeFaces())
     {
@@ -67,9 +79,8 @@ const labelList& primitiveMesh::edgeFaces(const label edgeI) const
 
         label i0 = 0;
         label i1 = 0;
-        label n = 0;
 
-        labels_.size() = allocSize_;
+        storage.clear();
 
         while (i0 < pFaces0.size() && i1 < pFaces1.size())
         {
@@ -84,26 +95,23 @@ const labelList& primitiveMesh::edgeFaces(const label edgeI) const
             else
             {
                 // Equal. Append.
-                if (n == allocSize_)
+                storage.append(pFaces0[i0]);
-                {
-                    // Have setSize copy contents so far
-                    labels_.size() = n;
-                    allocSize_ = allocSize_*2 + 1;
-                    labels_.setSize(allocSize_);
-                }
-                labels_[n++] = pFaces0[i0];
                 ++i0;
                 ++i1;
             }
         }
 
-        labels_.size() = n;
+        return storage;
-
-        return labels_;
     }
 }
 
 
+const labelList& primitiveMesh::edgeFaces(const label edgeI) const
+{
+    return edgeFaces(edgeI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C

@@ -441,11 +441,10 @@ void primitiveMesh::calcEdges(const bool doFaceEdges) const
         forAll(pe, pointI)
         {
             DynamicList<label>& pEdges = pe[pointI];
-            inplaceRenumber(oldToNew, pEdges);
             pEdges.shrink();
+            inplaceRenumber(oldToNew, pEdges);
             pointEdges[pointI].transfer(pEdges);
             Foam::sort(pointEdges[pointI]);
-            pEdges.clear();
         }
 
         // faceEdges
@@ -583,11 +582,15 @@ void primitiveMesh::clearOutEdges()
     deleteDemandDrivenData(pePtr_);
     deleteDemandDrivenData(fePtr_);
     labels_.clear();
-    allocSize_ = 0;
+    labelSet_.clear();
 }
 
 
-const labelList& primitiveMesh::faceEdges(const label faceI) const
+const labelList& primitiveMesh::faceEdges
+(
+    const label faceI,
+    DynamicList<label>& storage
+) const
 {
     if (hasFaceEdges())
     {
@@ -598,34 +601,40 @@ const labelList& primitiveMesh::faceEdges(const label faceI) const
         const labelListList& pointEs = pointEdges();
         const face& f = faces()[faceI];
 
-        labels_.size() = allocSize_;
+        storage.clear();
-
+        if (f.size() > storage.allocSize())
-        if (f.size() > allocSize_)
         {
-            labels_.clear();
+            storage.setSize(f.size());
-            allocSize_ = f.size();
-            labels_.setSize(allocSize_);
         }
 
-        label n = 0;
-
         forAll(f, fp)
         {
-            labels_[n++] = findFirstCommonElementFromSortedLists
+            storage.append
             (
-                pointEs[f[fp]],
+                findFirstCommonElementFromSortedLists
-                pointEs[f.nextLabel(fp)]
+                (
+                    pointEs[f[fp]],
+                    pointEs[f.nextLabel(fp)]
+                )
             );
         }
 
-        labels_.size() = n;
+        return storage;
-
-        return labels_;
     }
 }
 
 
-const labelList& primitiveMesh::cellEdges(const label cellI) const
+const labelList& primitiveMesh::faceEdges(const label faceI) const
+{
+    return faceEdges(faceI, labels_);
+}
+
+
+const labelList& primitiveMesh::cellEdges
+(
+    const label cellI,
+    DynamicList<label>& storage
+) const
 {
     if (hasCellEdges())
     {
@@ -647,29 +656,29 @@ const labelList& primitiveMesh::cellEdges(const label cellI) const
             }
         }
 
-        labels_.size() = allocSize_;
+        storage.clear();
 
-        if (labelSet_.size() > allocSize_)
+        if (labelSet_.size() > storage.allocSize())
         {
-            labels_.clear();
+            storage.setSize(labelSet_.size());
-            allocSize_ = labelSet_.size();
-            labels_.setSize(allocSize_);
         }
 
-        label n =0;
-
         forAllConstIter(labelHashSet, labelSet_, iter)
         {
-            labels_[n++] = iter.key();
+            storage.append(iter.key());
         }
 
-        labels_.size() = n;
+        return storage;
-
-        return labels_;
     }
 }
 
 
+const labelList& primitiveMesh::cellEdges(const label cellI) const
+{
+    return cellEdges(cellI, labels_);;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointCells.C

@@ -43,6 +43,14 @@ void primitiveMesh::calcPointCells() const
         Pout<< "primitiveMesh::calcPointCells() : "
             << "calculating pointCells"
             << endl;
+
+        if (debug == -1)
+        {
+            // For checking calls:abort so we can quickly hunt down
+            // origin of call
+            FatalErrorIn("primitiveMesh::calcPointCells()")
+                << abort(FatalError);
+        }
     }
 
     // It is an error to attempt to recalculate pointCells
@@ -114,7 +122,11 @@ const labelListList& primitiveMesh::pointCells() const
 }
 
 
-const labelList& primitiveMesh::pointCells(const label pointI) const
+const labelList& primitiveMesh::pointCells
+(
+    const label pointI,
+    DynamicList<label>& storage
+) const
 {
     if (hasPointCells())
     {
@@ -126,58 +138,48 @@ const labelList& primitiveMesh::pointCells(const label pointI) const
         const labelList& nei = faceNeighbour();
         const labelList& pFaces = pointFaces()[pointI];
 
-        labels_.size() = allocSize_;
+        storage.clear();
-
-        label n = 0;
 
         forAll(pFaces, i)
         {
             const label faceI = pFaces[i];
 
             // Append owner
-            if (n == allocSize_)
+            storage.append(own[faceI]);
-            {
-                labels_.size() = n;
-                allocSize_ = allocSize_*2 + 1;
-                labels_.setSize(allocSize_);
-            }
-            labels_[n++] = own[faceI];
 
             // Append neighbour
             if (faceI < nInternalFaces())
             {
-                if (n == allocSize_)
+                storage.append(nei[faceI]);
-                {
-                    labels_.size() = n;
-                    allocSize_ = allocSize_*2 + 1;
-                    labels_.setSize(allocSize_);
-                }
-                labels_[n++] = nei[faceI];
             }
         }
-        labels_.size() = n;
-
 
         // Filter duplicates
-        sort(labels_);
+        sort(storage);
 
-        n = 1;
+        label n = 1;
 
-        for (label i = 1; i < labels_.size(); i++)
+        for (label i = 1; i < storage.size(); i++)
         {
-            if (labels_[i] != labels_[i-1])
+            if (storage[i] != storage[i-1])
             {
-                labels_[n++] = labels_[i];
+                storage[n++] = storage[i];
             }
         }
 
-        labels_.size() = n;
+        storage.size() = n;
 
-        return labels_;
+        return storage;
     }
 }
 
 
+const labelList& primitiveMesh::pointCells(const label pointI) const
+{
+    return pointCells(pointI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointPoints.C

@@ -40,6 +40,14 @@ void primitiveMesh::calcPointPoints() const
         Pout<< "primitiveMesh::calcPointPoints() : "
             << "calculating pointPoints"
             << endl;
+
+        if (debug == -1)
+        {
+            // For checking calls:abort so we can quickly hunt down
+            // origin of call
+            FatalErrorIn("primitiveMesh::calcPointPoints()")
+                << abort(FatalError);
+        }
     }
 
     // It is an error to attempt to recalculate pointPoints
@@ -97,7 +105,11 @@ const labelListList& primitiveMesh::pointPoints() const
 }
 
 
-const labelList& primitiveMesh::pointPoints(const label pointI) const
+const labelList& primitiveMesh::pointPoints
+(
+    const label pointI,
+    DynamicList<label>& storage
+) const
 {
     if (hasPointPoints())
     {
@@ -108,30 +120,29 @@ const labelList& primitiveMesh::pointPoints(const label pointI) const
         const edgeList& edges = this->edges();
         const labelList& pEdges = pointEdges()[pointI];
 
-        labels_.size() = allocSize_;
+        storage.clear();
 
-        if (pEdges.size() > allocSize_)
+        if (pEdges.size() > storage.allocSize())
         {
-            // Set size() so memory allocation behaves as normal.
+            storage.setSize(pEdges.size());
-            labels_.clear();
-            allocSize_ = pEdges.size();
-            labels_.setSize(allocSize_);
         }
 
-        label n = 0;
-
         forAll(pEdges, i)
         {
-            labels_[n++] = edges[pEdges[i]].otherVertex(pointI);
+            storage.append(edges[pEdges[i]].otherVertex(pointI));
         }
 
-        labels_.size() = n;
+        return storage;
-
-        return labels_;
     }
 }
 
 
+const labelList& primitiveMesh::pointPoints(const label pointI) const
+{
+    return pointPoints(pointI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitivePatch/primitivePatch.H

@@ -26,7 +26,7 @@ Typedef
     Foam::primitivePatch
 
 Description
-    Foam::primitivePatch
+    Addressing for a faceList slice.
 
 \*---------------------------------------------------------------------------*/
 

src/OpenFOAM/meshes/primitiveShapes/line/linePointRef.H

@@ -26,6 +26,7 @@ Typedef
     Foam::linePointRef
 
 Description
+    Line using referred points
 
 \*---------------------------------------------------------------------------*/
 

src/OpenFOAM/primitives/contiguous/contiguous.H

@@ -38,6 +38,8 @@ Description
 #ifndef contiguous_H
 #define contiguous_H
 
+#include "label.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -45,54 +47,117 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Forward declaration of friend functions and operators
+template<class T, label Size> class FixedList;
+template<class T> class Pair;
+
+
 // Assume the data associated with type T is not contiguous
 template<class T>
 inline bool contiguous()                        {return false;}
 
 
-// Specify data associated with primitive types is contiguous
+
+// Specify data associated with primitive types (and simple fixed size
+// containers - only size 2 defined here) is contiguous
 
 template<>
-inline bool contiguous<bool>()                  {return true;}
+inline bool contiguous<bool>()                              {return true;}
+template<>
+inline bool contiguous<FixedList<bool, 2> >()               {return true;}
+template<>
+inline bool contiguous<Pair<bool> >()                       {return true;}
 
 template<>
-inline bool contiguous<char>()                  {return true;}
+inline bool contiguous<char>()                              {return true;}
+template<>
+inline bool contiguous<FixedList<char, 2> >()               {return true;}
+template<>
+inline bool contiguous<Pair<char> >()                       {return true;}
 
 template<>
-inline bool contiguous<unsigned char>()         {return true;}
+inline bool contiguous<unsigned char>()                     {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned char, 2> >()      {return true;}
+template<>
+inline bool contiguous<Pair<unsigned char> >()              {return true;}
 
 template<>
-inline bool contiguous<short>()                 {return true;}
+inline bool contiguous<short>()                             {return true;}
+template<>
+inline bool contiguous<FixedList<short, 2> >()              {return true;}
+template<>
+inline bool contiguous<Pair<short> >()                      {return true;}
 
 template<>
-inline bool contiguous<unsigned short>()        {return true;}
+inline bool contiguous<unsigned short>()                    {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned short, 2> >()     {return true;}
+template<>
+inline bool contiguous<Pair<unsigned short> >()             {return true;}
 
 template<>
-inline bool contiguous<int>()                   {return true;}
+inline bool contiguous<int>()                               {return true;}
+template<>
+inline bool contiguous<FixedList<int, 2> >()                {return true;}
+template<>
+inline bool contiguous<Pair<int> >()                        {return true;}
 
 template<>
-inline bool contiguous<unsigned int>()          {return true;}
+inline bool contiguous<unsigned int>()                      {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned int, 2> >()       {return true;}
+template<>
+inline bool contiguous<Pair<unsigned int> >()               {return true;}
 
 template<>
-inline bool contiguous<long>()                  {return true;}
+inline bool contiguous<long>()                              {return true;}
+template<>
+inline bool contiguous<FixedList<long, 2> >()               {return true;}
+template<>
+inline bool contiguous<Pair<long> >()                       {return true;}
 
 template<>
-inline bool contiguous<unsigned long>()         {return true;}
+inline bool contiguous<unsigned long>()                     {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned long, 2> >()      {return true;}
+template<>
+inline bool contiguous<Pair<unsigned long> >()              {return true;}
 
 template<>
-inline bool contiguous<long long>()             {return true;}
+inline bool contiguous<long long>()                         {return true;}
+template<>
+inline bool contiguous<FixedList<long long, 2> >()          {return true;}
+template<>
+inline bool contiguous<Pair<long long> >()                  {return true;}
 
 template<>
-inline bool contiguous<unsigned long long>()    {return true;}
+inline bool contiguous<unsigned long long>()                {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned long long, 2> >() {return true;}
+template<>
+inline bool contiguous<Pair<unsigned long long> >()         {return true;}
 
 template<>
-inline bool contiguous<float>()                 {return true;}
+inline bool contiguous<float>()                             {return true;}
+template<>
+inline bool contiguous<FixedList<float, 2> >()              {return true;}
+template<>
+inline bool contiguous<Pair<float> >()                       {return true;}
 
 template<>
-inline bool contiguous<double>()                {return true;}
+inline bool contiguous<double>()                            {return true;}
+template<>
+inline bool contiguous<FixedList<double, 2> >()             {return true;}
+template<>
+inline bool contiguous<Pair<double> >()                     {return true;}
 
 template<>
-inline bool contiguous<long double>()           {return true;}
+inline bool contiguous<long double>()                       {return true;}
+template<>
+inline bool contiguous<FixedList<long double, 2> >()        {return true;}
+template<>
+inline bool contiguous<Pair<long double> >()                {return true;}
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/Pstream/Allwmake

@@ -1,4 +1,5 @@
 #!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake libso dummy
@@ -18,3 +19,4 @@ GAMMA)
    ;;
 esac
 
+# ----------------------------------------------------------------- end-of-file

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C

@@ -25,6 +25,7 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "autoRefineDriver.H"
+#include "meshRefinement.H"
 #include "fvMesh.H"
 #include "Time.H"
 #include "boundBox.H"

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H

@@ -35,7 +35,7 @@ SourceFiles
 #ifndef autoRefineDriver_H
 #define autoRefineDriver_H
 
-#include "meshRefinement.H"
+#include "treeBoundBox.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,6 +45,9 @@ namespace Foam
 // Forward declaration of classes
 class featureEdgeMesh;
 class refinementParameters;
+class meshRefinement;
+class decompositionMethod;
+class fvMeshDistribute;
 
 /*---------------------------------------------------------------------------*\
                            Class autoRefineDriver Declaration

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C

@@ -1059,6 +1059,8 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
         // Surfaces with zone information
         const wordList& faceZoneNames = surfaces.faceZoneNames();
 
+        scalarField minSnapDist(snapDist);
+
         forAll(zonedSurfaces, i)
         {
             label zoneSurfI = zonedSurfaces[i];
@@ -1075,31 +1077,33 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
                 )
             );
 
-            pointField zonePoints(zonePointIndices.size());
-            forAll(zonePointIndices, i)
-            {
-                zonePoints[i] = localPoints[zonePointIndices[i]];
-            }
-
             // Find nearest for points both on faceZone and pp.
             List<pointIndexHit> hitInfo;
             labelList hitSurface;
             surfaces.findNearest
             (
                 labelList(1, zoneSurfI),
-                zonePoints,
+                pointField(localPoints, zonePointIndices),
-                sqr(4*snapDist),
+                sqr(4*scalarField(minSnapDist, zonePointIndices)),
                 hitSurface,
                 hitInfo
             );
 
-            forAll(hitInfo, pointI)
+            forAll(hitInfo, i)
             {
-                if (hitInfo[pointI].hit())
+                label pointI = zonePointIndices[i]; 
+
+                if (hitInfo[i].hit())
                 {
                     patchDisp[pointI] =
-                        hitInfo[pointI].hitPoint()
+                        hitInfo[i].hitPoint()
                       - localPoints[pointI];
+
+                    minSnapDist[pointI] = min
+                    (
+                        minSnapDist[pointI],
+                        mag(patchDisp[pointI])
+                    );
                 }
                 else
                 {
@@ -1107,7 +1111,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
                         << "For point:" << pointI
                         << " coordinate:" << localPoints[pointI]
                         << " did not find any surface within:"
-                        << 4*snapDist[pointI]
+                        << 4*minSnapDist[pointI]
                         << " meter." << endl;
                 }
             }

src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -51,6 +51,9 @@ License
 #include "meshTools.H"
 #include "OFstream.H"
 #include "geomDecomp.H"
+#include "Random.H"
+#include "searchableSurfaces.H"
+#include "treeBoundBox.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -122,13 +125,36 @@ void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
 {
     const pointField& cellCentres = mesh_.cellCentres();
 
-    label nTotEdges = returnReduce(surfaceIndex_.size(), sumOp<label>());
+    // Stats on edges to test. Count proc faces only once.
-    label nChangedFaces = returnReduce(changedFaces.size(), sumOp<label>());
+    PackedList<1> isMasterFace(syncTools::getMasterFaces(mesh_));
+
+    {
+        label nMasterFaces = 0;
+        forAll(isMasterFace, faceI)
+        {
+            if (isMasterFace.get(faceI) == 1)
+            {
+                nMasterFaces++;
+            }
+        }
+        reduce(nMasterFaces, sumOp<label>());
+
+        label nChangedFaces = 0;
+        forAll(changedFaces, i)
+        {
+            if (isMasterFace.get(changedFaces[i]) == 1)
+            {
+                nChangedFaces++;
+            }
+        }
+        reduce(nChangedFaces, sumOp<label>());
+
+        Info<< "Edge intersection testing:" << nl
+            << "    Number of edges             : " << nMasterFaces << nl
+            << "    Number of edges to retest   : " << nChangedFaces
+            << endl;
+    }
 
-    Info<< "Edge intersection testing:" << nl
-        << "    Number of edges             : " << nTotEdges << nl
-        << "    Number of edges to retest   : " << nChangedFaces
-        << endl;
 
     // Get boundary face centre and level. Coupled aware.
     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
@@ -838,11 +864,14 @@ Foam::meshRefinement::meshRefinement
 
 Foam::label Foam::meshRefinement::countHits() const
 {
+    // Stats on edges to test. Count proc faces only once.
+    PackedList<1> isMasterFace(syncTools::getMasterFaces(mesh_));
+
     label nHits = 0;
 
     forAll(surfaceIndex_, faceI)
     {
-        if (surfaceIndex_[faceI] >= 0)
+        if (surfaceIndex_[faceI] >= 0 && isMasterFace.get(faceI) == 1)
         {
             nHits++;
         }
@@ -996,11 +1025,6 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 
     if (Pstream::parRun())
     {
-        //Info<< nl
-        //    << "Doing final balancing" << nl
-        //    << "---------------------" << nl
-        //    << endl;
-        //
         //if (debug_)
         //{
         //    const_cast<Time&>(mesh_.time())++;
@@ -1014,17 +1038,12 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
             // Faces where owner and neighbour are not 'connected' so can
             // go to different processors.
             boolList blockedFace(mesh_.nFaces(), true);
+            label nUnblocked = 0;
             // Pairs of baffles
             List<labelPair> couples;
 
             if (keepZoneFaces)
             {
-                label nNamed = surfaces().getNamedSurfaces().size();
-
-                Info<< "Found " << nNamed << " surfaces with faceZones."
-                    << " Applying special decomposition to keep those together."
-                    << endl;
-
                 // Determine decomposition to keep/move surface zones
                 // on one processor. The reason is that snapping will make these
                 // into baffles, move and convert them back so if they were
@@ -1039,8 +1058,6 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                 // Get faces whose owner and neighbour should stay together,
                 // i.e. they are not 'blocked'.
 
-                label nZoned = 0;
-
                 forAll(fzNames, surfI)
                 {
                     if (fzNames[surfI].size() > 0)
@@ -1055,14 +1072,18 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                             if (blockedFace[fZone[i]])
                             {
                                 blockedFace[fZone[i]] = false;
-                                nZoned++;
+                                nUnblocked++;
                             }
                         }
                     }
                 }
-                Info<< "Found " << returnReduce(nZoned, sumOp<label>())
+            }
-                    << " zoned faces to keep together."
+            reduce(nUnblocked, sumOp<label>());
-                    << endl;
+
+            if (keepZoneFaces)
+            {
+                Info<< "Found " << nUnblocked
+                    << " zoned faces to keep together." << endl;
             }
 
             if (keepBaffles)
@@ -1073,29 +1094,43 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                     identity(mesh_.nFaces()-mesh_.nInternalFaces())
                    +mesh_.nInternalFaces()
                 );
+            }
+            label nCouples = returnReduce(couples.size(), sumOp<label>());
 
-                Info<< "Found " << returnReduce(couples.size(), sumOp<label>())
+            if (keepBaffles)
-                    << " baffles to keep together."
+            {
+                Info<< "Found " << nCouples << " baffles to keep together."
                     << endl;
             }
 
-            distribution = decomposeCombineRegions
+            if (nUnblocked > 0 || nCouples > 0)
-            (
-                blockedFace,
-                couples,
-                decomposer
-            );
-
-            labelList nProcCells(distributor.countCells(distribution));
-            Pstream::listCombineGather(nProcCells, plusEqOp<label>());
-            Pstream::listCombineScatter(nProcCells);
-
-            Info<< "Calculated decomposition:" << endl;
-            forAll(nProcCells, procI)
             {
-                Info<< "    " << procI << '\t' << nProcCells[procI] << endl;
+                Info<< "Applying special decomposition to keep baffles"
+                    << " and zoned faces together." << endl;
+
+                distribution = decomposeCombineRegions
+                (
+                    blockedFace,
+                    couples,
+                    decomposer
+                );
+
+                labelList nProcCells(distributor.countCells(distribution));
+                Pstream::listCombineGather(nProcCells, plusEqOp<label>());
+                Pstream::listCombineScatter(nProcCells);
+
+                Info<< "Calculated decomposition:" << endl;
+                forAll(nProcCells, procI)
+                {
+                    Info<< "    " << procI << '\t' << nProcCells[procI] << endl;
+                }
+                Info<< endl;
+            }
+            else
+            {
+                // Normal decomposition
+                distribution = decomposer.decompose(mesh_.cellCentres());
             }
-            Info<< endl;
         }
         else
         {
@@ -1715,6 +1750,36 @@ void Foam::meshRefinement::distribute(const mapDistributePolyMesh& map)
     {
         map.distributeFaceData(userFaceData_[i].second());
     }
+
+    // Redistribute surface and any fields on it.
+    {
+        Random rndGen(653213);
+
+        // Get local mesh bounding box. Single box for now.
+        List<treeBoundBox> meshBb(1);
+        treeBoundBox& bb = meshBb[0];
+        bb = boundBox(mesh_.points(), false);
+        bb = bb.extend(rndGen, 1E-4);
+
+        // Distribute all geometry (so refinementSurfaces and shellSurfaces)
+        searchableSurfaces& geometry =
+            const_cast<searchableSurfaces&>(surfaces_.geometry());
+
+        forAll(geometry, i)
+        {
+            autoPtr<mapDistribute> faceMap;
+            autoPtr<mapDistribute> pointMap;
+            geometry[i].distribute
+            (
+                meshBb,
+                false,          // do not keep outside triangles
+                faceMap,
+                pointMap
+            );
+            faceMap.clear();
+            pointMap.clear();
+        }
+    }
 }
 
 

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -386,7 +386,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
         FatalErrorIn
         (
             "meshRefinement::createBaffles"
-            "(const label, const labelList&, const labelList&)"
+            "(const labelList&, const labelList&)"
         )   << "Illegal size :"
             << " ownPatch:" << ownPatch.size()
             << " neiPatch:" << neiPatch.size()
@@ -412,7 +412,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
                 FatalErrorIn
                 (
                     "meshRefinement::createBaffles"
-                    "(const label, const labelList&, const labelList&)"
+                    "(const labelList&, const labelList&)"
                 )   << "Non synchronised at face:" << faceI
                     << " on patch:" << mesh_.boundaryMesh().whichPatch(faceI)
                     << " fc:" << mesh_.faceCentres()[faceI] << endl
@@ -461,7 +461,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
 
     //- Redo the intersections on the newly create baffle faces. Note that
     //  this changes also the cell centre positions.
-    labelHashSet baffledFacesSet(2*nBaffles);
+    faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles);
 
     const labelList& reverseFaceMap = map().reverseFaceMap();
     const labelList& faceMap = map().faceMap();
@@ -496,6 +496,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
             }
         }
     }
+    baffledFacesSet.sync(mesh_);
 
     updateMesh(map, baffledFacesSet.toc());
 
@@ -513,7 +514,7 @@ void Foam::meshRefinement::markBoundaryFace
 {
     isBoundaryFace[faceI] = true;
 
-    const labelList& fEdges = mesh_.faceEdges()[faceI];
+    const labelList& fEdges = mesh_.faceEdges(faceI);
 
     forAll(fEdges, fp)
     {
@@ -622,12 +623,16 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     // If so what is the remaining non-boundary anchor point?
     labelHashSet nonBoundaryAnchors(mesh_.nCells()/10000);
 
+    // On-the-fly addressing storage.
+    DynamicList<label> dynFEdges;
+    DynamicList<label> dynCPoints;
+
     // Count of faces marked for baffling
     label nBaffleFaces = 0;
 
     forAll(cellLevel, cellI)
     {
-        const labelList cPoints(meshCutter_.cellPoints(cellI));
+        const labelList& cPoints = mesh_.cellPoints(cellI, dynCPoints);
 
         // Get number of anchor points (pointLevel == cellLevel)
 
@@ -713,11 +718,14 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     // Loop over all points. If a point is connected to 4 or more cells
     // with 7 anchor points on the boundary set those cell's non-boundary faces
     // to baffles
+
+    DynamicList<label> dynPCells;
+
     forAllConstIter(labelHashSet, nonBoundaryAnchors, iter)
     {
         label pointI = iter.key();
 
-        const labelList& pCells = mesh_.pointCells()[pointI];
+        const labelList& pCells = mesh_.pointCells(pointI, dynPCells);
 
         // Count number of 'hasSevenBoundaryAnchorPoints' cells.
         label n = 0;
@@ -805,7 +813,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     {
         if (facePatch[faceI] == -1)
         {
-            const labelList& fEdges = mesh_.faceEdges()[faceI];
+            const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
             label nFaceBoundaryEdges = 0;
 
             forAll(fEdges, fe)
@@ -839,7 +847,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
             {
                 if (facePatch[faceI] == -1)
                 {
-                    const labelList& fEdges = mesh_.faceEdges()[faceI];
+                    const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
                     label nFaceBoundaryEdges = 0;
 
                     forAll(fEdges, fe)
@@ -1238,6 +1246,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
     labelList nBafflesPerEdge(mesh_.nEdges(), 0);
 
 
+
     // Count number of boundary faces per edge
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -1254,7 +1263,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
 
             forAll(pp, i)
             {
-                const labelList& fEdges = mesh_.faceEdges()[faceI];
+                const labelList& fEdges = mesh_.faceEdges(faceI);
 
                 forAll(fEdges, fEdgeI)
                 {
@@ -1266,19 +1275,23 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
     }
 
 
+    DynamicList<label> fe0;
+    DynamicList<label> fe1;
+
+
     // Count number of duplicate boundary faces per edge
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     forAll(couples, i)
     {
-        const labelList& fEdges0 = mesh_.faceEdges()[couples[i].first()];
+        const labelList& fEdges0 = mesh_.faceEdges(couples[i].first(), fe0);
 
         forAll(fEdges0, fEdgeI)
         {
             nBafflesPerEdge[fEdges0[fEdgeI]] += 1000000;
         }
 
-        const labelList& fEdges1 = mesh_.faceEdges()[couples[i].second()];
+        const labelList& fEdges1 = mesh_.faceEdges(couples[i].second(), fe1);
 
         forAll(fEdges1, fEdgeI)
         {
@@ -1313,7 +1326,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
          == patches.whichPatch(couple.second())
         )
         {
-            const labelList& fEdges = mesh_.faceEdges()[couples[i].first()];
+            const labelList& fEdges = mesh_.faceEdges(couples[i].first());
 
             forAll(fEdges, fEdgeI)
             {
@@ -1886,15 +1899,15 @@ void Foam::meshRefinement::baffleAndSplitMesh
 
         labelList facePatch
         (
-            //markFacesOnProblemCells
+            markFacesOnProblemCells
-            //(
-            //    globalToPatch
-            //)
-            markFacesOnProblemCellsGeometric
             (
-                motionDict,
                 globalToPatch
             )
+            //markFacesOnProblemCellsGeometric
+            //(
+            //    motionDict,
+            //    globalToPatch
+            //)
         );
         Info<< "Analyzed problem cells in = "
             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C

@@ -1091,8 +1091,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
         //    << "    local allowable refinement:" << nAllowRefine << nl
         //    << endl;
 
-        //- Disable refinement shortcut
+        //- Disable refinement shortcut. nAllowRefine is per processor limit.
-        label nAllowRefine = labelMax;
+        label nAllowRefine = labelMax / Pstream::nProcs();
 
         // Marked for refinement (>= 0) or not (-1). Actual value is the
         // index of the surface it intersects.

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -25,19 +25,13 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "refinementSurfaces.H"
-#include "orientedSurface.H"
 #include "Time.H"
 #include "searchableSurfaces.H"
 #include "shellSurfaces.H"
 #include "triSurfaceMesh.H"
 #include "labelPair.H"
 #include "searchableSurfacesQueries.H"
-
+#include "UPtrList.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -438,31 +432,6 @@ Foam::labelList Foam::refinementSurfaces::getClosedNamedSurfaces() const
 }
 
 
-// Orient surface (if they're closed) before any searching is done.
-void Foam::refinementSurfaces::orientSurface
-(
-    const point& keepPoint,
-    triSurfaceMesh& s
-)
-{
-    // Flip surface so keepPoint is outside.
-    bool anyFlipped = orientedSurface::orient(s, keepPoint, true);
-
-    if (anyFlipped)
-    {
-        // orientedSurface will have done a clearOut of the surface.
-        // we could do a clearout of the triSurfaceMeshes::trees()
-        // but these aren't affected by orientation (except for cached
-        // sideness which should not be set at this point. !!Should
-        // check!)
-
-        Info<< "orientSurfaces : Flipped orientation of surface "
-            << s.searchableSurface::name()
-            << " so point " << keepPoint << " is outside." << endl;
-    }
-}
-
-
 // Count number of triangles per surface region
 Foam::labelList Foam::refinementSurfaces::countRegions(const triSurface& s)
 {
@@ -485,7 +454,7 @@ void Foam::refinementSurfaces::setMinLevelFields
     const shellSurfaces& shells
 )
 {
-    minLevelFields_.setSize(surfaces_.size());
+    //minLevelFields_.setSize(surfaces_.size());
 
     forAll(surfaces_, surfI)
     {
@@ -495,26 +464,24 @@ void Foam::refinementSurfaces::setMinLevelFields
         {
             const triSurfaceMesh& triMesh = refCast<const triSurfaceMesh>(geom);
 
-            minLevelFields_.set
+            autoPtr<triSurfaceLabelField> minLevelFieldPtr
             (
-                surfI,
                 new triSurfaceLabelField
                 (
                     IOobject
                     (
-                        triMesh.searchableSurface::name(),
+                        "minLevel",
                         triMesh.objectRegistry::time().constant(),// directory
                         "triSurface",               // instance
                         triMesh,
                         IOobject::NO_READ,
-                        IOobject::AUTO_WRITE,
+                        IOobject::AUTO_WRITE
-                        false
                     ),
                     triMesh,
                     dimless
                 )
             );
-            triSurfaceLabelField& minLevelField = minLevelFields_[surfI];
+            triSurfaceLabelField& minLevelField = minLevelFieldPtr();
 
             const triSurface& s = static_cast<const triSurface&>(triMesh);
 
@@ -542,6 +509,9 @@ void Foam::refinementSurfaces::setMinLevelFields
                     shellLevel[triI]
                 );
             }
+
+            // Store field on triMesh
+            minLevelFieldPtr.ptr()->store();
         }
     }
 }
@@ -569,17 +539,38 @@ void Foam::refinementSurfaces::findHigherIntersection
         return;
     }
 
+
+    // Precalculate per surface whether it has a minlevelfield
+    UPtrList<triSurfaceLabelField> minLevelFields(surfaces_.size());
+    forAll(surfaces_, surfI)
+    {
+        const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
+
+        if (isA<triSurfaceMesh>(geom))
+        {
+            const triSurfaceMesh& triMesh = refCast<const triSurfaceMesh>(geom);
+            minLevelFields.set
+            (
+                surfI,
+               &const_cast<triSurfaceLabelField&>
+                (
+                    triMesh.lookupObject<triSurfaceLabelField>("minLevel")
+                )
+            );
+        }
+    }
+
     // Work arrays
     labelList hitMap(identity(start.size()));
     pointField p0(start);
     pointField p1(end);
     List<pointIndexHit> hitInfo(start.size());
 
-
     forAll(surfaces_, surfI)
     {
         allGeometry_[surfaces_[surfI]].findLineAny(p0, p1, hitInfo);
 
+
         // Copy all hits into arguments, continue with misses
         label newI = 0;
         forAll(hitInfo, hitI)
@@ -592,12 +583,12 @@ void Foam::refinementSurfaces::findHigherIntersection
                 // Check if minLevelField for this surface.
                 if
                 (
-                    minLevelFields_.set(surfI)
+                    minLevelFields.set(surfI)
-                 && minLevelFields_[surfI].size() > 0
+                 && minLevelFields[surfI].size() > 0
                 )
                 {
                     minLocalLevel =
-                        minLevelFields_[surfI][hitInfo[hitI].index()];
+                        minLevelFields[surfI][hitInfo[hitI].index()];
                 }
                 else
                 {
@@ -668,30 +659,57 @@ void Foam::refinementSurfaces::findAllHigherIntersections
     {
         allGeometry_[surfaces_[surfI]].findLineAll(start, end, hitInfo);
 
+        // Repack hits for surface into flat list
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // To avoid overhead of calling getRegion for every point
+
+        label n = 0;
+        forAll(hitInfo, pointI)
+        {
+            n += hitInfo[pointI].size();
+        }
+
+        List<pointIndexHit> surfInfo(n);
+        labelList pointMap(n);
+        n = 0;
+
         forAll(hitInfo, pointI)
         {
             const List<pointIndexHit>& pHits = hitInfo[pointI];
-            allGeometry_[surfaces_[surfI]].getRegion(pHits, pRegions);
-            allGeometry_[surfaces_[surfI]].getNormal(pHits, pNormals);
 
-            // Extract those hits that are on higher levelled surfaces.
+            forAll(pHits, i)
-            // Note: should move extraction of region, normal outside of loop
-            // below for if getRegion/getNormal have high overhead.
-
-            forAll(pHits, pHitI)
             {
-                label region = globalRegion(surfI, pRegions[pHitI]);
+                surfInfo[n] = pHits[i];
+                pointMap[n] = pointI;
+                n++;
+            }
+        }
 
-                if (maxLevel_[region] > currentLevel[pointI])
+        labelList surfRegion(n);
-                {
+        vectorField surfNormal(n);
-                    // Append to pointI info
+        allGeometry_[surfaces_[surfI]].getRegion(surfInfo, surfRegion);
-                    label sz = surfaceNormal[pointI].size();
+        allGeometry_[surfaces_[surfI]].getNormal(surfInfo, surfNormal);
-                    surfaceNormal[pointI].setSize(sz+1);
-                    surfaceNormal[pointI][sz] = pNormals[pHitI];
 
-                    surfaceLevel[pointI].setSize(sz+1);
+        surfInfo.clear();
-                    surfaceLevel[pointI][sz] = maxLevel_[region];
+
-                }
+
+        // Extract back into pointwise
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        forAll(surfRegion, i)
+        {
+            label region = globalRegion(surfI, surfRegion[i]);
+            label pointI = pointMap[i];
+
+            if (maxLevel_[region] > currentLevel[pointI])
+            {
+                // Append to pointI info
+                label sz = surfaceNormal[pointI].size();
+                surfaceNormal[pointI].setSize(sz+1);
+                surfaceNormal[pointI][sz] = surfNormal[i];
+
+                surfaceLevel[pointI].setSize(sz+1);
+                surfaceLevel[pointI][sz] = maxLevel_[region];
             }
         }
     }

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -90,9 +90,6 @@ class refinementSurfaces
         //- From global region number to refinement level
         labelList maxLevel_;
 
-        //- Per surface refinement level adapted for shells.
-        PtrList<triSurfaceLabelField> minLevelFields_;
-
 
     // Private Member Functions
 
@@ -213,13 +210,6 @@ public:
                 const shellSurfaces& shells
             );
 
-            //- Helper: orient (closed only) surfaces so keepPoint is outside.
-            static void orientSurface
-            (
-                const point& keepPoint,
-                triSurfaceMesh& surface
-            );
-
             //- Helper: count number of triangles per region
             static labelList countRegions(const triSurface&);
 

src/conversion/Allwmake

@@ -1,4 +1,7 @@
 #!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake libso
+
+# ----------------------------------------------------------------- end-of-file

src/conversion/meshTables/boundaryRegion.C

@@ -79,10 +79,15 @@ Foam::Map<Foam::word> Foam::boundaryRegion::names() const
 
     forAllConstIter(Map<dictionary>, *this, iter)
     {
-        word value = "boundaryRegion_" + Foam::name(iter.key());
+        lookup.insert
-        iter().readIfPresent("Label", value);
+        (
-
+            iter.key(),
-        lookup.insert(iter.key(), value);
+            iter().lookupOrDefault<word>
+            (
+                "Label",
+                "boundaryRegion_" + Foam::name(iter.key())
+            )
+        );
     }
 
     return lookup;
@@ -95,9 +100,11 @@ Foam::Map<Foam::word> Foam::boundaryRegion::boundaryTypes() const
 
     forAllConstIter(Map<dictionary>, *this, iter)
     {
-        word value = "patch";
+        lookup.insert
-        iter().readIfPresent("BoundaryType", value);
+        (
-        lookup.insert(iter.key(), value);
+            iter.key(),
+            iter().lookupOrDefault<word>("BoundaryType", "patch")
+        );
     }
 
     return lookup;
@@ -106,15 +113,16 @@ Foam::Map<Foam::word> Foam::boundaryRegion::boundaryTypes() const
 
 Foam::label Foam::boundaryRegion::findIndex(const word& name) const
 {
+    if (!name.size())
+    {
+        return -1;
+    }
+
     forAllConstIter(Map<dictionary>, *this, iter)
     {
-        word theName;
+        if (iter().lookupOrDefault<word>("Label", word::null) == name)
-        if (iter().readIfPresent("Label", theName))
         {
-            if (theName == name)
+            return iter.key();
-            {
-                return iter.key();
-            }
         }
     }
 
@@ -124,13 +132,12 @@ Foam::label Foam::boundaryRegion::findIndex(const word& name) const
 
 Foam::word Foam::boundaryRegion::boundaryType(const word& name) const
 {
-    word bndType = "patch";
+    word bndType("patch");
 
     label id = this->findIndex(name);
     if (id >= 0)
     {
-        const dictionary& dict = operator[](id);
+        operator[](id).readIfPresent<word>("BoundaryType", bndType);
-        dict.readIfPresent("BoundaryType", bndType);
     }
 
     return bndType;
@@ -192,7 +199,8 @@ void Foam::boundaryRegion::writeDict
         )
     );
 
-    ioObj.note() = "persistent data for thirdParty mesh <-> OpenFOAM translation";
+    ioObj.note() =
+        "persistent data for thirdParty mesh <-> OpenFOAM translation";
 
     Info<< "Writing " << ioObj.name() << " to " << ioObj.objectPath() << endl;
 
@@ -218,16 +226,19 @@ void Foam::boundaryRegion::operator=(const Map<dictionary>& rhs)
 
 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 
-void Foam::boundaryRegion::rename(const dictionary& dict)
+void Foam::boundaryRegion::rename(const dictionary& mapDict)
 {
-    if (!dict.size())
+    if (!mapDict.size())
     {
         return;
     }
 
-    Map<word> mapping;
+    // Use 1st pass to collect all the regions to be changed
+    // and 2nd pass to relabel regions.
+    // This avoid re-matching any renamed regions
 
-    forAllConstIter(dictionary, dict, iter)
+    Map<word> mapping;
+    forAllConstIter(dictionary, mapDict, iter)
     {
         word oldName(iter().stream());
 
@@ -238,21 +249,14 @@ void Foam::boundaryRegion::rename(const dictionary& dict)
         }
     }
 
-    if (mapping.size())
+    forAllConstIter(Map<word>, mapping, iter)
     {
-        forAllConstIter(Map<word>, mapping, iter)
+        dictionary& dict = operator[](iter.key());
-        {
-            label id = iter.key();
-            word oldName(operator[](id).lookup("Label"));
-            word newName(iter());
 
-            dictionary newDict(operator[](id));
+        Info<< "rename patch: " << iter()
-            newDict.remove("Label");
+            << " <- " << word(dict.lookup("Label")) << nl;
-            newDict.add("Label", newName);
 
-            this->set(id, newDict);
+        dict.set("Label", iter());
-            Info<< "rename patch: " << newName << " <- " << oldName << endl;
-        }
     }
 }
 

src/conversion/meshTables/cellTable.C

@@ -37,26 +37,17 @@ const char* const Foam::cellTable::defaultMaterial_ = "fluid";
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::labelList Foam::cellTable::zoneMap() const
+Foam::Map<Foam::label> Foam::cellTable::zoneMap() const
 {
-    label maxId = 0;
+    Map<label> lookup;
-    forAllConstIter(Map<dictionary>, *this, iter)
-    {
-        if (maxId < iter.key())
-        {
-            maxId = iter.key();
-        }
-    }
 
     label zoneI = 0;
-    labelList list(maxId+1, -1);
-
     forAllConstIter(Map<dictionary>, *this, iter)
     {
-        list[iter.key()] = zoneI++;
+        lookup.insert(iter.key(), zoneI++);
     }
 
-    return list;
+    return lookup;
 }
 
 
@@ -65,14 +56,12 @@ Foam::wordList Foam::cellTable::namesList() const
     Map<word> lookup = names();
     wordList lst(lookup.size());
 
-    label n = 0;
+    label zoneI = 0;
     forAllConstIter(Map<word>, lookup, iter)
     {
-        lst[n] = iter();
+        lst[zoneI++] = iter();
     }
 
-    lst.setSize(n);
-
     return lst;
 }
 
@@ -162,11 +151,15 @@ Foam::Map<Foam::word> Foam::cellTable::names() const
 
     forAllConstIter(Map<dictionary>, *this, iter)
     {
-        word theName = "cellTable_" + Foam::name(iter.key());
+        lookup.insert
-
+        (
-        iter().readIfPresent("Label", theName);
+            iter.key(),
-
+            iter().lookupOrDefault<word>
-        lookup.insert(iter.key(), theName);
+            (
+                "Label",
+                "cellTable_" + Foam::name(iter.key())
+            )
+        );
     }
 
     return lookup;
@@ -175,9 +168,9 @@ Foam::Map<Foam::word> Foam::cellTable::names() const
 
 Foam::word Foam::cellTable::name(const label& id) const
 {
-    const_iterator iter = find(id);
+    word theName("cellTable_" + Foam::name(id));
 
-    word theName = "cellTable_" + Foam::name(id);
+    const_iterator iter = find(id);
     if (iter != end())
     {
         iter().readIfPresent("Label", theName);
@@ -189,15 +182,16 @@ Foam::word Foam::cellTable::name(const label& id) const
 
 Foam::label Foam::cellTable::findIndex(const word& name) const
 {
+    if (!name.size())
+    {
+        return -1;
+    }
+
     forAllConstIter(Map<dictionary>, *this, iter)
     {
-        word theName;
+        if (iter().lookupOrDefault<word>("Label", word::null) == name)
-        if (iter().readIfPresent("Label", theName))
         {
-            if (theName == name)
+            return iter.key();
-            {
-                return iter.key();
-            }
         }
     }
 
@@ -205,26 +199,44 @@ Foam::label Foam::cellTable::findIndex(const word& name) const
 }
 
 
-Foam::Map<Foam::word> Foam::cellTable::selectType
+Foam::Map<Foam::word> Foam::cellTable::materialTypes() const
-(
-    const word& materialType
-) const
 {
     Map<word> lookup;
 
     forAllConstIter(Map<dictionary>, *this, iter)
     {
-        word matl(defaultMaterial_);
+        lookup.insert
+        (
+            iter.key(),
+            iter().lookupOrDefault<word>("MaterialType", defaultMaterial_)
+        );
+    }
 
-        iter().readIfPresent("MaterialType", matl);
+    return lookup;
+}
 
-        if (matl == materialType)
+
+Foam::Map<Foam::word> Foam::cellTable::selectType(const word& matl) const
+{
+    Map<word> lookup;
+
+    forAllConstIter(Map<dictionary>, *this, iter)
+    {
+        if
+        (
+            matl
+         == iter().lookupOrDefault<word>("MaterialType", defaultMaterial_)
+        )
         {
-            word theName = "cellTable_" + Foam::name(iter.key());
+            lookup.insert
-
+            (
-            iter().readIfPresent("Label", theName);
+                iter.key(),
-
+                iter().lookupOrDefault<word>
-            lookup.insert(iter.key(), theName);
+                (
+                    "Label",
+                    "cellTable_" + Foam::name(iter.key())
+                )
+            );
         }
     }
 
@@ -250,22 +262,6 @@ Foam::Map<Foam::word> Foam::cellTable::shells() const
 }
 
 
-Foam::Map<Foam::word> Foam::cellTable::materialTypes() const
-{
-    Map<word> lookup;
-
-    forAllConstIter(Map<dictionary>, *this, iter)
-    {
-        word matlType(defaultMaterial_);
-
-        iter().readIfPresent("MaterialType", matlType);
-
-        lookup.insert(iter.key(), matlType);
-    }
-
-    return lookup;
-}
-
 
 void Foam::cellTable::setMaterial(const label& id, const word& matlType)
 {
@@ -346,7 +342,8 @@ void Foam::cellTable::writeDict
         )
     );
 
-    ioObj.note() = "persistent data for thirdParty mesh <-> OpenFOAM translation";
+    ioObj.note() =
+        "persistent data for thirdParty mesh <-> OpenFOAM translation";
 
     Info<< "Writing " << ioObj.name() << " to " << ioObj.objectPath() << endl;
 
@@ -406,10 +403,7 @@ void Foam::cellTable::operator=(const polyMesh& mesh)
         zoneDict.insert
         (
             unZonedType,
-            dictionary
+            dictionary(IStringStream("Label cells;")())
-            (
-                IStringStream("Label cells;")()
-            )
         );
     }
 
@@ -426,22 +420,21 @@ void Foam::cellTable::addCellZones
     const labelList& tableIds
 ) const
 {
-    labelList typeToZone = zoneMap();
+    Map<label> typeToZone = zoneMap();
-    wordList  zoneNames = namesList();
-
     List<DynamicList<label> > zoneCells(size());
 
     forAll(tableIds, cellI)
     {
-        label zoneI = typeToZone[tableIds[cellI]];
+        Map<label>::const_iterator iter = typeToZone.find(tableIds[cellI]);
-        if (zoneI >= 0)
+        if (iter != typeToZone.end())
         {
-            zoneCells[zoneI].append(cellI);
+            zoneCells[iter()].append(cellI);
         }
     }
 
-    // avoid empty zones
+    // track which zones were actually used
     labelList zoneUsed(zoneCells.size());
+    wordList  zoneNames(namesList());
 
     label nZone = 0;
     forAll(zoneCells, zoneI)
@@ -454,23 +447,25 @@ void Foam::cellTable::addCellZones
     }
     zoneUsed.setSize(nZone);
 
-    mesh.cellZones().clear();
+    cellZoneMesh& czMesh = mesh.cellZones();
+
+    czMesh.clear();
     if (nZone <= 1)
     {
         Info<< "cellZones not used" << endl;
         return;
     }
-    mesh.cellZones().setSize(nZone);
+    czMesh.setSize(nZone);
 
     forAll(zoneUsed, zoneI)
     {
         const label origZoneI = zoneUsed[zoneI];
 
         Info<< "cellZone " << zoneI
-            << " (size: " << zoneCells[origZoneI].size() << ") name: "
+            << " (size: "  << zoneCells[origZoneI].size()
-            << zoneNames[origZoneI] << endl;
+            << ") name: "  << zoneNames[origZoneI] << endl;
 
-        mesh.cellZones().set
+        czMesh.set
         (
             zoneI,
             new cellZone
@@ -478,17 +473,17 @@ void Foam::cellTable::addCellZones
                 zoneNames[origZoneI],
                 zoneCells[origZoneI],
                 zoneI,
-                mesh.cellZones()
+                czMesh
             )
         );
     }
-    mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
+    czMesh.writeOpt() = IOobject::AUTO_WRITE;
 }
 
 
-void Foam::cellTable::combine(const dictionary& dict, labelList& tableIds)
+void Foam::cellTable::combine(const dictionary& mapDict, labelList& tableIds)
 {
-    if (!dict.size())
+    if (!mapDict.size())
     {
         return;
     }
@@ -496,13 +491,13 @@ void Foam::cellTable::combine(const dictionary& dict, labelList& tableIds)
     bool remap = false;
     labelList mapping(identity(max(this->toc()) + 1));
 
-    forAllConstIter (dictionary, dict, iter)
+    forAllConstIter (dictionary, mapDict, iter)
     {
         wordList  zoneNames(iter().stream());
         labelList zoneIndex(zoneNames.size());
 
         label nElem = 0;
-        forAll (zoneNames, zoneI)
+        forAll(zoneNames, zoneI)
         {
             zoneIndex[nElem] = this->findIndex(zoneNames[zoneI]);
             if (zoneIndex[nElem] >= 0)
@@ -532,7 +527,7 @@ void Foam::cellTable::combine(const dictionary& dict, labelList& tableIds)
             {
                 Info<< " = (";
             }
-            forAll (zoneNames, zoneI)
+            forAll(zoneNames, zoneI)
             {
                 Info<< " " << zoneNames[zoneI];
             }
@@ -542,14 +537,10 @@ void Foam::cellTable::combine(const dictionary& dict, labelList& tableIds)
             if (targetId < 0)
             {
                 targetId = min(zoneIndex);
-                dictionary newDict(operator[](targetId));
+                operator[](targetId).set("Label", iter().keyword());
-
-                newDict.remove("Label");
-                newDict.add("Label", iter().keyword());
-                this->set(targetId, newDict);
             }
 
-            forAll (zoneIndex, zoneI)
+            forAll(zoneIndex, zoneI)
             {
                 label idx = zoneIndex[zoneI];
                 if (idx != targetId && idx >= 0)
@@ -567,5 +558,4 @@ void Foam::cellTable::combine(const dictionary& dict, labelList& tableIds)
     }
 }
 
-
 // ************************************************************************* //

src/conversion/meshTables/cellTable.H

@@ -86,8 +86,8 @@ class cellTable
 
     // Private Member Functions
 
-        //- Map from cellTable ID => zone number (unmapped = -1)
+        //- Map from cellTable ID => zone number
-        labelList zoneMap() const;
+        Map<label> zoneMap() const;
 
         //- A contiguous list of cellTable names
         List<word> namesList() const;

src/decompositionAgglomeration/Allwmake

@@ -1,4 +1,5 @@
 #!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake libso decompositionMethods

src/dynamicMesh/meshCut/cellLooper/hexCellLooper.C

@@ -29,7 +29,6 @@ License
 #include "polyMesh.H"
 #include "cellModeller.H"
 #include "mathematicalConstants.H"
-#include "DynamicList.H"
 #include "plane.H"
 #include "ListOps.H"
 #include "meshTools.H"

src/dynamicMesh/meshCut/meshModifiers/meshCutAndRemove/meshCutAndRemove.H

@@ -40,7 +40,6 @@ SourceFiles
 #include "boolList.H"
 #include "labelList.H"
 #include "typeInfo.H"
-#include "DynamicList.H"
 #include "Map.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.H

@@ -113,10 +113,8 @@ SourceFiles
 #define meshCutter_H
 
 #include "edgeVertex.H"
-#include "boolList.H"
 #include "labelList.H"
 #include "typeInfo.H"
-#include "DynamicList.H"
 #include "Map.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/dynamicMesh/meshCut/meshModifiers/refinementIterator/refinementIterator.C

@@ -26,12 +26,10 @@ License
 
 #include "refinementIterator.H"
 #include "polyMesh.H"
-#include "polyTopoChanger.H"
 #include "Time.H"
 #include "refineCell.H"
 #include "undoableMeshCutter.H"
 #include "polyTopoChange.H"
-#include "DynamicList.H"
 #include "mapPolyMesh.H"
 #include "cellCuts.H"
 #include "OFstream.H"

src/dynamicMesh/motionSmoother/motionSmoother.C

@@ -343,7 +343,7 @@ void Foam::motionSmoother::getAffectedFacesAndPoints
 
         forAllConstIter(pointSet, nbrPoints, iter)
         {
-            const labelList& pCells = mesh_.pointCells()[iter.key()];  
+            const labelList& pCells = mesh_.pointCells(iter.key());
 
             forAll(pCells, pCellI)
             {

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C

@@ -330,15 +330,14 @@ Foam::label Foam::addPatchCellLayer::addSideFace
     const label meshEdgeI,              // corresponding mesh edge
     const label layerI,                 // layer
     const label numEdgeFaces,           // number of layers for edge
+    const labelList& meshFaces,         // precalculated edgeFaces
     polyTopoChange& meshMod
 ) const
 {
     // Edge to 'inflate' from
     label inflateEdgeI = -1;
 
-    // Mesh faces using edge
+    // Check mesh faces using edge
-    const labelList& meshFaces = mesh_.edgeFaces()[meshEdgeI];
-
     forAll(meshFaces, i)
     {
         if (mesh_.isInternalFace(meshFaces[i]))
@@ -620,6 +619,9 @@ void Foam::addPatchCellLayer::setRefinement
 
     const labelList& meshPoints = pp.meshPoints();
 
+    // Some storage for edge-face-addressing.
+    DynamicList<label> ef;
+
     // Precalculate mesh edges for pp.edges.
     labelList meshEdges(calcMeshEdges(mesh_, pp));
 
@@ -777,7 +779,9 @@ void Foam::addPatchCellLayer::setRefinement
                 label meshEdgeI = meshEdges[edgeI];
 
                 // Mesh faces using edge
-                const labelList& meshFaces = mesh_.edgeFaces()[meshEdgeI];
+
+                // Mesh faces using edge
+                const labelList& meshFaces = mesh_.edgeFaces(meshEdgeI, ef);
 
                 // Check that there is only one patchface using edge.
                 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
@@ -1353,6 +1357,12 @@ void Foam::addPatchCellLayer::setRefinement
                             patchFaceI
                         );
 
+                        const labelList& meshFaces = mesh_.edgeFaces
+                        (
+                            meshEdgeI,
+                            ef
+                        );
+
                         addSideFace
                         (
                             pp,
@@ -1365,6 +1375,7 @@ void Foam::addPatchCellLayer::setRefinement
                             meshEdgeI,      // corresponding mesh edge
                             i,
                             numEdgeSideFaces,
+                            meshFaces,
                             meshMod
                         );
                     }

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H

@@ -232,6 +232,7 @@ class addPatchCellLayer
             const label meshEdgeI,
             const label layerI,
             const label numEdgeFaces,
+            const labelList& meshFaces,
             polyTopoChange&
         ) const;
 

src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C

@@ -125,11 +125,11 @@ void Foam::combineFaces::regioniseFaces
 (
     const scalar minCos,
     const label cellI,
+    const labelList& cEdges,
     Map<label>& faceRegion
 ) const
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-    const labelList& cEdges = mesh_.cellEdges()[cellI];
 
     forAll(cEdges, i)
     {
@@ -220,9 +220,10 @@ bool Foam::combineFaces::faceNeighboursValid
         return true;
     }
 
-    const labelListList& faceEdges = mesh_.faceEdges();
     const cell& cFaces = mesh_.cells()[cellI];
 
+    DynamicList<label> storage;
+
     // Test for face collapsing to edge since too many neighbours merged.
     forAll(cFaces, cFaceI)
     {
@@ -230,7 +231,7 @@ bool Foam::combineFaces::faceNeighboursValid
 
         if (!faceRegion.found(faceI))
         {
-            const labelList& fEdges = faceEdges[faceI];
+            const labelList& fEdges = mesh_.faceEdges(faceI, storage);
 
             // Count number of remaining faces neighbouring faceI. This has
             // to be 3 or more.
@@ -299,6 +300,8 @@ Foam::labelListList Foam::combineFaces::getMergeSets
 {
     // Lists of faces that can be merged.
     DynamicList<labelList> allFaceSets(boundaryCells.size() / 10);
+    // Storage for on-the-fly cell-edge addressing.
+    DynamicList<label> storage;
 
     // On all cells regionise the faces
     forAllConstIter(labelHashSet, boundaryCells, iter)
@@ -307,9 +310,11 @@ Foam::labelListList Foam::combineFaces::getMergeSets
 
         const cell& cFaces = mesh_.cells()[cellI];
 
+        const labelList& cEdges = mesh_.cellEdges(cellI, storage);
+
         // Region per face
         Map<label> faceRegion(cFaces.size());
-        regioniseFaces(featureCos, cellI, faceRegion);
+        regioniseFaces(featureCos, cellI, cEdges, faceRegion);
 
         // Now we have in faceRegion for every face the region with planar
         // face sharing the same region. We now check whether the resulting

src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.H

@@ -103,6 +103,7 @@ class combineFaces
         (
             const scalar minCos,
             const label cellI,
+            const labelList& cEdges,
             Map<label>& faceRegion
         ) const;
 

src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C

@@ -595,6 +595,7 @@ void Foam::edgeCollapser::updateMesh(const mapPolyMesh& map)
 {
     pointRegion_.setSize(mesh_.nPoints());
     pointRegion_ = -1;
+    // Reset count, do not remove underlying storage
     pointRegionMaster_.clear();
     freeRegions_.clear();
 }

src/dynamicMesh/polyTopoChange/polyTopoChange/faceCollapser.H

@@ -54,7 +54,6 @@ SourceFiles
 #define faceCollapser_H
 
 #include "labelList.H"
-#include "DynamicList.H"
 #include "point.H"
 #include "Map.H"
 #include "labelHashSet.H"

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C

@@ -372,7 +372,7 @@ Foam::scalar Foam::hexRef8::getLevel0EdgeLength() const
         {
             const label cLevel = cellLevel_[cellI];
 
-            const labelList& cEdges = mesh_.cellEdges()[cellI];
+            const labelList& cEdges = mesh_.cellEdges(cellI);
 
             forAll(cEdges, i)
             {
@@ -447,7 +447,7 @@ Foam::scalar Foam::hexRef8::getLevel0EdgeLength() const
     {
         const label cLevel = cellLevel_[cellI];
 
-        const labelList& cEdges = mesh_.cellEdges()[cellI];
+        const labelList& cEdges = mesh_.cellEdges(cellI);
 
         forAll(cEdges, i)
         {
@@ -1190,6 +1190,10 @@ void Foam::hexRef8::createInternalFaces
     // From edge mid to face mids
     Map<edge> midPointToFaceMids(24);
 
+    // Storage for on-the-fly addressing
+    DynamicList<label> storage;
+
+
     // Running count of number of internal faces added so far.
     label nFacesAdded = 0;
 
@@ -1198,7 +1202,7 @@ void Foam::hexRef8::createInternalFaces
         label faceI = cFaces[i];
 
         const face& f = mesh_.faces()[faceI];
-        const labelList& fEdges = mesh_.faceEdges()[faceI];
+        const labelList& fEdges = mesh_.faceEdges(faceI, storage);
 
         // We are on the cellI side of face f. The face will have 1 or 4
         // cLevel points and lots of higher numbered ones.
@@ -1299,7 +1303,7 @@ void Foam::hexRef8::createInternalFaces
                     {
                         dumpCell(cellI);
 
-                        const labelList cPoints(cellPoints(cellI));
+                        const labelList& cPoints = mesh_.cellPoints(cellI);
 
                         FatalErrorIn("createInternalFaces(..)")
                             << "cell:" << cellI << " cLevel:" << cLevel
@@ -1372,7 +1376,7 @@ void Foam::hexRef8::createInternalFaces
                     {
                         dumpCell(cellI);
 
-                        const labelList cPoints(cellPoints(cellI));
+                        const labelList& cPoints = mesh_.cellPoints(cellI);
 
                         FatalErrorIn("createInternalFaces(..)")
                             << "cell:" << cellI << " cLevel:" << cLevel
@@ -1454,7 +1458,7 @@ void Foam::hexRef8::walkFaceToMid
 ) const
 {
     const face& f = mesh_.faces()[faceI];
-    const labelList& fEdges = mesh_.faceEdges()[faceI];
+    const labelList& fEdges = mesh_.faceEdges(faceI);
 
     label fp = startFp;
 
@@ -1503,7 +1507,7 @@ void Foam::hexRef8::walkFaceFromMid
 ) const
 {
     const face& f = mesh_.faces()[faceI];
-    const labelList& fEdges = mesh_.faceEdges()[faceI];
+    const labelList& fEdges = mesh_.faceEdges(faceI);
 
     label fp = f.rcIndex(startFp);
 
@@ -2013,27 +2017,6 @@ Foam::hexRef8::hexRef8
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-//- Get points of a cell (without using cellPoints addressing)
-Foam::labelList Foam::hexRef8::cellPoints(const label cellI) const
-{
-    // Pick up points of the cell
-    const cell& cFaces = mesh_.cells()[cellI];
-
-    labelHashSet cPoints(4*cFaces.size());
-
-    forAll(cFaces, i)
-    {
-        const face& f = mesh_.faces()[cFaces[i]];
-
-        forAll(f, fp)
-        {
-            cPoints.insert(f[fp]);
-        }
-    }
-    return cPoints.toc();
-}
-
-
 Foam::labelList Foam::hexRef8::consistentRefinement
 (
     const labelList& cellsToRefine,
@@ -2358,13 +2341,11 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
         // as cell level purely for ease)
         labelList maxPointCount(mesh_.nPoints(), 0);
 
-        const labelListList& pointCells = mesh_.pointCells();
+        forAll(maxPointCount, pointI)
-
-        forAll(pointCells, pointI)
         {
             label& pLevel = maxPointCount[pointI];
 
-            const labelList& pCells = pointCells[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             forAll(pCells, i)
             {
@@ -2395,7 +2376,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
             // Loop over all cells using the point and check whether their
             // refinement level is much less than the maximum.
 
-            const labelList& pCells = pointCells[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             forAll(pCells, pCellI)
             {
@@ -3021,7 +3002,9 @@ Foam::labelListList Foam::hexRef8::setRefinement
             << " Checking initial mesh just to make sure" << endl;
 
         checkMesh();
-        checkRefinementLevels(-1, labelList(0));
+        // Cannot call checkRefinementlevels since hanging points might
+        // get triggered by the mesher after subsetting.
+        //checkRefinementLevels(-1, labelList(0));
     }
 
     // Clear any saved point/cell data.
@@ -3119,7 +3102,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
     {
         if (cellMidPoint[cellI] >= 0)
         {
-            const labelList& cEdges = mesh_.cellEdges()[cellI];
+            const labelList& cEdges = mesh_.cellEdges(cellI);
 
             forAll(cEdges, i)
             {
@@ -3456,7 +3439,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
 
         forAll(pointLevel_, pointI)
         {
-            const labelList& pCells = mesh_.pointCells()[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             forAll(pCells, pCellI)
             {
@@ -3496,7 +3479,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
                 {
                     dumpCell(cellI);
 
-                    const labelList cPoints(cellPoints(cellI));
+                    const labelList& cPoints = mesh_.cellPoints(cellI);
 
                     FatalErrorIn
                     (
@@ -3608,7 +3591,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
         {
             if (edgeMidPoint[edgeI] >= 0)
             {
-                const labelList& eFaces = mesh_.edgeFaces()[edgeI];
+                const labelList& eFaces = mesh_.edgeFaces(edgeI);
 
                 forAll(eFaces, i)
                 {
@@ -3766,13 +3749,16 @@ Foam::labelListList Foam::hexRef8::setRefinement
             << endl;
     }
 
+    DynamicList<label> eFacesStorage;
+    DynamicList<label> fEdgesStorage;
+
     forAll(edgeMidPoint, edgeI)
     {
         if (edgeMidPoint[edgeI] >= 0)
         {
             // Split edge. Check that face not already handled above.
 
-            const labelList& eFaces = mesh_.edgeFaces()[edgeI];
+            const labelList& eFaces = mesh_.edgeFaces(edgeI, eFacesStorage);
 
             forAll(eFaces, i)
             {
@@ -3783,7 +3769,11 @@ Foam::labelListList Foam::hexRef8::setRefinement
                     // Unsplit face. Add edge splits to face.
 
                     const face& f = mesh_.faces()[faceI];
-                    const labelList& fEdges = mesh_.faceEdges()[faceI];
+                    const labelList& fEdges = mesh_.faceEdges
+                    (
+                        faceI,
+                        fEdgesStorage
+                    );
 
                     DynamicList<label> newFaceVerts(f.size());
 
@@ -4713,14 +4703,12 @@ void Foam::hexRef8::checkRefinementLevels
     // Check 2:1 across points (instead of faces)
     if (maxPointDiff != -1)
     {
-        const labelListList& pointCells = mesh_.pointCells();
-
         // Determine per point the max cell level.
         labelList maxPointLevel(mesh_.nPoints(), 0);
 
-        forAll(pointCells, pointI)
+        forAll(maxPointLevel, pointI)
         {
-            const labelList& pCells = pointCells[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             label& pLevel = maxPointLevel[pointI];
 
@@ -4745,7 +4733,7 @@ void Foam::hexRef8::checkRefinementLevels
         {
             label pointI = pointsToCheck[i];
 
-            const labelList& pCells = pointCells[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             forAll(pCells, i)
             {
@@ -4813,7 +4801,7 @@ void Foam::hexRef8::checkRefinementLevels
             false                   // no separation
         );
 
-        OFstream str(mesh_.time().path()/"hangingPoints.obj");
+        //OFstream str(mesh_.time().path()/"hangingPoints.obj");
 
         label nHanging = 0;
 
@@ -4826,8 +4814,7 @@ void Foam::hexRef8::checkRefinementLevels
                 Pout<< "Hanging boundary point " << pointI
                     << " at " << mesh_.points()[pointI]
                     << endl;
-
+                //meshTools::writeOBJ(str, mesh_.points()[pointI]);
-                meshTools::writeOBJ(str, mesh_.points()[pointI]);
             }
         }
 
@@ -4880,11 +4867,11 @@ Foam::labelList Foam::hexRef8::getSplitPoints() const
     labelList splitMasterLevel(mesh_.nPoints(), 0);
 
     // Unmark all with not 8 cells
-    const labelListList& pointCells = mesh_.pointCells();
+    //const labelListList& pointCells = mesh_.pointCells();
 
-    forAll(pointCells, pointI)
+    for (label pointI = 0; pointI < mesh_.nPoints(); pointI++)
     {
-        const labelList& pCells = pointCells[pointI];
+        const labelList& pCells = mesh_.pointCells(pointI);
 
         if (pCells.size() != 8)
         {
@@ -4897,8 +4884,7 @@ Foam::labelList Foam::hexRef8::getSplitPoints() const
 
     forAll(visibleCells, cellI)
     {
-        //const labelList& cPoints = mesh_.cellPoints()[cellI];
+        const labelList& cPoints = mesh_.cellPoints(cellI);
-        const labelList cPoints(cellPoints(cellI));
 
         if (visibleCells[cellI] != -1 && history_.parentIndex(cellI) >= 0)
         {
@@ -5103,7 +5089,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
         {
             if (unrefinePoint.get(pointI) == 1)
             {
-                const labelList& pCells = mesh_.pointCells()[pointI];
+                const labelList& pCells = mesh_.pointCells(pointI);
 
                 forAll(pCells, j)
                 {
@@ -5243,7 +5229,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
         {
             if (unrefinePoint.get(pointI) == 1)
             {
-                const labelList& pCells = mesh_.pointCells()[pointI];
+                const labelList& pCells = mesh_.pointCells(pointI);
 
                 forAll(pCells, j)
                 {
@@ -5328,7 +5314,7 @@ void Foam::hexRef8::setUnrefinement
 
         forAll(splitPointLabels, i)
         {
-            const labelList& pCells = mesh_.pointCells()[splitPointLabels[i]];
+            const labelList& pCells = mesh_.pointCells(splitPointLabels[i]);
 
             forAll(pCells, j)
             {
@@ -5394,7 +5380,7 @@ void Foam::hexRef8::setUnrefinement
 
         // Get original cell label
 
-        const labelList& pCells = mesh_.pointCells()[pointI];
+        const labelList& pCells = mesh_.pointCells(pointI);
 
         // Check
         if (pCells.size() != 8)
@@ -5462,7 +5448,7 @@ void Foam::hexRef8::setUnrefinement
     {
         label pointI = splitPointLabels[i];
 
-        const labelList& pCells = mesh_.pointCells()[pointI];
+        const labelList& pCells = mesh_.pointCells(pointI);
 
         label masterCellI = min(pCells);
 

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H

@@ -370,9 +370,6 @@ public:
 
         // Refinement
 
-            //- Helper:get points of a cell without using cellPoints addressing
-            labelList cellPoints(const label cellI) const;
-
             //- Given valid mesh and current cell level and proposed
             //  cells to refine calculate any clashes (due to 2:1) and return
             //  ok list of cells to refine.

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChangeTemplates.C

@@ -96,13 +96,5 @@ void Foam::polyTopoChange::renumberKey
     elems.transfer(newElems);
 }
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
 
 // ************************************************************************* //

src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C

@@ -83,6 +83,7 @@ Foam::label Foam::removeFaces::changeFaceRegion
     const labelList& nFacesPerEdge,
     const label faceI,
     const label newRegion,
+    const labelList& fEdges,
     labelList& faceRegion
 ) const
 {
@@ -94,27 +95,33 @@ Foam::label Foam::removeFaces::changeFaceRegion
 
         nChanged = 1;
 
+        // Storage for on-the-fly addressing
+        DynamicList<label> fe;
+        DynamicList<label> ef;
+
         // Step to neighbouring faces across edges that will get removed
-
-        const labelList& fEdges = mesh_.faceEdges()[faceI];
-
         forAll(fEdges, i)
         {
             label edgeI = fEdges[i];
 
             if (nFacesPerEdge[edgeI] >= 0 && nFacesPerEdge[edgeI] <= 2)
             {
-                const labelList& eFaces = mesh_.edgeFaces()[edgeI];
+                const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
 
                 forAll(eFaces, j)
                 {
+                    label nbrFaceI = eFaces[j];
+
+                    const labelList& fEdges1 = mesh_.faceEdges(nbrFaceI, fe);
+
                     nChanged += changeFaceRegion
                     (
                         cellRegion,
                         removedFace,
                         nFacesPerEdge,
-                        eFaces[j],
+                        nbrFaceI,
                         newRegion,
+                        fEdges1,
                         faceRegion
                     );
                 }
@@ -166,7 +173,7 @@ Foam::boolList Foam::removeFaces::getFacesAffected
     //  Mark faces affected by removal of edges
     forAllConstIter(labelHashSet, edgesToRemove, iter)
     {
-        const labelList& eFaces = mesh_.edgeFaces()[iter.key()];
+        const labelList& eFaces = mesh_.edgeFaces(iter.key());
 
         forAll(eFaces, eFaceI)
         {
@@ -814,6 +821,10 @@ void Foam::removeFaces::setRefinement
     // Number of connected face regions
     label nRegions = 0;
 
+    // Storage for on-the-fly addressing
+    DynamicList<label> fe;
+    DynamicList<label> ef;
+
 
     {
         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
@@ -827,7 +838,7 @@ void Foam::removeFaces::setRefinement
         {
             label faceI = faceLabels[i];
 
-            const labelList& fEdges = mesh_.faceEdges()[faceI];
+            const labelList& fEdges = mesh_.faceEdges(faceI, fe);
 
             forAll(fEdges, i)
             {
@@ -835,8 +846,7 @@ void Foam::removeFaces::setRefinement
 
                 if (nFacesPerEdge[edgeI] == -1)
                 {
-                    nFacesPerEdge[edgeI] =
+                    nFacesPerEdge[edgeI] = mesh_.edgeFaces(edgeI, ef).size()-1;
-                        mesh_.edgeFaces()[edgeI].size()-1;
                 }
                 else
                 {
@@ -849,16 +859,15 @@ void Foam::removeFaces::setRefinement
         // Note that this only needs to be done for possibly coupled edges
         // so we could choose to loop only over boundary faces and use faceEdges
         // of those.
-        const labelListList& edgeFaces = mesh_.edgeFaces();
 
-        forAll(edgeFaces, edgeI)
+        forAll(mesh_.edges(), edgeI)
         {
             if (nFacesPerEdge[edgeI] == -1)
             {
                 // Edge not yet handled in loop above so is not used by any
                 // face to be removed.
 
-                const labelList& eFaces = edgeFaces[edgeI];
+                const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
 
                 if (eFaces.size() > 2)
                 {
@@ -922,7 +931,7 @@ void Foam::removeFaces::setRefinement
                 label f0 = -1;
                 label f1 = -1;
 
-                const labelList& eFaces = mesh_.edgeFaces()[edgeI];
+                const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
 
                 forAll(eFaces, i)
                 {
@@ -1152,6 +1161,7 @@ void Foam::removeFaces::setRefinement
                 nFacesPerEdge,
                 startFaceI,
                 nRegions,
+                mesh_.faceEdges(startFaceI, fe),
                 faceRegion
             );
 

src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.H

@@ -93,6 +93,7 @@ class removeFaces
             const labelList& nFacesPerEdge,
             const label faceI,
             const label newRegion,
+            const labelList& fEdges,
             labelList& faceRegion
         ) const;
 

src/edgeMesh/edgeMesh.C

@@ -156,7 +156,7 @@ label edgeMesh::regions(labelList& edgeRegion) const
                 }
             }
 
-            edgesToVisit.transfer(newEdgesToVisit.shrink());
+            edgesToVisit.transfer(newEdgesToVisit);
         }
 
         currentRegion++;

src/finiteVolume/Make/files

@@ -69,9 +69,11 @@ $(constraintFvPatchFields)/wedge/wedgeFvPatchScalarField.C
 derivedFvPatchFields = $(fvPatchFields)/derived
 $(derivedFvPatchFields)/advective/advectiveFvPatchFields.C
 $(derivedFvPatchFields)/directMappedFixedValue/directMappedFixedValueFvPatchFields.C
+$(derivedFvPatchFields)/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.C
 $(derivedFvPatchFields)/fan/fanFvPatchFields.C
 $(derivedFvPatchFields)/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
+$(derivedFvPatchFields)/fixedInternalValueFvPatchField/fixedInternalValueFvPatchFields.C
 $(derivedFvPatchFields)/fixedNormalSlip/fixedNormalSlipFvPatchFields.C
 $(derivedFvPatchFields)/fluxCorrectedVelocity/fluxCorrectedVelocityFvPatchVectorField.C
 $(derivedFvPatchFields)/freestream/freestreamFvPatchFields.C
@@ -263,6 +265,7 @@ cfdTools/general/porousMedia/porousZone.C
 cfdTools/general/porousMedia/porousZones.C
 cfdTools/general/MRF/MRFZone.C
 cfdTools/general/MRF/MRFZones.C
+cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.C
 cfdTools/general/fieldSources/timeActivatedExplicitSource/timeActivatedExplicitSource.C
 cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C
 cfdTools/general/SRF/SRFModel/SRFModel/newSRFModel.C

src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.C

@@ -0,0 +1,170 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "pressureGradientExplicitSource.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::pressureGradientExplicitSource::pressureGradientExplicitSource
+(
+    const word& sourceName,
+    const fvMesh& mesh,
+    volVectorField& U
+)
+:
+    dict_
+    (
+        IOobject
+        (
+            sourceName + "Properties",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    mesh_(mesh),
+    U_(U),
+    Ubar_(dict_.lookup("Ubar")),
+    gradPini_(readScalar(dict_.lookup("gradPini"))),
+    gradP_(gradPini_),
+    flowDir_(Ubar_/mag(Ubar_)),
+    cellSource_(dict_.lookup("cellSource")),
+    cellSelector_
+    (
+        topoSetSource::New
+        (
+            cellSource_,
+            mesh,
+            dict_.subDict(cellSource_ + "Coeffs")
+        )
+    ),
+    selectedCellSet_
+    (
+        mesh,
+        "pressureGradientExplicitSourceCellSet",
+        mesh.nCells()/10 + 1  // Reasonable size estimate.
+    )
+{
+    // Create the cell set
+    cellSelector_->applyToSet
+    (
+        topoSetSource::NEW,
+        selectedCellSet_
+    );
+
+    // Give some feedback
+    Info<< "pressureGradientExplicitSource(" << sourceName << ")" << nl
+        << "Selected " << returnReduce(selectedCellSet_.size(), sumOp<label>())
+        << " cells." << endl;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh> >
+Foam::pressureGradientExplicitSource::Su() const
+{
+    tmp<DimensionedField<vector, volMesh> > tSource
+    (
+        new DimensionedField<vector, volMesh>
+        (
+            IOobject
+            (
+                "pressureGradientExplicitSource",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedVector("zero", dimVelocity/dimTime, vector::zero)
+        )
+    );
+
+    DimensionedField<vector, volMesh>& sourceField = tSource();
+
+    forAllConstIter(cellSet, selectedCellSet_, iter)
+    {
+        label cellI = iter.key();
+
+        sourceField[cellI] = flowDir_*gradP_;
+    }
+
+    return tSource;
+}
+
+
+void Foam::pressureGradientExplicitSource::update()
+{
+    const volScalarField& rUA =
+        mesh_.lookupObject<volScalarField>("(1|A(" + U_.name() + "))");
+
+    // Integrate flow variables over cell set
+    scalar volTot = 0.0;
+    scalar magUbarAve = 0.0;
+    scalar rUAave = 0.0;
+    forAllConstIter(cellSet, selectedCellSet_, iter)
+    {
+        label cellI = iter.key();
+
+        scalar volCell = mesh_.V()[cellI];
+        volTot += volCell;
+
+        magUbarAve += (flowDir_ & U_[cellI])*volCell;
+        rUAave += rUA[cellI]*volCell;
+    }
+
+    // Collect across all processors
+    reduce(volTot, sumOp<scalar>());
+    reduce(magUbarAve, sumOp<scalar>());
+    reduce(rUAave, sumOp<scalar>());
+
+    // Volume averages
+    magUbarAve /= volTot;
+    rUAave /= volTot;
+
+    // Calculate the pressure gradient increment needed to adjust the average
+    // flow-rate to the desired value
+    scalar gradPplus = (mag(Ubar_) - magUbarAve)/rUAave;
+
+    // Apply correction to velocity field
+    forAllConstIter(cellSet, selectedCellSet_, iter)
+    {
+        label cellI = iter.key();
+        U_[cellI] += flowDir_*rUA[cellI]*gradPplus;
+    }
+
+    // Update pressure gradient
+    gradP_ += gradPplus;
+
+    Info<< "Uncorrected Ubar = " << magUbarAve << tab
+        << "Pressure gradient = " << gradP_ << endl;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.H

@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::pressureGradientExplicitSource
+
+Description
+    Creates a cell set pressure gradient source
+
+    Note: Currently only handles kinematic pressure
+
+SourceFiles
+    pressureGradientExplicitSource.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pressureGradientExplicitSource_H
+#define pressureGradientExplicitSource_H
+
+#include "autoPtr.H"
+#include "topoSetSource.H"
+#include "cellSet.H"
+#include "fvMesh.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+              Class pressureGradientExplicitSource Declaration
+\*---------------------------------------------------------------------------*/
+
+class pressureGradientExplicitSource
+{
+    // Private data
+
+        //- Properties dictionary
+        IOdictionary dict_;
+
+        //- Reference to the mesh
+        const fvMesh& mesh_;
+
+        //- Reference to the velocity field
+        volVectorField& U_;
+
+        //- Average velocity
+        vector Ubar_;
+
+        //- Initial pressure gradient
+        scalar gradPini_;
+
+        //- Pressure gradient
+        scalar gradP_;
+
+        //- Flow direction
+        vector flowDir_;
+
+        //- Name of cell source
+        word cellSource_;
+
+        //- The method by which the cells will be selecetd
+        autoPtr<topoSetSource> cellSelector_;
+
+        //- The set of selected cells
+        cellSet selectedCellSet_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        pressureGradientExplicitSource(const pressureGradientExplicitSource&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const pressureGradientExplicitSource&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from explicit source name and mesh
+        pressureGradientExplicitSource
+        (
+            const word& sourceName,
+            const fvMesh& mesh,
+            volVectorField& U
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return a tmp field of the source
+            tmp<DimensionedField<vector, volMesh> > Su() const;
+
+
+        //- Correct driving force for a constant mass flow rate
+        void update();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -3,20 +3,11 @@
     int nOuterCorr(readInt(pimple.lookup("nOuterCorrectors")));
     int nCorr(readInt(pimple.lookup("nCorrectors")));
 
-    int nNonOrthCorr = 0;
+    int nNonOrthCorr =
-    if (pimple.found("nNonOrthogonalCorrectors"))
+        pimple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-    {
-        nNonOrthCorr = readInt(pimple.lookup("nNonOrthogonalCorrectors"));
-    }
 
-    bool momentumPredictor = true;
+    bool momentumPredictor =
-    if (pimple.found("momentumPredictor"))
+        pimple.lookupOrDefault<Switch>("momentumPredictor", true);
-    {
-        momentumPredictor = Switch(pimple.lookup("momentumPredictor"));
-    }
 
-    bool transonic = false;
+    bool transonic =
-    if (pimple.found("transonic"))
+        pimple.lookupOrDefault<Switch>("transonic", false);
-    {
-        transonic = Switch(pimple.lookup("transonic"));
-    }

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

@@ -2,32 +2,18 @@
 
     int nCorr(readInt(piso.lookup("nCorrectors")));
 
-    int nNonOrthCorr = 0;
+    int nNonOrthCorr =
-    if (piso.found("nNonOrthogonalCorrectors"))
+        piso.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-    {
-        nNonOrthCorr = readInt(piso.lookup("nNonOrthogonalCorrectors"));
-    }
 
-    bool momentumPredictor = true;
+    bool momentumPredictor =
-    if (piso.found("momentumPredictor"))
+        piso.lookupOrDefault<Switch>("momentumPredictor", true);
-    {
-        momentumPredictor = Switch(piso.lookup("momentumPredictor"));
-    }
 
-    bool transonic = false;
+    bool transonic =
-    if (piso.found("transonic"))
+        piso.lookupOrDefault<Switch>("transonic", false);
-    {
-        transonic = Switch(piso.lookup("transonic"));
-    }
 
-    int nOuterCorr = 1;
+    int nOuterCorr =
-    if (piso.found("nOuterCorrectors"))
+        piso.lookupOrDefault<int>("nOuterCorrectors", 1);
-    {
+
-        nOuterCorr = readInt(piso.lookup("nOuterCorrectors"));
+    bool ddtPhiCorr =
-    }
+        piso.lookupOrDefault<Switch>("ddtPhiCorr", false);
 
-    bool ddtPhiCorr = false;
-    if (piso.found("ddtPhiCorr"))
-    {
-        ddtPhiCorr = Switch(piso.lookup("ddtPhiCorr"));
-    }

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

@@ -1,25 +1,14 @@
     dictionary simple = mesh.solutionDict().subDict("SIMPLE");
 
-    int nNonOrthCorr = 0;
+    int nNonOrthCorr =
-    if (simple.found("nNonOrthogonalCorrectors"))
+        simple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-    {
-        nNonOrthCorr = readInt(simple.lookup("nNonOrthogonalCorrectors"));
-    }
 
-    bool momentumPredictor = true;
+    bool momentumPredictor =
-    if (simple.found("momentumPredictor"))
+        simple.lookupOrDefault<Switch>("momentumPredictor", true);
-    {
-        momentumPredictor = Switch(simple.lookup("momentumPredictor"));
-    }
 
-    bool fluxGradp = false;
+    bool fluxGradp =
-    if (simple.found("fluxGradp"))
+        simple.lookupOrDefault<Switch>("fluxGradp", false);
-    {
+
-        fluxGradp = Switch(simple.lookup("fluxGradp"));
+    bool transonic =
-    }
+        simple.lookupOrDefault<Switch>("transonic", false);
 
-    bool transonic = false;
-    if (simple.found("transonic"))
-    {
-        transonic = Switch(simple.lookup("transonic"));
-    }

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

@@ -40,11 +40,7 @@ scalar maxCo
     readScalar(runTime.controlDict().lookup("maxCo"))
 );
 
-scalar maxDeltaT = GREAT;
+scalar maxDeltaT =
-
+    runTime.controlDict().lookupOrDefault<scalar>("maxDeltaT", GREAT);
-if (runTime.controlDict().found("maxDeltaT"))
-(
-    maxDeltaT = readScalar(runTime.controlDict().lookup("maxDeltaT"))
-);
 
 // ************************************************************************* //