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

src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.C

@@ -37,6 +37,18 @@ void Foam::OutputFilterFunctionObject<OutputFilter>::readDict()
     dict_.readIfPresent("dictionary", dictName_);
     dict_.readIfPresent("enabled", enabled_);
     dict_.readIfPresent("storeFilter", storeFilter_);
+    dict_.readIfPresent("timeStart", timeStart_);
+    dict_.readIfPresent("timeEnd", timeEnd_);
+}
+
+
+template<class OutputFilter>
+bool Foam::OutputFilterFunctionObject<OutputFilter>::active() const
+{
+    return
+        enabled_
+     && time_.value() >= timeStart_
+     && time_.value() <= timeEnd_;
 }
 
 
@@ -94,6 +106,8 @@ Foam::OutputFilterFunctionObject<OutputFilter>::OutputFilterFunctionObject
     dictName_(),
     enabled_(true),
     storeFilter_(true),
+    timeStart_(-VGREAT),
+    timeEnd_(VGREAT),
     outputControl_(t, dict)
 {
     readDict();
@@ -121,7 +135,7 @@ bool Foam::OutputFilterFunctionObject<OutputFilter>::start()
 {
     readDict();
 
-    if (enabled_&&storeFilter_)
+    if (enabled_ && storeFilter_)
     {
         allocateFilter();
     }
@@ -136,7 +150,7 @@ bool Foam::OutputFilterFunctionObject<OutputFilter>::execute
     const bool forceWrite
 )
 {
-    if (enabled_)
+    if (active())
     {
         if (!storeFilter_)
         {
@@ -163,7 +177,7 @@ bool Foam::OutputFilterFunctionObject<OutputFilter>::execute
 template<class OutputFilter>
 bool Foam::OutputFilterFunctionObject<OutputFilter>::end()
 {
-    if (enabled_)
+    if (active())
     {
         if (!storeFilter_)
         {

src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.H

@@ -68,18 +68,28 @@ class OutputFilterFunctionObject
         //- Input dictionary
         dictionary dict_;
 
-        //- Name of region
-        word regionName_;
 
-        //- Optional dictionary name to supply required inputs
-        word dictName_;
+        // Optional user inputs
 
-        //- Switch for the execution of the functionObject
-        bool enabled_;
+            //- Name of region - defaults to name of polyMesh::defaultRegion
+            word regionName_;
+
+            //- Dictionary name to supply required inputs
+            word dictName_;
+
+            //- Switch for the execution - defaults to 'yes/on'
+            bool enabled_;
+
+            //- Switch to store filter in between writes or use on-the-fly
+            //  construction - defaults to true
+            bool storeFilter_;
+
+            //- Activation time - defaults to -VGREAT
+            scalar timeStart_;
+
+            //- De-activation time - defaults to VGREAT
+            scalar timeEnd_;
 
-        //- Switch to store filter in between writes or use on-the-fly
-        //  construction
-        bool storeFilter_;
 
         //- Output controls
         outputFilterOutputControl outputControl_;
@@ -99,6 +109,9 @@ class OutputFilterFunctionObject
         //- Destroys most of the data associated with this object.
         void destroyFilter();
 
+        //- Returns true if active (enabled and within time bounds)
+        bool active() const;
+
         //- Disallow default bitwise copy construct
         OutputFilterFunctionObject(const OutputFilterFunctionObject&);
 

src/OpenFOAM/interpolations/interpolateSplineXY/interpolateSplineXY.C

@@ -40,7 +40,7 @@ Foam::Field<Type> Foam::interpolateSplineXY
 
     forAll(xNew, i)
     {
-        yNew[i] = interpolateSmoothXY(xNew[i], xOld, yOld);
+        yNew[i] = interpolateSplineXY(xNew[i], xOld, yOld);
     }
 
     return yNew;

src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C

@@ -514,7 +514,8 @@ Foam::label Foam::polyBoundaryMesh::whichPatch(const label faceIndex) const
         FatalErrorIn
         (
             "polyBoundaryMesh::whichPatch(const label faceIndex) const"
-        )   << "given label greater than the number of geometric faces"
+        )   << "given label " << faceIndex
+            << " greater than the number of geometric faces " << mesh().nFaces()
             << abort(FatalError);
     }
 

src/conversion/meshReader/starcd/STARCDMeshReader.C

@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "STARCDMeshReader.H"
-#include "cyclicPolyPatch.H"
+#include "oldCyclicPolyPatch.H"
 #include "emptyPolyPatch.H"
 #include "wallPolyPatch.H"
 #include "symmetryPolyPatch.H"
@@ -950,7 +950,7 @@ void Foam::meshReaders::STARCD::readBoundary(const fileName& inputName)
         {
             // incorrect. should be cyclicPatch but this
             // requires info on connected faces.
-            patchTypes_[patchI] = cyclicPolyPatch::typeName;
+            patchTypes_[patchI] = oldCyclicPolyPatch::typeName;
             patchPhysicalTypes_[patchI] = patchTypes_[patchI];
         }
         else if (origType == "baffle")

src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C

@@ -186,7 +186,16 @@ void Foam::dynamicRefineFvMesh::readDict()
         ).subDict(typeName + "Coeffs")
     );
 
-    correctFluxes_ = List<Pair<word> >(refineDict.lookup("correctFluxes"));
+    List<Pair<word> > fluxVelocities = List<Pair<word> >
+    (
+        refineDict.lookup("correctFluxes")
+    );
+    // Rework into hashtable.
+    correctFluxes_.resize(fluxVelocities.size());
+    forAll(fluxVelocities, i)
+    {
+        correctFluxes_.insert(fluxVelocities[i][0], fluxVelocities[i][1]);
+    }
 
     dumpLevel_ = Switch(refineDict.lookup("dumpLevel"));
 }
@@ -289,23 +298,46 @@ Foam::dynamicRefineFvMesh::refine
                 << " split faces " << endl;
         }
 
-        forAll(correctFluxes_, i)
+        HashTable<const surfaceScalarField*> fluxes
+        (
+            lookupClass<surfaceScalarField>()
+        );
+        forAllConstIter(HashTable<const surfaceScalarField*>, fluxes, iter)
         {
+            if (!correctFluxes_.found(iter.key()))
+            {
+                WarningIn("dynamicRefineFvMesh::refine(const labelList&)")
+                    << "Cannot find surfaceScalarField " << iter.key()
+                    << " in user-provided flux mapping table "
+                    << correctFluxes_ << endl
+                    << "    The flux mapping table is used to recreate the"
+                    << " flux on newly created faces." << endl
+                    << "    Either add the entry if it is a flux or use ("
+                    << iter.key() << " none) to suppress this warning."
+                    << endl;
+                continue;
+            }
+
+            const word& UName = correctFluxes_[iter.key()];
+
+            if (UName == "none")
+            {
+                continue;
+            }
+
             if (debug)
             {
-                Info<< "Mapping flux " << correctFluxes_[i][0]
-                    << " using interpolated flux " << correctFluxes_[i][1]
+                Info<< "Mapping flux " << iter.key()
+                    << " using interpolated flux " << UName
                     << endl;
             }
-            surfaceScalarField& phi = const_cast<surfaceScalarField&>
-            (
-                lookupObject<surfaceScalarField>(correctFluxes_[i][0])
-            );
+
+            surfaceScalarField& phi = const_cast<surfaceScalarField&>(*iter());
             const surfaceScalarField phiU
             (
                 fvc::interpolate
                 (
-                    lookupObject<volVectorField>(correctFluxes_[i][1])
+                    lookupObject<volVectorField>(UName)
                 )
               & Sf()
             );
@@ -482,27 +514,51 @@ Foam::dynamicRefineFvMesh::unrefine
         const labelList& reversePointMap = map().reversePointMap();
         const labelList& reverseFaceMap = map().reverseFaceMap();
 
-        forAll(correctFluxes_, i)
+        HashTable<const surfaceScalarField*> fluxes
+        (
+            lookupClass<surfaceScalarField>()
+        );
+        forAllConstIter(HashTable<const surfaceScalarField*>, fluxes, iter)
         {
+            if (!correctFluxes_.found(iter.key()))
+            {
+                WarningIn("dynamicRefineFvMesh::refine(const labelList&)")
+                    << "Cannot find surfaceScalarField " << iter.key()
+                    << " in user-provided flux mapping table "
+                    << correctFluxes_ << endl
+                    << "    The flux mapping table is used to recreate the"
+                    << " flux on newly created faces." << endl
+                    << "    Either add the entry if it is a flux or use ("
+                    << iter.key() << " none) to suppress this warning."
+                    << endl;
+                continue;
+            }
+
+            const word& UName = correctFluxes_[iter.key()];
+
+            if (UName == "none")
+            {
+                continue;
+            }
+
             if (debug)
             {
-                Info<< "Mapping flux " << correctFluxes_[i][0]
-                    << " using interpolated flux " << correctFluxes_[i][1]
+                Info<< "Mapping flux " << iter.key()
+                    << " using interpolated flux " << UName
                     << endl;
             }
-            surfaceScalarField& phi = const_cast<surfaceScalarField&>
-            (
-                lookupObject<surfaceScalarField>(correctFluxes_[i][0])
-            );
-            surfaceScalarField phiU
+
+            surfaceScalarField& phi = const_cast<surfaceScalarField&>(*iter());
+            const surfaceScalarField phiU
             (
                 fvc::interpolate
                 (
-                    lookupObject<volVectorField>(correctFluxes_[i][1])
+                    lookupObject<volVectorField>(UName)
                 )
               & Sf()
             );
 
+
             forAllConstIter(Map<label>, faceToSplitPoint, iter)
             {
                 label oldFaceI = iter.key();

src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H

@@ -64,7 +64,7 @@ protected:
         Switch dumpLevel_;
 
         //- Fluxes to map
-        List<Pair<word> > correctFluxes_;
+        HashTable<word> correctFluxes_;
 
         //- Number of refinement/unrefinement steps done so far.
         label nRefinementIterations_;

src/finiteVolume/Make/files

@@ -127,7 +127,6 @@ $(derivedFvPatchFields)/fan/fanFvPatchFields.C
 $(derivedFvPatchFields)/fanPressure/fanPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/buoyantPressure/buoyantPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
-$(derivedFvPatchFields)/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/fixedInternalValueFvPatchField/fixedInternalValueFvPatchFields.C
 $(derivedFvPatchFields)/fixedNormalSlip/fixedNormalSlipFvPatchFields.C
 $(derivedFvPatchFields)/fixedPressureCompressibleDensity/fixedPressureCompressibleDensityFvPatchScalarField.C

src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C

@@ -107,7 +107,7 @@ bool Foam::adjustPhi
         {
             massCorr = (massIn - fixedMassOut)/adjustableMassOut;
         }
-        else if (mag(fixedMassOut - massIn)/totalFlux > 1e-10)
+        else if (mag(fixedMassOut - massIn)/totalFlux > 1e-8)
         {
             FatalErrorIn
             (

src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C

@@ -59,7 +59,7 @@ addRadialActuationDiskAxialInertialResistance
     const Field<scalar> zoneCellVolumes(mesh().cellVolumes(), cells);
 
     const vector avgCentre = gSum(zoneCellVolumes*zoneCellCentres)/V();
-    const scalar maxR = mag(max(zoneCellCentres - avgCentre));
+    const scalar maxR = gMax(mag(zoneCellCentres - avgCentre));
 
     scalar intCoeffs =
         coeffs_[0]

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C

@@ -59,36 +59,37 @@ bool Foam::pimpleControl::criteriaSatisfied()
     bool firstIter = corr_ == 1;
 
     bool achieved = true;
-    const dictionary& solverDict = mesh_.solverPerformanceDict();
+    bool checked = false;    // safety that some checks were indeed performed
 
+    const dictionary& solverDict = mesh_.solverPerformanceDict();
     forAllConstIter(dictionary, solverDict, iter)
     {
         const word& variableName = iter().keyword();
-        label fieldI = applyToField(variableName);
+        const label fieldI = applyToField(variableName);
         if (fieldI != -1)
         {
             const List<lduMatrix::solverPerformance> sp(iter().stream());
             const scalar residual = sp.last().initialResidual();
 
+            checked = true;
+
             if (firstIter)
             {
                 residualControl_[fieldI].initialResidual =
                     sp.first().initialResidual();
             }
 
-            bool absCheck = residual < residualControl_[fieldI].absTol;
-
+            const bool absCheck = residual < residualControl_[fieldI].absTol;
             bool relCheck = false;
 
             scalar relative = 0.0;
             if (!firstIter)
             {
-                scalar iniRes =
+                const scalar iniRes =
                     residualControl_[fieldI].initialResidual
                   + ROOTVSMALL;
 
                 relative = residual/iniRes;
-
                 relCheck = relative < residualControl_[fieldI].relTol;
             }
 
@@ -110,7 +111,7 @@ bool Foam::pimpleControl::criteriaSatisfied()
         }
     }
 
-    return achieved;
+    return checked && achieved;
 }
 
 
@@ -129,7 +130,13 @@ Foam::pimpleControl::pimpleControl(fvMesh& mesh)
     if (nOuterCorr_ > 1)
     {
         Info<< nl;
-        if (!residualControl_.empty())
+        if (residualControl_.empty())
+        {
+            Info<< algorithmName_ << ": no residual control data found. "
+                << "Calculations will employ " << nOuterCorr_
+                << " corrector loops" << nl << endl;
+        }
+        else
         {
             Info<< algorithmName_ << ": max iterations = " << nOuterCorr_
                 << endl;
@@ -142,12 +149,6 @@ Foam::pimpleControl::pimpleControl(fvMesh& mesh)
             }
             Info<< endl;
         }
-        else
-        {
-            Info<< algorithmName_ << ": no residual control data found. " << nl
-                << "Calculations will employ " << nOuterCorr_
-                << " corrector loops" << nl << endl;
-        }
     }
     else
     {

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H

@@ -69,10 +69,10 @@ protected:
 
     // Protected Member Functions
 
-        //- Read constrols from fvSolution dictionary
+        //- Read controls from fvSolution dictionary
         virtual void read();
 
-        //- Return true if all convergence checks are satified
+        //- Return true if all convergence checks are satisfied
         virtual bool criteriaSatisfied();
 
         //- Disallow default bitwise copy construct

src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C

@@ -50,17 +50,20 @@ bool Foam::simpleControl::criteriaSatisfied()
     }
 
     bool achieved = true;
-    const dictionary& solverDict = mesh_.solverPerformanceDict();
+    bool checked = false;    // safety that some checks were indeed performed
 
+    const dictionary& solverDict = mesh_.solverPerformanceDict();
     forAllConstIter(dictionary, solverDict, iter)
     {
         const word& variableName = iter().keyword();
-        label fieldI = applyToField(variableName);
+        const label fieldI = applyToField(variableName);
         if (fieldI != -1)
         {
             const List<lduMatrix::solverPerformance> sp(iter().stream());
             const scalar residual = sp.first().initialResidual();
 
+            checked = true;
+
             bool absCheck = residual < residualControl_[fieldI].absTol;
             achieved = achieved && absCheck;
 
@@ -75,7 +78,7 @@ bool Foam::simpleControl::criteriaSatisfied()
         }
     }
 
-    return achieved;
+    return checked && achieved;
 }
 
 
@@ -90,7 +93,13 @@ Foam::simpleControl::simpleControl(fvMesh& mesh)
 
     Info<< nl;
 
-    if (residualControl_.size() > 0)
+    if (residualControl_.empty())
+    {
+        Info<< algorithmName_ << ": no convergence criteria found. "
+            << "Calculations will run for " << mesh_.time().endTime().value()
+            << " steps." << nl << endl;
+    }
+    else
     {
         Info<< algorithmName_ << ": convergence criteria" << nl;
         forAll(residualControl_, i)
@@ -101,12 +110,6 @@ Foam::simpleControl::simpleControl(fvMesh& mesh)
         }
         Info<< endl;
     }
-    else
-    {
-        Info<< algorithmName_ << ": no convergence criteria found. "
-            << "Calculations will run for " << mesh_.time().endTime().value()
-            << " steps." << nl << endl;
-    }
 }
 
 

src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.H

@@ -59,10 +59,10 @@ protected:
 
     // Protected Member Functions
 
-        //- Read constrols from fvSolution dictionary
+        //- Read controls from fvSolution dictionary
         void read();
 
-        //- Return true if all convergence checks are satified
+        //- Return true if all convergence checks are satisfied
         bool criteriaSatisfied();
 
         //- Disallow default bitwise copy construct

src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C

@@ -76,7 +76,8 @@ void Foam::solutionControl::read(const bool absTolOnly)
                 {
                     FatalErrorIn("bool Foam::solutionControl::read()")
                         << "Residual data for " << iter().keyword()
-                        << " must be specified as a dictionary";
+                        << " must be specified as a dictionary"
+                        << exit(FatalError);
                 }
             }
 

src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H

@@ -40,7 +40,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                  Class solutionControl Declaration
+                       Class solutionControl Declaration
 \*---------------------------------------------------------------------------*/
 
 class solutionControl
@@ -78,19 +78,19 @@ protected:
             //- Flag to indicate to solve for momentum
             bool momentumPredictor_;
 
-            //- Flag to indictae to solve using transonic algorithm
+            //- Flag to indicate to solve using transonic algorithm
             bool transonic_;
 
 
     // Protected Member Functions
 
-        //- Read constrols from fvSolution dictionary
+        //- Read controls from fvSolution dictionary
         virtual void read(const bool absTolOnly);
 
         //- Return index of field in residualControl_ if present
         virtual label applyToField(const word& fieldName) const;
 
-        //- Return true if all convergence checks are satified
+        //- Return true if all convergence checks are satisfied
         virtual bool criteriaSatisfied() = 0;
 
         //- Store previous iteration fields
@@ -142,7 +142,7 @@ public:
             //- Flag to indicate to solve for momentum
             inline bool momentumPredictor() const;
 
-            //- Flag to indictae to solve using transonic algorithm
+            //- Flag to indicate to solve using transonic algorithm
             inline bool transonic() const;
 };
 

src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.C

@@ -204,14 +204,13 @@ void mappedFieldFvPatchField<Type>::updateCoeffs()
     {
         case NEARESTCELL:
         {
-            const mapDistribute& distMap = mappedPatchBase::map();
             newValues = sampleField();
 
-            distMap.distribute(newValues);
+            this->distribute(newValues);
 
             break;
         }
-        case NEARESTPATCHFACE:
+        case NEARESTPATCHFACE: case NEARESTPATCHFACEAMI:
         {
             const label nbrPatchID =
                 nbrMesh.boundaryMesh().findPatchID(samplePatch());
@@ -228,36 +227,8 @@ void mappedFieldFvPatchField<Type>::updateCoeffs()
 
             const fieldType& nbrField = sampleField();
 
-            const mapDistribute& distMap = mappedPatchBase::map();
             newValues = nbrField.boundaryField()[nbrPatchID];
-            distMap.distribute(newValues);
-
-            break;
-        }
-        case mappedPatchBase::NEARESTPATCHFACEAMI:
-        {
-            const label nbrPatchID =
-                nbrMesh.boundaryMesh().findPatchID(samplePatch());
-
-            if (nbrPatchID < 0)
-            {
-                FatalErrorIn
-                (
-                    "void mappedFixedValueFvPatchField<Type>::updateCoeffs()"
-                )<< "Unable to find sample patch " << samplePatch()
-                 << " in region " << sampleRegion()
-                 << " for patch " << this->patch().name() << nl
-                 << abort(FatalError);
-            }
-
-//            const fieldType& nbrField = sampleField();
-//            newValues = mpp.AMI().interpolateToSource(nbrField);
-
-            notImplemented
-            (
-                "void mappedFieldFvPatchField<Type>::updateCoeffs() "
-                "with mappedPatchBase::NEARESTPATCHFACEAMI"
-            );
+            this->distribute(newValues);
 
             break;
         }
@@ -279,9 +250,7 @@ void mappedFieldFvPatchField<Type>::updateCoeffs()
                 }
             }
 
-            const mapDistribute& distMap = mappedPatchBase::map();
-            distMap.distribute(allValues);
-
+            this->distribute(allValues);
             newValues.transfer(allValues);
 
             break;

src/finiteVolume/fields/fvPatchFields/derived/mappedFixedInternalValue/mappedFixedInternalValueFvPatchField.C

@@ -114,42 +114,65 @@ void Foam::mappedFixedInternalValueFvPatchField<Type>::updateCoeffs()
     const mappedPatchBase& mpp =
         refCast<const mappedPatchBase>(this->patch().patch());
     const fvMesh& nbrMesh = refCast<const fvMesh>(mpp.sampleMesh());
-    const label samplePatchI = mpp.samplePolyPatch().index();
-    const fvPatch& nbrPatch = nbrMesh.boundary()[samplePatchI];
 
-
-    // Retrieve the neighbour field
-    const fvPatchField<Type>& nbrField =
-        nbrPatch.template lookupPatchField<FieldType, Type>
-        (
-            this->dimensionedInternalField().name()
-        );
-
-    // Retrieve the neighbour patch internal field
-    Field<Type> nbrIntFld(nbrField.patchInternalField());
+    Field<Type> nbrIntFld;
 
     switch (mpp.mode())
     {
-        case (mappedPatchBase::NEARESTPATCHFACEAMI):
+        case mappedPatchBase::NEARESTCELL:
         {
-            // Retrieve the neighbour patch internal field
-            nbrIntFld = mpp.AMI().interpolateToSource(nbrIntFld);
+            FatalErrorIn
+            (
+                "void mappedFixedValueFvPatchField<Type>::updateCoeffs()"
+            )   << "Cannot apply "
+                << mappedPatchBase::sampleModeNames_
+                   [
+                       mappedPatchBase::NEARESTCELL
+                   ]
+                << " mapping mode for patch " << this->patch().name()
+                << exit(FatalError);
+
+            break;
+        }
+        case mappedPatchBase::NEARESTPATCHFACE:
+        {
+            const label samplePatchI = mpp.samplePolyPatch().index();
+            const fvPatchField<Type>& nbrPatchField =
+                this->sampleField().boundaryField()[samplePatchI];
+            nbrIntFld = nbrPatchField.patchInternalField();
+            mpp.distribute(nbrIntFld);
+
+            break;
+        }
+        case mappedPatchBase::NEARESTFACE:
+        {
+            Field<Type> allValues(nbrMesh.nFaces(), pTraits<Type>::zero);
+
+            const FieldType& nbrField = this->sampleField();
+
+            forAll(nbrField.boundaryField(), patchI)
+            {
+                const fvPatchField<Type>& pf = nbrField.boundaryField()[patchI];
+                const Field<Type> pif(pf.patchInternalField());
+
+                label faceStart = pf.patch().start();
+
+                forAll(pf, faceI)
+                {
+                    allValues[faceStart++] = pif[faceI];
+                }
+            }
+
+            mpp.distribute(allValues);
+            nbrIntFld.transfer(allValues);
+
             break;
         }
         default:
         {
-            const mapDistribute& distMap = mpp.map();
-
-            mapDistribute::distribute
-            (
-                Pstream::defaultCommsType,
-                distMap.schedule(),
-                distMap.constructSize(),
-                distMap.subMap(),           // what to send
-                distMap.constructMap(),     // what to receive
-                nbrIntFld
-            );
-            break;
+            FatalErrorIn("mappedFixedValueFvPatchField<Type>::updateCoeffs()")
+                << "Unknown sampling mode: " << mpp.mode()
+                << abort(FatalError);
         }
     }
 

src/finiteVolume/fields/fvPatchFields/derived/mappedFixedValue/mappedFixedValueFvPatchField.C

@@ -283,30 +283,6 @@ void mappedFixedValueFvPatchField<Type>::updateCoeffs()
             break;
         }
         case mappedPatchBase::NEARESTPATCHFACE:
-        {
-            const mapDistribute& distMap = mpp.map();
-
-            const label nbrPatchID =
-                nbrMesh.boundaryMesh().findPatchID(mpp.samplePatch());
-
-            if (nbrPatchID < 0)
-            {
-                FatalErrorIn
-                (
-                    "void mappedFixedValueFvPatchField<Type>::updateCoeffs()"
-                )<< "Unable to find sample patch " << mpp.samplePatch()
-                 << " in region " << mpp.sampleRegion()
-                 << " for patch " << this->patch().name() << nl
-                 << abort(FatalError);
-            }
-
-            const fieldType& nbrField = sampleField();
-
-            newValues = nbrField.boundaryField()[nbrPatchID];
-            distMap.distribute(newValues);
-
-            break;
-        }
         case mappedPatchBase::NEARESTPATCHFACEAMI:
         {
             const label nbrPatchID =
@@ -324,9 +300,9 @@ void mappedFixedValueFvPatchField<Type>::updateCoeffs()
             }
 
             const fieldType& nbrField = sampleField();
-            newValues = nbrField.boundaryField()[nbrPatchID];
 
-            newValues = mpp.AMI().interpolateToSource(newValues);
+            newValues = nbrField.boundaryField()[nbrPatchID];
+            mpp.distribute(newValues);
 
             break;
         }

src/finiteVolume/fields/fvPatchFields/derived/mappedFixedValue/mappedFixedValueFvPatchField.H

@@ -71,7 +71,10 @@ class mappedFixedValueFvPatchField
 :
     public fixedValueFvPatchField<Type>
 {
-    // Private data
+
+protected:
+
+    // Protected data
 
         //- Name of field to sample - defaults to field associated with this
         //  patchField if not specified
@@ -90,7 +93,7 @@ class mappedFixedValueFvPatchField
         mutable autoPtr<interpolation<Type> > interpolator_;
 
 
-    // Private Member Functions
+    // Protected Member Functions
 
         //- Field to sample. Either on my or nbr mesh
         const GeometricField<Type, fvPatchField, volMesh>& sampleField() const;

src/finiteVolume/fields/fvPatchFields/derived/mappedFlowRate/mappedFlowRateFvPatchVectorField.C

@@ -28,7 +28,6 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "fvPatchFieldMapper.H"
 #include "mappedPatchBase.H"
-#include "mapDistribute.H"
 #include "surfaceFields.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -128,7 +127,7 @@ void Foam::mappedFlowRateFvPatchVectorField::updateCoeffs()
     scalarList phi =
         nbrPatch.lookupPatchField<surfaceScalarField, scalar>(nbrPhiName_);
 
-    mpp.map().distribute(phi);
+    mpp.distribute(phi);
 
 
     const surfaceScalarField& phiName =

src/finiteVolume/fields/fvPatchFields/derived/mappedVelocityFluxFixedValue/mappedVelocityFluxFixedValueFvPatchField.C

@@ -176,13 +176,10 @@ void Foam::mappedVelocityFluxFixedValueFvPatchField::updateCoeffs()
     (
         mappedVelocityFluxFixedValueFvPatchField::patch().patch()
     );
-    const mapDistribute& distMap = mpp.map();
     const fvMesh& nbrMesh = refCast<const fvMesh>(mpp.sampleMesh());
     const word& fieldName = dimensionedInternalField().name();
-    const volVectorField& UField = nbrMesh.lookupObject<volVectorField>
-    (
-        fieldName
-    );
+    const volVectorField& UField =
+        nbrMesh.lookupObject<volVectorField>(fieldName);
 
     surfaceScalarField& phiField = const_cast<surfaceScalarField&>
     (
@@ -213,26 +210,25 @@ void Foam::mappedVelocityFluxFixedValueFvPatchField::updateCoeffs()
                 }
             }
 
-            distMap.distribute(allUValues);
+            mpp.distribute(allUValues);
             newUValues.transfer(allUValues);
 
-            distMap.distribute(allPhiValues);
+            mpp.distribute(allPhiValues);
             newPhiValues.transfer(allPhiValues);
 
             break;
         }
         case mappedPolyPatch::NEARESTPATCHFACE:
+        case mappedPolyPatch::NEARESTPATCHFACEAMI:
         {
-            const label nbrPatchID = nbrMesh.boundaryMesh().findPatchID
-            (
-                mpp.samplePatch()
-            );
+            const label nbrPatchID =
+                nbrMesh.boundaryMesh().findPatchID(mpp.samplePatch());
 
             newUValues = UField.boundaryField()[nbrPatchID];
-            distMap.distribute(newUValues);
+            mpp.distribute(newUValues);
 
             newPhiValues = phiField.boundaryField()[nbrPatchID];
-            distMap.distribute(newPhiValues);
+            mpp.distribute(newPhiValues);
 
             break;
         }
@@ -242,8 +238,9 @@ void Foam::mappedVelocityFluxFixedValueFvPatchField::updateCoeffs()
             (
                 "mappedVelocityFluxFixedValueFvPatchField::"
                 "updateCoeffs()"
-            )   << "patch can only be used in NEARESTPATCHFACE or NEARESTFACE "
-                << "mode" << nl << abort(FatalError);
+            )   << "patch can only be used in NEARESTPATCHFACE, "
+                << "NEARESTPATCHFACEAMI or NEARESTFACE mode" << nl
+                << abort(FatalError);
         }
     }
 

src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C

@@ -1,177 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "multiphaseFixedFluxPressureFvPatchScalarField.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-#include "surfaceFields.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(p, iF),
-    UName_("U"),
-    phiName_("phi"),
-    rhoName_("rho")
-{}
-
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
-    const multiphaseFixedFluxPressureFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
-    UName_(ptf.UName_),
-    phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_)
-{}
-
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedGradientFvPatchScalarField(p, iF),
-    UName_(dict.lookupOrDefault<word>("U", "U")),
-    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
-{
-    if (dict.found("gradient"))
-    {
-        gradient() = scalarField("gradient", dict, p.size());
-        fixedGradientFvPatchScalarField::updateCoeffs();
-        fixedGradientFvPatchScalarField::evaluate();
-    }
-    else
-    {
-        fvPatchField<scalar>::operator=(patchInternalField());
-        gradient() = 0.0;
-    }
-}
-
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
-    const multiphaseFixedFluxPressureFvPatchScalarField& wbppsf
-)
-:
-    fixedGradientFvPatchScalarField(wbppsf),
-    UName_(wbppsf.UName_),
-    phiName_(wbppsf.phiName_),
-    rhoName_(wbppsf.rhoName_)
-{}
-
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
-    const multiphaseFixedFluxPressureFvPatchScalarField& wbppsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(wbppsf, iF),
-    UName_(wbppsf.UName_),
-    phiName_(wbppsf.phiName_),
-    rhoName_(wbppsf.rhoName_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::multiphaseFixedFluxPressureFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const fvPatchField<vector>& Up =
-        patch().lookupPatchField<volVectorField, vector>(UName_);
-
-    const surfaceScalarField& phi =
-        db().lookupObject<surfaceScalarField>(phiName_);
-
-    fvsPatchField<scalar> phip =
-        patch().patchField<surfaceScalarField, scalar>(phi);
-
-    if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
-    {
-        const fvPatchField<scalar>& rhop =
-            patch().lookupPatchField<volScalarField, scalar>(rhoName_);
-
-        phip /= rhop;
-    }
-
-    const fvsPatchField<scalar>& Dpp =
-        patch().lookupPatchField<surfaceScalarField, scalar>("Dp");
-
-    gradient() = (phip - (patch().Sf() & Up))/patch().magSf()/Dpp;
-
-    fixedGradientFvPatchScalarField::updateCoeffs();
-}
-
-
-void Foam::multiphaseFixedFluxPressureFvPatchScalarField::write
-(
-    Ostream& os
-) const
-{
-    fvPatchScalarField::write(os);
-    writeEntryIfDifferent<word>(os, "U", "U", UName_);
-    writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
-    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
-    gradient().writeEntry("gradient", os);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    makePatchTypeField
-    (
-        fvPatchScalarField,
-        multiphaseFixedFluxPressureFvPatchScalarField
-    );
-}
-
-// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H

@@ -1,154 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::multiphaseFixedFluxPressureFvPatchScalarField
-
-Description
-    Foam::multiphaseFixedFluxPressureFvPatchScalarField
-
-SourceFiles
-    multiphaseFixedFluxPressureFvPatchScalarField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef multiphaseFixedFluxPressureFvPatchScalarFields_H
-#define multiphaseFixedFluxPressureFvPatchScalarFields_H
-
-#include "fvPatchFields.H"
-#include "fixedGradientFvPatchFields.H"
-#include "Switch.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-              Class multiphaseFixedFluxPressureFvPatch Declaration
-\*---------------------------------------------------------------------------*/
-
-class multiphaseFixedFluxPressureFvPatchScalarField
-:
-    public fixedGradientFvPatchScalarField
-{
-    // Private data
-
-        //- Name of the velocity field
-        word UName_;
-
-        //- Name of the flux transporting the field
-        word phiName_;
-
-        //- Name of the density field used to normalise the mass flux
-        //  if neccessary
-        word rhoName_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("multiphaseFixedFluxPressure");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        multiphaseFixedFluxPressureFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        multiphaseFixedFluxPressureFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given
-        //  multiphaseFixedFluxPressureFvPatchScalarField onto a new patch
-        multiphaseFixedFluxPressureFvPatchScalarField
-        (
-            const multiphaseFixedFluxPressureFvPatchScalarField&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct as copy
-        multiphaseFixedFluxPressureFvPatchScalarField
-        (
-            const multiphaseFixedFluxPressureFvPatchScalarField&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new multiphaseFixedFluxPressureFvPatchScalarField(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        multiphaseFixedFluxPressureFvPatchScalarField
-        (
-            const multiphaseFixedFluxPressureFvPatchScalarField&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new multiphaseFixedFluxPressureFvPatchScalarField(*this, iF)
-            );
-        }
-
-
-    // Member functions
-
-        //- Update the coefficients associated with the patch field
-        virtual void updateCoeffs();
-
-        //- Write
-        virtual void write(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/selfContainedMapped/selfContainedMappedFixedValueFvPatchField.C

@@ -229,7 +229,6 @@ void selfContainedMappedFixedValueFvPatchField<Type>::updateCoeffs()
 
     const fvMesh& thisMesh = this->patch().boundaryMesh().mesh();
     const fvMesh& nbrMesh = refCast<const fvMesh>(sampleMesh());
-    const mapDistribute& distMap = mappedPatchBase::map();
 
     // Result of obtaining remote values
     Field<Type> newValues;
@@ -238,6 +237,8 @@ void selfContainedMappedFixedValueFvPatchField<Type>::updateCoeffs()
     {
         case NEARESTCELL:
         {
+            const mapDistribute& distMap = mappedPatchBase::map();
+
             if (interpolationScheme_ != interpolationCell<Type>::typeName)
             {
                 // Need to do interpolation so need cells to sample.
@@ -275,12 +276,10 @@ void selfContainedMappedFixedValueFvPatchField<Type>::updateCoeffs()
 
             break;
         }
-        case NEARESTPATCHFACE:
+        case NEARESTPATCHFACE: case NEARESTPATCHFACEAMI:
         {
-            const label nbrPatchID = nbrMesh.boundaryMesh().findPatchID
-            (
-                samplePatch()
-            );
+            const label nbrPatchID =
+                nbrMesh.boundaryMesh().findPatchID(samplePatch());
             if (nbrPatchID < 0)
             {
                 FatalErrorIn
@@ -297,7 +296,7 @@ void selfContainedMappedFixedValueFvPatchField<Type>::updateCoeffs()
             const fieldType& nbrField = sampleField();
 
             newValues = nbrField.boundaryField()[nbrPatchID];
-            distMap.distribute(newValues);
+            this->distribute(newValues);
 
             break;
         }
@@ -319,7 +318,7 @@ void selfContainedMappedFixedValueFvPatchField<Type>::updateCoeffs()
                 }
             }
 
-            distMap.distribute(allValues);
+            this->distribute(allValues);
 
             newValues.transfer(allValues);
 

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C

@@ -324,6 +324,8 @@ bool Foam::KinematicParcel<ParcelType>::move
         }
 
         p.age() += dt;
+
+        td.cloud().functions().postMove(p, cellI, dt);
     }
 
     return td.keepParticle;

src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.H

@@ -31,6 +31,7 @@ License
 #include "FacePostProcessing.H"
 #include "ParticleTracks.H"
 #include "PatchPostProcessing.H"
+#include "VoidFraction.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -40,7 +41,8 @@ License
                                                                               \
     makeCloudFunctionObjectType(FacePostProcessing, CloudType);               \
     makeCloudFunctionObjectType(ParticleTracks, CloudType);                   \
-    makeCloudFunctionObjectType(PatchPostProcessing, CloudType);
+    makeCloudFunctionObjectType(PatchPostProcessing, CloudType);              \
+    makeCloudFunctionObjectType(VoidFraction, CloudType);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.C

@@ -91,6 +91,18 @@ void Foam::CloudFunctionObject<CloudType>::postEvolve()
 }
 
 
+template<class CloudType>
+void Foam::CloudFunctionObject<CloudType>::postMove
+(
+    const typename CloudType::parcelType&,
+    const label,
+    const scalar
+)
+{
+    // do nothing
+}
+
+
 template<class CloudType>
 void Foam::CloudFunctionObject<CloudType>::postPatch
 (
@@ -98,14 +110,7 @@ void Foam::CloudFunctionObject<CloudType>::postPatch
     const label
 )
 {
-    notImplemented
-    (
-        "void Foam::CloudFunctionObject<CloudType>::postPatch"
-        "("
-            "const typename CloudType::parcelType&,"
-            "const label"
-        ")"
-    );
+    // do nothing
 }
 
 
@@ -115,13 +120,7 @@ void Foam::CloudFunctionObject<CloudType>::postFace
     const typename CloudType::parcelType&
 )
 {
-    notImplemented
-    (
-        "void Foam::CloudFunctionObject<CloudType>::postFace"
-        "("
-            "const typename CloudType::parcelType&"
-        ")"
-    );
+    // do nothing
 }
 
 

src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.H

@@ -129,6 +129,14 @@ public:
             //- Post-evolve hook
             virtual void postEvolve();
 
+            //- Post-move hook
+            virtual void postMove
+            (
+                const typename CloudType::parcelType& p,
+                const label cellI,
+                const scalar dt
+            );
+
             //- Post-patch hook
             virtual void postPatch
             (

src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.C

@@ -127,6 +127,21 @@ void Foam::CloudFunctionObjectList<CloudType>::postEvolve()
 }
 
 
+template<class CloudType>
+void Foam::CloudFunctionObjectList<CloudType>::postMove
+(
+    const typename CloudType::parcelType& p,
+    const label cellI,
+    const scalar dt
+)
+{
+    forAll(*this, i)
+    {
+        this->operator[](i).postMove(p, cellI, dt);
+    }
+}
+
+
 template<class CloudType>
 void Foam::CloudFunctionObjectList<CloudType>::postPatch
 (

src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.H

@@ -109,6 +109,14 @@ public:
             //- Post-evolve hook
             virtual void postEvolve();
 
+            //- Post-move hook
+            virtual void postMove
+            (
+                const typename CloudType::parcelType& p,
+                const label cellI,
+                const scalar dt
+            );
+
             //- Post-patch hook
             virtual void postPatch
             (

src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.C

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "VoidFraction.H"
+#
+// * * * * * * * * * * * * * Protectd Member Functions * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::VoidFraction<CloudType>::write()
+{
+    if (thetaPtr_.valid())
+    {
+        thetaPtr_->write();
+    }
+    else
+    {
+        FatalErrorIn("void Foam::VoidFraction<CloudType>::write()")
+            << "thetaPtr not valid" << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::VoidFraction<CloudType>::VoidFraction
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    CloudFunctionObject<CloudType>(owner),
+    thetaPtr_(NULL)
+{}
+
+
+template<class CloudType>
+Foam::VoidFraction<CloudType>::VoidFraction
+(
+    const VoidFraction<CloudType>& vf
+)
+:
+    CloudFunctionObject<CloudType>(vf),
+    thetaPtr_(NULL)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::VoidFraction<CloudType>::~VoidFraction()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::VoidFraction<CloudType>::preEvolve()
+{
+    if (thetaPtr_.valid())
+    {
+        thetaPtr_->internalField() = 0.0;
+    }
+    else
+    {
+        const fvMesh& mesh = this->owner().mesh();
+
+        thetaPtr_.reset
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    this->owner().name() + "Theta",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedScalar("zero", dimless, 0.0)
+            )
+        );
+    }
+}
+
+
+template<class CloudType>
+void Foam::VoidFraction<CloudType>::postEvolve()
+{
+    volScalarField& theta = thetaPtr_();
+
+    const fvMesh& mesh = this->owner().mesh();
+
+    theta.internalField() /= mesh.time().deltaTValue()*mesh.V();
+
+    CloudFunctionObject<CloudType>::postEvolve();
+}
+
+
+template<class CloudType>
+void Foam::VoidFraction<CloudType>::postMove
+(
+    const parcelType& p,
+    const label cellI,
+    const scalar dt
+)
+{
+    volScalarField& theta = thetaPtr_();
+
+    theta[cellI] += dt*p.nParticle()*p.volume();
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.H

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::VoidFraction
+
+Description
+    Creates particle void fraction field on carrier phase
+
+SourceFiles
+    VoidFraction.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef VoidFraction_H
+#define VoidFraction_H
+
+#include "CloudFunctionObject.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class VoidFraction Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class VoidFraction
+:
+    public CloudFunctionObject<CloudType>
+{
+    // Private Data
+
+        // Typedefs
+
+            //- Convenience typedef for parcel type
+            typedef typename CloudType::parcelType parcelType;
+
+
+        //- Void fraction field
+        autoPtr<volScalarField> thetaPtr_;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Write post-processing info
+        virtual void write();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("voidFraction");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        VoidFraction(const dictionary& dict, CloudType& owner);
+
+        //- Construct copy
+        VoidFraction(const VoidFraction<CloudType>& vf);
+
+        //- Construct and return a clone
+        virtual autoPtr<CloudFunctionObject<CloudType> > clone() const
+        {
+            return autoPtr<CloudFunctionObject<CloudType> >
+            (
+                new VoidFraction<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~VoidFraction();
+
+
+    // Member Functions
+
+        // Evaluation
+
+            //- Pre-evolve hook
+            virtual void preEvolve();
+
+            //- Post-evolve hook
+            virtual void postEvolve();
+
+            //- Post-move hook
+            virtual void postMove
+            (   
+                const parcelType& p,
+                const label cellI,
+                const scalar dt
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "VoidFraction.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C

@@ -220,11 +220,14 @@ Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
 :
     InjectionModel<CloudType>(im),
     injectionMethod_(im.injectionMethod_),
+    flowType_(im.flowType_),
     outerDiameter_(im.outerDiameter_),
     innerDiameter_(im.innerDiameter_),
     duration_(im.duration_),
     position_(im.position_),
     injectorCell_(im.injectorCell_),
+    tetFaceI_(im.tetFaceI_),
+    tetPtI_(im.tetPtI_),
     direction_(im.direction_),
     parcelsPerSecond_(im.parcelsPerSecond_),
     flowRateProfile_(im.flowRateProfile_().clone().ptr()),
@@ -235,9 +238,18 @@ Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
     tanVec2_(im.tanVec1_),
     normal_(im.normal_),
     UMag_(im.UMag_),
-    Cd_(im.Cd_().clone().ptr()),
-    Pinj_(im.Pinj_().clone().ptr())
-{}
+    Cd_(NULL),
+    Pinj_(NULL)
+{
+    if (im.Cd_.valid())
+    {
+        Cd_.reset(im.Cd_().clone().ptr());
+    }
+    if (im.Pinj_.valid())
+    {
+        Pinj_.reset(im.Pinj_().clone().ptr());
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C

@@ -114,8 +114,8 @@ bool Foam::SurfaceFilmModel<CloudType>::transferParcel
         "bool Foam::SurfaceFilmModel<CloudType>::transferParcel"
         "("
             "parcelType&, "
-            "const label, "
-            "const bool&"
+            "const polyPatch&, "
+            "bool&"
         ")"
     );
 
@@ -156,11 +156,9 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
         const label filmPatchI = filmPatches[i];
         const label primaryPatchI = primaryPatches[i];
 
-        const mappedPatchBase& mapPatch = filmModel.mappedPatches()[filmPatchI];
-
         const labelList& injectorCellsPatch = pbm[primaryPatchI].faceCells();
 
-        cacheFilmFields(filmPatchI, primaryPatchI, mapPatch, filmModel);
+        cacheFilmFields(filmPatchI, primaryPatchI, filmModel);
 
         const vectorField& Cf = mesh.C().boundaryField()[primaryPatchI];
         const vectorField& Sf = mesh.Sf().boundaryField()[primaryPatchI];
@@ -172,13 +170,11 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
             {
                 const label cellI = injectorCellsPatch[j];
 
-                // The position is at the cell centre, which could be
-                // in any tet of the decomposed cell, so arbitrarily
-                // choose the first face of the cell as the tetFace
-                // and the first point on the face after the base
-                // point as the tetPt.  The tracking will
-                // pick the cell consistent with the motion in the
-                // first tracking step.
+                // The position could bein any tet of the decomposed cell,
+                // so arbitrarily choose the first face of the cell as the
+                // tetFace and the first point on the face after the base
+                // point as the tetPt.  The tracking will pick the cell
+                // consistent with the motion in the first tracking step.
                 const label tetFaceI = this->owner().mesh().cells()[cellI][0];
                 const label tetPtI = 1;
 
@@ -208,14 +204,22 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
 
                 setParcelProperties(*pPtr, j);
 
-                // Check new parcel properties
-//                td.cloud().checkParcelProperties(*pPtr, 0.0, true);
-                td.cloud().checkParcelProperties(*pPtr, 0.0, false);
+                if (pPtr->nParticle() > 0.001)
+                {
+                    // Check new parcel properties
+    //                td.cloud().checkParcelProperties(*pPtr, 0.0, true);
+                    td.cloud().checkParcelProperties(*pPtr, 0.0, false);
 
-                // Add the new parcel to the cloud
-                td.cloud().addParticle(pPtr);
+                    // Add the new parcel to the cloud
+                    td.cloud().addParticle(pPtr);
 
-                nParcelsInjected_++;
+                    nParcelsInjected_++;
+                }
+                else
+                {
+                    // TODO: cache mass and re-distribute?
+                    delete pPtr;
+                }
             }
         }
     }
@@ -227,26 +231,25 @@ void Foam::SurfaceFilmModel<CloudType>::cacheFilmFields
 (
     const label filmPatchI,
     const label primaryPatchI,
-    const mappedPatchBase& mapPatch,
     const regionModels::surfaceFilmModels::surfaceFilmModel& filmModel
 )
 {
     massParcelPatch_ = filmModel.cloudMassTrans().boundaryField()[filmPatchI];
-    mapPatch.distribute(massParcelPatch_);
+    filmModel.toPrimary(filmPatchI, massParcelPatch_);
 
     diameterParcelPatch_ =
         filmModel.cloudDiameterTrans().boundaryField()[filmPatchI];
-    mapPatch.distribute(diameterParcelPatch_);
+    filmModel.toPrimary(filmPatchI, diameterParcelPatch_);
 
     UFilmPatch_ = filmModel.Us().boundaryField()[filmPatchI];
-    mapPatch.distribute(UFilmPatch_);
+    filmModel.toPrimary(filmPatchI, UFilmPatch_);
 
     rhoFilmPatch_ = filmModel.rho().boundaryField()[filmPatchI];
-    mapPatch.distribute(rhoFilmPatch_);
+    filmModel.toPrimary(filmPatchI, rhoFilmPatch_);
 
     deltaFilmPatch_[primaryPatchI] =
         filmModel.delta().boundaryField()[filmPatchI];
-    mapPatch.distribute(deltaFilmPatch_[primaryPatchI]);
+    filmModel.toPrimary(filmPatchI, deltaFilmPatch_[primaryPatchI]);
 }
 
 

src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H

@@ -116,7 +116,6 @@ protected:
         (
             const label filmPatchI,
             const label primaryPatchI,
-            const mappedPatchBase& mapPatch,
             const regionModels::surfaceFilmModels::surfaceFilmModel& filmModel
         );
 

src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C

@@ -643,7 +643,6 @@ void Foam::ThermoSurfaceFilm<CloudType>::cacheFilmFields
 (
     const label filmPatchI,
     const label primaryPatchI,
-    const mappedPatchBase& mapPatch,
     const regionModels::surfaceFilmModels::surfaceFilmModel& filmModel
 )
 {
@@ -651,15 +650,14 @@ void Foam::ThermoSurfaceFilm<CloudType>::cacheFilmFields
     (
         filmPatchI,
         primaryPatchI,
-        mapPatch,
         filmModel
     );
 
     TFilmPatch_ = filmModel.Ts().boundaryField()[filmPatchI];
-    mapPatch.distribute(TFilmPatch_);
+    filmModel.toPrimary(filmPatchI, TFilmPatch_);
 
     CpFilmPatch_ = filmModel.Cp().boundaryField()[filmPatchI];
-    mapPatch.distribute(CpFilmPatch_);
+    filmModel.toPrimary(filmPatchI, CpFilmPatch_);
 }
 
 

src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H

@@ -227,7 +227,6 @@ protected:
             (
                 const label filmPatchI,
                 const label primaryPatchI,
-                const mappedPatchBase& distMap,
                 const regionModels::surfaceFilmModels::surfaceFilmModel&
                     filmModel
             );

src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.C

@@ -110,6 +110,11 @@ Foam::SprayCloud<CloudType>::SprayCloud
 
         Info << "Average parcel mass: " << averageParcelMass_ << endl;
     }
+
+    if (this->solution().resetSourcesOnStartup())
+    {
+        CloudType::resetSourceTerms();
+    }
 }
 
 

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

@@ -46,6 +46,7 @@ Description
 #include "combineFaces.H"
 #include "IOmanip.H"
 #include "globalIndex.H"
+#include "DynamicField.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1990,6 +1991,7 @@ Foam::label Foam::autoLayerDriver::checkAndUnmark
 (
     const addPatchCellLayer& addLayer,
     const dictionary& meshQualityDict,
+    const bool additionalReporting,
     const List<labelPair>& baffles,
     const indirectPrimitivePatch& pp,
     const fvMesh& newMesh,
@@ -2032,6 +2034,12 @@ Foam::label Foam::autoLayerDriver::checkAndUnmark
     );
 
     // Check if any of the faces in error uses any face of an added cell
+    // - if additionalReporting print the few remaining areas for ease of
+    //   finding out where the problems are.
+
+    const label nReportMax = 10;
+    DynamicField<point> disabledFaceCentres(nReportMax);
+
     forAll(addedCells, oldPatchFaceI)
     {
         // Get the cells (in newMesh labels) per old patch face (in mesh
@@ -2052,12 +2060,58 @@ Foam::label Foam::autoLayerDriver::checkAndUnmark
                 )
             )
             {
+                if (additionalReporting && (nChanged < nReportMax))
+                {
+                    disabledFaceCentres.append
+                    (
+                        pp.faceCentres()[oldPatchFaceI]
+                    );
+                }
+
                 nChanged++;
             }
         }
     }
 
-    return returnReduce(nChanged, sumOp<label>());
+
+    label nChangedTotal = returnReduce(nChanged, sumOp<label>());
+
+    if (additionalReporting)
+    {
+        // Limit the number of points to be printed so that
+        // not too many points are reported when running in parallel
+        // Not accurate, i.e. not always nReportMax points are written,
+        // but this estimation avoid some communication here.
+        // The important thing, however, is that when only a few faces
+        // are disabled, their coordinates are printed, and this should be
+        // the case
+        label nReportLocal =
+            min
+            (
+                max(nChangedTotal / Pstream::nProcs(), 1),
+                min
+                (
+                    nChanged,
+                    max(nReportMax / Pstream::nProcs(), 1)
+                )
+            );
+
+        Pout<< "Checked mesh with layers. Disabled extrusion at " << endl;
+        for (label i=0; i < nReportLocal; i++)
+        {
+            Pout<< "    " << disabledFaceCentres[i] << endl;
+        }
+
+        label nReportTotal = returnReduce(nReportLocal, sumOp<label>());
+
+        if (nReportTotal < nChangedTotal)
+        {
+            Info<< "Suppressed disabled extrusion message for other "
+                << nChangedTotal - nReportTotal << " faces." << endl;
+        }
+    }
+
+    return nChangedTotal;
 }
 
 
@@ -2858,6 +2912,7 @@ void Foam::autoLayerDriver::addLayers
         (
             addLayer,
             meshQualityDict,
+            layerParams.additionalReporting(),
             newMeshBaffles,
             pp(),
             newMesh,

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

@@ -331,6 +331,7 @@ class autoLayerDriver
                 (
                     const addPatchCellLayer& addLayer,
                     const dictionary& motionDict,
+                    const bool additionalReporting,
                     const List<labelPair>& baffles,
                     const indirectPrimitivePatch& pp,
                     const fvMesh&,

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

@@ -137,98 +137,98 @@ Foam::labelList Foam::layerParameters::readNumLayers
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from dictionary
-Foam::layerParameters::layerParameters
-(
-    const PtrList<dictionary>& surfaceDicts,
-    const refinementSurfaces& refineSurfaces,
-    const labelList& globalToPatch,
-    const dictionary& dict,
-    const polyBoundaryMesh& boundaryMesh
-)
-:
-    numLayers_
-    (
-        readNumLayers
-        (
-            surfaceDicts,
-            refineSurfaces,
-            globalToPatch,
-            boundaryMesh
-        )
-    ),
-    expansionRatio_
-    (
-        numLayers_.size(),
-        readScalar(dict.lookup("expansionRatio"))
-    ),
-    relativeSizes_(false),
-    finalLayerThickness_
-    (
-        numLayers_.size(),
-        readScalar(dict.lookup("finalLayerRatio"))
-    ),
-    minThickness_
-    (
-        numLayers_.size(),
-        readScalar(dict.lookup("minThickness"))
-    ),
-    featureAngle_(readScalar(dict.lookup("featureAngle"))),
-    concaveAngle_
-    (
-        dict.lookupOrDefault("concaveAngle", defaultConcaveAngle)
-    ),
-    nGrow_(readLabel(dict.lookup("nGrow"))),
-    nSmoothSurfaceNormals_
-    (
-        readLabel(dict.lookup("nSmoothSurfaceNormals"))
-    ),
-    nSmoothNormals_(readLabel(dict.lookup("nSmoothNormals"))),
-    nSmoothThickness_(readLabel(dict.lookup("nSmoothThickness"))),
-    maxFaceThicknessRatio_
-    (
-        readScalar(dict.lookup("maxFaceThicknessRatio"))
-    ),
-    layerTerminationCos_
-    (
-        Foam::cos(degToRad(0.5*featureAngle_))
-    ),
-    maxThicknessToMedialRatio_
-    (
-        readScalar(dict.lookup("maxThicknessToMedialRatio"))
-    ),
-    minMedianAxisAngleCos_
-    (
-        Foam::cos(degToRad(readScalar(dict.lookup("minMedianAxisAngle"))))
-    ),
-    nBufferCellsNoExtrude_
-    (
-        readLabel(dict.lookup("nBufferCellsNoExtrude"))
-    ),
-    nSnap_(readLabel(dict.lookup("nSnap"))),
-    nLayerIter_(readLabel(dict.lookup("nLayerIter"))),
-    nRelaxedIter_(labelMax)
-{
-    if (nGrow_ > 0)
-    {
-        WarningIn("layerParameters::layerParameters(..)")
-            << "The nGrow parameter effect has changed with respect to 1.6.x."
-            << endl
-            << "Please set nGrow=0 for 1.6.x behaviour."
-            << endl;
-    }
-
-    dict.readIfPresent("nRelaxedIter", nRelaxedIter_);
-
-    if (nLayerIter_ < 0 || nRelaxedIter_ < 0)
-    {
-        FatalErrorIn("layerParameters::layerParameters(..)")
-            << "Layer iterations should be >= 0." << endl
-            << "nLayerIter:" << nLayerIter_
-            << " nRelaxedIter:" << nRelaxedIter_
-            << exit(FatalError);
-    }
-}
+//// Construct from dictionary
+//Foam::layerParameters::layerParameters
+//(
+//    const PtrList<dictionary>& surfaceDicts,
+//    const refinementSurfaces& refineSurfaces,
+//    const labelList& globalToPatch,
+//    const dictionary& dict,
+//    const polyBoundaryMesh& boundaryMesh
+//)
+//:
+//    numLayers_
+//    (
+//        readNumLayers
+//        (
+//            surfaceDicts,
+//            refineSurfaces,
+//            globalToPatch,
+//            boundaryMesh
+//        )
+//    ),
+//    expansionRatio_
+//    (
+//        numLayers_.size(),
+//        readScalar(dict.lookup("expansionRatio"))
+//    ),
+//    relativeSizes_(false),
+//    finalLayerThickness_
+//    (
+//        numLayers_.size(),
+//        readScalar(dict.lookup("finalLayerRatio"))
+//    ),
+//    minThickness_
+//    (
+//        numLayers_.size(),
+//        readScalar(dict.lookup("minThickness"))
+//    ),
+//    featureAngle_(readScalar(dict.lookup("featureAngle"))),
+//    concaveAngle_
+//    (
+//        dict.lookupOrDefault("concaveAngle", defaultConcaveAngle)
+//    ),
+//    nGrow_(readLabel(dict.lookup("nGrow"))),
+//    nSmoothSurfaceNormals_
+//    (
+//        readLabel(dict.lookup("nSmoothSurfaceNormals"))
+//    ),
+//    nSmoothNormals_(readLabel(dict.lookup("nSmoothNormals"))),
+//    nSmoothThickness_(readLabel(dict.lookup("nSmoothThickness"))),
+//    maxFaceThicknessRatio_
+//    (
+//        readScalar(dict.lookup("maxFaceThicknessRatio"))
+//    ),
+//    layerTerminationCos_
+//    (
+//        Foam::cos(degToRad(0.5*featureAngle_))
+//    ),
+//    maxThicknessToMedialRatio_
+//    (
+//        readScalar(dict.lookup("maxThicknessToMedialRatio"))
+//    ),
+//    minMedianAxisAngleCos_
+//    (
+//        Foam::cos(degToRad(readScalar(dict.lookup("minMedianAxisAngle"))))
+//    ),
+//    nBufferCellsNoExtrude_
+//    (
+//        readLabel(dict.lookup("nBufferCellsNoExtrude"))
+//    ),
+//    nSnap_(readLabel(dict.lookup("nSnap"))),
+//    nLayerIter_(readLabel(dict.lookup("nLayerIter"))),
+//    nRelaxedIter_(labelMax)
+//{
+//    if (nGrow_ > 0)
+//    {
+//        WarningIn("layerParameters::layerParameters(..)")
+//            << "The nGrow parameter effect has changed with respect to 1.6.x."
+//            << endl
+//            << "Please set nGrow=0 for 1.6.x behaviour."
+//            << endl;
+//    }
+//
+//    dict.readIfPresent("nRelaxedIter", nRelaxedIter_);
+//
+//    if (nLayerIter_ < 0 || nRelaxedIter_ < 0)
+//    {
+//        FatalErrorIn("layerParameters::layerParameters(..)")
+//            << "Layer iterations should be >= 0." << endl
+//            << "nLayerIter:" << nLayerIter_
+//            << " nRelaxedIter:" << nRelaxedIter_
+//            << exit(FatalError);
+//    }
+//}
 
 
 // Construct from dictionary
@@ -289,7 +289,8 @@ Foam::layerParameters::layerParameters
     ),
     nSnap_(readLabel(dict.lookup("nRelaxIter"))),
     nLayerIter_(readLabel(dict.lookup("nLayerIter"))),
-    nRelaxedIter_(labelMax)
+    nRelaxedIter_(labelMax),
+    additionalReporting_(dict.lookupOrDefault("additionalReporting", false))
 {
     if (nGrow_ > 0)
     {

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

@@ -105,6 +105,7 @@ class layerParameters
 
         label nRelaxedIter_;
 
+        Switch additionalReporting_;
 
     // Private Member Functions
 
@@ -128,15 +129,15 @@ public:
 
     // Constructors
 
-        //- Construct from dictionary - old syntax
-        layerParameters
-        (
-            const PtrList<dictionary>& surfaceDicts,
-            const refinementSurfaces& refineSurfaces,
-            const labelList& globalToPatch,
-            const dictionary& dict,
-            const polyBoundaryMesh& boundaryMesh
-        );
+        ////- Construct from dictionary - old syntax
+        //layerParameters
+        //(
+        //    const PtrList<dictionary>& surfaceDicts,
+        //    const refinementSurfaces& refineSurfaces,
+        //    const labelList& globalToPatch,
+        //    const dictionary& dict,
+        //    const polyBoundaryMesh& boundaryMesh
+        //);
 
         //- Construct from dictionary - new syntax
         layerParameters(const dictionary& dict, const polyBoundaryMesh&);
@@ -259,6 +260,12 @@ public:
                 return nSnap_;
             }
 
+            const Switch& additionalReporting() const
+            {
+                return additionalReporting_;
+            }
+
+
             // Overall
 
                 //- Number of overall layer addition iterations

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

@@ -97,9 +97,10 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::checkPatches
     boundBox bbSrc(srcPatch.points(), srcPatch.meshPoints());
     boundBox bbTgt(tgtPatch.points(), tgtPatch.meshPoints());
 
-    bbTgt.inflate(maxBoundsError);
+    boundBox bbTgtInf(bbTgt);
+    bbTgtInf.inflate(maxBoundsError);
 
-    if (!bbTgt.contains(bbSrc))
+    if (!bbTgtInf.contains(bbSrc))
     {
         WarningIn
         (
@@ -109,14 +110,46 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::checkPatches
                 "const primitivePatch&"
             ")"
         )   << "Source and target patch bounding boxes are not similar" << nl
-            << "    src span : " << bbSrc.span() << nl
-            << "    tgt span : " << bbTgt.span() << nl
-            << "    source: " << bbSrc << nl
-            << "    target: " << bbTgt << endl;
+            << "    source box span     : " << bbSrc.span() << nl
+            << "    target box span     : " << bbTgt.span() << nl
+            << "    source box          : " << bbSrc << nl
+            << "    target box          : " << bbTgt << nl
+            << "    inflated target box : " << bbTgtInf << endl;
     }
 }
 
 
+template<class SourcePatch, class TargetPatch>
+bool Foam::AMIInterpolation<SourcePatch, TargetPatch>::distributed
+(
+    const primitivePatch& srcPatch,
+    const primitivePatch& tgtPatch
+)
+{
+    if (Pstream::parRun())
+    {
+        List<label> facesPresentOnProc(Pstream::nProcs(), 0);
+        if ((srcPatch.size() > 0) || (tgtPatch.size() > 0))
+        {
+            facesPresentOnProc[Pstream::myProcNo()] = 1;
+        }
+        else
+        {
+            facesPresentOnProc[Pstream::myProcNo()] = 0;
+        }
+
+        Pstream::gatherList(facesPresentOnProc);
+        Pstream::scatterList(facesPresentOnProc);
+        if (sum(facesPresentOnProc) > 1)
+        {
+            return true;
+        }
+    }
+    
+    return false;
+}
+
+
 template<class SourcePatch, class TargetPatch>
 Foam::label
 Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcOverlappingProcs
@@ -1158,7 +1191,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
 {
     static label patchI = 0;
 
-    if (Pstream::parRun())
+    if (Pstream::parRun() && distributed(srcPatch, tgtPatch))
     {
         // convert local addressing to global addressing
         globalIndex globalSrcFaces(srcPatch.size());

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H

@@ -167,6 +167,13 @@ class AMIInterpolation
 
         // Parallel functionality
 
+            //- Return true if faces are spread over multiple domains
+            bool distributed
+            (
+                const primitivePatch& srcPatch,
+                const primitivePatch& tgtPatch
+            );
+
             label calcOverlappingProcs
             (
                 const List<treeBoundBoxList>& procBb,

src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H

@@ -112,6 +112,10 @@ public:
             NORMAL              // use face normal + distance
         };
 
+        static const NamedEnum<sampleMode, 4> sampleModeNames_;
+
+        static const NamedEnum<offsetMode, 3> offsetModeNames_;
+
 
     //- Helper class for finding nearest
     // Nearest:
@@ -142,13 +146,9 @@ public:
     };
 
 
-private:
+protected:
 
-    // Private data
-
-        static const NamedEnum<sampleMode, 4> sampleModeNames_;
-
-        static const NamedEnum<offsetMode, 3> offsetModeNames_;
+    // Protected data
 
         //- Patch to sample
         const polyPatch& patch_;
@@ -199,7 +199,7 @@ private:
             dictionary surfDict_;
 
 
-    // Private Member Functions
+    // Protected Member Functions
 
         //- Collect single list of samples and originating processor+face.
         void collectSamples
@@ -339,6 +339,26 @@ public:
         }
 
 
+        //- Wrapper around map/interpolate data distribution
+        template<class Type>
+        void reverseDistribute(List<Type>& lst) const
+        {
+            switch (mode_)
+            {
+                case NEARESTPATCHFACEAMI:
+                {
+                    lst = AMI().interpolateToTarget(Field<Type>(lst.xfer()));
+                    break;
+                }
+                default:
+                {
+                    label cSize = patch_.size();
+                    map().reverseDistribute(cSize, lst);
+                }
+            }
+        }
+
+
         //- Return reference to the parallel distribution map
         const mapDistribute& map() const
         {

src/parallel/decompose/Allwmake

@@ -53,4 +53,6 @@ fi
 
 wmake $makeType decompositionMethods
 
+wmake $makeType decompose
+
 # ----------------------------------------------------------------- end-of-file

src/parallel/decompose/decompose/Make/files

@@ -0,0 +1,3 @@
+fvFieldDecomposer.C
+
+LIB = $(FOAM_LIBBIN)/libdecompose

src/parallel/decompose/decompose/Make/options

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -llagrangian

src/parallel/decompose/decompose/fvFieldDecomposer.C

@@ -0,0 +1,216 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvFieldDecomposer.H"
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fvFieldDecomposer::patchFieldDecomposer::patchFieldDecomposer
+(
+    const labelUList& addressingSlice,
+    const label addressingOffset
+)
+:
+    directAddressing_(addressingSlice)
+{
+    forAll(directAddressing_, i)
+    {
+        // Subtract one to align addressing.
+        directAddressing_[i] -= addressingOffset + 1;
+    }
+}
+
+
+Foam::fvFieldDecomposer::processorVolPatchFieldDecomposer::
+processorVolPatchFieldDecomposer
+(
+    const fvMesh& mesh,
+    const labelUList& addressingSlice
+)
+:
+    directAddressing_(addressingSlice.size())
+{
+    const labelList& own = mesh.faceOwner();
+    const labelList& neighb = mesh.faceNeighbour();
+
+    forAll(directAddressing_, i)
+    {
+        // Subtract one to align addressing.
+        label ai = mag(addressingSlice[i]) - 1;
+
+        if (ai < neighb.size())
+        {
+            // This is a regular face. it has been an internal face
+            // of the original mesh and now it has become a face
+            // on the parallel boundary.
+            // Give face the value of the neighbour.
+
+            if (addressingSlice[i] >= 0)
+            {
+                // I have the owner so use the neighbour value
+                directAddressing_[i] = neighb[ai];
+            }
+            else
+            {
+                directAddressing_[i] = own[ai];
+            }
+        }
+        else
+        {
+            // This is a face that used to be on a cyclic boundary
+            // but has now become a parallel patch face. I cannot
+            // do the interpolation properly (I would need to look
+            // up the different (face) list of data), so I will
+            // just grab the value from the owner cell
+
+            directAddressing_[i] = own[ai];
+        }
+    }
+}
+
+
+Foam::fvFieldDecomposer::processorSurfacePatchFieldDecomposer::
+processorSurfacePatchFieldDecomposer
+(
+    const labelUList& addressingSlice
+)
+:
+    addressing_(addressingSlice.size()),
+    weights_(addressingSlice.size())
+{
+    forAll(addressing_, i)
+    {
+        addressing_[i].setSize(1);
+        weights_[i].setSize(1);
+
+        addressing_[i][0] = mag(addressingSlice[i]) - 1;
+        weights_[i][0] = sign(addressingSlice[i]);
+    }
+}
+
+
+Foam::fvFieldDecomposer::fvFieldDecomposer
+(
+    const fvMesh& completeMesh,
+    const fvMesh& procMesh,
+    const labelList& faceAddressing,
+    const labelList& cellAddressing,
+    const labelList& boundaryAddressing
+)
+:
+    completeMesh_(completeMesh),
+    procMesh_(procMesh),
+    faceAddressing_(faceAddressing),
+    cellAddressing_(cellAddressing),
+    boundaryAddressing_(boundaryAddressing),
+    patchFieldDecomposerPtrs_
+    (
+        procMesh_.boundary().size(),
+        static_cast<patchFieldDecomposer*>(NULL)
+    ),
+    processorVolPatchFieldDecomposerPtrs_
+    (
+        procMesh_.boundary().size(),
+        static_cast<processorVolPatchFieldDecomposer*>(NULL)
+    ),
+    processorSurfacePatchFieldDecomposerPtrs_
+    (
+        procMesh_.boundary().size(),
+        static_cast<processorSurfacePatchFieldDecomposer*>(NULL)
+    )
+{
+    forAll(boundaryAddressing_, patchi)
+    {
+        if
+        (
+            boundaryAddressing_[patchi] >= 0
+        && !isA<processorLduInterface>(procMesh.boundary()[patchi])
+        )
+        {
+            patchFieldDecomposerPtrs_[patchi] = new patchFieldDecomposer
+            (
+                procMesh_.boundary()[patchi].patchSlice(faceAddressing_),
+                completeMesh_.boundaryMesh()
+                [
+                    boundaryAddressing_[patchi]
+                ].start()
+            );
+        }
+        else
+        {
+            processorVolPatchFieldDecomposerPtrs_[patchi] =
+                new processorVolPatchFieldDecomposer
+                (
+                    completeMesh_,
+                    procMesh_.boundary()[patchi].patchSlice(faceAddressing_)
+                );
+
+            processorSurfacePatchFieldDecomposerPtrs_[patchi] =
+                new processorSurfacePatchFieldDecomposer
+                (
+                    static_cast<const labelUList&>
+                    (
+                        procMesh_.boundary()[patchi].patchSlice
+                        (
+                            faceAddressing_
+                        )
+                    )
+                );
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::fvFieldDecomposer::~fvFieldDecomposer()
+{
+    forAll(patchFieldDecomposerPtrs_, patchi)
+    {
+        if (patchFieldDecomposerPtrs_[patchi])
+        {
+            delete patchFieldDecomposerPtrs_[patchi];
+        }
+    }
+
+    forAll(processorVolPatchFieldDecomposerPtrs_, patchi)
+    {
+        if (processorVolPatchFieldDecomposerPtrs_[patchi])
+        {
+            delete processorVolPatchFieldDecomposerPtrs_[patchi];
+        }
+    }
+
+    forAll(processorSurfacePatchFieldDecomposerPtrs_, patchi)
+    {
+        if (processorSurfacePatchFieldDecomposerPtrs_[patchi])
+        {
+            delete processorSurfacePatchFieldDecomposerPtrs_[patchi];
+        }
+    }
+}
+
+// ************************************************************************* //

src/parallel/decompose/decompose/fvFieldDecomposer.H

@@ -0,0 +1,275 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fvFieldDecomposer
+
+Description
+    Finite Volume volume and surface field decomposer.
+
+SourceFiles
+    fvFieldDecomposer.C
+    fvFieldDecomposerDecomposeFields.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fvFieldDecomposer_H
+#define fvFieldDecomposer_H
+
+#include "fvMesh.H"
+#include "fvPatchFieldMapper.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class IOobjectList;
+
+/*---------------------------------------------------------------------------*\
+                    Class fvFieldDecomposer Declaration
+\*---------------------------------------------------------------------------*/
+
+class fvFieldDecomposer
+{
+public:
+
+        //- Patch field decomposer class
+        class patchFieldDecomposer
+        :
+            public fvPatchFieldMapper
+        {
+            // Private data
+
+                labelList directAddressing_;
+
+        public:
+
+            // Constructors
+
+                //- Construct given addressing
+                patchFieldDecomposer
+                (
+                    const labelUList& addressingSlice,
+                    const label addressingOffset
+                );
+
+
+            // Member functions
+
+                label size() const
+                {
+                    return directAddressing_.size();
+                }
+
+                bool direct() const
+                {
+                    return true;
+                }
+
+                const labelUList& directAddressing() const
+                {
+                    return directAddressing_;
+                }
+        };
+
+
+        //- Processor patch field decomposer class. Maps either owner or
+        //  neighbour data (no interpolate anymore - processorFvPatchField
+        //  holds neighbour data)
+        class processorVolPatchFieldDecomposer
+        :
+            public fvPatchFieldMapper
+        {
+            // Private data
+
+                labelList directAddressing_;
+
+        public:
+
+            //- Construct given addressing
+            processorVolPatchFieldDecomposer
+            (
+                const fvMesh& mesh,
+                const labelUList& addressingSlice
+            );
+
+
+            // Member functions
+
+                label size() const
+                {
+                    return directAddressing_.size();
+                }
+
+                bool direct() const
+                {
+                    return true;
+                }
+
+                const labelUList& directAddressing() const
+                {
+                    return directAddressing_;
+                }
+        };
+
+
+        //- Processor patch field decomposer class. Surface field is assumed
+        //  to have direction (so manipulates sign when mapping)
+        class processorSurfacePatchFieldDecomposer
+        :
+            public fvPatchFieldMapper
+        {
+            labelListList addressing_;
+            scalarListList weights_;
+
+        public:
+
+            //- Construct given addressing
+            processorSurfacePatchFieldDecomposer
+            (
+                const labelUList& addressingSlice
+            );
+
+
+            // Member functions
+
+                label size() const
+                {
+                    return addressing_.size();
+                }
+
+                bool direct() const
+                {
+                    return false;
+                }
+
+                const labelListList& addressing() const
+                {
+                    return addressing_;
+                }
+
+                const scalarListList& weights() const
+                {
+                    return weights_;
+                }
+        };
+
+
+private:
+
+    // Private data
+
+        //- Reference to complete mesh
+        const fvMesh& completeMesh_;
+
+        //- Reference to processor mesh
+        const fvMesh& procMesh_;
+
+        //- Reference to face addressing
+        const labelList& faceAddressing_;
+
+        //- Reference to cell addressing
+        const labelList& cellAddressing_;
+
+        //- Reference to boundary addressing
+        const labelList& boundaryAddressing_;
+
+        //- List of patch field decomposers
+        List<patchFieldDecomposer*> patchFieldDecomposerPtrs_;
+
+        List<processorVolPatchFieldDecomposer*>
+            processorVolPatchFieldDecomposerPtrs_;
+
+        List<processorSurfacePatchFieldDecomposer*>
+            processorSurfacePatchFieldDecomposerPtrs_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        fvFieldDecomposer(const fvFieldDecomposer&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const fvFieldDecomposer&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        fvFieldDecomposer
+        (
+            const fvMesh& completeMesh,
+            const fvMesh& procMesh,
+            const labelList& faceAddressing,
+            const labelList& cellAddressing,
+            const labelList& boundaryAddressing
+        );
+
+
+    //- Destructor
+    ~fvFieldDecomposer();
+
+
+    // Member Functions
+
+        //- Decompose volume field
+        template<class Type>
+        tmp<GeometricField<Type, fvPatchField, volMesh> >
+        decomposeField
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& field,
+            const bool allowUnknownPatchFields = false
+        ) const;
+
+        //- Decompose surface field
+        template<class Type>
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+        decomposeField
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& field
+        ) const;
+
+        template<class GeoField>
+        void decomposeFields(const PtrList<GeoField>& fields) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "fvFieldDecomposerDecomposeFields.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/parallel/decompose/decompose/fvFieldDecomposerDecomposeFields.C

@@ -0,0 +1,286 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvFieldDecomposer.H"
+#include "processorFvPatchField.H"
+#include "processorFvsPatchField.H"
+#include "processorCyclicFvPatchField.H"
+#include "processorCyclicFvsPatchField.H"
+#include "emptyFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::fvFieldDecomposer::decomposeField
+(
+    const GeometricField<Type, fvPatchField, volMesh>& field,
+    const bool allowUnknownPatchFields
+) const
+{
+    // Create and map the internal field values
+    Field<Type> internalField(field.internalField(), cellAddressing_);
+
+    // Create and map the patch field values
+    PtrList<fvPatchField<Type> > patchFields(boundaryAddressing_.size());
+
+    forAll(boundaryAddressing_, patchi)
+    {
+        if (patchFieldDecomposerPtrs_[patchi])
+        {
+            patchFields.set
+            (
+                patchi,
+                fvPatchField<Type>::New
+                (
+                    field.boundaryField()[boundaryAddressing_[patchi]],
+                    procMesh_.boundary()[patchi],
+                    DimensionedField<Type, volMesh>::null(),
+                    *patchFieldDecomposerPtrs_[patchi]
+                )
+            );
+        }
+        else if (isA<processorCyclicFvPatch>(procMesh_.boundary()[patchi]))
+        {
+            patchFields.set
+            (
+                patchi,
+                new processorCyclicFvPatchField<Type>
+                (
+                    procMesh_.boundary()[patchi],
+                    DimensionedField<Type, volMesh>::null(),
+                    Field<Type>
+                    (
+                        field.internalField(),
+                        *processorVolPatchFieldDecomposerPtrs_[patchi]
+                    )
+                )
+            );
+        }
+        else if (isA<processorFvPatch>(procMesh_.boundary()[patchi]))
+        {
+            patchFields.set
+            (
+                patchi,
+                new processorFvPatchField<Type>
+                (
+                    procMesh_.boundary()[patchi],
+                    DimensionedField<Type, volMesh>::null(),
+                    Field<Type>
+                    (
+                        field.internalField(),
+                        *processorVolPatchFieldDecomposerPtrs_[patchi]
+                    )
+                )
+            );
+        }
+        else if (allowUnknownPatchFields)
+        {
+            patchFields.set
+            (
+                patchi,
+                new emptyFvPatchField<Type>
+                (
+                    procMesh_.boundary()[patchi],
+                    DimensionedField<Type, volMesh>::null()
+                )
+            );
+        }
+        else
+        {
+            FatalErrorIn("fvFieldDecomposer::decomposeField()")
+                << "Unknown type." << abort(FatalError);
+        }
+    }
+
+    // Create the field for the processor
+    return tmp<GeometricField<Type, fvPatchField, volMesh> >
+    (
+        new GeometricField<Type, fvPatchField, volMesh>
+        (
+            IOobject
+            (
+                field.name(),
+                procMesh_.time().timeName(),
+                procMesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            procMesh_,
+            field.dimensions(),
+            internalField,
+            patchFields
+        )
+    );
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
+Foam::fvFieldDecomposer::decomposeField
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& field
+) const
+{
+    labelList mapAddr
+    (
+        labelList::subList
+        (
+            faceAddressing_,
+            procMesh_.nInternalFaces()
+        )
+    );
+    forAll(mapAddr, i)
+    {
+        mapAddr[i] -= 1;
+    }
+
+    // Create and map the internal field values
+    Field<Type> internalField
+    (
+        field.internalField(),
+        mapAddr
+    );
+
+    // Problem with addressing when a processor patch picks up both internal
+    // faces and faces from cyclic boundaries. This is a bit of a hack, but
+    // I cannot find a better solution without making the internal storage
+    // mechanism for surfaceFields correspond to the one of faces in polyMesh
+    // (i.e. using slices)
+    Field<Type> allFaceField(field.mesh().nFaces());
+
+    forAll(field.internalField(), i)
+    {
+        allFaceField[i] = field.internalField()[i];
+    }
+
+    forAll(field.boundaryField(), patchi)
+    {
+        const Field<Type> & p = field.boundaryField()[patchi];
+
+        const label patchStart = field.mesh().boundaryMesh()[patchi].start();
+
+        forAll(p, i)
+        {
+            allFaceField[patchStart + i] = p[i];
+        }
+    }
+
+    // Create and map the patch field values
+    PtrList<fvsPatchField<Type> > patchFields(boundaryAddressing_.size());
+
+    forAll(boundaryAddressing_, patchi)
+    {
+        if (patchFieldDecomposerPtrs_[patchi])
+        {
+            patchFields.set
+            (
+                patchi,
+                fvsPatchField<Type>::New
+                (
+                    field.boundaryField()[boundaryAddressing_[patchi]],
+                    procMesh_.boundary()[patchi],
+                    DimensionedField<Type, surfaceMesh>::null(),
+                    *patchFieldDecomposerPtrs_[patchi]
+                )
+            );
+        }
+        else if (isA<processorCyclicFvPatch>(procMesh_.boundary()[patchi]))
+        {
+            patchFields.set
+            (
+                patchi,
+                new processorCyclicFvsPatchField<Type>
+                (
+                    procMesh_.boundary()[patchi],
+                    DimensionedField<Type, surfaceMesh>::null(),
+                    Field<Type>
+                    (
+                        allFaceField,
+                        *processorSurfacePatchFieldDecomposerPtrs_[patchi]
+                    )
+                )
+            );
+        }
+        else if (isA<processorFvPatch>(procMesh_.boundary()[patchi]))
+        {
+            patchFields.set
+            (
+                patchi,
+                new processorFvsPatchField<Type>
+                (
+                    procMesh_.boundary()[patchi],
+                    DimensionedField<Type, surfaceMesh>::null(),
+                    Field<Type>
+                    (
+                        allFaceField,
+                        *processorSurfacePatchFieldDecomposerPtrs_[patchi]
+                    )
+                )
+            );
+        }
+        else
+        {
+            FatalErrorIn("fvFieldDecomposer::decomposeField()")
+                << "Unknown type." << abort(FatalError);
+        }
+    }
+
+    // Create the field for the processor
+    return tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            IOobject
+            (
+                field.name(),
+                procMesh_.time().timeName(),
+                procMesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            procMesh_,
+            field.dimensions(),
+            internalField,
+            patchFields
+        )
+    );
+}
+
+
+template<class GeoField>
+void Foam::fvFieldDecomposer::decomposeFields
+(
+    const PtrList<GeoField>& fields
+) const
+{
+    forAll(fields, fieldI)
+    {
+        decomposeField(fields[fieldI])().write();
+    }
+}
+
+
+// ************************************************************************* //

src/parallel/decompose/ptscotchDecomp/Make/options

@@ -3,9 +3,10 @@ sinclude $(RULES)/mplib$(WM_MPLIB)
 
 EXE_INC = \
     $(PFLAGS) $(PINC) \
+    -I$(SCOTCH_ROOT)/include \
     -I$(SCOTCH_ARCH_PATH)/include/$(FOAM_MPI) \
     -I/usr/include/scotch \
     -I../decompositionMethods/lnInclude
 
 LIB_LIBS = \
-    -L$(FOAM_EXT_LIBBIN)/$(FOAM_MPI) -lptscotch -lptscotcherrexit -lrt
+    -L$(SCOTCH_ROOT)/lib -L$(FOAM_EXT_LIBBIN)/$(FOAM_MPI) -lptscotch -lptscotcherrexit -lrt

src/parallel/decompose/scotchDecomp/Make/options

@@ -7,9 +7,10 @@ sinclude $(RULES)/mplib$(WM_MPLIB)
 
 EXE_INC = \
     $(PFLAGS) $(PINC) \
+    -I$(SCOTCH_ROOT)/include \
     -I$(SCOTCH_ARCH_PATH)/include \
     -I/usr/include/scotch \
     -I../decompositionMethods/lnInclude
 
 LIB_LIBS = \
-    -L$(FOAM_EXT_LIBBIN) -lscotch -lscotcherrexit -lrt
+    -L$(SCOTCH_ROOT)/lib -L$(FOAM_EXT_LIBBIN) -lscotch -lscotcherrexit -lrt

src/parallel/reconstruct/reconstruct/fvFieldReconstructor.C

@@ -42,7 +42,40 @@ Foam::fvFieldReconstructor::fvFieldReconstructor
     cellProcAddressing_(cellProcAddressing),
     boundaryProcAddressing_(boundaryProcAddressing),
     nReconstructed_(0)
-{}
+{
+    forAll(procMeshes_, procI)
+    {
+        const fvMesh& procMesh = procMeshes_[procI];
+        if
+        (
+            faceProcAddressing[procI].size() != procMesh.nFaces()
+         || cellProcAddressing[procI].size() != procMesh.nCells()
+         || boundaryProcAddressing[procI].size() != procMesh.boundary().size()
+        )
+        {
+            FatalErrorIn
+            (
+                "fvFieldReconstructor::fvFieldReconstructor\n"
+                "(\n"
+                "   fvMesh&,\n"
+                "   const PtrList<fvMesh>&,\n"
+                "   const PtrList<labelIOList>&,\n"
+                "   const PtrList<labelIOList>&,\n"
+                "   const PtrList<labelIOList>&\n"
+                ")"
+            )   << "Size of maps does not correspond to size of mesh"
+                << " for processor " << procI << endl
+                << "faceProcAddressing : " << faceProcAddressing[procI].size()
+                << " nFaces : " << procMesh.nFaces() << endl
+                << "cellProcAddressing : " << cellProcAddressing[procI].size()
+                << " nCell : " << procMesh.nCells() << endl
+                << "boundaryProcAddressing : "
+                << boundaryProcAddressing[procI].size()
+                << " nFaces : " << procMesh.boundary().size()
+                << exit(FatalError);
+        }
+    }
+}
 
 
 // ************************************************************************* //

src/parallel/reconstruct/reconstruct/fvFieldReconstructor.H

@@ -86,6 +86,7 @@ class fvFieldReconstructor
 
 public:
 
+        //- Mapper for sizing only - does not do any actual mapping.
         class fvPatchFieldReconstructor
         :
             public fvPatchFieldMapper
@@ -146,19 +147,48 @@ public:
         //- Reconstruct volume internal field
         template<class Type>
         tmp<DimensionedField<Type, volMesh> >
-        reconstructFvVolumeInternalField(const IOobject& fieldIoObject);
+        reconstructFvVolumeInternalField
+        (
+            const IOobject& fieldIoObject,
+            const PtrList<DimensionedField<Type, volMesh> >& procFields
+        ) const;
+
+        //- Read and reconstruct volume internal field
+        template<class Type>
+        tmp<DimensionedField<Type, volMesh> >
+        reconstructFvVolumeInternalField(const IOobject& fieldIoObject) const;
+
 
         //- Reconstruct volume field
         template<class Type>
         tmp<GeometricField<Type, fvPatchField, volMesh> >
-        reconstructFvVolumeField(const IOobject& fieldIoObject);
+        reconstructFvVolumeField
+        (
+            const IOobject& fieldIoObject,
+            const PtrList<GeometricField<Type, fvPatchField, volMesh> >&
+        ) const;
+
+        //- Read and reconstruct volume field
+        template<class Type>
+        tmp<GeometricField<Type, fvPatchField, volMesh> >
+        reconstructFvVolumeField(const IOobject& fieldIoObject) const;
+
 
         //- Reconstruct surface field
         template<class Type>
         tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
-        reconstructFvSurfaceField(const IOobject& fieldIoObject);
+        reconstructFvSurfaceField
+        (
+            const IOobject& fieldIoObject,
+            const PtrList<GeometricField<Type, fvsPatchField, surfaceMesh> >&
+        ) const;
 
-        //- Reconstruct and write all/selected volume internal fields
+        //- Read and reconstruct surface field
+        template<class Type>
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+        reconstructFvSurfaceField(const IOobject& fieldIoObject) const;
+
+        //- Read, reconstruct and write all/selected volume internal fields
         template<class Type>
         void reconstructFvVolumeInternalFields
         (
@@ -166,7 +196,7 @@ public:
             const HashSet<word>& selectedFields
         );
 
-        //- Reconstruct and write all/selected volume fields
+        //- Read, reconstruct and write all/selected volume fields
         template<class Type>
         void reconstructFvVolumeFields
         (
@@ -174,7 +204,7 @@ public:
             const HashSet<word>& selectedFields
         );
 
-        //- Reconstruct and write all/selected surface fields
+        //- Read, reconstruct and write all/selected surface fields
         template<class Type>
         void reconstructFvSurfaceFields
         (

src/parallel/reconstruct/reconstruct/fvFieldReconstructorReconstructFields.C

@@ -33,12 +33,48 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+Foam::tmp<Foam::DimensionedField<Type, Foam::volMesh> >
+Foam::fvFieldReconstructor::reconstructFvVolumeInternalField
+(
+    const IOobject& fieldIoObject,
+    const PtrList<DimensionedField<Type, volMesh> >& procFields
+) const
+{
+    // Create the internalField
+    Field<Type> internalField(mesh_.nCells());
+
+    forAll(procMeshes_, procI)
+    {
+        const DimensionedField<Type, volMesh>& procField = procFields[procI];
+
+        // Set the cell values in the reconstructed field
+        internalField.rmap
+        (
+            procField.field(),
+            cellProcAddressing_[procI]
+        );
+    }
+
+    return tmp<DimensionedField<Type, volMesh> >
+    (
+        new DimensionedField<Type, volMesh>
+        (
+            fieldIoObject,
+            mesh_,
+            procFields[0].dimensions(),
+            internalField
+        )
+    );
+}
+
+
 template<class Type>
 Foam::tmp<Foam::DimensionedField<Type, Foam::volMesh> >
 Foam::fvFieldReconstructor::reconstructFvVolumeInternalField
 (
     const IOobject& fieldIoObject
-)
+) const
 {
     // Read the field for all the processors
     PtrList<DimensionedField<Type, volMesh> > procFields
@@ -67,37 +103,17 @@ Foam::fvFieldReconstructor::reconstructFvVolumeInternalField
     }
 
 
-    // Create the internalField
-    Field<Type> internalField(mesh_.nCells());
-
-    forAll(procMeshes_, procI)
-    {
-        const DimensionedField<Type, volMesh>& procField = procFields[procI];
-
-        // Set the cell values in the reconstructed field
-        internalField.rmap
-        (
-            procField.field(),
-            cellProcAddressing_[procI]
-        );
-    }
-
-    return tmp<DimensionedField<Type, volMesh> >
+    return reconstructFvVolumeInternalField
     (
-        new DimensionedField<Type, volMesh>
+        IOobject
         (
-            IOobject
-            (
-                fieldIoObject.name(),
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
+            fieldIoObject.name(),
+            mesh_.time().timeName(),
             mesh_,
-            procFields[0].dimensions(),
-            internalField
-        )
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        procFields
     );
 }
 
@@ -106,44 +122,17 @@ template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
 Foam::fvFieldReconstructor::reconstructFvVolumeField
 (
-    const IOobject& fieldIoObject
-)
+    const IOobject& fieldIoObject,
+    const PtrList<GeometricField<Type, fvPatchField, volMesh> >& procFields
+) const
 {
-    // Read the field for all the processors
-    PtrList<GeometricField<Type, fvPatchField, volMesh> > procFields
-    (
-        procMeshes_.size()
-    );
-
-    forAll(procMeshes_, procI)
-    {
-        procFields.set
-        (
-            procI,
-            new GeometricField<Type, fvPatchField, volMesh>
-            (
-                IOobject
-                (
-                    fieldIoObject.name(),
-                    procMeshes_[procI].time().timeName(),
-                    procMeshes_[procI],
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                ),
-                procMeshes_[procI]
-            )
-        );
-    }
-
-
     // Create the internalField
     Field<Type> internalField(mesh_.nCells());
 
     // Create the patch fields
     PtrList<fvPatchField<Type> > patchFields(mesh_.boundary().size());
 
-
-    forAll(procMeshes_, procI)
+    forAll(procFields, procI)
     {
         const GeometricField<Type, fvPatchField, volMesh>& procField =
             procFields[procI];
@@ -163,7 +152,7 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
 
             // Get addressing slice for this patch
             const labelList::subList cp =
-                procMeshes_[procI].boundary()[patchI].patchSlice
+                procField.mesh().boundary()[patchI].patchSlice
                 (
                     faceProcAddressing_[procI]
                 );
@@ -198,11 +187,32 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
 
                 forAll(cp, faceI)
                 {
+                    // Check
+                    if (cp[faceI] <= 0)
+                    {
+                        FatalErrorIn
+                        (
+                            "fvFieldReconstructor::reconstructFvVolumeField\n"
+                            "(\n"
+                            "    const IOobject&,\n"
+                            "    const PtrList<GeometricField<Type,"
+                            " fvPatchField, volMesh> >&\n"
+                            ") const\n"
+                        )   << "Processor " << procI
+                            << " patch "
+                            << procField.mesh().boundary()[patchI].name()
+                            << " face " << faceI
+                            << " originates from reversed face since "
+                            << cp[faceI]
+                            << exit(FatalError);
+                    }
+
                     // Subtract one to take into account offsets for
                     // face direction.
                     reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart;
                 }
 
+
                 patchFields[curBPatch].rmap
                 (
                     procField.boundaryField()[patchI],
@@ -283,14 +293,7 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
     (
         new GeometricField<Type, fvPatchField, volMesh>
         (
-            IOobject
-            (
-                fieldIoObject.name(),
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
+            fieldIoObject,
             mesh_,
             procFields[0].dimensions(),
             internalField,
@@ -301,14 +304,14 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
 
 
 template<class Type>
-Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
-Foam::fvFieldReconstructor::reconstructFvSurfaceField
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::fvFieldReconstructor::reconstructFvVolumeField
 (
     const IOobject& fieldIoObject
-)
+) const
 {
     // Read the field for all the processors
-    PtrList<GeometricField<Type, fvsPatchField, surfaceMesh> > procFields
+    PtrList<GeometricField<Type, fvPatchField, volMesh> > procFields
     (
         procMeshes_.size()
     );
@@ -318,7 +321,7 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField
         procFields.set
         (
             procI,
-            new GeometricField<Type, fvsPatchField, surfaceMesh>
+            new GeometricField<Type, fvPatchField, volMesh>
             (
                 IOobject
                 (
@@ -333,7 +336,29 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField
         );
     }
 
+    return reconstructFvVolumeField
+    (
+        IOobject
+        (
+            fieldIoObject.name(),
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        procFields
+    );
+}
 
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
+Foam::fvFieldReconstructor::reconstructFvSurfaceField
+(
+    const IOobject& fieldIoObject,
+    const PtrList<GeometricField<Type, fvsPatchField, surfaceMesh> >& procFields
+) const
+{
     // Create the internalField
     Field<Type> internalField(mesh_.nInternalFaces());
 
@@ -503,14 +528,7 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField
     (
         new GeometricField<Type, fvsPatchField, surfaceMesh>
         (
-            IOobject
-            (
-                fieldIoObject.name(),
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
+            fieldIoObject,
             mesh_,
             procFields[0].dimensions(),
             internalField,
@@ -520,6 +538,54 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField
 }
 
 
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
+Foam::fvFieldReconstructor::reconstructFvSurfaceField
+(
+    const IOobject& fieldIoObject
+) const
+{
+    // Read the field for all the processors
+    PtrList<GeometricField<Type, fvsPatchField, surfaceMesh> > procFields
+    (
+        procMeshes_.size()
+    );
+
+    forAll(procMeshes_, procI)
+    {
+        procFields.set
+        (
+            procI,
+            new GeometricField<Type, fvsPatchField, surfaceMesh>
+            (
+                IOobject
+                (
+                    fieldIoObject.name(),
+                    procMeshes_[procI].time().timeName(),
+                    procMeshes_[procI],
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                ),
+                procMeshes_[procI]
+            )
+        );
+    }
+
+    return reconstructFvSurfaceField
+    (
+        IOobject
+        (
+            fieldIoObject.name(),
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        procFields
+    );
+}
+
+
 template<class Type>
 void Foam::fvFieldReconstructor::reconstructFvVolumeInternalFields
 (

src/postProcessing/functionObjects/forces/forces/forces.C

@@ -190,6 +190,7 @@ Foam::forces::forces
     rhoRef_(VGREAT),
     pRef_(0),
     coordSys_(),
+    localSystem_(false),
     forcesFilePtr_(NULL)
 {
     // Check if the available mesh is an fvMesh otherise deactivate
@@ -239,6 +240,7 @@ Foam::forces::forces
     rhoRef_(rhoInf),
     pRef_(pRef),
     coordSys_(coordSys),
+    localSystem_(false),
     forcesFilePtr_(NULL)
 {}
 
@@ -330,6 +332,7 @@ void Foam::forces::read(const dictionary& dict)
         if (!dict.readIfPresent<point>("CofR", coordSys_.origin()))
         {
             coordSys_ = coordinateSystem(dict, obr_);
+            localSystem_ = true;
         }
     }
 }
@@ -382,10 +385,17 @@ void Foam::forces::writeFileHeader()
     {
         forcesFilePtr_()
             << "# Time" << tab
-            << "forces(pressure, viscous) moment(pressure, viscous)"
-            << tab
-            << "local forces(pressure, viscous) local moment(pressure, viscous)"
-            << endl;
+            << "forces(pressure, viscous) moment(pressure, viscous)";
+
+        if (localSystem_)
+        {
+            forcesFilePtr_()
+                << tab
+                << "local forces(pressure, viscous) "
+                << "local moment(pressure, viscous)";
+        }
+
+        forcesFilePtr_()<< endl;
     }
 }
 
@@ -413,33 +423,49 @@ void Foam::forces::write()
 
         if (Pstream::master())
         {
-            forcesMoments fmLocal;
-
-            fmLocal.first().first() =
-                coordSys_.localVector(fm.first().first());
-
-            fmLocal.first().second() =
-                coordSys_.localVector(fm.first().second());
-
-            fmLocal.second().first() =
-                coordSys_.localVector(fm.second().first());
-
-            fmLocal.second().second() =
-                coordSys_.localVector(fm.second().second());
-
-            forcesFilePtr_() << obr_.time().value()
-                << tab << fm
-                << tab << fmLocal << endl;
-
             if (log_)
             {
                 Info<< "forces output:" << nl
                     << "    forces(pressure, viscous)" << fm.first() << nl
-                    << "    moment(pressure, viscous)" << fm.second() << nl
-                    << "  local:" << nl
-                    << "    forces(pressure, viscous)" << fmLocal.first() << nl
-                    << "    moment(pressure, viscous)" << fmLocal.second() << nl
-                    << endl;
+                    << "    moment(pressure, viscous)" << fm.second() << nl;
+            }
+
+            forcesFilePtr_() << obr_.time().value() << tab << fm;
+
+            if (localSystem_)
+            {
+                forcesMoments fmLocal;
+
+                fmLocal.first().first() =
+                    coordSys_.localVector(fm.first().first());
+
+                fmLocal.first().second() =
+                    coordSys_.localVector(fm.first().second());
+
+                fmLocal.second().first() =
+                    coordSys_.localVector(fm.second().first());
+
+                fmLocal.second().second() =
+                    coordSys_.localVector(fm.second().second());
+
+                forcesFilePtr_() << tab << fmLocal;
+
+                if (log_)
+                {
+                    Info<< "  local:" << nl
+                        << "    forces(pressure, viscous)" << fmLocal.first()
+                        << nl
+                        << "    moment(pressure, viscous)" << fmLocal.second()
+                        << nl;
+                }
+            }
+
+
+            forcesFilePtr_() << endl;
+
+            if (log_)
+            {
+                Info<< endl;
             }
         }
     }

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

@@ -165,6 +165,9 @@ protected:
             //- Coordinate system used when evaluting forces/moments
             coordinateSystem coordSys_;
 
+            //- Flag to indicate whether we are using a local co-ordinate sys
+            bool localSystem_;
+
 
         //- Forces/moment file ptr
         autoPtr<OFstream> forcesFilePtr_;

src/postProcessing/postCalc/postCalc.C

@@ -57,11 +57,17 @@ namespace Foam
 int main(int argc, char *argv[])
 {
     Foam::timeSelector::addOptions();
+#   include "addRegionOption.H"
     Foam::argList::addBoolOption
     (
         "noWrite",
         "suppress writing results"
     );
+    Foam::argList::addBoolOption
+    (
+        "noFlow",
+        "suppress creating flow models (execFlowFunctionObjects only)"
+    );
     Foam::argList::addOption
     (
         "dict",
@@ -72,7 +78,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     Foam::instantList timeDirs = Foam::timeSelector::select0(runTime, args);
-    #include "createMesh.H"
+    #include "createNamedMesh.H"
 
     forAll(timeDirs, timeI)
     {

src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisTemperatureCoupled/filmPyrolysisTemperatureCoupledFvPatchScalarField.C

@@ -151,13 +151,11 @@ void Foam::filmPyrolysisTemperatureCoupledFvPatchScalarField::updateCoeffs()
 
     const label filmPatchI = filmModel.regionPatchID(patchI);
 
-    const mappedPatchBase& filmMap = filmModel.mappedPatches()[filmPatchI];
-
     scalarField deltaFilm = filmModel.delta().boundaryField()[filmPatchI];
-    filmMap.distribute(deltaFilm);
+    filmModel.toPrimary(filmPatchI, deltaFilm);
 
     scalarField TFilm = filmModel.Ts().boundaryField()[filmPatchI];
-    filmMap.distribute(TFilm);
+    filmModel.toPrimary(filmPatchI, TFilm);
 
 
     // Retrieve pyrolysis model
@@ -166,10 +164,8 @@ void Foam::filmPyrolysisTemperatureCoupledFvPatchScalarField::updateCoeffs()
 
     const label pyrPatchI = pyrModel.regionPatchID(patchI);
 
-    const mappedPatchBase& pyrMap = pyrModel.mappedPatches()[pyrPatchI];
-
     scalarField TPyr = pyrModel.T().boundaryField()[pyrPatchI];
-    pyrMap.distribute(TPyr);
+    pyrModel.toPrimary(pyrPatchI, TPyr);
 
 
     forAll(deltaFilm, i)

src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.C

@@ -154,13 +154,11 @@ void Foam::filmPyrolysisVelocityCoupledFvPatchVectorField::updateCoeffs()
 
     const label filmPatchI = filmModel.regionPatchID(patchI);
 
-    const mappedPatchBase& filmMap = filmModel.mappedPatches()[filmPatchI];
-
     scalarField deltaFilm = filmModel.delta().boundaryField()[filmPatchI];
-    filmMap.distribute(deltaFilm);
+    filmModel.toPrimary(filmPatchI, deltaFilm);
 
     vectorField UFilm = filmModel.Us().boundaryField()[filmPatchI];
-    filmMap.distribute(UFilm);
+    filmModel.toPrimary(filmPatchI, UFilm);
 
 
     // Retrieve pyrolysis model
@@ -172,10 +170,8 @@ void Foam::filmPyrolysisVelocityCoupledFvPatchVectorField::updateCoeffs()
 
     const label pyrPatchI = pyrModel.regionPatchID(patchI);
 
-    const mappedPatchBase& pyrMap = pyrModel.mappedPatches()[pyrPatchI];
-
     scalarField phiPyr = pyrModel.phiGas().boundaryField()[pyrPatchI];
-    pyrMap.distribute(phiPyr);
+    pyrModel.toPrimary(pyrPatchI, phiPyr);
 
 
     const surfaceScalarField& phi =

src/regionModels/regionModel/regionModel/regionModel.C

@@ -93,8 +93,6 @@ void Foam::regionModels::regionModel::initialise()
     DynamicList<label> primaryPatchIDs;
     DynamicList<label> intCoupledPatchIDs;
     const polyBoundaryMesh& rbm = regionMesh().boundaryMesh();
-    const polyBoundaryMesh& pbm = primaryMesh().boundaryMesh();
-    mappedPatches_.setSize(rbm.size());
 
     forAll(rbm, patchI)
     {
@@ -116,19 +114,6 @@ void Foam::regionModels::regionModel::initialise()
 
             const label primaryPatchI = mapPatch.samplePolyPatch().index();
             primaryPatchIDs.append(primaryPatchI);
-
-            mappedPatches_.set
-            (
-                patchI,
-                new mappedPatchBase
-                (
-                    pbm[primaryPatchI],
-                    regionMesh().name(),
-                    mapPatch.mode(),
-                    regionPatch.name(),
-                    vector::zero
-                )
-            );
         }
     }
 
@@ -212,8 +197,7 @@ Foam::regionModels::regionModel::regionModel(const fvMesh& mesh)
     regionMeshPtr_(NULL),
     coeffs_(dictionary::null),
     primaryPatchIDs_(),
-    intCoupledPatchIDs_(),
-    mappedPatches_()
+    intCoupledPatchIDs_()
 {}
 
 
@@ -244,8 +228,7 @@ Foam::regionModels::regionModel::regionModel
     regionMeshPtr_(NULL),
     coeffs_(subOrEmptyDict(modelName + "Coeffs")),
     primaryPatchIDs_(),
-    intCoupledPatchIDs_(),
-    mappedPatches_()
+    intCoupledPatchIDs_()
 {
     if (active_)
     {
@@ -290,8 +273,7 @@ Foam::regionModels::regionModel::regionModel
     regionMeshPtr_(NULL),
     coeffs_(dict.subOrEmptyDict(modelName + "Coeffs")),
     primaryPatchIDs_(),
-    intCoupledPatchIDs_(),
-    mappedPatches_()
+    intCoupledPatchIDs_()
 {
     if (active_)
     {

src/regionModels/regionModel/regionModel/regionModel.H

@@ -117,9 +117,6 @@ protected:
             //- List of patch IDs internally coupled with the primary region
             labelList intCoupledPatchIDs_;
 
-            //- List of patch map info
-            PtrList<mappedPatchBase> mappedPatches_;
-
 
     // Protected member functions
 
@@ -212,13 +209,29 @@ public:
                 //  primary region
                 inline const labelList& intCoupledPatchIDs() const;
 
-                //- Return the list of patch map info
-                inline const PtrList<mappedPatchBase>& mappedPatches() const;
-
                 //- Return region ID corresponding to primaryPatchID
                 inline label regionPatchID(const label primaryPatchID) const;
 
 
+        // Helper
+
+            //- Convert a local region field to the primary region
+            template<class Type>
+            void toPrimary
+            (
+                const label regionPatchI,
+                List<Type>& regionField
+            ) const;
+
+            //- Convert a primary region field to the local region
+            template<class Type>
+            void toRegion
+            (
+                const label regionPatchI,
+                List<Type>& primaryFieldField
+            ) const;
+
+
         // Evolution
 
             //- Pre-evolve region
@@ -249,6 +262,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+    #include "regionModelTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/regionModels/regionModel/regionModel/regionModelI.H

@@ -132,13 +132,6 @@ Foam::regionModels::regionModel::intCoupledPatchIDs() const
 }
 
 
-inline const Foam::PtrList<Foam::mappedPatchBase>&
-Foam::regionModels::regionModel::mappedPatches() const
-{
-    return mappedPatches_;
-}
-
-
 inline Foam::label Foam::regionModels::regionModel::regionPatchID
 (
     const label primaryPatchID

src/regionModels/regionModel/regionModel/regionModelTemplates.C

@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+void Foam::regionModels::regionModel::toPrimary
+(
+    const label regionPatchI,
+    List<Type>& regionField
+) const
+{
+    forAll(intCoupledPatchIDs_, i)
+    {
+        if (intCoupledPatchIDs_[i] == regionPatchI)
+        {
+            const mappedPatchBase& mpb =
+                refCast<const mappedPatchBase>
+                (
+                    regionMesh().boundaryMesh()[regionPatchI]
+                );
+            mpb.reverseDistribute(regionField);
+            return;
+        }
+    }
+
+    FatalErrorIn("const void toPrimary(const label, List<Type>&) const")
+        << "Region patch ID " << regionPatchI << " not found in region mesh"
+        << abort(FatalError);
+}
+
+
+template<class Type>
+void Foam::regionModels::regionModel::toRegion
+(
+    const label regionPatchI,
+    List<Type>& primaryField
+) const
+{
+    forAll(intCoupledPatchIDs_, i)
+    {
+        if (intCoupledPatchIDs_[i] == regionPatchI)
+        {
+            const mappedPatchBase& mpb =
+                refCast<const mappedPatchBase>
+                (
+                    regionMesh().boundaryMesh()[regionPatchI]
+                );
+            mpb.distribute(primaryField);
+            return;
+        }
+    }
+
+    FatalErrorIn("const void toRegion(const label, List<Type>&) const")
+        << "Region patch ID " << regionPatchI << " not found in region mesh"
+        << abort(FatalError);
+}
+
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/alphatFilmWallFunction/alphatFilmWallFunctionFvPatchScalarField.C

@@ -158,11 +158,9 @@ void alphatFilmWallFunctionFvPatchScalarField::updateCoeffs()
 
     const label filmPatchI = filmModel.regionPatchID(patchI);
 
-    const mappedPatchBase& filmMap = filmModel.mappedPatches()[filmPatchI];
-
     tmp<volScalarField> mDotFilm(filmModel.primaryMassTrans());
     scalarField mDotFilmp = mDotFilm().boundaryField()[filmPatchI];
-    filmMap.distribute(mDotFilmp);
+    filmModel.toPrimary(filmPatchI, mDotFilmp);
 
     // Retrieve RAS turbulence model
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");

src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/mutkFilmWallFunction/mutkFilmWallFunctionFvPatchScalarField.C

@@ -70,11 +70,9 @@ tmp<scalarField> mutkFilmWallFunctionFvPatchScalarField::calcUTau
 
     const label filmPatchI = filmModel.regionPatchID(patchI);
 
-    const mappedPatchBase& filmMap = filmModel.mappedPatches()[filmPatchI];
-
     tmp<volScalarField> mDotFilm(filmModel.primaryMassTrans());
     scalarField mDotFilmp = mDotFilm().boundaryField()[filmPatchI];
-    filmMap.distribute(mDotFilmp);
+    filmModel.toPrimary(filmPatchI, mDotFilmp);
 
 
     // Retrieve RAS turbulence model

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C

@@ -827,7 +827,7 @@ void kinematicSingleLayer::preEvolveRegion()
 //    availableMass_ = mass();
     availableMass_ = netMass();
     cloudMassTrans_ == dimensionedScalar("zero", dimMass, 0.0);
-    cloudDiameterTrans_ == dimensionedScalar("zero", dimLength, -1.0);
+    cloudDiameterTrans_ == dimensionedScalar("zero", dimLength, 0.0);
 }
 
 

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerTemplates.C

@@ -1,77 +1,77 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "kinematicSingleLayer.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace regionModels
-{
-namespace surfaceFilmModels
-{
-
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class Type>
-void kinematicSingleLayer::constrainFilmField
-(
-    Type& field,
-    const typename Type::cmptType& value
-)
-{
-    forAll(intCoupledPatchIDs_, i)
-    {
-        label patchI = intCoupledPatchIDs_[i];
-        field.boundaryField()[patchI] = value;
-        if (debug)
-        {
-            Info<< "Constraining " << field.name()
-                << " boundary " << field.boundaryField()[patchI].patch().name()
-                << " to " << value << endl;
-        }
-    }
-    forAll(passivePatchIDs_, i)
-    {
-        label patchI = passivePatchIDs_[i];
-        field.boundaryField()[patchI] = value;
-        if (debug)
-        {
-            Info<< "Constraining " << field.name()
-                << " boundary " << field.boundaryField()[patchI].patch().name()
-                << " to " << value << endl;
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // end namespace Foam
-} // end namespace regionModels
-} // end namespace surfaceFilmModels
-
-// ************************************************************************* //
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "kinematicSingleLayer.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void kinematicSingleLayer::constrainFilmField
+(
+    Type& field,
+    const typename Type::cmptType& value
+)
+{
+    forAll(intCoupledPatchIDs_, i)
+    {
+        label patchI = intCoupledPatchIDs_[i];
+        field.boundaryField()[patchI] = value;
+        if (debug)
+        {
+            Info<< "Constraining " << field.name()
+                << " boundary " << field.boundaryField()[patchI].patch().name()
+                << " to " << value << endl;
+        }
+    }
+    forAll(passivePatchIDs_, i)
+    {
+        label patchI = passivePatchIDs_[i];
+        field.boundaryField()[patchI] = value;
+        if (debug)
+        {
+            Info<< "Constraining " << field.name()
+                << " boundary " << field.boundaryField()[patchI].patch().name()
+                << " to " << value << endl;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // end namespace Foam
+} // end namespace regionModels
+} // end namespace surfaceFilmModels
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/drippingInjection/drippingInjection.C

@@ -139,7 +139,7 @@ void drippingInjection::correct
         {
             // Mass below minimum threshold - cannot be injected
             massToInject[cellI] = 0.0;
-            diameterToInject[cellI] = -1.0;
+            diameterToInject[cellI] = 0.0;
         }
     }
 }

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C

@@ -705,12 +705,11 @@ tmp<DimensionedField<scalar, volMesh> > thermoSingleLayer::Srho() const
     forAll(intCoupledPatchIDs(), i)
     {
         const label filmPatchI = intCoupledPatchIDs()[i];
-        const mappedPatchBase& filmMap = mappedPatches_[filmPatchI];
 
         scalarField patchMass =
             primaryMassPCTrans_.boundaryField()[filmPatchI];
 
-        filmMap.distribute(patchMass);
+        toPrimary(filmPatchI, patchMass);
 
         const label primaryPatchI = primaryPatchIDs()[i];
         const unallocLabelList& cells =
@@ -761,12 +760,11 @@ tmp<DimensionedField<scalar, volMesh> > thermoSingleLayer::Srho
         forAll(intCoupledPatchIDs_, i)
         {
             const label filmPatchI = intCoupledPatchIDs_[i];
-            const mappedPatchBase& filmMap = mappedPatches_[filmPatchI];
 
             scalarField patchMass =
                 primaryMassPCTrans_.boundaryField()[filmPatchI];
 
-            filmMap.distribute(patchMass);
+            toPrimary(filmPatchI, patchMass);
 
             const label primaryPatchI = primaryPatchIDs()[i];
             const unallocLabelList& cells =
@@ -812,12 +810,11 @@ tmp<DimensionedField<scalar, volMesh> > thermoSingleLayer::Sh() const
     forAll(intCoupledPatchIDs_, i)
     {
         const label filmPatchI = intCoupledPatchIDs_[i];
-        const mappedPatchBase& filmMap = mappedPatches_[filmPatchI];
 
         scalarField patchEnergy =
             primaryEnergyPCTrans_.boundaryField()[filmPatchI];
 
-        filmMap.distribute(patchEnergy);
+        toPrimary(filmPatchI, patchEnergy);
 
         const label primaryPatchI = primaryPatchIDs()[i];
         const unallocLabelList& cells =

src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C

@@ -287,16 +287,16 @@ Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface
 
         if (shared[0] != -1)
         {
-            //vector n0 = tri0.normal(localPoints);
-            //n0 /= mag(n0);
-            //vector n1 = tri1.normal(localPoints);
-            //n1 /= mag(n1);
-            //
-            //if ((n0 & n1) < 0)
-            //{
-            //    // Too big an angle. Do not simplify.
-            //}
-            //else
+            vector n0 = tri0.normal(localPoints);
+            n0 /= mag(n0);
+            vector n1 = tri1.normal(localPoints);
+            n1 /= mag(n1);
+
+            if ((n0 & n1) < 0)
+            {
+                // Too big an angle. Do not simplify.
+            }
+            else
             {
                 info.setPoint
                 (
@@ -332,18 +332,18 @@ Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface
         if (nZones == 1)
         {
             // Check that all normals make a decent angle
-            //scalar minCos = GREAT;
-            //const vector& n0 = surf.faceNormals()[0];
-            //for (label i = 1; i < surf.size(); i++)
-            //{
-            //    scalar cosAngle = (n0 & surf.faceNormals()[i]);
-            //    if (cosAngle < minCos)
-            //    {
-            //        minCos = cosAngle;
-            //    }
-            //}
-            //
-            //if (minCos > 0)
+            scalar minCos = GREAT;
+            const vector& n0 = surf.faceNormals()[0];
+            for (label i = 1; i < surf.size(); i++)
+            {
+                scalar cosAngle = (n0 & surf.faceNormals()[i]);
+                if (cosAngle < minCos)
+                {
+                    minCos = cosAngle;
+                }
+            }
+
+            if (minCos > 0)
             {
                 info.setPoint(calcCentre(surf));
                 info.setHit();
@@ -443,7 +443,7 @@ void Foam::isoSurfaceCell::calcSnappedCc
                     label faceI = cFaces[cFaceI];
                     const face& f = mesh_.faces()[faceI];
 
-                    // Do a tetrahedrisation. Each face to cc becomes pyr.
+                    // Do a tetrahedralisation. Each face to cc becomes pyr.
                     // Each pyr gets split into tets by diagonalisation
                     // of face.
 
@@ -467,14 +467,18 @@ void Foam::isoSurfaceCell::calcSnappedCc
                          && (s[2] >= 0.0 && s[2] <= 0.5)
                         )
                         {
-                            if (mesh_.faceOwner()[faceI] == cellI)
+                            if
+                            (
+                                (mesh_.faceOwner()[faceI] == cellI)
+                             == (cVal >= pVals[tri[0]])
+                            )
                             {
                                 localTris.append
                                 (
                                     labelledTri
                                     (
-                                        pointToLocal[tri[0]],
                                         pointToLocal[tri[1]],
+                                        pointToLocal[tri[0]],
                                         pointToLocal[tri[2]],
                                         0
                                     )
@@ -486,8 +490,8 @@ void Foam::isoSurfaceCell::calcSnappedCc
                                 (
                                     labelledTri
                                     (
-                                        pointToLocal[tri[1]],
                                         pointToLocal[tri[0]],
+                                        pointToLocal[tri[1]],
                                         pointToLocal[tri[2]],
                                         0
                                     )
@@ -575,7 +579,11 @@ void Foam::isoSurfaceCell::genPointTris
             point p1 = (1.0-s[1])*pts[pointI] + s[1]*pts[c];
             point p2 = (1.0-s[2])*pts[pointI] + s[2]*cc[cellI];
 
-            if (mesh_.faceOwner()[faceI] == cellI)
+            if
+            (
+                (mesh_.faceOwner()[faceI] == cellI)
+             == (pointValues[pointI] > cellValues[cellI])
+            )
             {
                 localTriPoints.append(p0);
                 localTriPoints.append(p1);
@@ -754,7 +762,7 @@ void Foam::isoSurfaceCell::calcSnappedPoint
         }
 
 
-        // Loop over all pointCells (by using pointFaces)
+        // Do a pointCells walk (by using pointFaces)
 
         localPointCells.clear();
         localTriPoints.clear();
@@ -766,6 +774,8 @@ void Foam::isoSurfaceCell::calcSnappedPoint
 
             if (isTet.get(own) == 1)
             {
+                // Since tets have no cell centre to include make sure
+                // we only generate a triangle once per point.
                 if (localPointCells.insert(own))
                 {
                     genPointTris(pVals, pointI, faceI, own, localTriPoints);
@@ -848,20 +858,19 @@ void Foam::isoSurfaceCell::calcSnappedPoint
 
             if (nZones == 1)
             {
-                //// Check that all normals make a decent angle
-                //scalar minCos = GREAT;
-                //const vector& n0 = surf.faceNormals()[0];
-                //for (label i = 1; i < surf.size(); i++)
-                //{
-                //    const vector& n = surf.faceNormals()[i];
-                //    scalar cosAngle = (n0 & n);
-                //    if (cosAngle < minCos)
-                //    {
-                //        minCos = cosAngle;
-                //    }
-                //}
-                //
-                //if (minCos > 0)
+                // Check that all normals make a decent angle
+                scalar minCos = GREAT;
+                const vector& n0 = surf.faceNormals()[0];
+                for (label i = 1; i < surf.size(); i++)
+                {
+                    const vector& n = surf.faceNormals()[i];
+                    scalar cosAngle = (n0 & n);
+                    if (cosAngle < minCos)
+                    {
+                        minCos = cosAngle;
+                    }
+                }
+                if (minCos > 0)
                 {
                     collapsedPoint[pointI] = calcCentre(surf);
                 }
@@ -882,7 +891,9 @@ void Foam::isoSurfaceCell::calcSnappedPoint
 
     forAll(collapsedPoint, pointI)
     {
-        if (collapsedPoint[pointI] != greatPoint)
+        // Cannot do == comparison since might be transformed so have
+        // truncation errors.
+        if (magSqr(collapsedPoint[pointI]) < 0.5*magSqr(greatPoint))
         {
             snappedPoint[pointI] = snappedPoints.size();
             snappedPoints.append(collapsedPoint[pointI]);
@@ -1208,142 +1219,142 @@ void Foam::isoSurfaceCell::calcAddressing
 }
 
 
-void Foam::isoSurfaceCell::walkOrientation
-(
-    const triSurface& surf,
-    const List<FixedList<label, 3> >& faceEdges,
-    const labelList& edgeFace0,
-    const labelList& edgeFace1,
-    const label seedTriI,
-    labelList& flipState
-)
-{
-    // Do walk for consistent orientation.
-    DynamicList<label> changedFaces(surf.size());
-
-    changedFaces.append(seedTriI);
-
-    while (changedFaces.size())
-    {
-        DynamicList<label> newChangedFaces(changedFaces.size());
-
-        forAll(changedFaces, i)
-        {
-            label triI = changedFaces[i];
-            const labelledTri& tri = surf[triI];
-            const FixedList<label, 3>& fEdges = faceEdges[triI];
-
-            forAll(fEdges, fp)
-            {
-                label edgeI = fEdges[fp];
-
-                // my points:
-                label p0 = tri[fp];
-                label p1 = tri[tri.fcIndex(fp)];
-
-                label nbrI =
-                (
-                    edgeFace0[edgeI] != triI
-                  ? edgeFace0[edgeI]
-                  : edgeFace1[edgeI]
-                );
-
-                if (nbrI != -1 && flipState[nbrI] == -1)
-                {
-                    const labelledTri& nbrTri = surf[nbrI];
-
-                    // nbr points
-                    label nbrFp = findIndex(nbrTri, p0);
-                    label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
-
-                    bool sameOrientation = (p1 == nbrP1);
-
-                    if (flipState[triI] == 0)
-                    {
-                        flipState[nbrI] = (sameOrientation ? 0 : 1);
-                    }
-                    else
-                    {
-                        flipState[nbrI] = (sameOrientation ? 1 : 0);
-                    }
-                    newChangedFaces.append(nbrI);
-                }
-            }
-        }
-
-        changedFaces.transfer(newChangedFaces);
-    }
-}
-
-
-void Foam::isoSurfaceCell::orientSurface
-(
-    triSurface& surf,
-    const List<FixedList<label, 3> >& faceEdges,
-    const labelList& edgeFace0,
-    const labelList& edgeFace1,
-    const Map<labelList>& edgeFacesRest
-)
-{
-    // -1 : unvisited
-    //  0 : leave as is
-    //  1 : flip
-    labelList flipState(surf.size(), -1);
-
-    label seedTriI = 0;
-
-    while (true)
-    {
-        // Find first unvisited triangle
-        for
-        (
-            ;
-            seedTriI < surf.size() && flipState[seedTriI] != -1;
-            seedTriI++
-        )
-        {}
-
-        if (seedTriI == surf.size())
-        {
-            break;
-        }
-
-        // Note: Determine orientation of seedTriI?
-        // for now assume it is ok
-        flipState[seedTriI] = 0;
-
-        walkOrientation
-        (
-            surf,
-            faceEdges,
-            edgeFace0,
-            edgeFace1,
-            seedTriI,
-            flipState
-        );
-    }
-
-    // Do actual flipping
-    surf.clearOut();
-    forAll(surf, triI)
-    {
-        if (flipState[triI] == 1)
-        {
-            labelledTri tri(surf[triI]);
-
-            surf[triI][0] = tri[0];
-            surf[triI][1] = tri[2];
-            surf[triI][2] = tri[1];
-        }
-        else if (flipState[triI] == -1)
-        {
-            FatalErrorIn
-            (
-                "isoSurfaceCell::orientSurface(triSurface&, const label)"
-            )   << "problem" << abort(FatalError);
-        }
-    }
-}
+//void Foam::isoSurfaceCell::walkOrientation
+//(
+//    const triSurface& surf,
+//    const List<FixedList<label, 3> >& faceEdges,
+//    const labelList& edgeFace0,
+//    const labelList& edgeFace1,
+//    const label seedTriI,
+//    labelList& flipState
+//)
+//{
+//    // Do walk for consistent orientation.
+//    DynamicList<label> changedFaces(surf.size());
+//
+//    changedFaces.append(seedTriI);
+//
+//    while (changedFaces.size())
+//    {
+//        DynamicList<label> newChangedFaces(changedFaces.size());
+//
+//        forAll(changedFaces, i)
+//        {
+//            label triI = changedFaces[i];
+//            const labelledTri& tri = surf[triI];
+//            const FixedList<label, 3>& fEdges = faceEdges[triI];
+//
+//            forAll(fEdges, fp)
+//            {
+//                label edgeI = fEdges[fp];
+//
+//                // my points:
+//                label p0 = tri[fp];
+//                label p1 = tri[tri.fcIndex(fp)];
+//
+//                label nbrI =
+//                (
+//                    edgeFace0[edgeI] != triI
+//                  ? edgeFace0[edgeI]
+//                  : edgeFace1[edgeI]
+//                );
+//
+//                if (nbrI != -1 && flipState[nbrI] == -1)
+//                {
+//                    const labelledTri& nbrTri = surf[nbrI];
+//
+//                    // nbr points
+//                    label nbrFp = findIndex(nbrTri, p0);
+//                    label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
+//
+//                    bool sameOrientation = (p1 == nbrP1);
+//
+//                    if (flipState[triI] == 0)
+//                    {
+//                        flipState[nbrI] = (sameOrientation ? 0 : 1);
+//                    }
+//                    else
+//                    {
+//                        flipState[nbrI] = (sameOrientation ? 1 : 0);
+//                    }
+//                    newChangedFaces.append(nbrI);
+//                }
+//            }
+//        }
+//
+//        changedFaces.transfer(newChangedFaces);
+//    }
+//}
+//
+//
+//void Foam::isoSurfaceCell::orientSurface
+//(
+//    triSurface& surf,
+//    const List<FixedList<label, 3> >& faceEdges,
+//    const labelList& edgeFace0,
+//    const labelList& edgeFace1,
+//    const Map<labelList>& edgeFacesRest
+//)
+//{
+//    // -1 : unvisited
+//    //  0 : leave as is
+//    //  1 : flip
+//    labelList flipState(surf.size(), -1);
+//
+//    label seedTriI = 0;
+//
+//    while (true)
+//    {
+//        // Find first unvisited triangle
+//        for
+//        (
+//            ;
+//            seedTriI < surf.size() && flipState[seedTriI] != -1;
+//            seedTriI++
+//        )
+//        {}
+//
+//        if (seedTriI == surf.size())
+//        {
+//            break;
+//        }
+//
+//        // Note: Determine orientation of seedTriI?
+//        // for now assume it is ok
+//        flipState[seedTriI] = 0;
+//
+//        walkOrientation
+//        (
+//            surf,
+//            faceEdges,
+//            edgeFace0,
+//            edgeFace1,
+//            seedTriI,
+//            flipState
+//        );
+//    }
+//
+//    // Do actual flipping
+//    surf.clearOut();
+//    forAll(surf, triI)
+//    {
+//        if (flipState[triI] == 1)
+//        {
+//            labelledTri tri(surf[triI]);
+//
+//            surf[triI][0] = tri[0];
+//            surf[triI][1] = tri[2];
+//            surf[triI][2] = tri[1];
+//        }
+//        else if (flipState[triI] == -1)
+//        {
+//            FatalErrorIn
+//            (
+//                "isoSurfaceCell::orientSurface(triSurface&, const label)"
+//            )   << "problem" << abort(FatalError);
+//        }
+//    }
+//}
 
 
 // Checks if triangle is connected through edgeI only.
@@ -1625,7 +1636,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
     (
         stitchTriPoints
         (
-            true,               // check for duplicate tris
+            regularise,         // check for duplicate tris
             triPoints,
             triPointMergeMap_,  // unmerged to merged point
             triMap
@@ -1657,7 +1668,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
     }
 
 
-    if (false)
+    if (regularise)
     {
         List<FixedList<label, 3> > faceEdges;
         labelList edgeFace0, edgeFace1;
@@ -1717,7 +1728,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
             inplaceRenumber(reversePointMap, triPointMergeMap_);
         }
 
-        orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
+        //orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
     }
 }
 

src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H

@@ -214,15 +214,19 @@ class isoSurfaceCell
         void generateTriPoints
         (
             const DynamicList<Type>& snapped,
+            const scalar isoVal0,
             const scalar s0,
             const Type& p0,
             const label p0Index,
+            const scalar isoVal1,
             const scalar s1,
             const Type& p1,
             const label p1Index,
+            const scalar isoVal2,
             const scalar s2,
             const Type& p2,
             const label p2Index,
+            const scalar isoVal3,
             const scalar s3,
             const Type& p3,
             const label p3Index,
@@ -267,26 +271,26 @@ class isoSurfaceCell
             Map<labelList>& edgeFacesRest
         ) const;
 
-        //- Determine orientation
-        static void walkOrientation
-        (
-            const triSurface& surf,
-            const List<FixedList<label, 3> >& faceEdges,
-            const labelList& edgeFace0,
-            const labelList& edgeFace1,
-            const label seedTriI,
-            labelList& flipState
-        );
+        ////- Determine orientation
+        //static void walkOrientation
+        //(
+        //    const triSurface& surf,
+        //    const List<FixedList<label, 3> >& faceEdges,
+        //    const labelList& edgeFace0,
+        //    const labelList& edgeFace1,
+        //    const label seedTriI,
+        //    labelList& flipState
+        //);
 
-        //- Orient surface
-        static void orientSurface
-        (
-            triSurface&,
-            const List<FixedList<label, 3> >& faceEdges,
-            const labelList& edgeFace0,
-            const labelList& edgeFace1,
-            const Map<labelList>& edgeFacesRest
-        );
+        ////- Orient surface
+        //static void orientSurface
+        //(
+        //    triSurface&,
+        //    const List<FixedList<label, 3> >& faceEdges,
+        //    const labelList& edgeFace0,
+        //    const labelList& edgeFace1,
+        //    const Map<labelList>& edgeFacesRest
+        //);
 
         //- Is triangle (given by 3 edges) not fully connected?
         static bool danglingTriangle

src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C

@@ -76,23 +76,27 @@ void Foam::isoSurfaceCell::generateTriPoints
 (
     const DynamicList<Type>& snapped,
 
+    const scalar isoVal0,
     const scalar s0,
     const Type& p0,
     const label p0Index,
 
+    const scalar isoVal1,
     const scalar s1,
     const Type& p1,
     const label p1Index,
 
+    const scalar isoVal2,
     const scalar s2,
     const Type& p2,
     const label p2Index,
 
+    const scalar isoVal3,
     const scalar s3,
     const Type& p3,
     const label p3Index,
 
-    DynamicList<Type>& points
+    DynamicList<Type>& pts
 ) const
 {
     int triIndex = 0;
@@ -122,77 +126,160 @@ void Foam::isoSurfaceCell::generateTriPoints
 
         case 0x0E:
         case 0x01:
-            points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
-            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+        {
+            // 0 is common point. Orient such that normal points in positive
+            // gradient direction
+            if (isoVal0 >= isoVal1)
+            {
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            }
+            else
+            {
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            }
+        }
         break;
 
         case 0x0D:
         case 0x02:
-            points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+        {
+            // 1 is common point
+            if (isoVal1 >= isoVal0)
+            {
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            }
+            else
+            {
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            }
+        }
         break;
 
         case 0x0C:
         case 0x03:
         {
-            Type tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
-            Type tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
+            Type s02 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
+            Type s13 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
 
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-            points.append(tp1);
-            points.append(tp2);
-            points.append(tp2);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-            points.append(tp1);
+            if (isoVal0 >= isoVal3)
+            {
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+                pts.append(s02);
+                pts.append(s13);
+                pts.append(s13);
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+                pts.append(s02);
+            }
+            else
+            {
+                pts.append(s02);
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+                pts.append(s13);
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+                pts.append(s13);
+                pts.append(s02);
+            }
         }
         break;
 
         case 0x0B:
         case 0x04:
         {
-            points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
-            points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+            // 2 is common point
+            if (isoVal2 >= isoVal0)
+            {
+                pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
+                pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
+                pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+            }
+            else
+            {
+                pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
+                pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
+                pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+            }
         }
         break;
 
         case 0x0A:
         case 0x05:
         {
-            Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+            Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
 
-            points.append(tp0);
-            points.append(tp1);
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-            points.append(tp0);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-            points.append(tp1);
+            if (isoVal3 >= isoVal0)
+            {
+                pts.append(s01);
+                pts.append(s23);
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+                pts.append(s01);
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+                pts.append(s23);
+            }
+            else
+            {
+                pts.append(s23);
+                pts.append(s01);
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+                pts.append(s01);
+                pts.append(s23);
+            }
         }
         break;
 
         case 0x09:
         case 0x06:
         {
-            Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+            Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
 
-            points.append(tp0);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-            points.append(tp1);
-            points.append(tp0);
-            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-            points.append(tp1);
+            if (isoVal3 >= isoVal1)
+            {
+                pts.append(s01);
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+                pts.append(s23);
+                pts.append(s01);
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+                pts.append(s23);
+            }
+            else
+            {
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+                pts.append(s01);
+                pts.append(s23);
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+                pts.append(s01);
+                pts.append(s23);
+            }
         }
         break;
 
         case 0x07:
         case 0x08:
-            points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
-            points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+        {
+            // 3 is common point
+            if (isoVal3 >= isoVal0)
+            {
+                pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
+                pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
+                pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+            }
+            else
+            {
+                pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
+                pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
+                pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+            }
+        }
         break;
     }
 }
@@ -252,18 +339,23 @@ void Foam::isoSurfaceCell::generateTriPoints
                     generateTriPoints
                     (
                         snappedPoints,
+
+                        pVals_[f0[1]],
                         pVals[f0[1]],
                         pCoords[f0[1]],
                         snappedPoint[f0[1]],
 
+                        pVals_[f0[0]],
                         pVals[f0[0]],
                         pCoords[f0[0]],
                         snappedPoint[f0[0]],
 
+                        pVals_[f0[2]],
                         pVals[f0[2]],
                         pCoords[f0[2]],
                         snappedPoint[f0[2]],
 
+                        pVals_[oppositeI],
                         pVals[oppositeI],
                         pCoords[oppositeI],
                         snappedPoint[oppositeI],
@@ -276,18 +368,23 @@ void Foam::isoSurfaceCell::generateTriPoints
                     generateTriPoints
                     (
                         snappedPoints,
+
+                        pVals_[f0[0]],
                         pVals[f0[0]],
                         pCoords[f0[0]],
                         snappedPoint[f0[0]],
 
+                        pVals_[f0[1]],
                         pVals[f0[1]],
                         pCoords[f0[1]],
                         snappedPoint[f0[1]],
 
+                        pVals_[f0[2]],
                         pVals[f0[2]],
                         pCoords[f0[2]],
                         snappedPoint[f0[2]],
 
+                        pVals_[oppositeI],
                         pVals[oppositeI],
                         pCoords[oppositeI],
                         snappedPoint[oppositeI],
@@ -321,18 +418,22 @@ void Foam::isoSurfaceCell::generateTriPoints
                             (
                                 snappedPoints,
 
+                                pVals_[tri[1]],
                                 pVals[tri[1]],
                                 pCoords[tri[1]],
                                 snappedPoint[tri[1]],
 
+                                pVals_[tri[0]],
                                 pVals[tri[0]],
                                 pCoords[tri[0]],
                                 snappedPoint[tri[0]],
 
+                                pVals_[tri[2]],
                                 pVals[tri[2]],
                                 pCoords[tri[2]],
                                 snappedPoint[tri[2]],
 
+                                cVals_[cellI],
                                 cVals[cellI],
                                 cCoords[cellI],
                                 snappedCc[cellI],
@@ -346,18 +447,22 @@ void Foam::isoSurfaceCell::generateTriPoints
                             (
                                 snappedPoints,
 
+                                pVals_[tri[0]],
                                 pVals[tri[0]],
                                 pCoords[tri[0]],
                                 snappedPoint[tri[0]],
 
+                                pVals_[tri[1]],
                                 pVals[tri[1]],
                                 pCoords[tri[1]],
                                 snappedPoint[tri[1]],
 
+                                pVals_[tri[2]],
                                 pVals[tri[2]],
                                 pCoords[tri[2]],
                                 snappedPoint[tri[2]],
 
+                                cVals_[cellI],
                                 cVals[cellI],
                                 cCoords[cellI],
                                 snappedCc[cellI],

src/sampling/sampledSurface/sampledPatchInternalField/sampledPatchInternalField.C

@@ -58,24 +58,80 @@ Foam::sampledPatchInternalField::sampledPatchInternalField
     sampledPatch(name, mesh, dict),
     mappers_(patchIDs().size())
 {
-    const scalar distance = readScalar(dict.lookup("distance"));
-
-    forAll(patchIDs(), i)
+    mappedPatchBase::offsetMode mode = mappedPatchBase::NORMAL;
+    if (dict.found("offsetMode"))
     {
-        label patchI = patchIDs()[i];
-        mappers_.set
+        mode = mappedPatchBase::offsetModeNames_.read
         (
-            i,
-            new mappedPatchBase
-            (
-                mesh.boundaryMesh()[patchI],
-                mesh.name(),                        // sampleRegion
-                mappedPatchBase::NEARESTCELL, // sampleMode
-                word::null,                         // samplePatch
-                -distance                           // sample inside my domain
-            )
+            dict.lookup("offsetMode")
         );
     }
+
+    switch (mode)
+    {
+        case mappedPatchBase::NORMAL:
+        {
+            const scalar distance = readScalar(dict.lookup("distance"));
+            forAll(patchIDs(), i)
+            {
+                mappers_.set
+                (
+                    i,
+                    new mappedPatchBase
+                    (
+                        mesh.boundaryMesh()[patchIDs()[i]],
+                        mesh.name(),                        // sampleRegion
+                        mappedPatchBase::NEARESTCELL,       // sampleMode
+                        word::null,                         // samplePatch
+                        -distance                  // sample inside my domain
+                    )
+                );
+            }
+        }
+        break;
+
+        case mappedPatchBase::UNIFORM:
+        {
+            const point offset(dict.lookup("offset"));
+            forAll(patchIDs(), i)
+            {
+                mappers_.set
+                (
+                    i,
+                    new mappedPatchBase
+                    (
+                        mesh.boundaryMesh()[patchIDs()[i]],
+                        mesh.name(),                        // sampleRegion
+                        mappedPatchBase::NEARESTCELL,       // sampleMode
+                        word::null,                         // samplePatch
+                        offset                  // sample inside my domain
+                    )
+                );
+            }
+        }
+        break;
+
+        case mappedPatchBase::NONUNIFORM:
+        {
+            const pointField offsets(dict.lookup("offsets"));
+            forAll(patchIDs(), i)
+            {
+                mappers_.set
+                (
+                    i,
+                    new mappedPatchBase
+                    (
+                        mesh.boundaryMesh()[patchIDs()[i]],
+                        mesh.name(),                        // sampleRegion
+                        mappedPatchBase::NEARESTCELL,       // sampleMode
+                        word::null,                         // samplePatch
+                        offsets                  // sample inside my domain
+                    )
+                );
+            }
+        }
+        break;
+    }
 }
 
 

src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/temperatureThermoBaffle1D/temperatureThermoBaffle1DFvPatchScalarField.H

@@ -25,7 +25,7 @@ Class
     Foam::temperatureThermoBaffle1DFvPatchScalarField
 
 Description
-    Bounday which solves the 1D steady state heat transfer equation
+    Boundary which solves the 1D steady state heat transfer equation
     through a baffle.
 
 SourceFiles

src/turbulenceModels/incompressible/LES/derivedFvPatchFields/wallFunctions/nuSgsWallFunctions/nuSgsUSpaldingWallFunction/nuSgsUSpaldingWallFunctionFvPatchScalarField.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "nuSgsUSpaldingWallFunctionFvPatchScalarField.H"
+#include "LESModel.H"
 #include "fvPatchFieldMapper.H"
 #include "volFields.H"
 #include "addToRunTimeSelectionTable.H"
@@ -47,8 +48,6 @@ nuSgsUSpaldingWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(p, iF),
-    UName_("U"),
-    nuName_("nu"),
     kappa_(0.41),
     E_(9.8)
 {}
@@ -64,8 +63,6 @@ nuSgsUSpaldingWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(ptf, p, iF, mapper),
-    UName_(ptf.UName_),
-    nuName_(ptf.nuName_),
     kappa_(ptf.kappa_),
     E_(ptf.E_)
 {}
@@ -80,8 +77,6 @@ nuSgsUSpaldingWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(p, iF, dict),
-    UName_(dict.lookupOrDefault<word>("U", "U")),
-    nuName_(dict.lookupOrDefault<word>("nu", "nu")),
     kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
     E_(dict.lookupOrDefault<scalar>("E", 9.8))
 {}
@@ -94,8 +89,6 @@ nuSgsUSpaldingWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(nwfpsf),
-    UName_(nwfpsf.UName_),
-    nuName_(nwfpsf.nuName_),
     kappa_(nwfpsf.kappa_),
     E_(nwfpsf.E_)
 {}
@@ -109,8 +102,6 @@ nuSgsUSpaldingWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(nwfpsf, iF),
-    UName_(nwfpsf.UName_),
-    nuName_(nwfpsf.nuName_),
     kappa_(nwfpsf.kappa_),
     E_(nwfpsf.E_)
 {}
@@ -123,16 +114,15 @@ void nuSgsUSpaldingWallFunctionFvPatchScalarField::evaluate
     const Pstream::commsTypes
 )
 {
+    const LESModel& lesModel = db().lookupObject<LESModel>("LESProperties");
+    const label patchi = patch().index();
+    const fvPatchVectorField& U = lesModel.U().boundaryField()[patchi];
+    const scalarField nuw = lesModel.nu()().boundaryField()[patchi];
+
     const scalarField& ry = patch().deltaCoeffs();
 
-    const fvPatchVectorField& U =
-        patch().lookupPatchField<volVectorField, vector>(UName_);
-
     const scalarField magUp(mag(U.patchInternalField() - U));
 
-    const scalarField& nuw =
-        patch().lookupPatchField<volScalarField, scalar>(nuName_);
-
     scalarField& nuSgsw = *this;
 
     const scalarField magFaceGradU(mag(U.snGrad()));
@@ -185,8 +175,6 @@ void nuSgsUSpaldingWallFunctionFvPatchScalarField::evaluate
 void nuSgsUSpaldingWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchField<scalar>::write(os);
-    writeEntryIfDifferent<word>(os, "U", "U", UName_);
-    writeEntryIfDifferent<word>(os, "nu", "nu", nuName_);
     os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
     os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
     writeEntry("value", os);
@@ -201,6 +189,7 @@ makePatchTypeField
     nuSgsUSpaldingWallFunctionFvPatchScalarField
 );
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace LESModels

src/turbulenceModels/incompressible/LES/derivedFvPatchFields/wallFunctions/nuSgsWallFunctions/nuSgsUSpaldingWallFunction/nuSgsUSpaldingWallFunctionFvPatchScalarField.H

@@ -58,12 +58,6 @@ class nuSgsUSpaldingWallFunctionFvPatchScalarField
 {
     // Private data
 
-        //- Name of velocity field
-        word UName_;
-
-        //- Name of laminar viscosity field
-        word nuName_;
-
         //- Von Karman constant
         scalar kappa_;