Merge branch 'master' into particleInteractions

applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C

@@ -270,7 +270,7 @@ void PDRkEpsilon::correct()
     }
 
     tmp<volTensorField> tgradU = fvc::grad(U_);
-    volScalarField G = 2*mut_*(tgradU() && dev(symm(tgradU())));
+    volScalarField G("RASModel::G", mut_*(tgradU() && dev(twoSymm(tgradU()))));
     tgradU.clear();
 
     // Update espsilon and G at the wall

applications/solvers/combustion/fireFoam/combustionModels/Make/options

@@ -1,6 +1,4 @@
 EXE_INC = \
-    -I../sensibleEnthalpyCombustionThermophysicalModels/basic/lnInclude \
-    -I../sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \

applications/solvers/electromagnetics/mhdFoam/createFields.H

@@ -61,7 +61,7 @@
         mesh
     );
 
-#   include "createPhi.H"
+    #include "createPhi.H"
 
     Info<< "Reading field pB\n" << endl;
     volScalarField pB
@@ -93,15 +93,15 @@
     );
 
 
-#   include "createPhiB.H"
+    #include "createPhiB.H"
 
-dimensionedScalar DB = 1.0/(mu*sigma);
-DB.name() = "DB";
+    dimensionedScalar DB = 1.0/(mu*sigma);
+    DB.name() = "DB";
 
-dimensionedScalar DBU = 1.0/(2.0*mu*rho);
-DBU.name() = "DBU";
+    dimensionedScalar DBU = 1.0/(2.0*mu*rho);
+    DBU.name() = "DBU";
 
 
-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);

applications/solvers/electromagnetics/mhdFoam/mhdFoam.C

@@ -49,8 +49,6 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "string.H"
-#include "Time.H"
 #include "fvCFD.H"
 #include "OSspecific.H"
 

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C

@@ -38,7 +38,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(p, iF),
+    fixedGradientFvPatchScalarField(p, iF),
     q_(p.size(), 0.0),
     KName_("undefined-K")
 {}
@@ -53,7 +53,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
     q_(ptf.q_, mapper),
     KName_(ptf.KName_)
 {}
@@ -67,7 +67,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const dictionary& dict
 )
 :
-    fixedValueFvPatchScalarField(p, iF, dict),
+    fixedGradientFvPatchScalarField(p, iF, dict),
     q_("q", dict, p.size()),
     KName_(dict.lookup("K"))
 {}
@@ -79,7 +79,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const solidWallHeatFluxTemperatureFvPatchScalarField& tppsf
 )
 :
-    fixedValueFvPatchScalarField(tppsf),
+    fixedGradientFvPatchScalarField(tppsf),
     q_(tppsf.q_),
     KName_(tppsf.KName_)
 {}
@@ -92,7 +92,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(tppsf, iF),
+    fixedGradientFvPatchScalarField(tppsf, iF),
     q_(tppsf.q_),
     KName_(tppsf.KName_)
 {}
@@ -105,7 +105,7 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::autoMap
     const fvPatchFieldMapper& m
 )
 {
-    fixedValueFvPatchScalarField::autoMap(m);
+    fixedGradientFvPatchScalarField::autoMap(m);
     q_.autoMap(m);
 }
 
@@ -116,7 +116,7 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::rmap
     const labelList& addr
 )
 {
-    fixedValueFvPatchScalarField::rmap(ptf, addr);
+    fixedGradientFvPatchScalarField::rmap(ptf, addr);
 
     const solidWallHeatFluxTemperatureFvPatchScalarField& hfptf =
         refCast<const solidWallHeatFluxTemperatureFvPatchScalarField>(ptf);
@@ -132,14 +132,14 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    const scalarField& Kw =
-        patch().lookupPatchField<volScalarField, scalar>(KName_);
+    const scalarField& Kw = patch().lookupPatchField<volScalarField, scalar>
+    (
+        KName_
+    );
 
-    const fvPatchScalarField& Tw = *this;
+    gradient() = q_/Kw;
 
-    operator==(q_/(patch().deltaCoeffs()*Kw) + Tw.patchInternalField());
-
-    fixedValueFvPatchScalarField::updateCoeffs();
+    fixedGradientFvPatchScalarField::updateCoeffs();
 }
 
 
@@ -148,9 +148,10 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::write
     Ostream& os
 ) const
 {
-    fixedValueFvPatchScalarField::write(os);
+    fixedGradientFvPatchScalarField::write(os);
     q_.writeEntry("q", os);
     os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
+    this->writeEntry("value", os);
 }
 
 

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H

@@ -46,7 +46,7 @@ SourceFiles
 #ifndef solidWallHeatFluxTemperatureFvPatchScalarField_H
 #define solidWallHeatFluxTemperatureFvPatchScalarField_H
 
-#include "fixedValueFvPatchFields.H"
+#include "fixedGradientFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,7 +59,7 @@ namespace Foam
 
 class solidWallHeatFluxTemperatureFvPatchScalarField
 :
-    public fixedValueFvPatchScalarField
+    public fixedGradientFvPatchScalarField
 {
     // Private data
 

applications/solvers/incompressible/ajointShapeOptimizationFoam/Make/files

@@ -0,0 +1,5 @@
+adjointOutletPressure/adjointOutletPressureFvPatchScalarField.C
+adjointOutletVelocity/adjointOutletVelocityFvPatchVectorField.C
+adjointShapeOptimizationFoam.C
+
+EXE = $(FOAM_APPBIN)/adjointShapeOptimizationFoam

applications/solvers/incompressible/ajointShapeOptimizationFoam/Make/options

@@ -0,0 +1,11 @@
+EXE_INC = \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lincompressibleRASModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume

applications/solvers/incompressible/ajointShapeOptimizationFoam/adjointContinuityErrs.H

@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Global
+    continuityErrs
+
+Description
+    Calculates and prints the continuity errors.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    scalar sumLocalContErr = runTime.deltaT().value()*
+        mag(fvc::div(phia))().weightedAverage(mesh.V()).value();
+
+    scalar globalContErr = runTime.deltaT().value()*
+        fvc::div(phia)().weightedAverage(mesh.V()).value();
+    cumulativeContErr += globalContErr;
+
+    Info<< "Adjoint continuity errors : sum local = " << sumLocalContErr
+        << ", global = " << globalContErr
+        << ", cumulative = " << cumulativeContErr
+        << endl;
+}
+
+// ************************************************************************* //

applications/solvers/incompressible/ajointShapeOptimizationFoam/adjointOutletPressure/adjointOutletPressureFvPatchScalarField.C

@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "adjointOutletPressureFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchMapper.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::adjointOutletPressureFvPatchScalarField::
+adjointOutletPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF)
+{}
+
+
+Foam::adjointOutletPressureFvPatchScalarField::
+adjointOutletPressureFvPatchScalarField
+(
+    const adjointOutletPressureFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper)
+{}
+
+
+Foam::adjointOutletPressureFvPatchScalarField::
+adjointOutletPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF)
+{
+    fvPatchField<scalar>::operator=
+    (
+        scalarField("value", dict, p.size())
+    );
+}
+
+
+Foam::adjointOutletPressureFvPatchScalarField::
+adjointOutletPressureFvPatchScalarField
+(
+    const adjointOutletPressureFvPatchScalarField& tppsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(tppsf, iF)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::adjointOutletPressureFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const fvsPatchField<scalar>& phip =
+        patch().lookupPatchField<surfaceScalarField, scalar>("phi");
+
+    const fvsPatchField<scalar>& phiap =
+        patch().lookupPatchField<surfaceScalarField, scalar>("phia");
+
+    const fvPatchField<vector>& Up =
+        patch().lookupPatchField<volVectorField, vector>("U");
+
+    const fvPatchField<vector>& Uap =
+        patch().lookupPatchField<volVectorField, vector>("Ua");
+
+    operator==((phiap/patch().magSf() - 1.0)*phip/patch().magSf() + (Up & Uap));
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::adjointOutletPressureFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        adjointOutletPressureFvPatchScalarField
+    );
+}
+
+// ************************************************************************* //

applications/solvers/incompressible/ajointShapeOptimizationFoam/adjointOutletPressure/adjointOutletPressureFvPatchScalarField.H

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    adjointOutletPressureFvPatchScalarField
+
+Description
+
+SourceFiles
+    adjointOutletPressureFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef adjointOutletPressureFvPatchScalarField_H
+#define adjointOutletPressureFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class adjointOutletPressureFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class adjointOutletPressureFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("adjointOutletPressure");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        adjointOutletPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        adjointOutletPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given adjointOutletPressureFvPatchScalarField
+        //  onto a new patch
+        adjointOutletPressureFvPatchScalarField
+        (
+            const adjointOutletPressureFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new adjointOutletPressureFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        adjointOutletPressureFvPatchScalarField
+        (
+            const adjointOutletPressureFvPatchScalarField&,
+            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 adjointOutletPressureFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/incompressible/ajointShapeOptimizationFoam/adjointOutletVelocity/adjointOutletVelocityFvPatchVectorField.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "adjointOutletVelocityFvPatchVectorField.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "surfaceFields.H"
+#include "fvPatchFieldMapper.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::adjointOutletVelocityFvPatchVectorField::
+adjointOutletVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchVectorField(p, iF)
+{}
+
+
+Foam::adjointOutletVelocityFvPatchVectorField::
+adjointOutletVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchVectorField(p, iF)
+{
+    fvPatchVectorField::operator=(vectorField("value", dict, p.size()));
+}
+
+
+Foam::adjointOutletVelocityFvPatchVectorField::
+adjointOutletVelocityFvPatchVectorField
+(
+    const adjointOutletVelocityFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchVectorField(ptf, p, iF, mapper)
+{}
+
+
+Foam::adjointOutletVelocityFvPatchVectorField::
+adjointOutletVelocityFvPatchVectorField
+(
+    const adjointOutletVelocityFvPatchVectorField& pivpvf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchVectorField(pivpvf, iF)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+// Update the coefficients associated with the patch field
+void Foam::adjointOutletVelocityFvPatchVectorField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const fvsPatchField<scalar>& phiap =
+        patch().lookupPatchField<surfaceScalarField, scalar>("phia");
+
+    const fvPatchField<vector>& Up =
+        patch().lookupPatchField<volVectorField, vector>("U");
+
+    scalarField Un = mag(patch().nf() & Up);
+    vectorField UtHat = (Up - patch().nf()*Un)/(Un + SMALL);
+
+    vectorField Uan = patch().nf()*(patch().nf() & patchInternalField());
+
+    vectorField::operator=(phiap*patch().Sf()/sqr(patch().magSf()) + UtHat);
+    //vectorField::operator=(Uan + UtHat);
+
+    fixedValueFvPatchVectorField::updateCoeffs();
+}
+
+
+void Foam::adjointOutletVelocityFvPatchVectorField::write(Ostream& os) const
+{
+    fvPatchVectorField::write(os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchVectorField,
+        adjointOutletVelocityFvPatchVectorField
+    );
+}
+
+
+// ************************************************************************* //

applications/solvers/incompressible/ajointShapeOptimizationFoam/adjointOutletVelocity/adjointOutletVelocityFvPatchVectorField.H

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    adjointOutletVelocityFvPatchVectorField
+
+Description
+
+SourceFiles
+    adjointOutletVelocityFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef adjointOutletVelocityFvPatchVectorField_H
+#define adjointOutletVelocityFvPatchVectorField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class adjointOutletVelocityFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class adjointOutletVelocityFvPatchVectorField
+:
+    public fixedValueFvPatchVectorField
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("adjointOutletVelocity");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        adjointOutletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        adjointOutletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given adjointOutletVelocityFvPatchVectorField
+        //  onto a new patch
+        adjointOutletVelocityFvPatchVectorField
+        (
+            const adjointOutletVelocityFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new adjointOutletVelocityFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        adjointOutletVelocityFvPatchVectorField
+        (
+            const adjointOutletVelocityFvPatchVectorField&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchVectorField> clone(const DimensionedField<vector, volMesh>& iF) const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new adjointOutletVelocityFvPatchVectorField(*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
+
+// ************************************************************************* //

applications/solvers/incompressible/ajointShapeOptimizationFoam/adjointShapeOptimizationFoam.C

@@ -0,0 +1,226 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    ajointShapeOptimizationFoam
+
+Description
+    Steady-state solver for incompressible, turbulent flow of non-Newtonian
+    fluids with optimisation of duct shape by applying "blockage" in regions
+    causing pressure loss as estimated using an adjoint formulation.
+
+    References:
+    @verbatim
+        "Implementation of a continuous adjoint for topology optimization of
+         ducted flows"
+        C. Othmer,
+        E. de Villiers,
+        H.G. Weller
+        AIAA-2007-3947
+        http://pdf.aiaa.org/preview/CDReadyMCFD07_1379/PV2007_3947.pdf
+    @endverbatim
+
+    Note that this solver optimises for total pressure loss whereas the
+    above paper describes the method for optimising power-loss.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "RASModel.H"
+
+template<class Type>
+void zeroCells
+(
+    GeometricField<Type, fvPatchField, volMesh>& vf,
+    const labelList& cells
+)
+{
+    forAll(cells, i)
+    {
+        vf[cells[i]] = pTraits<Type>::zero;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    for (runTime++; !runTime.end(); runTime++)
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "readSIMPLEControls.H"
+
+        p.storePrevIter();
+
+        laminarTransport.lookup("lambda") >> lambda;
+
+        //alpha +=
+        //    mesh.relaxationFactor("alpha")
+        //   *(lambda*max(Ua & U, zeroSensitivity) - alpha);
+        alpha +=
+            mesh.relaxationFactor("alpha")
+           *(min(max(alpha + lambda*(Ua & U), zeroAlpha), alphaMax) - alpha);
+
+        zeroCells(alpha, inletCells);
+        //zeroCells(alpha, outletCells);
+
+        // Pressure-velocity SIMPLE corrector
+        {
+            // Momentum predictor
+
+            tmp<fvVectorMatrix> UEqn
+            (
+                fvm::div(phi, U)
+              + turbulence->divDevReff(U)
+              + fvm::Sp(alpha, U)
+            );
+
+            UEqn().relax();
+
+            solve(UEqn() == -fvc::grad(p));
+
+            p.boundaryField().updateCoeffs();
+            volScalarField rAU = 1.0/UEqn().A();
+            U = rAU*UEqn().H();
+            UEqn.clear();
+            phi = fvc::interpolate(U) & mesh.Sf();
+            adjustPhi(phi, U, p);
+
+            // Non-orthogonal pressure corrector loop
+            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            {
+                fvScalarMatrix pEqn
+                (
+                    fvm::laplacian(rAU, p) == fvc::div(phi)
+                );
+
+                pEqn.setReference(pRefCell, pRefValue);
+                pEqn.solve();
+
+                if (nonOrth == nNonOrthCorr)
+                {
+                    phi -= pEqn.flux();
+                }
+            }
+
+            #include "continuityErrs.H"
+
+            // Explicitly relax pressure for momentum corrector
+            p.relax();
+
+            // Momentum corrector
+            U -= rAU*fvc::grad(p);
+            U.correctBoundaryConditions();
+        }
+
+        pa.storePrevIter();
+
+        // Adjoint Pressure-velocity SIMPLE corrector
+        {
+            // Adjoint Momentum predictor
+
+            volVectorField adjointTransposeConvection = (fvc::grad(Ua) & U);
+            //volVectorField adjointTransposeConvection = fvc::reconstruct
+            //(
+            //    mesh.magSf()*(fvc::snGrad(Ua) & fvc::interpolate(U))
+            //);
+
+            zeroCells(adjointTransposeConvection, inletCells);
+
+            tmp<fvVectorMatrix> UaEqn
+            (
+                fvm::div(-phi, Ua)
+              - adjointTransposeConvection
+              + turbulence->divDevReff(Ua)
+              + fvm::Sp(alpha, Ua)
+            );
+
+            UaEqn().relax();
+
+            solve(UaEqn() == -fvc::grad(pa));
+
+            pa.boundaryField().updateCoeffs();
+            volScalarField rAUa = 1.0/UaEqn().A();
+            Ua = rAUa*UaEqn().H();
+            UaEqn.clear();
+            phia = fvc::interpolate(Ua) & mesh.Sf();
+            adjustPhi(phia, Ua, pa);
+
+            // Non-orthogonal pressure corrector loop
+            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            {
+                fvScalarMatrix paEqn
+                (
+                    fvm::laplacian(rAUa, pa) == fvc::div(phia)
+                );
+
+                paEqn.setReference(paRefCell, paRefValue);
+                paEqn.solve();
+
+                if (nonOrth == nNonOrthCorr)
+                {
+                    phia -= paEqn.flux();
+                }
+            }
+
+            #include "adjointContinuityErrs.H"
+
+            // Explicitly relax pressure for adjoint momentum corrector
+            pa.relax();
+
+            // Adjoint momentum corrector
+            Ua -= rAUa*fvc::grad(pa);
+            Ua.correctBoundaryConditions();
+        }
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = "
+            << runTime.elapsedCpuTime()
+            << " s\n\n" << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/incompressible/ajointShapeOptimizationFoam/createFields.H

@@ -0,0 +1,105 @@
+    Info << "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info << "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "createPhi.H"
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
+
+
+    Info << "Reading field pa\n" << endl;
+    volScalarField pa
+    (
+        IOobject
+        (
+            "pa",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info << "Reading field Ua\n" << endl;
+    volVectorField Ua
+    (
+        IOobject
+        (
+            "Ua",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "createPhia.H"
+
+
+    label paRefCell = 0;
+    scalar paRefValue = 0.0;
+    setRefCell(pa, mesh.solutionDict().subDict("SIMPLE"), paRefCell, paRefValue);
+
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::RASModel> turbulence
+    (
+        incompressible::RASModel::New(U, phi, laminarTransport)
+    );
+
+
+    dimensionedScalar zeroSensitivity("0", dimVelocity*dimVelocity, 0.0);
+    dimensionedScalar zeroAlpha("0", dimless/dimTime, 0.0);
+
+    dimensionedScalar lambda(laminarTransport.lookup("lambda"));
+    dimensionedScalar alphaMax(laminarTransport.lookup("alphaMax"));
+
+    const labelList& inletCells =
+        mesh.boundary()[mesh.boundaryMesh().findPatchID("inlet")].faceCells();
+    //const labelList& outletCells =
+    //    mesh.boundary()[mesh.boundaryMesh().findPatchID("outlet")].faceCells();
+
+    volScalarField alpha
+    (
+        IOobject
+        (
+            "alpha",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        lambda*max(Ua & U, zeroSensitivity)
+    );
+    zeroCells(alpha, inletCells);
+    //zeroCells(alpha, outletCells);

applications/solvers/incompressible/ajointShapeOptimizationFoam/createPhia.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,70 +22,33 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Class
-    Foam::MCLimiter
+Global
+    createPhia
 
 Description
-    Class with limiter function which returns the limiter for the
-    monotonised centred differencing scheme based on r obtained from
-    the LimiterFunc class.
-
-    Used in conjunction with the template class LimitedScheme.
-
-SourceFiles
-    MC.C
+    Creates and initialises the face-flux field phia.
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef MC_H
-#define MC_H
-
-#include "vector.H"
+#ifndef createPhia_H
+#define createPhia_H
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
+Info<< "Reading/calculating face flux field phia\n" << endl;
 
-/*---------------------------------------------------------------------------*\
-                       Class MCLimiter Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class LimiterFunc>
-class MCLimiter
-:
-    public LimiterFunc
-{
-
-public:
-
-    MCLimiter(Istream&)
-    {}
-
-    scalar limiter
+surfaceScalarField phia
+(
+    IOobject
     (
-        const scalar cdWeight,
-        const scalar faceFlux,
-        const typename LimiterFunc::phiType& phiP,
-        const typename LimiterFunc::phiType& phiN,
-        const typename LimiterFunc::gradPhiType& gradcP,
-        const typename LimiterFunc::gradPhiType& gradcN,
-        const vector& d
-    ) const
-    {
-        scalar r = LimiterFunc::r
-        (
-            faceFlux, phiP, phiN, gradcP, gradcN, d
-        );
-
-        return max(min(min(2*r, 0.5*(1 + r)), 2), 0);
-    }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
+        "phia",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    linearInterpolate(Ua) & mesh.Sf()
+);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -71,7 +71,7 @@ Foam::label Foam::checkTopology
             {
                 Info<< " ***FaceZone " << mesh.faceZones()[zoneI].name()
                     << " is not correctly synchronised"
-                    << " acrosss coupled boundaries."
+                    << " across coupled boundaries."
                     << " (coupled faces both"
                     << " present in set but with opposite flipmap)" << endl;
                 noFailedChecks++;

applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C

@@ -110,40 +110,30 @@ void writeAllData
 {
     if (nPrims)
     {
-        PstreamBuffers pBufs(Pstream::nonBlocking);
-
-        if (!Pstream::master())
-        {
-            UOPstream toMaster(Pstream::masterNo(), pBufs);
-            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
-            {
-                toMaster<< map(vf, prims, cmpt);
-            }
-        }
-
-        pBufs.finishedSends();
-
         if (Pstream::master())
         {
             ensightFile << key << nl;
 
-            PtrList<UIPstream> fromSlaves(Pstream::nProcs());
-            for (int slave=1; slave<Pstream::nProcs(); slave++)
-            {
-                fromSlaves.set(slave, new UIPstream(slave, pBufs));
-            }
-
             for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
             {
                 writeData(map(vf, prims, cmpt), ensightFile);
 
                 for (int slave=1; slave<Pstream::nProcs(); slave++)
                 {
-                    scalarField data(fromSlaves[slave]);
+                    IPstream fromSlave(Pstream::scheduled, slave);
+                    scalarField data(fromSlave);
                     writeData(data, ensightFile);
                 }
             }
         }
+        else
+        {
+            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+            {
+                OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+                toMaster<< map(vf, prims, cmpt);
+            }
+        }
     }
 }
 
@@ -160,40 +150,30 @@ void writeAllDataBinary
 {
     if (nPrims)
     {
-        PstreamBuffers pBufs(Pstream::nonBlocking);
-
-        if (!Pstream::master())
-        {
-            UOPstream toMaster(Pstream::masterNo(), pBufs);
-            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
-            {
-                toMaster<< map(vf, prims, cmpt);
-            }
-        }
-
-        pBufs.finishedSends();
-
         if (Pstream::master())
         {
             writeEnsDataBinary(key,ensightFile);
 
-            PtrList<UIPstream> fromSlaves(Pstream::nProcs());
-            for (int slave=1; slave<Pstream::nProcs(); slave++)
-            {
-                fromSlaves.set(slave, new UIPstream(slave, pBufs));
-            }
-
             for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
             {
                 writeEnsDataBinary(map(vf, prims, cmpt), ensightFile);
 
                 for (int slave=1; slave<Pstream::nProcs(); slave++)
                 {
-                    scalarField data(fromSlaves[slave]);
+                    IPstream fromSlave(Pstream::scheduled, slave);
+                    scalarField data(fromSlave);
                     writeEnsDataBinary(data, ensightFile);
                 }
             }
         }
+        else
+        {
+            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+            {
+                OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+                toMaster<< map(vf, prims, cmpt);
+            }
+        }
     }
 }
 
@@ -210,40 +190,31 @@ void writeAllFaceData
 {
     if (nPrims)
     {
-        PstreamBuffers pBufs(Pstream::nonBlocking);
-
-        if (!Pstream::master())
-        {
-            UOPstream toMaster(Pstream::masterNo(), pBufs);
-            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
-            {
-                toMaster<< map(pf, prims, cmpt);
-            }
-        }
-
-        pBufs.finishedSends();
-
         if (Pstream::master())
         {
             ensightFile << key << nl;
 
-            PtrList<UIPstream> fromSlaves(Pstream::nProcs());
-            for (int slave=1; slave<Pstream::nProcs(); slave++)
-            {
-                fromSlaves.set(slave, new UIPstream(slave, pBufs));
-            }
-
             for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
             {
                 writeData(map(pf, prims, cmpt), ensightFile);
 
                 for (int slave=1; slave<Pstream::nProcs(); slave++)
                 {
-                    scalarField pf(fromSlaves[slave]);
+                    IPstream fromSlave(Pstream::scheduled, slave);
+                    scalarField pf(fromSlave);
+
                     writeData(pf, ensightFile);
                 }
             }
         }
+        else
+        {
+            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+            {
+                OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+                toMaster<< map(pf, prims, cmpt);
+            }
+        }
     }
 }
 
@@ -260,40 +231,31 @@ void writeAllFaceDataBinary
 {
     if (nPrims)
     {
-        PstreamBuffers pBufs(Pstream::nonBlocking);
-
-        if (!Pstream::master())
-        {
-            UOPstream toMaster(Pstream::masterNo(), pBufs);
-            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
-            {
-                toMaster<< map(pf, prims, cmpt);
-            }
-        }
-
-        pBufs.finishedSends();
-
         if (Pstream::master())
         {
             writeEnsDataBinary(key,ensightFile);
 
-            PtrList<UIPstream> fromSlaves(Pstream::nProcs());
-            for (int slave=1; slave<Pstream::nProcs(); slave++)
-            {
-                fromSlaves.set(slave, new UIPstream(slave, pBufs));
-            }
-
             for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
             {
                 writeEnsDataBinary(map(pf, prims, cmpt), ensightFile);
 
                 for (int slave=1; slave<Pstream::nProcs(); slave++)
                 {
-                    scalarField pf(fromSlaves[slave]);
+                    IPstream fromSlave(Pstream::scheduled, slave);
+                    scalarField pf(fromSlave);
+
                     writeEnsDataBinary(pf, ensightFile);
                 }
             }
         }
+        else
+        {
+            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+            {
+                OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+                toMaster<< map(pf, prims, cmpt);
+            }
+        }
     }
 }
 
@@ -539,6 +501,10 @@ void ensightFieldAscii
     const wordHashSet& patchNames = eMesh.patchNames();
     const HashTable<ensightMesh::nFacePrimitives>&
         nPatchPrims = eMesh.nPatchPrims();
+    const List<faceSets>& faceZoneFaceSets = eMesh.faceZoneFaceSets();
+    const wordHashSet& faceZoneNames = eMesh.faceZoneNames();
+    const HashTable<ensightMesh::nFacePrimitives>&
+        nFaceZonePrims = eMesh.nFaceZonePrims();
 
     const labelList& tets = meshCellSets.tets;
     const labelList& pyrs = meshCellSets.pyrs;
@@ -557,7 +523,7 @@ void ensightFieldAscii
             postProcPath/ensightFileName,
             runTime.writeFormat(),
             runTime.writeVersion(),
-            runTime.writeCompression()
+            IOstream::UNCOMPRESSED
         );
     }
 
@@ -567,6 +533,8 @@ void ensightFieldAscii
 
     if (patchNames.empty())
     {
+        eMesh.barrier();
+
         if (Pstream::master())
         {
             if (timeIndex == 0)
@@ -592,29 +560,10 @@ void ensightFieldAscii
 
         if (meshCellSets.nHexesWedges)
         {
-            PstreamBuffers pBufs(Pstream::nonBlocking);
-
-            if (!Pstream::master())
-            {
-                UOPstream toMaster(Pstream::masterNo(), pBufs);
-                for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
-                {
-                    toMaster<< map(vf, hexes, wedges, cmpt);
-                }
-            }
-
-            pBufs.finishedSends();
-
             if (Pstream::master())
             {
                 ensightFile << "hexa8" << nl;
 
-                PtrList<UIPstream> fromSlaves(Pstream::nProcs());
-                for (int slave=1; slave<Pstream::nProcs(); slave++)
-                {
-                    fromSlaves.set(slave, new UIPstream(slave, pBufs));
-                }
-
                 for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
                 {
                     writeData
@@ -625,17 +574,57 @@ void ensightFieldAscii
 
                     for (int slave=1; slave<Pstream::nProcs(); slave++)
                     {
-                        scalarField data(fromSlaves[slave]);
+                        IPstream fromSlave(Pstream::scheduled, slave);
+                        scalarField data(fromSlave);
                         writeData(data, ensightFile);
                     }
                 }
             }
+            else
+            {
+                for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+                {
+                    OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+                    toMaster<< map(vf, hexes, wedges, cmpt);
+                }
+            }
         }
 
-        writeAllData("penta6", vf, prisms, meshCellSets.nPrisms, ensightFile);
-        writeAllData("pyramid5", vf, pyrs, meshCellSets.nPyrs, ensightFile);
-        writeAllData("tetra4", vf, tets, meshCellSets.nTets, ensightFile);
-        writeAllData("nfaced", vf, polys, meshCellSets.nPolys, ensightFile);
+        writeAllData
+        (
+            "penta6",
+            vf,
+            prisms,
+            meshCellSets.nPrisms,
+            ensightFile
+        );
+
+        writeAllData
+        (
+            "pyramid5",
+            vf,
+            pyrs,
+            meshCellSets.nPyrs,
+            ensightFile
+        );
+
+        writeAllData
+        (
+            "tetra4",
+            vf,
+            tets,
+            meshCellSets.nTets,
+            ensightFile
+        );
+
+        writeAllData
+        (
+            "nfaced",
+            vf,
+            polys,
+            meshCellSets.nPolys,
+            ensightFile
+        );
     }
 
     label ensightPatchI = eMesh.patchPartOffset();
@@ -644,6 +633,8 @@ void ensightFieldAscii
     {
         const word& patchName = allPatchNames[patchi];
 
+        eMesh.barrier();
+
         if (patchNames.empty() || patchNames.found(patchName))
         {
             if
@@ -664,6 +655,80 @@ void ensightFieldAscii
         }
     }
 
+    // write faceZones, if requested
+    if (faceZoneNames.size())
+    {
+        // Interpolates cell values to faces - needed only when exporting
+        // faceZones...
+        GeometricField<Type, fvsPatchField, surfaceMesh> sf
+        (
+            linearInterpolate(vf)
+        );
+
+        forAllConstIter(wordHashSet, faceZoneNames, iter)
+        {
+            const word& faceZoneName = iter.key();
+
+            eMesh.barrier();
+
+            label zoneID = mesh.faceZones().findZoneID(faceZoneName);
+
+            const faceZone& fz = mesh.faceZones()[zoneID];
+
+            // Prepare data to write
+            label nIncluded = 0;
+            forAll(fz, i)
+            {
+                if (eMesh.faceToBeIncluded(fz[i]))
+                {
+                    ++nIncluded;
+                }
+            }
+
+            Field<Type> values(nIncluded);
+
+            // Loop on the faceZone and store the needed field values
+            label j = 0;
+            forAll(fz, i)
+            {
+                label faceI = fz[i];
+                if (mesh.isInternalFace(faceI))
+                {
+                    values[j] = sf[faceI];
+                    ++j;
+                }
+                else
+                {
+                    if (eMesh.faceToBeIncluded(faceI))
+                    {
+                        label patchI = mesh.boundaryMesh().whichPatch(faceI);
+                        const polyPatch& pp = mesh.boundaryMesh()[patchI];
+                        label patchFaceI = pp.whichFace(faceI);
+                        Type value = sf.boundaryField()[patchI][patchFaceI];
+                        values[j] = value;
+                        ++j;
+                    }
+                }
+            }
+
+            if
+            (
+                writePatchField
+                (
+                    values,
+                    zoneID,
+                    ensightPatchI,
+                    faceZoneFaceSets[zoneID],
+                    nFaceZonePrims.find(faceZoneName)(),
+                    ensightFile
+                )
+            )
+            {
+                ensightPatchI++;
+            }
+        }
+    }
+
     if (Pstream::master())
     {
         delete ensightFilePtr;
@@ -695,6 +760,10 @@ void ensightFieldBinary
     const wordHashSet& patchNames = eMesh.patchNames();
     const HashTable<ensightMesh::nFacePrimitives>&
         nPatchPrims = eMesh.nPatchPrims();
+    const List<faceSets>& faceZoneFaceSets = eMesh.faceZoneFaceSets();
+    const wordHashSet& faceZoneNames = eMesh.faceZoneNames();
+    const HashTable<ensightMesh::nFacePrimitives>&
+        nFaceZonePrims = eMesh.nFaceZonePrims();
 
     const labelList& tets = meshCellSets.tets;
     const labelList& pyrs = meshCellSets.pyrs;
@@ -722,6 +791,8 @@ void ensightFieldBinary
 
     if (patchNames.empty())
     {
+        eMesh.barrier();
+
         if (Pstream::master())
         {
             if (timeIndex == 0)
@@ -743,29 +814,10 @@ void ensightFieldBinary
 
         if (meshCellSets.nHexesWedges)
         {
-            PstreamBuffers pBufs(Pstream::nonBlocking);
-
-            if (!Pstream::master())
-            {
-                UOPstream toMaster(Pstream::masterNo(), pBufs);
-                for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
-                {
-                    toMaster<< map(vf, hexes, wedges, cmpt);
-                }
-            }
-
-            pBufs.finishedSends();
-
             if (Pstream::master())
             {
                 writeEnsDataBinary("hexa8",ensightFile);
 
-                PtrList<UIPstream> fromSlaves(Pstream::nProcs());
-                for (int slave=1; slave<Pstream::nProcs(); slave++)
-                {
-                    fromSlaves.set(slave, new UIPstream(slave, pBufs));
-                }
-
                 for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
                 {
                     writeEnsDataBinary
@@ -776,11 +828,20 @@ void ensightFieldBinary
 
                     for (int slave=1; slave<Pstream::nProcs(); slave++)
                     {
-                        scalarField data(fromSlaves[slave]);
+                        IPstream fromSlave(Pstream::scheduled, slave);
+                        scalarField data(fromSlave);
                         writeEnsDataBinary(data, ensightFile);
                     }
                 }
             }
+            else
+            {
+                for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+                {
+                    OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+                    toMaster<< map(vf, hexes, wedges, cmpt);
+                }
+            }
         }
 
         writeAllDataBinary
@@ -826,6 +887,8 @@ void ensightFieldBinary
     {
         const word& patchName = allPatchNames[patchi];
 
+        eMesh.barrier();
+
         if (patchNames.empty() || patchNames.found(patchName))
         {
             if
@@ -844,6 +907,81 @@ void ensightFieldBinary
                 ensightPatchI++;
             }
         }
+
+    }
+
+    // write faceZones, if requested
+    if (faceZoneNames.size())
+    {
+        // Interpolates cell values to faces - needed only when exporting
+        // faceZones...
+        GeometricField<Type, fvsPatchField, surfaceMesh> sf
+        (
+            linearInterpolate(vf)
+        );
+
+        forAllConstIter(wordHashSet, faceZoneNames, iter)
+        {
+            const word& faceZoneName = iter.key();
+
+            eMesh.barrier();
+
+            label zoneID = mesh.faceZones().findZoneID(faceZoneName);
+
+            const faceZone& fz = mesh.faceZones()[zoneID];
+
+            // Prepare data to write
+            label nIncluded = 0;
+            forAll(fz, i)
+            {
+                if (eMesh.faceToBeIncluded(fz[i]))
+                {
+                    ++nIncluded;
+                }
+            }
+
+            Field<Type> values(nIncluded);
+
+            // Loop on the faceZone and store the needed field values
+            label j = 0;
+            forAll(fz, i)
+            {
+                label faceI = fz[i];
+                if (mesh.isInternalFace(faceI))
+                {
+                    values[j] = sf[faceI];
+                    ++j;
+                }
+                else
+                {
+                    if (eMesh.faceToBeIncluded(faceI))
+                    {
+                        label patchI = mesh.boundaryMesh().whichPatch(faceI);
+                        const polyPatch& pp = mesh.boundaryMesh()[patchI];
+                        label patchFaceI = pp.whichFace(faceI);
+                        Type value = sf.boundaryField()[patchI][patchFaceI];
+                        values[j] = value;
+                        ++j;
+                    }
+                }
+            }
+
+            if
+            (
+                writePatchFieldBinary
+                (
+                    values,
+                    zoneID,
+                    ensightPatchI,
+                    faceZoneFaceSets[zoneID],
+                    nFaceZonePrims.find(faceZoneName)(),
+                    ensightFile
+                )
+            )
+            {
+                ensightPatchI++;
+            }
+        }
     }
 
     if (Pstream::master())

applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.H

@@ -42,7 +42,7 @@ SourceFiles
 #include "fvMesh.H"
 #include "OFstream.H"
 #include <fstream>
-#include "globalIndex.H"
+#include "PackedBoolList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,6 +51,7 @@ namespace Foam
 
 class fvMesh;
 class argList;
+class globalIndex;
 
 /*---------------------------------------------------------------------------*\
                            Class ensightMesh Declaration
@@ -98,6 +99,16 @@ private:
 
         HashTable<nFacePrimitives> nPatchPrims_;
 
+        // faceZone - related variables
+        List<faceSets> faceZoneFaceSets_;
+
+        wordHashSet faceZoneNames_;
+
+        HashTable<nFacePrimitives> nFaceZonePrims_;
+
+        //- Per boundary face whether to include or not
+        PackedBoolList boundaryFaceToBeIncluded_;
+
 
     // Private Member Functions
 
@@ -206,6 +217,38 @@ private:
             OFstream& ensightGeometryFile
         ) const;
 
+        void writeAllInternalPoints
+        (
+            const pointField& uniquePoints,
+            const label nPoints,
+            OFstream& ensightGeometryFile
+        ) const;
+
+        void writeAllPatchPoints
+        (
+            label ensightPatchI,
+            const word& patchName,
+            const pointField& uniquePoints,
+            const label nPoints,
+            OFstream& ensightGeometryFile
+        ) const;
+
+        void writeAllInternalPointsBinary
+        (
+            const pointField& uniquePoints,
+            const label nPoints,
+            std::ofstream& ensightGeometryFile
+        ) const;
+
+        void writeAllPatchPointsBinary
+        (
+            label ensightPatchI,
+            const word& patchName,
+            const pointField& uniquePoints,
+            const label nPoints,
+            std::ofstream& ensightGeometryFile
+        ) const;
+
         void writeAscii
         (
             const fileName& postProcPath,
@@ -311,7 +354,12 @@ public:
     // Constructors
 
         //- Construct from fvMesh
-        ensightMesh(const fvMesh&, const argList& args, const bool binary);
+        ensightMesh
+        (
+            const fvMesh&,
+            const argList& args,
+            const bool binary
+        );
 
 
     //- Destructor
@@ -352,12 +400,36 @@ public:
                 return nPatchPrims_;
             }
 
+            const List<faceSets>& faceZoneFaceSets() const
+            {
+                return faceZoneFaceSets_;
+            }
+
+            const wordHashSet& faceZoneNames() const
+            {
+                return faceZoneNames_;
+            }
+
+            const HashTable<nFacePrimitives>& nFaceZonePrims() const
+            {
+                return nFaceZonePrims_;
+            }
+
             //- The ensight part id for the first patch
             label patchPartOffset() const
             {
                 return patchPartOffset_;
             }
 
+    // Other
+
+        //- When exporting faceZones, check if a given face has to be included
+        //  or not (i.e. faces on processor boundaries)
+        bool faceToBeIncluded(const label faceI) const;
+
+        //- Helper to cause barrier. Necessary on Quadrics.
+        static void barrier();
+
 
     // I-O
 
@@ -368,6 +440,7 @@ public:
             const label timeIndex,
             Ostream& ensightCaseFile
         ) const;
+
 };
 
 

applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C

@@ -41,6 +41,9 @@ Usage
     @param -noPatches \n
     Suppress writing any patches.
 
+    @param -faceZones zoneList \n
+    Specify faceZones to write, with wildcards
+
 Note
     Parallel support for cloud data is not supported
     - writes to @a EnSight directory to avoid collisions with foamToEnsightParts
@@ -113,6 +116,12 @@ int main(int argc, char *argv[])
         "specify particular patches to write - eg '(inlet outlet)'. "
         "An empty list suppresses writing the internalMesh."
     );
+    argList::addOption
+    (
+        "faceZones",
+        "wordList",
+        "specify faceZones to write, with wildcards - eg '(mfp-.*)'. "
+    );
 
 #   include "setRootCase.H"
 

applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C

@@ -1019,6 +1019,7 @@ int main(int argc, char *argv[])
             }
         }
         writer.writeEnd();
+        Info<< endl;
 
 
 
@@ -1067,7 +1068,7 @@ int main(int argc, char *argv[])
             writer.writeInit
             (
                 runTime.caseName(),
-                cellVarNames,
+                allVarNames,
                 patchFileName,
                 DataFileType_Full
             );
@@ -1163,6 +1164,8 @@ int main(int argc, char *argv[])
                         << nl << endl;
                 }
             }
+            writer.writeEnd();
+            Info<< endl;
         }
 
 

src/OpenFOAM/db/IOstreams/Pstreams/IPstream.C

@@ -38,7 +38,17 @@ Foam::IPstream::IPstream
 )
 :
     Pstream(commsType, bufSize),
-    UIPstream(commsType, fromProcNo, buf_, externalBufPosition_),
+    UIPstream
+    (
+        commsType,
+        fromProcNo,
+        buf_,
+        externalBufPosition_,
+        UPstream::msgType(),        // tag
+        false,                      // do not clear buf_ if at end
+        format,
+        version
+    ),
     externalBufPosition_(0)
 {}
 

src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.C

@@ -77,6 +77,21 @@ inline void Foam::UIPstream::readFromBuffer
 }
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::UIPstream::~UIPstream()
+{
+    if (clearAtEnd_ && eof())
+    {
+        if (debug)
+        {
+            Pout<< "UIPstream::~UIPstream() : clearing externalBuf_ of size "
+                << externalBuf_.size() << endl;
+        }
+        externalBuf_.clearStorage();
+    }
+}
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::Istream& Foam::UIPstream::read(token& t)

src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.H

@@ -67,6 +67,8 @@ class UIPstream
 
         const int tag_;
 
+        const bool clearAtEnd_;
+
         label messageSize_;
 
 
@@ -96,6 +98,7 @@ public:
             DynamicList<char>& externalBuf,
             label& externalBufPosition,
             const int tag = UPstream::msgType(),
+            const bool clearAtEnd = false,   // destroy externalBuf if at end
             streamFormat format=BINARY,
             versionNumber version=currentVersion
         );
@@ -104,6 +107,11 @@ public:
         UIPstream(const int fromProcNo, PstreamBuffers&);
 
 
+    // Destructor
+
+        ~UIPstream();
+
+
     // Member functions
 
         // Inquiry

src/OpenFOAM/db/Time/Time.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -625,21 +625,25 @@ void Foam::Time::setDeltaT(const scalar deltaT)
 Foam::TimeState Foam::Time::subCycle(const label nSubCycles)
 {
     subCycling_ = true;
+    prevTimeState_.set(new TimeState(*this));
 
-    TimeState ts = *this;
     setTime(*this - deltaT(), (timeIndex() - 1)*nSubCycles);
     deltaT_ /= nSubCycles;
     deltaT0_ /= nSubCycles;
     deltaTSave_ = deltaT0_;
 
-    return ts;
+    return prevTimeState();
 }
 
 
-void Foam::Time::endSubCycle(const TimeState& ts)
+void Foam::Time::endSubCycle()
 {
-    subCycling_ = false;
-    TimeState::operator=(ts);
+    if (subCycling_)
+    {
+        subCycling_ = false;
+        TimeState::operator=(prevTimeState());
+        prevTimeState_.clear();
+    }
 }
 
 

src/OpenFOAM/db/Time/Time.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -127,6 +127,9 @@ protected:
         //- Is the time currently being sub-cycled?
         bool subCycling_;
 
+        //- If time is being sub-cycled this is the previous TimeState
+        autoPtr<TimeState> prevTimeState_;
+
         //- Time directory name format
         static fmtflags format_;
 
@@ -347,6 +350,18 @@ public:
                 return functionObjects_;
             }
 
+            //- Return true if time currently being sub-cycled, otherwise false
+            bool subCycling() const
+            {
+                return subCycling_;
+            }
+
+            //- Return previous TimeState if time is being sub-cycled
+            const TimeState& prevTimeState() const
+            {
+                return prevTimeState_();
+            }
+
 
         // Check
 
@@ -427,8 +442,8 @@ public:
             //- Set time to sub-cycle for the given number of steps
             virtual TimeState subCycle(const label nSubCycles);
 
-            //- Reset time after sub-cycling back to given TimeState
-            virtual void endSubCycle(const TimeState&);
+            //- Reset time after sub-cycling back to previous TimeState
+            virtual void endSubCycle();
 
             //- Return non-const access to the list of function objects
             functionObjectList& functionObjects()

src/OpenFOAM/db/Time/subCycleTime.C

@@ -32,9 +32,10 @@ Foam::subCycleTime::subCycleTime(Time& t, const label nSubCycles)
 :
     time_(t),
     nSubCycles_(nSubCycles),
-    subCycleIndex_(0),
-    initialTimeState_(time_.subCycle(nSubCycles_))
-{}
+    subCycleIndex_(0)
+{
+    time_.subCycle(nSubCycles_);
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -55,7 +56,7 @@ bool Foam::subCycleTime::end() const
 
 void Foam::subCycleTime::endSubCycle()
 {
-    time_.endSubCycle(initialTimeState_);
+    time_.endSubCycle();
 }
 
 

src/OpenFOAM/db/Time/subCycleTime.H

@@ -55,7 +55,6 @@ class subCycleTime
 
         label nSubCycles_;
         label subCycleIndex_;
-        TimeState initialTimeState_;
 
 public:
 

src/Pstream/dummy/UIPread.C

@@ -38,6 +38,7 @@ Foam::UIPstream::UIPstream
     DynamicList<char>& externalBuf,
     label& externalBufPosition,
     const int tag,
+    const bool clearAtEnd,
     streamFormat format,
     versionNumber version
 )
@@ -48,6 +49,7 @@ Foam::UIPstream::UIPstream
     externalBuf_(externalBuf),
     externalBufPosition_(externalBufPosition),
     tag_(tag),
+    clearAtEnd_(clearAtEnd),
     messageSize_(0)
 {
     notImplemented
@@ -59,6 +61,7 @@ Foam::UIPstream::UIPstream
             "DynamicList<char>&,"
             "label&,"
             "const int tag,"
+            "const bool,"
             "streamFormat, versionNumber"
         ")"
     );

src/Pstream/mpi/UIPread.C

@@ -42,6 +42,7 @@ Foam::UIPstream::UIPstream
     DynamicList<char>& externalBuf,
     label& externalBufPosition,
     const int tag,
+    const bool clearAtEnd,
     streamFormat format,
     versionNumber version
 )
@@ -52,6 +53,7 @@ Foam::UIPstream::UIPstream
     externalBuf_(externalBuf),
     externalBufPosition_(externalBufPosition),
     tag_(tag),
+    clearAtEnd_(clearAtEnd),
     messageSize_(0)
 {
     setOpened();
@@ -125,6 +127,7 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers)
     externalBuf_(buffers.recvBuf_[fromProcNo]),
     externalBufPosition_(buffers.recvBufPos_[fromProcNo]),
     tag_(buffers.tag_),
+    clearAtEnd_(true),
     messageSize_(0)
 {
     if (commsType() != UPstream::scheduled && !buffers.finishedSendsCalled_)

src/finiteVolume/Make/files

@@ -242,7 +242,6 @@ $(limitedSchemes)/SuperBee/SuperBee.C
 $(limitedSchemes)/QUICK/QUICK.C
 $(limitedSchemes)/MUSCL/MUSCL.C
 $(limitedSchemes)/UMIST/UMIST.C
-$(limitedSchemes)/MC/MC.C
 $(limitedSchemes)/Phi/Phi.C
 $(limitedSchemes)/filteredLinear/filteredLinear.C
 $(limitedSchemes)/filteredLinear2/filteredLinear2.C

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,7 @@ License
 #include "fvcSurfaceIntegrate.H"
 #include "slicedSurfaceFields.H"
 #include "syncTools.H"
+
 #include "fvm.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -107,7 +108,7 @@ void Foam::MULES::explicitSolve
     {
         psiIf =
         (
-            mesh.V0()*rho.oldTime()*psi0/(deltaT*mesh.V())
+            mesh.Vsc0()*rho.oldTime()*psi0/(deltaT*mesh.Vsc())
           + Su.field()
           - psiIf
         )/(rho/deltaT - Sp.field());
@@ -325,7 +326,8 @@ void Foam::MULES::limiter
 
     const unallocLabelList& owner = mesh.owner();
     const unallocLabelList& neighb = mesh.neighbour();
-    const scalarField& V = mesh.V();
+    tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
+    const scalarField& V = tVsc();
     const scalar deltaT = mesh.time().deltaTValue();
 
     const scalarField& phiBDIf = phiBD;
@@ -452,14 +454,16 @@ void Foam::MULES::limiter
 
     if (mesh.moving())
     {
+        tmp<volScalarField::DimensionedInternalField> V0 = mesh.Vsc0();
+
         psiMaxn =
             V*((rho/deltaT - Sp)*psiMaxn - Su)
-          - (mesh.V0()/deltaT)*rho.oldTime()*psi0
+          - (V0()/deltaT)*rho.oldTime()*psi0
           + sumPhiBD;
 
         psiMinn =
             V*(Su - (rho/deltaT - Sp)*psiMinn)
-          + (mesh.V0()/deltaT)*rho.oldTime()*psi0
+          + (V0/deltaT)*rho.oldTime()*psi0
           - sumPhiBD;
     }
     else

src/finiteVolume/fvMesh/fvMesh.H

@@ -269,6 +269,12 @@ public:
             //- Return old-old-time cell volumes
             const DimensionedField<scalar, volMesh>& V00() const;
 
+            //- Return sub-cycle cell volumes
+            tmp<DimensionedField<scalar, volMesh> > Vsc() const;
+
+            //- Return sub-cycl old-time cell volumes
+            tmp<DimensionedField<scalar, volMesh> > Vsc0() const;
+
             //- Return cell face area vectors
             const surfaceVectorField& Sf() const;
 

src/finiteVolume/fvMesh/fvMeshGeometry.C

@@ -285,6 +285,63 @@ const volScalarField::DimensionedInternalField& fvMesh::V00() const
 }
 
 
+tmp<volScalarField::DimensionedInternalField> fvMesh::Vsc() const
+{
+    if (moving() && time().subCycling())
+    {
+        const TimeState& ts = time();
+        const TimeState& ts0 = time().prevTimeState();
+
+        scalar tFrac =
+        (
+            ts.value() - (ts0.value() - ts0.deltaTValue())
+        )/ts0.deltaTValue();
+
+        if (tFrac < (1 - SMALL))
+        {
+            return V0() + tFrac*(V() - V0());
+        }
+        else
+        {
+            return V();
+        }
+    }
+    else
+    {
+        return V();
+    }
+}
+
+
+tmp<volScalarField::DimensionedInternalField> fvMesh::Vsc0() const
+{
+    if (moving() && time().subCycling())
+    {
+        const TimeState& ts = time();
+        const TimeState& ts0 = time().prevTimeState();
+
+        scalar t0Frac =
+        (
+            (ts.value() - ts.deltaTValue())
+          - (ts0.value() - ts0.deltaTValue())
+        )/ts0.deltaTValue();
+
+        if (t0Frac > SMALL)
+        {
+            return V0() + t0Frac*(V() - V0());
+        }
+        else
+        {
+            return V0();
+        }
+    }
+    else
+    {
+        return V0();
+    }
+}
+
+
 const surfaceVectorField& fvMesh::Sf() const
 {
     if (!SfPtr_)

src/fvMotionSolver/Make/files

@@ -1,6 +1,8 @@
 fvMotionSolvers/fvMotionSolver/fvMotionSolver.C
 fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
 fvMotionSolvers/displacement/displacementFvMotionSolver/displacementFvMotionSolver.C
+fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.C
+fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.C
 fvMotionSolvers/displacement/interpolation/displacementInterpolationFvMotionSolver.C
 fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
 fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C

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

@@ -0,0 +1,567 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "displacementLayeredMotionFvMotionSolver.H"
+#include "addToRunTimeSelectionTable.H"
+#include "pointEdgeStructuredWalk.H"
+#include "pointFields.H"
+#include "PointEdgeWave.H"
+#include "syncTools.H"
+#include "interpolationTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(displacementLayeredMotionFvMotionSolver, 0);
+
+    addToRunTimeSelectionTable
+    (
+        fvMotionSolver,
+        displacementLayeredMotionFvMotionSolver,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::displacementLayeredMotionFvMotionSolver::calcZoneMask
+(
+    const label cellZoneI, 
+    PackedBoolList& isZonePoint,
+    PackedBoolList& isZoneEdge
+) const
+{
+    if (cellZoneI == -1)
+    {
+        isZonePoint.setSize(mesh().nPoints());
+        isZonePoint = 1;
+
+        isZoneEdge.setSize(mesh().nEdges());
+        isZoneEdge = 1;
+    }
+    else
+    {
+        const cellZone& cz = mesh().cellZones()[cellZoneI];
+
+        label nPoints = 0;
+        forAll(cz, i)
+        {
+            const labelList& cPoints = mesh().cellPoints(cz[i]);
+            forAll(cPoints, cPointI)
+            {
+                if (!isZonePoint[cPoints[cPointI]])
+                {
+                    isZonePoint[cPoints[cPointI]] = 1;
+                    nPoints++;
+                }
+            }
+        }
+        syncTools::syncPointList
+        (
+            mesh(),
+            isZonePoint,
+            orEqOp<unsigned int>(),
+            0
+        );
+
+
+        // Mark edge inside cellZone
+        label nEdges = 0;
+        forAll(cz, i)
+        {
+            const labelList& cEdges = mesh().cellEdges(cz[i]);
+            forAll(cEdges, cEdgeI)
+            {
+                if (!isZoneEdge[cEdges[cEdgeI]])
+                {
+                    isZoneEdge[cEdges[cEdgeI]] = 1;
+                    nEdges++;
+                }
+            }
+        }
+        syncTools::syncEdgeList
+        (
+            mesh(),
+            isZoneEdge,
+            orEqOp<unsigned int>(),
+            0
+        );
+
+
+        Info<< "On cellZone " << cz.name()
+            << " marked " << returnReduce(nPoints, sumOp<label>())
+            << " points and " << returnReduce(nEdges, sumOp<label>())
+            << " edges." << endl;
+    }
+}
+
+
+// Find distance to starting point
+void Foam::displacementLayeredMotionFvMotionSolver::walkStructured
+(
+    const label cellZoneI,
+    const PackedBoolList& isZonePoint,
+    const PackedBoolList& isZoneEdge,
+    const labelList& seedPoints,
+    const vectorField& seedData,
+    scalarField& distance,
+    vectorField& data
+) const
+{
+    const pointField& points = mesh().points();
+
+    List<pointEdgeStructuredWalk> seedInfo(seedPoints.size());
+
+    forAll(seedPoints, i)
+    {
+        seedInfo[i] = pointEdgeStructuredWalk
+        (
+            true,
+            points[seedPoints[i]],
+            0.0,
+            seedData[i]
+        );
+    }
+
+    // Current info on points
+    List<pointEdgeStructuredWalk> allPointInfo(mesh().nPoints());
+    // Mark points inside cellZone
+    forAll(isZonePoint, pointI)
+    {
+        if (isZonePoint[pointI])
+        {
+            allPointInfo[pointI].inZone() = true;
+        }
+    }
+
+    // Current info on edges
+    List<pointEdgeStructuredWalk> allEdgeInfo(mesh().nEdges());
+    // Mark edges inside cellZone
+    forAll(isZoneEdge, edgeI)
+    {
+        if (isZoneEdge[edgeI])
+        {
+            allEdgeInfo[edgeI].inZone() = true;
+        }
+    }
+
+    // Walk
+    PointEdgeWave<pointEdgeStructuredWalk> wallCalc
+    (
+        mesh(),
+        seedPoints,
+        seedInfo,
+
+        allPointInfo,
+        allEdgeInfo,
+        mesh().globalData().nTotalPoints()  // max iterations
+    );
+
+    // Extract distance and passive data
+    forAll(allPointInfo, pointI)
+    {
+        if (isZonePoint[pointI])
+        {
+            distance[pointI] = allPointInfo[pointI].dist();
+            data[pointI] = allPointInfo[pointI].data();
+        }
+    }
+}
+
+
+// Evaluate faceZone patch
+Foam::tmp<Foam::vectorField>
+Foam::displacementLayeredMotionFvMotionSolver::faceZoneEvaluate
+(
+    const faceZone& fz,
+    const labelList& meshPoints,
+    const dictionary& dict,
+    const PtrList<pointVectorField>& patchDisp,
+    const label patchI
+) const
+{
+    tmp<vectorField> tfld(new vectorField(meshPoints.size()));
+    vectorField& fld = tfld();
+
+    const word& type = dict.lookup("type");
+
+    if (type == "fixedValue")
+    {
+        fld = vectorField("value", dict, meshPoints.size());
+    }
+    else if (type == "timeVaryingUniformFixedValue")
+    {
+        interpolationTable<vector> timeSeries(dict);
+
+        fld = timeSeries(mesh().time().timeOutputValue());
+    }
+    else if (type == "slip")
+    {
+        if ((patchI % 2) != 1)
+        {
+            FatalIOErrorIn
+            (
+                "displacementLayeredMotionFvMotionSolver::faceZoneEvaluate(..)",
+                *this
+            )   << "slip can only be used on second faceZonePatch of pair."
+                << "FaceZone:" << fz.name()
+                << exit(FatalIOError);
+        }
+        // Use field set by previous bc
+        fld = vectorField(patchDisp[patchI-1], meshPoints);
+    }
+    else if (type == "follow")
+    {
+        // Only on boundary faces - follow boundary conditions
+        fld = vectorField(pointDisplacement_, meshPoints);
+    }
+    else
+    {
+        FatalIOErrorIn
+        (
+            "displacementLayeredMotionFvMotionSolver::faceZoneEvaluate(..)",
+            *this
+        )   << "Unknown faceZonePatch type " << type << " for faceZone "
+            << fz.name() << exit(FatalIOError);
+    }
+    return tfld;
+}
+
+
+void Foam::displacementLayeredMotionFvMotionSolver::cellZoneSolve
+(
+    const label cellZoneI,
+    const dictionary& zoneDict
+)
+{
+    PackedBoolList isZonePoint(mesh().nPoints());
+    PackedBoolList isZoneEdge(mesh().nEdges());
+    calcZoneMask(cellZoneI, isZonePoint, isZoneEdge);
+
+    const dictionary& patchesDict = zoneDict.subDict("boundaryField");
+
+    if (patchesDict.size() != 2)
+    {
+        FatalIOErrorIn
+        (
+            "displacementLayeredMotionFvMotionSolver::"
+            "correctBoundaryConditions(..)",
+            *this
+        )   << "Can only handle 2 faceZones (= patches) per cellZone. "
+            << " cellZone:" << cellZoneI
+            << " patches:" << patchesDict.toc()
+            << exit(FatalIOError);
+    }
+
+    PtrList<scalarField> patchDist(patchesDict.size());
+    PtrList<pointVectorField> patchDisp(patchesDict.size());
+
+    // Allocate the fields
+    label patchI = 0;
+    forAllConstIter(dictionary, patchesDict, patchIter)
+    {
+        const word& faceZoneName = patchIter().keyword();
+        label zoneI = mesh().faceZones().findZoneID(faceZoneName);
+        if (zoneI == -1)
+        {
+            FatalIOErrorIn
+            (
+                "displacementLayeredMotionFvMotionSolver::"
+                "correctBoundaryConditions(..)",
+                *this
+            )   << "Cannot find faceZone " << faceZoneName
+                << endl << "Valid zones are " << mesh().faceZones().names()
+                << exit(FatalIOError);
+        }
+
+        // Determine the points of the faceZone within the cellZone
+        const faceZone& fz = mesh().faceZones()[zoneI];
+
+        patchDist.set(patchI, new scalarField(mesh().nPoints()));
+        patchDisp.set
+        (
+            patchI,
+            new pointVectorField
+            (
+                IOobject
+                (
+                    mesh().cellZones()[cellZoneI].name() + "_" + fz.name(),
+                    mesh().time().timeName(),
+                    mesh(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                pointDisplacement_  // to inherit the boundary conditions
+            )
+        );
+
+        patchI++;
+    }
+
+
+
+    // 'correctBoundaryConditions'
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Loops over all the faceZones and walks their boundary values
+
+    // Make sure we can pick up bc values from field
+    pointDisplacement_.correctBoundaryConditions();
+
+    patchI = 0;
+    forAllConstIter(dictionary, patchesDict, patchIter)
+    {
+        const word& faceZoneName = patchIter().keyword();
+        const dictionary& faceZoneDict = patchIter().dict();
+
+        // Determine the points of the faceZone within the cellZone
+        label zoneI = mesh().faceZones().findZoneID(faceZoneName);
+        const faceZone& fz = mesh().faceZones()[zoneI];
+        const labelList& fzMeshPoints = fz().meshPoints();
+        DynamicList<label> meshPoints(fzMeshPoints.size());
+        forAll(fzMeshPoints, i)
+        {
+            if (isZonePoint[fzMeshPoints[i]])
+            {
+                meshPoints.append(fzMeshPoints[i]);
+            }
+        }
+
+        // Get initial value for all the faceZone points
+        tmp<vectorField> tseed = faceZoneEvaluate
+        (
+            fz,
+            meshPoints,
+            faceZoneDict,
+            patchDisp,
+            patchI
+        );
+
+
+Info<< "For cellZone:" << cellZoneI
+    << " for faceZone:" << fz.name() << " nPoints:" << tseed().size()
+    << " have patchField:"
+    << " max:" << gMax(tseed())
+    << " min:" << gMin(tseed())
+    << " avg:" << gAverage(tseed())
+    << endl;
+
+        // Set distance and transported value
+        walkStructured
+        (
+            cellZoneI,
+            isZonePoint,
+            isZoneEdge,
+
+            meshPoints,
+            tseed,
+            patchDist[patchI],
+            patchDisp[patchI]
+        );
+
+        // Implement real bc.
+        patchDisp[patchI].correctBoundaryConditions();
+
+
+//Info<< "Writing displacement for faceZone " << fz.name()
+//    << " to " << patchDisp[patchI].name() << endl;
+//patchDisp[patchI].write();
+
+//        // Copy into pointDisplacement for other fields to use
+//        forAll(isZonePoint, pointI)
+//        {
+//            if (isZonePoint[pointI])
+//            {
+//                pointDisplacement_[pointI] = patchDisp[patchI][pointI];
+//            }
+//        }
+//        pointDisplacement_.correctBoundaryConditions();
+
+
+        patchI++;
+    }
+
+
+    // Solve
+    // ~~~~~
+    // solving the interior is just interpolating
+
+//    // Get normalised distance
+//    pointScalarField distance
+//    (
+//        IOobject
+//        (
+//            "distance",
+//            mesh().time().timeName(),
+//            mesh(),
+//            IOobject::NO_READ,
+//            IOobject::NO_WRITE,
+//            false
+//        ),
+//        pointMesh::New(mesh()),
+//        dimensionedScalar("distance", dimLength, 0.0)
+//    );
+//    forAll(distance, pointI)
+//    {
+//        if (isZonePoint[pointI])
+//        {
+//            scalar d1 = patchDist[0][pointI];
+//            scalar d2 = patchDist[1][pointI];
+//            if (d1+d2 > SMALL)
+//            {
+//                scalar s = d1/(d1+d2);
+//                distance[pointI] = s;
+//            }
+//        }
+//    }
+//    Info<< "Writing distance pointScalarField to " << mesh().time().timeName()
+//        << endl;
+//    distance.write();
+
+    // Average
+    forAll(pointDisplacement_, pointI)
+    {
+        if (isZonePoint[pointI])
+        {
+            scalar d1 = patchDist[0][pointI];
+            scalar d2 = patchDist[1][pointI];
+
+            scalar s = d1/(d1+d2+VSMALL);
+
+            pointDisplacement_[pointI] =
+                (1-s)*patchDisp[0][pointI]
+              + s*patchDisp[1][pointI];
+        }
+    }
+    pointDisplacement_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::displacementLayeredMotionFvMotionSolver::
+displacementLayeredMotionFvMotionSolver
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    displacementFvMotionSolver(mesh, is),
+    pointDisplacement_
+    (
+        IOobject
+        (
+            "pointDisplacement",
+            fvMesh_.time().timeName(),
+            fvMesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pointMesh::New(fvMesh_)
+    )
+{
+    pointDisplacement_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::displacementLayeredMotionFvMotionSolver::
+~displacementLayeredMotionFvMotionSolver()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::pointField>
+Foam::displacementLayeredMotionFvMotionSolver::curPoints() const
+{
+    tmp<pointField> tcurPoints
+    (
+        points0() + pointDisplacement_.internalField()
+    );
+
+    twoDCorrectPoints(tcurPoints());
+
+//    const pointField& pts = tcurPoints();
+//    forAll(pts, pointI)
+//    {
+//        Info<< "    from:" << mesh().points()[pointI]
+//            << " to:" << pts[pointI]
+//            << endl;
+//    }
+
+
+    return tcurPoints;
+}
+
+
+void Foam::displacementLayeredMotionFvMotionSolver::solve()
+{
+    const dictionary& ms = mesh().lookupObject<motionSolver>("dynamicMeshDict");
+    const dictionary& solverDict = ms.subDict(typeName + "Coeffs");
+
+    // Apply all regions (=cellZones)
+
+    const dictionary& regionDicts = solverDict.subDict("regions");
+    forAllConstIter(dictionary, regionDicts, regionIter)
+    {
+        const word& cellZoneName = regionIter().keyword();
+        const dictionary& regionDict = regionIter().dict();
+
+        label zoneI = mesh().cellZones().findZoneID(cellZoneName);
+
+        Info<< "solve : zone:" << cellZoneName << " index:" << zoneI
+            << endl;
+
+        if (zoneI == -1)
+        {
+            FatalIOErrorIn
+            (
+                "displacementLayeredMotionFvMotionSolver::solve(..)",
+                *this
+            )   << "Cannot find cellZone " << cellZoneName
+                << endl << "Valid zones are " << mesh().cellZones().names()
+                << exit(FatalIOError);
+        }
+
+        cellZoneSolve(zoneI, regionDict);
+    }
+}
+
+
+void Foam::displacementLayeredMotionFvMotionSolver::updateMesh
+(
+    const mapPolyMesh& mpm
+)
+{
+    displacementFvMotionSolver::updateMesh(mpm);
+}
+
+
+// ************************************************************************* //

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.H

@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::displacementLayeredMotionFvMotionSolver
+
+Description
+    Mesh motion solver for an (multi-block) extruded fvMesh. Gets given the
+    structure of the mesh blocks and boundary conditions on these blocks.
+
+    Note: should not be an fvMotionSolver but just a motionSolver. Only here
+    so we can reuse displacementFvMotionSolver functionality (e.g. surface
+    following boundary conditions)
+
+
+    The displacementLayeredMotionCoeffs subdict of dynamicMeshDict specifies
+    per region (=cellZone) the boundary conditions on two opposing patches
+    (=faceZones). It then interpolates the boundary values using topological
+    walking so requires the cellZone to be a layered mesh.
+
+    The boundary conditions on faceZones are currently:
+
+    follow: the faceZone follows the overall mesh displacement.
+            Use this for faceZones on boundary faces (so it uses the
+            proper boundary conditions on the pointDisplacement).
+
+    fixedValue: fixed value.
+
+    timeVaryingUniformFixedValue: table-driven fixed value.
+
+    slip: the second of a pair of faceZones follows the tangential movement
+          specified by the first faceZone. (= removes the normal component).
+
+SourceFiles
+    displacementLayeredMotionFvMotionSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef displacementLayeredMotionFvMotionSolver_H
+#define displacementLayeredMotionFvMotionSolver_H
+
+#include "displacementFvMotionSolver.H"
+#include "PackedBoolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward class declarations
+
+/*---------------------------------------------------------------------------*\
+             Class displacementLayeredMotionFvMotionSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class displacementLayeredMotionFvMotionSolver
+:
+    public displacementFvMotionSolver
+{
+    // Private data
+
+        //- Point motion field
+        mutable pointVectorField pointDisplacement_;
+
+
+    // Private Member Functions
+
+        void calcZoneMask
+        (
+            const label cellZoneI, 
+            PackedBoolList& isZonePoint,
+            PackedBoolList& isZoneEdge
+        ) const;
+
+        void walkStructured
+        (
+            const label cellZoneI,
+            const PackedBoolList& isZonePoint,
+            const PackedBoolList& isZoneEdge,
+            const labelList& seedPoints,
+            const vectorField& seedData,
+            scalarField& distance,
+            vectorField& data
+        ) const;
+
+        tmp<vectorField> faceZoneEvaluate
+        (
+            const faceZone& fz,
+            const labelList& meshPoints,
+            const dictionary& dict,
+            const PtrList<pointVectorField>& patchDisp,
+            const label patchI
+        ) const;
+
+        void cellZoneSolve
+        (
+            const label cellZoneI,
+            const dictionary& zoneDict
+        );
+
+
+        //- Disallow default bitwise copy construct
+        displacementLayeredMotionFvMotionSolver
+        (
+            const displacementLayeredMotionFvMotionSolver&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const displacementLayeredMotionFvMotionSolver&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("displacementLayeredMotion");
+
+
+    // Constructors
+
+        //- Construct from polyMesh and data stream
+        displacementLayeredMotionFvMotionSolver
+        (
+            const polyMesh&,
+            Istream& msDataUnused
+        );
+
+
+    // Destructor
+
+        ~displacementLayeredMotionFvMotionSolver();
+
+
+    // Member Functions
+
+        //- Return reference to the point motion displacement field
+        pointVectorField& pointDisplacement()
+        {
+            return pointDisplacement_;
+        }
+
+        //- Return const reference to the point motion displacement field
+        const pointVectorField& pointDisplacement() const
+        {
+            return pointDisplacement_;
+        }
+
+        //- Return point location obtained from the current motion field
+        virtual tmp<pointField> curPoints() const;
+
+        //- Solve for motion
+        virtual void solve();
+
+        //- Update topology
+        virtual void updateMesh(const mapPolyMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -24,15 +24,31 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "LimitedScheme.H"
-#include "MC.H"
+#include "pointEdgeStructuredWalk.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
-namespace Foam
+Foam::Ostream& Foam::operator<<
+(
+    Foam::Ostream& os,
+    const Foam::pointEdgeStructuredWalk& wDist
+)
 {
-    makeLimitedSurfaceInterpolationScheme(MC, MCLimiter)
-    makeLimitedVSurfaceInterpolationScheme(MCV, MCLimiter)
+    return os
+        << wDist.inZone_ << wDist.previousPoint_
+        << wDist.dist_ << wDist.data_;
 }
 
+Foam::Istream& Foam::operator>>
+(
+    Foam::Istream& is,
+    Foam::pointEdgeStructuredWalk& wDist
+)
+{
+    return is
+        >> wDist.inZone_ >> wDist.previousPoint_
+        >> wDist.dist_ >> wDist.data_;
+}
+
+
 // ************************************************************************* //

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.H

@@ -0,0 +1,224 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::pointEdgeStructuredWalk
+
+Description
+    Determines length of string of edges walked to point.
+
+SourceFiles
+    pointEdgeStructuredWalkI.H
+    pointEdgeStructuredWalk.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointEdgeStructuredWalk_H
+#define pointEdgeStructuredWalk_H
+
+#include "point.H"
+#include "label.H"
+#include "scalar.H"
+#include "tensor.H"
+#include "pTraits.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyPatch;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class pointEdgeStructuredWalk Declaration
+\*---------------------------------------------------------------------------*/
+
+class pointEdgeStructuredWalk
+{
+    // Private data
+
+        //- Is this point/edge inside the walk zone
+        bool inZone_;
+
+        //- Previuos point
+        point previousPoint_;
+
+        //- Sum of distance
+        scalar dist_;
+
+        //- Passive data
+        vector data_;
+
+    // Private Member Functions
+
+        //- Evaluate distance to point.
+        inline bool update
+        (
+            const point&,
+            const pointEdgeStructuredWalk& w2,
+            const scalar tol
+        );
+
+        //- Combine current with w2. Update distSqr, origin if w2 has smaller
+        //  quantities and returns true.
+        inline bool update
+        (
+            const pointEdgeStructuredWalk& w2,
+            const scalar tol
+        );
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline pointEdgeStructuredWalk();
+
+        //- Construct from components
+        inline pointEdgeStructuredWalk
+        (
+            const bool,
+            const point&,
+            const scalar,
+            const vector&
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            inline bool inZone() const;
+
+            inline bool& inZone();
+
+            inline const point& previousPoint() const;
+
+            inline scalar dist() const;
+
+            inline const vector& data() const;
+
+
+        // Needed by meshWave
+
+            //- Check whether still contains original (invalid) value.
+            inline bool valid() const;
+
+            //- Check for identical geometrical data. Used for cyclics checking.
+            inline bool sameGeometry
+            (
+                const pointEdgeStructuredWalk&,
+                const scalar tol
+            ) const;
+
+            //- Convert origin to relative vector to leaving point
+            //  (= point coordinate)
+            inline void leaveDomain
+            (
+                const polyPatch& patch,
+                const label patchPointI,
+                const point& pos
+            );
+
+            //- Convert relative origin to absolute by adding entering point
+            inline void enterDomain
+            (
+                const polyPatch& patch,
+                const label patchPointI,
+                const point& pos
+            );
+
+            //- Apply rotation matrix to origin
+            inline void transform(const tensor& rotTensor);
+
+            //- Influence of edge on point
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const label edgeI,
+                const pointEdgeStructuredWalk& edgeInfo,
+                const scalar tol
+            );
+
+            //- Influence of different value on same point.
+            //  Merge new and old info.
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const pointEdgeStructuredWalk& newPointInfo,
+                const scalar tol
+            );
+
+            //- Influence of different value on same point.
+            //  No information about current position whatsoever.
+            inline bool updatePoint
+            (
+                const pointEdgeStructuredWalk& newPointInfo,
+                const scalar tol
+            );
+
+            //- Influence of point on edge.
+            inline bool updateEdge
+            (
+                const polyMesh& mesh,
+                const label edgeI,
+                const label pointI,
+                const pointEdgeStructuredWalk& pointInfo,
+                const scalar tol
+            );
+
+
+    // Member Operators
+
+        //Note: Used to determine whether to call update.
+        inline bool operator==(const pointEdgeStructuredWalk&) const;
+
+        inline bool operator!=(const pointEdgeStructuredWalk&) const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const pointEdgeStructuredWalk&);
+        friend Istream& operator>>(Istream&, pointEdgeStructuredWalk&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "pointEdgeStructuredWalkI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalkI.H

@@ -0,0 +1,302 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyMesh.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Update this with w2.
+inline bool Foam::pointEdgeStructuredWalk::update
+(
+    const point& pt,
+    const pointEdgeStructuredWalk& w2,
+    const scalar tol
+)
+{
+    if (!valid())
+    {
+        // current not yet set. Walked from w2 to here (=pt)
+        dist_ = w2.dist_ + mag(pt-w2.previousPoint_);
+        previousPoint_ = pt;
+        data_ = w2.data_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+    
+// Update this with w2 (information on same point)
+inline bool Foam::pointEdgeStructuredWalk::update
+(
+    const pointEdgeStructuredWalk& w2,
+    const scalar tol
+)
+{
+    if (!valid())
+    {
+        // current not yet set so use neighbour
+        dist_ = w2.dist_;
+        previousPoint_ = w2.previousPoint_;
+        data_ = w2.data_;
+
+        return true;
+    }
+    else
+    {
+        // already nearer to pt
+        return false;
+    }
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Null constructor
+inline Foam::pointEdgeStructuredWalk::pointEdgeStructuredWalk()
+:
+    inZone_(false),
+    previousPoint_(vector::max),
+    dist_(0),
+    data_(vector::zero)
+{}
+
+
+// Construct from origin, distance
+inline Foam::pointEdgeStructuredWalk::pointEdgeStructuredWalk
+(
+    const bool inZone,
+    const point& previousPoint, 
+    const scalar dist,
+    const vector& data
+)
+:
+    inZone_(inZone),
+    previousPoint_(previousPoint),
+    dist_(dist),
+    data_(data)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline bool Foam::pointEdgeStructuredWalk::inZone() const
+{
+    return inZone_;
+}
+
+
+inline bool& Foam::pointEdgeStructuredWalk::inZone()
+{
+    return inZone_;
+}
+
+
+inline const Foam::point& Foam::pointEdgeStructuredWalk::previousPoint() const
+{
+    return previousPoint_;
+}
+
+
+inline Foam::scalar Foam::pointEdgeStructuredWalk::dist() const
+{
+    return dist_;
+}
+
+
+inline const Foam::vector& Foam::pointEdgeStructuredWalk::data() const
+{
+    return data_;
+}
+
+
+inline bool Foam::pointEdgeStructuredWalk::valid() const
+{
+    return previousPoint_ != vector::max;
+}
+
+
+// Checks for cyclic points
+inline bool Foam::pointEdgeStructuredWalk::sameGeometry
+(
+    const pointEdgeStructuredWalk& w2,
+    const scalar tol
+) const
+{
+    scalar diff = Foam::mag(dist() - w2.dist());
+
+    if (diff < SMALL)
+    {
+        return true;
+    }
+    else
+    {
+        if ((dist() > SMALL) && ((diff/dist()) < tol))
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+}
+
+
+inline void Foam::pointEdgeStructuredWalk::leaveDomain
+(
+    const polyPatch& patch,
+    const label patchPointI,
+    const point& coord
+)
+{
+    previousPoint_ -= coord;
+}
+
+
+inline void Foam::pointEdgeStructuredWalk::transform(const tensor& rotTensor)
+{
+    previousPoint_ = Foam::transform(rotTensor, previousPoint_);
+}
+
+
+// Update absolute geometric quantities. Note that distance (dist_)
+// is not affected by leaving/entering domain.
+inline void Foam::pointEdgeStructuredWalk::enterDomain
+(
+    const polyPatch& patch,
+    const label patchPointI,
+    const point& coord
+)
+{
+    // back to absolute form
+    previousPoint_ += coord;
+}
+
+
+// Update this with information from connected edge
+inline bool Foam::pointEdgeStructuredWalk::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const label edgeI,
+    const pointEdgeStructuredWalk& edgeInfo,
+    const scalar tol
+)
+{
+    if (inZone_)
+    {
+        return update(mesh.points()[pointI], edgeInfo, tol);
+    }
+    else
+    {
+        return false;
+    }
+}    
+
+
+// Update this with new information on same point
+inline bool Foam::pointEdgeStructuredWalk::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const pointEdgeStructuredWalk& newPointInfo,
+    const scalar tol
+)
+{
+    if (inZone_)
+    {
+        return update(mesh.points()[pointI], newPointInfo, tol);
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// Update this with new information on same point. No extra information.
+inline bool Foam::pointEdgeStructuredWalk::updatePoint
+(
+    const pointEdgeStructuredWalk& newPointInfo,
+    const scalar tol
+)
+{
+    return update(newPointInfo, tol);
+}    
+
+
+// Update this with information from connected point
+inline bool Foam::pointEdgeStructuredWalk::updateEdge
+(
+    const polyMesh& mesh,
+    const label edgeI,
+    const label pointI,
+    const pointEdgeStructuredWalk& pointInfo,
+    const scalar tol
+)
+{
+    if (inZone_)
+    {
+        const pointField& points = mesh.points();
+
+        const edge& e = mesh.edges()[edgeI];
+
+        const point edgeMid(0.5*(points[e[0]] + points[e[1]]));
+
+        return update(edgeMid, pointInfo, tol);
+    }
+    else
+    {
+        return false;
+    }
+}    
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::pointEdgeStructuredWalk::operator==
+(
+    const Foam::pointEdgeStructuredWalk& rhs
+) const
+{
+    return previousPoint() == rhs.previousPoint();
+}
+
+
+inline bool Foam::pointEdgeStructuredWalk::operator!=
+(
+    const Foam::pointEdgeStructuredWalk& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //

src/sampling/sampledSurface/distanceSurface/distanceSurface.C

@@ -30,8 +30,6 @@ License
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
-#include "isoSurface.H"
-// #include "isoSurfaceCell.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -244,6 +242,14 @@ void Foam::distanceSurface::createGeometry()
             distance_,
             regularise_
         )
+        //new isoSurfaceCell
+        //(
+        //    fvm,
+        //    cellDistance,
+        //    pointDistance_,
+        //    distance_,
+        //    regularise_
+        //)
     );
 
     if (debug)
@@ -284,6 +290,7 @@ Foam::distanceSurface::distanceSurface
     distance_(readScalar(dict.lookup("distance"))),
     signed_(readBool(dict.lookup("signed"))),
     regularise_(dict.lookupOrDefault("regularise", true)),
+    average_(dict.lookupOrDefault("average", false)),
     zoneName_(word::null),
     needsUpdate_(true),
     isoSurfPtr_(NULL),

src/sampling/sampledSurface/distanceSurface/distanceSurface.H

@@ -68,6 +68,9 @@ class distanceSurface
         //- Whether to coarsen
         const Switch regularise_;
 
+        //- Whether to recalculate cell values as average of point values
+        const Switch average_;
+
         //- zone name (if restricted to zones)
         word zoneName_;
 
@@ -82,6 +85,7 @@ class distanceSurface
         scalarField pointDistance_;
 
         //- Constructed iso surface
+        //autoPtr<isoSurfaceCell> isoSurfPtr_;
         autoPtr<isoSurface> isoSurfPtr_;
 
         //- triangles converted to faceList
@@ -166,6 +170,7 @@ public:
         }
 
 
+        //const isoSurfaceCell& surface() const
         const isoSurface& surface() const
         {
             return isoSurfPtr_();

src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C

@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "distanceSurface.H"
-#include "isoSurface.H"
 #include "volFieldsFwd.H"
 #include "pointFields.H"
 #include "volPointInterpolation.H"
@@ -62,7 +61,15 @@ Foam::distanceSurface::interpolateField
     );
 
     // Sample.
-    return surface().interpolate(volFld, pointFld());
+    return surface().interpolate
+    (
+        (
+            average_
+          ? pointAverage(pointFld())()
+          : volFld
+        ),
+        pointFld()
+    );
 }
 
 

src/sampling/sampledSurface/isoSurface/isoSurface.C

@@ -2054,61 +2054,68 @@ Foam::isoSurface::isoSurface
     }
 
 
-    //if (false)
-    //{
-    List<FixedList<label, 3> > faceEdges;
-    labelList edgeFace0, edgeFace1;
-    Map<labelList> edgeFacesRest;
-
-
-    while (true)
+    if (false)
     {
-        // Calculate addressing
-        calcAddressing(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
+        List<FixedList<label, 3> > faceEdges;
+        labelList edgeFace0, edgeFace1;
+        Map<labelList> edgeFacesRest;
 
-        // See if any dangling triangles
-        boolList keepTriangles;
-        label nDangling = markDanglingTriangles
-        (
-            faceEdges,
-            edgeFace0,
-            edgeFace1,
-            edgeFacesRest,
-            keepTriangles
-        );
 
-        if (debug)
+        while (true)
         {
-            Pout<< "isoSurface : detected " << nDangling
-                << " dangling triangles." << endl;
-        }
-
-        if (nDangling == 0)
-        {
-            break;
-        }
-
-        // Create face map (new to old)
-        labelList subsetTriMap(findIndices(keepTriangles, true));
-
-        labelList subsetPointMap;
-        labelList reversePointMap;
-        triSurface::operator=
-        (
-            subsetMesh
+            // Calculate addressing
+            calcAddressing
             (
                 *this,
-                subsetTriMap,
-                reversePointMap,
-                subsetPointMap
-            )
-        );
-        meshCells_ = labelField(meshCells_, subsetTriMap);
-        inplaceRenumber(reversePointMap, triPointMergeMap_);
-    }
+                faceEdges,
+                edgeFace0,
+                edgeFace1,
+                edgeFacesRest
+            );
 
-    orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
-    //}
+            // See if any dangling triangles
+            boolList keepTriangles;
+            label nDangling = markDanglingTriangles
+            (
+                faceEdges,
+                edgeFace0,
+                edgeFace1,
+                edgeFacesRest,
+                keepTriangles
+            );
+
+            if (debug)
+            {
+                Pout<< "isoSurface : detected " << nDangling
+                    << " dangling triangles." << endl;
+            }
+
+            if (nDangling == 0)
+            {
+                break;
+            }
+
+            // Create face map (new to old)
+            labelList subsetTriMap(findIndices(keepTriangles, true));
+
+            labelList subsetPointMap;
+            labelList reversePointMap;
+            triSurface::operator=
+            (
+                subsetMesh
+                (
+                    *this,
+                    subsetTriMap,
+                    reversePointMap,
+                    subsetPointMap
+                )
+            );
+            meshCells_ = labelField(meshCells_, subsetTriMap);
+            inplaceRenumber(reversePointMap, triPointMergeMap_);
+        }
+
+        orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
+    }
 
 
     if (debug)

src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C

@@ -1398,6 +1398,8 @@ Foam::isoSurfaceCell::isoSurfaceCell
 )
 :
     mesh_(mesh),
+    cVals_(cVals),
+    pVals_(pVals),
     iso_(iso),
     mergeDistance_(mergeTol*mesh.bounds().mag())
 {
@@ -1563,59 +1565,68 @@ Foam::isoSurfaceCell::isoSurfaceCell
     }
 
 
-    List<FixedList<label, 3> > faceEdges;
-    labelList edgeFace0, edgeFace1;
-    Map<labelList> edgeFacesRest;
-
-
-    while (true)
+    if (false)
     {
-        // Calculate addressing
-        calcAddressing(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
+        List<FixedList<label, 3> > faceEdges;
+        labelList edgeFace0, edgeFace1;
+        Map<labelList> edgeFacesRest;
 
-        // See if any dangling triangles
-        boolList keepTriangles;
-        label nDangling = markDanglingTriangles
-        (
-            faceEdges,
-            edgeFace0,
-            edgeFace1,
-            edgeFacesRest,
-            keepTriangles
-        );
 
-        if (debug)
+        while (true)
         {
-            Pout<< "isoSurfaceCell : detected " << nDangling
-                << " dangling triangles." << endl;
-        }
-
-        if (nDangling == 0)
-        {
-            break;
-        }
-
-        // Create face map (new to old)
-        labelList subsetTriMap(findIndices(keepTriangles, true));
-
-        labelList subsetPointMap;
-        labelList reversePointMap;
-        triSurface::operator=
-        (
-            subsetMesh
+            // Calculate addressing
+            calcAddressing
             (
                 *this,
-                subsetTriMap,
-                reversePointMap,
-                subsetPointMap
-            )
-        );
-        meshCells_ = labelField(meshCells_, subsetTriMap);
-        inplaceRenumber(reversePointMap, triPointMergeMap_);
+                faceEdges,
+                edgeFace0,
+                edgeFace1,
+                edgeFacesRest
+            );
+
+            // See if any dangling triangles
+            boolList keepTriangles;
+            label nDangling = markDanglingTriangles
+            (
+                faceEdges,
+                edgeFace0,
+                edgeFace1,
+                edgeFacesRest,
+                keepTriangles
+            );
+
+            if (debug)
+            {
+                Pout<< "isoSurfaceCell : detected " << nDangling
+                    << " dangling triangles." << endl;
+            }
+
+            if (nDangling == 0)
+            {
+                break;
+            }
+
+            // Create face map (new to old)
+            labelList subsetTriMap(findIndices(keepTriangles, true));
+
+            labelList subsetPointMap;
+            labelList reversePointMap;
+            triSurface::operator=
+            (
+                subsetMesh
+                (
+                    *this,
+                    subsetTriMap,
+                    reversePointMap,
+                    subsetPointMap
+                )
+            );
+            meshCells_ = labelField(meshCells_, subsetTriMap);
+            inplaceRenumber(reversePointMap, triPointMergeMap_);
+        }
+
+        orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
     }
-
-    orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
-
     //combineCellTriangles();
 }
 

src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H

@@ -85,6 +85,10 @@ class isoSurfaceCell
         //- Reference to mesh
         const polyMesh& mesh_;
 
+        const scalarField& cVals_;
+
+        const scalarField& pVals_;
+
         //- isoSurfaceCell value
         const scalar iso_;
 
@@ -352,6 +356,14 @@ public:
             const Field<Type>& cCoords,
             const Field<Type>& pCoords
         ) const;
+
+        //- Interpolates cCoords,pCoords.
+        template <class Type>
+        tmp<Field<Type> > interpolate
+        (
+            const Field<Type>& cCoords,
+            const Field<Type>& pCoords
+        ) const;
 };
 
 

src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C

@@ -381,4 +381,56 @@ Foam::isoSurfaceCell::interpolate
 }
 
 
+template <class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::isoSurfaceCell::interpolate
+(
+    const Field<Type>& cCoords,
+    const Field<Type>& pCoords
+) const
+{
+    DynamicList<Type> triPoints(nCutCells_);
+    DynamicList<label> triMeshCells(nCutCells_);
+
+    // Dummy snap data
+    DynamicList<Type> snappedPoints;
+    labelList snappedCc(mesh_.nCells(), -1);
+    labelList snappedPoint(mesh_.nPoints(), -1);
+
+
+    generateTriPoints
+    (
+        cVals_,
+        pVals_,
+
+        cCoords,
+        pCoords,
+
+        snappedPoints,
+        snappedCc,
+        snappedPoint,
+
+        triPoints,
+        triMeshCells
+    );
+
+
+    // One value per point
+    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
+    Field<Type>& values = tvalues();
+
+    forAll(triPoints, i)
+    {
+        label mergedPointI = triPointMergeMap_[i];
+
+        if (mergedPointI >= 0)
+        {
+            values[mergedPointI] = triPoints[i];
+        }
+    }
+
+    return tvalues;
+}
+
+
 // ************************************************************************* //

src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C

@@ -164,7 +164,10 @@ void Foam::sampledIsoSurface::getIsoFields() const
         // point field.
         if (average_)
         {
-            storedVolFieldPtr_.reset(average(fvm, *pointFieldPtr_).ptr());
+            storedVolFieldPtr_.reset
+            (
+                pointAverage(*pointFieldPtr_).ptr()
+            );
             volFieldPtr_ = storedVolFieldPtr_.operator->();
         }
 
@@ -265,7 +268,7 @@ void Foam::sampledIsoSurface::getIsoFields() const
         {
             storedVolSubFieldPtr_.reset
             (
-                average(subFvm, *pointSubFieldPtr_).ptr()
+                pointAverage(*pointSubFieldPtr_).ptr()
             );
             volSubFieldPtr_ = storedVolSubFieldPtr_.operator->();
         }
@@ -286,99 +289,6 @@ void Foam::sampledIsoSurface::getIsoFields() const
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::sampledIsoSurface::average
-(
-    const fvMesh& mesh,
-    const pointScalarField& pfld
-) const
-{
-    tmp<volScalarField> tcellAvg
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "cellAvg",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh,
-            dimensionedScalar("zero", dimless, scalar(0.0))
-        )
-    );
-    volScalarField& cellAvg = tcellAvg();
-
-    labelField nPointCells(mesh.nCells(), 0);
-    {
-        for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
-        {
-            const labelList& pCells = mesh.pointCells(pointI);
-
-            forAll(pCells, i)
-            {
-                label cellI = pCells[i];
-
-                cellAvg[cellI] += pfld[pointI];
-                nPointCells[cellI]++;
-            }
-        }
-    }
-    forAll(cellAvg, cellI)
-    {
-        cellAvg[cellI] /= nPointCells[cellI];
-    }
-    // Give value to calculatedFvPatchFields
-    cellAvg.correctBoundaryConditions();
-
-    return tcellAvg;
-}
-
-
-Foam::tmp<Foam::pointScalarField> Foam::sampledIsoSurface::average
-(
-    const pointMesh& pMesh,
-    const volScalarField& fld
-) const
-{
-    tmp<pointScalarField> tpointAvg
-    (
-        new pointScalarField
-        (
-            IOobject
-            (
-                "pointAvg",
-                fld.time().timeName(),
-                fld.db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            pMesh,
-            dimensionedScalar("zero", dimless, scalar(0.0))
-        )
-    );
-    pointScalarField& pointAvg = tpointAvg();
-
-    for (label pointI = 0; pointI < fld.mesh().nPoints(); pointI++)
-    {
-        const labelList& pCells = fld.mesh().pointCells(pointI);
-
-        forAll(pCells, i)
-        {
-            pointAvg[pointI] += fld[pCells[i]];
-        }
-        pointAvg[pointI] /= pCells.size();
-    }
-    // Give value to calculatedFvPatchFields
-    pointAvg.correctBoundaryConditions();
-
-    return tpointAvg;
-}
-
-
 bool Foam::sampledIsoSurface::updateGeometry() const
 {
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());

src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H

@@ -118,18 +118,6 @@ class sampledIsoSurface
         //- Get fields needed to recreate iso surface.
         void getIsoFields() const;
 
-        tmp<volScalarField> average
-        (
-            const fvMesh&,
-            const pointScalarField&
-        ) const;
-
-        tmp<pointScalarField> average
-        (
-            const pointMesh&,
-            const volScalarField& fld
-        ) const;
-
         //- Create iso surface (if time has changed)
         //  Do nothing (and return false) if no update was needed
         bool updateGeometry() const;

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C

@@ -71,7 +71,15 @@ Foam::sampledIsoSurface::interpolateField
             volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
 
         // Sample.
-        return surface().interpolate(volSubFld, tpointSubFld());
+        return surface().interpolate
+        (
+            (
+                average_
+              ? pointAverage(tpointSubFld())()
+              : volSubFld
+            ),
+            tpointSubFld()
+        );
     }
     else
     {
@@ -79,7 +87,15 @@ Foam::sampledIsoSurface::interpolateField
             volPointInterpolation::New(volFld.mesh()).interpolate(volFld);
 
         // Sample.
-        return surface().interpolate(volFld, tpointFld());
+        return surface().interpolate
+        (
+            (
+                average_
+              ? pointAverage(tpointFld())()
+              : volFld
+            ),
+            tpointFld()
+        );
     }
 }
 

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C

@@ -30,7 +30,6 @@ License
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
-#include "isoSurface.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -173,7 +172,7 @@ void Foam::sampledCuttingPlane::createGeometry()
 
                 forAll(fld, i)
                 {
-                    fld[i] =  (cc[i] - plane_.refPoint()) & plane_.normal();
+                    fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal();
                 }
             }
         }
@@ -231,6 +230,14 @@ void Foam::sampledCuttingPlane::createGeometry()
             0.0,
             regularise_
         )
+        //new isoSurfaceCell
+        //(
+        //    fvm,
+        //    cellDistance,
+        //    pointDistance_,
+        //    0.0,
+        //    regularise_
+        //)
     );
 
     if (debug)
@@ -254,6 +261,7 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
     plane_(dict),
     mergeTol_(dict.lookupOrDefault("mergeTol", 1E-6)),
     regularise_(dict.lookupOrDefault("regularise", true)),
+    average_(dict.lookupOrDefault("average", false)),
     zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
     exposedPatchName_(word::null),
     needsUpdate_(true),

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H

@@ -38,6 +38,7 @@ SourceFiles
 
 #include "sampledSurface.H"
 #include "isoSurface.H"
+//#include "isoSurfaceCell.H"
 #include "plane.H"
 #include "ZoneIDs.H"
 #include "fvMeshSubset.H"
@@ -66,6 +67,9 @@ class sampledCuttingPlane
         //- Whether to coarsen
         const Switch regularise_;
 
+        //- Whether to recalculate cell values as average of point values
+        const Switch average_;
+
         //- zone name/index (if restricted to zones)
         mutable cellZoneID zoneID_;
 
@@ -86,6 +90,7 @@ class sampledCuttingPlane
         scalarField pointDistance_;
 
         //- Constructed iso surface
+        //autoPtr<isoSurfaceCell> isoSurfPtr_;
         autoPtr<isoSurface> isoSurfPtr_;
 
         //- triangles converted to faceList
@@ -170,6 +175,7 @@ public:
         }
 
 
+        //const isoSurfaceCell& surface() const
         const isoSurface& surface() const
         {
             return isoSurfPtr_();

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C

@@ -26,7 +26,6 @@ License
 
 #include "volPointInterpolation.H"
 #include "sampledCuttingPlane.H"
-#include "isoSurface.H"
 #include "volFieldsFwd.H"
 #include "pointFields.H"
 #include "volPointInterpolation.H"
@@ -67,7 +66,15 @@ Foam::sampledCuttingPlane::interpolateField
             volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
 
         // Sample.
-        return surface().interpolate(volSubFld, tpointSubFld());
+        return surface().interpolate
+        (
+            (
+                average_
+              ? pointAverage(tpointSubFld())()
+              : volSubFld
+            ),
+            tpointSubFld()
+        );
     }
     else
     {
@@ -77,7 +84,15 @@ Foam::sampledCuttingPlane::interpolateField
         );
 
         // Sample.
-        return surface().interpolate(volFld, tpointFld());
+        return surface().interpolate
+        (
+            (
+                average_
+              ? pointAverage(tpointFld())()
+              : volFld
+            ),
+            tpointFld()
+        );
     }
 }
 

src/sampling/sampledSurface/sampledSurface/sampledSurface.H

@@ -312,6 +312,13 @@ public:
         //- Project field onto surface
         tmp<Field<vector> > project(const Field<tensor>&) const;
 
+        //- Interpolate from points to cell centre
+        template<class Type>
+        tmp<GeometricField<Type, fvPatchField, volMesh> > pointAverage
+        (
+            const GeometricField<Type, pointPatchField, pointMesh>& pfld
+        ) const;
+
         //- Sample field on surface
         virtual tmp<scalarField> sample
         (

src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C

@@ -155,4 +155,59 @@ Foam::sampledSurface::project
 }
 
 
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::sampledSurface::pointAverage
+(
+    const GeometricField<Type, pointPatchField, pointMesh>& pfld
+) const
+{
+    const fvMesh& mesh = dynamic_cast<const fvMesh&>(pfld.mesh()());
+
+    tmp<GeometricField<Type, fvPatchField, volMesh> > tcellAvg
+    (
+        new GeometricField<Type, fvPatchField, volMesh>
+        (
+            IOobject
+            (
+                "cellAvg",
+                mesh.time().timeName(),
+                pfld.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimensioned<Type>("zero", dimless, pTraits<Type>::zero)
+        )
+    );
+    GeometricField<Type, fvPatchField, volMesh>& cellAvg = tcellAvg();
+
+    labelField nPointCells(mesh.nCells(), 0);
+    {
+        for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
+        {
+            const labelList& pCells = mesh.pointCells(pointI);
+
+            forAll(pCells, i)
+            {
+                label cellI = pCells[i];
+
+                cellAvg[cellI] += pfld[pointI];
+                nPointCells[cellI]++;
+            }
+        }
+    }
+    forAll(cellAvg, cellI)
+    {
+        cellAvg[cellI] /= nPointCells[cellI];
+    }
+    // Give value to calculatedFvPatchFields
+    cellAvg.correctBoundaryConditions();
+
+    return tcellAvg;
+}
+
+
+
 // ************************************************************************* //

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,11 +25,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sampledTriSurfaceMesh.H"
-#include "dictionary.H"
-#include "polyMesh.H"
-#include "polyPatch.H"
-#include "volFields.H"
+#include "treeDataPoint.H"
 #include "meshSearch.H"
+#include "Tuple2.H"
+#include "globalIndex.H"
 
 #include "addToRunTimeSelectionTable.H"
 
@@ -44,6 +43,25 @@ namespace Foam
         sampledTriSurfaceMesh,
         word
     );
+
+    //- Private class for finding nearest
+    //  - global index
+    //  - sqr(distance)
+    typedef Tuple2<scalar, label> nearInfo;
+
+    class nearestEqOp
+    {
+
+    public:
+
+        void operator()(nearInfo& x, const nearInfo& y) const
+        {
+            if (y.first() < x.first())
+            {
+                x = y;
+            }
+        }
+    };
 }
 
 
@@ -63,16 +81,16 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
         (
             name,
             mesh.time().constant(), // instance
-            "triSurface",       // local
-            mesh,               // registry
+            "triSurface",           // local
+            mesh,                   // registry
             IOobject::MUST_READ,
-            IOobject::NO_WRITE
+            IOobject::NO_WRITE,
+            false
         )
     ),
     needsUpdate_(true),
     cellLabels_(0),
-    pointMap_(0),
-    faceMap_(0)
+    pointToFace_(0)
 {}
 
 
@@ -90,16 +108,16 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
         (
             dict.lookup("surface"),
             mesh.time().constant(), // instance
-            "triSurface",       // local
-            mesh,               // registry
+            "triSurface",           // local
+            mesh,                   // registry
             IOobject::MUST_READ,
-            IOobject::NO_WRITE
+            IOobject::NO_WRITE,
+            false
         )
     ),
     needsUpdate_(true),
     cellLabels_(0),
-    pointMap_(0),
-    faceMap_(0)
+    pointToFace_(0)
 {}
 
 
@@ -128,8 +146,7 @@ bool Foam::sampledTriSurfaceMesh::expire()
     sampledSurface::clearGeom();
     MeshStorage::clear();
     cellLabels_.clear();
-    pointMap_.clear();
-    faceMap_.clear();
+    pointToFace_.clear();
 
     needsUpdate_ = true;
     return true;
@@ -146,95 +163,139 @@ bool Foam::sampledTriSurfaceMesh::update()
 
     // Find the cells the triangles of the surface are in.
     // Does approximation by looking at the face centres only
-    const pointField fc = surface_.faceCentres();
-
-    cellLabels_.setSize(fc.size());
+    const pointField& fc = surface_.faceCentres();
 
     meshSearch meshSearcher(mesh(), false);
 
+    const indexedOctree<treeDataPoint>& cellCentreTree =
+        meshSearcher.cellCentreTree();
+
+
+    // Global numbering for cells - only used to uniquely identify local cells.
+    globalIndex globalCells(mesh().nCells());
+    List<nearInfo> nearest(fc.size());
+    forAll(nearest, i)
+    {
+        nearest[i].first() = GREAT;
+        nearest[i].second() = labelMax;
+    }
+
+    // Search triangles using nearest on local mesh
     forAll(fc, triI)
     {
-        cellLabels_[triI] = meshSearcher.findCell
+        pointIndexHit nearInfo = cellCentreTree.findNearest
         (
             fc[triI],
-            -1,                 // seed
-            true                // use tree
+            sqr(GREAT)
         );
-    }
-
-
-    boolList include(surface_.size(), true);
-
-    label nNotFound = 0;
-
-    forAll(cellLabels_, triI)
-    {
-        if (cellLabels_[triI] == -1)
+        if (nearInfo.hit())
         {
-            include[triI] = false;
-            nNotFound++;
+            nearest[triI].first() = magSqr(nearInfo.hitPoint()-fc[triI]);
+            nearest[triI].second() = globalCells.toGlobal(nearInfo.index());
         }
     }
-    label nTotalNotFound = returnReduce(nNotFound, sumOp<label>());
+
+    // See which processor has the nearest.
+    Pstream::listCombineGather(nearest, nearestEqOp());
+    Pstream::listCombineScatter(nearest);
+
+    boolList include(surface_.size(), false);
+
+    cellLabels_.setSize(fc.size());
+    cellLabels_ = -1;
+
+    label nFound = 0;
+    forAll(nearest, triI)
+    {
+        if (nearest[triI].second() == labelMax)
+        {
+            // Not found on any processor. How to map?
+        }
+        else if (globalCells.isLocal(nearest[triI].second()))
+        {
+            cellLabels_[triI] = globalCells.toLocal(nearest[triI].second());
+
+            include[triI] = true;
+            nFound++;
+        }
+    }
+
 
     if (debug)
     {
-        Pout<< "surface:" << surface_.size()
-            << "  included:" << surface_.size()-nNotFound
-            << "  total not found" << nTotalNotFound << endl;
+        Pout<< "Local out of faces:" << cellLabels_.size()
+            << " keeping:" << nFound << endl;
     }
 
+    // Now subset the surface. Do not use triSurface::subsetMesh since requires
+    // original surface to be in compact numbering.
+
+    const triSurface& s = surface_;
+
+    // Compact to original triangle
+    labelList faceMap(s.size());
+    // Compact to original points
+    labelList pointMap(s.points().size());
+    // From original point to compact points
+    labelList reversePointMap(s.points().size(), -1);
 
-    // Add to master all triangles are outside all meshes.
     {
-        boolList onAnyProc(surface_.size(), false);
-        Pstream::listCombineGather(onAnyProc, orEqOp<bool>());
+        label newPointI = 0;
+        label newTriI = 0;
 
-        if (Pstream::master())
+        forAll(s, triI)
         {
-            label nAdded = 0;
-            forAll(onAnyProc, triI)
+            if (include[triI])
             {
-                if (!onAnyProc[triI])
+                faceMap[newTriI++] = triI;
+
+                const labelledTri& f = s[triI];
+                forAll(f, fp)
                 {
-                    include[triI] = true;
-                    nAdded++;
+                    if (reversePointMap[f[fp]] == -1)
+                    {
+                        pointMap[newPointI] = f[fp];
+                        reversePointMap[f[fp]] = newPointI++;
+                    }
                 }
             }
+        }
+        faceMap.setSize(newTriI);
+        pointMap.setSize(newPointI);
+    }
 
-            if (debug)
-            {
-                Pout<< "nAdded to master:" << nAdded << endl;
-            }
+    // Subset cellLabels
+    cellLabels_ = UIndirectList<label>(cellLabels_, faceMap)();
+
+    // Store any face per point
+    pointToFace_.setSize(pointMap.size());
+
+    // Create faces and points for subsetted surface
+    faceList& faces = this->storedFaces();
+    faces.setSize(faceMap.size());
+    forAll(faceMap, i)
+    {
+        const triFace& f = s[faceMap[i]];
+        triFace newF
+        (
+            reversePointMap[f[0]],
+            reversePointMap[f[1]],
+            reversePointMap[f[2]]
+        );
+        faces[i] = newF.triFaceFace();
+
+        forAll(newF, fp)
+        {
+            pointToFace_[newF[fp]] = i;
         }
     }
 
-
-    // Now subset the surface
-    triSurface localSurface
-    (
-        surface_.subsetMesh
-        (
-            include,
-            pointMap_,
-            faceMap_
-        )
-    );
-
-    // And convert into faces.
-    faceList& faces = this->storedFaces();
-    faces.setSize(localSurface.size());
-    forAll(localSurface, i)
-    {
-        faces[i] = localSurface[i].triFaceFace();
-    }
-
-    this->storedPoints() = localSurface.points();
+    this->storedPoints() = pointField(s.points(), pointMap);
 
     if (debug)
     {
         print(Pout);
-        Pout << endl;
+        Pout<< endl;
     }
 
     needsUpdate_ = false;

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,15 +29,13 @@ Description
     A sampledSurface from a triSurfaceMesh. It samples on the points/triangles
     of the triSurface.
 
+    It samples using the cell nearest to the triangle centre so does
+    not check the cell the centre is actually in ...
+
     In parallel every processor just operates on the part of the surface
     where the face centres are inside the mesh. It is then up to the
     caller to stitch the partial surfaces together.
 
-    No check is done for triangle centres being on multiple processors
-    (should never happen but truncation errors ...)
-
-    Any value on a triangle outside the mesh will be set to 0.
-
 SourceFiles
     sampledTriSurfaceMesh.C
 
@@ -75,14 +73,11 @@ class sampledTriSurfaceMesh
         //- Track if the surface needs an update
         mutable bool needsUpdate_;
 
-        //- Local cell labels
+        //- From local surface triangle to mesh cell.
         labelList cellLabels_;
 
         //- From local surface back to surface_
-        labelList pointMap_;
-
-        //- From local surface back to surface_
-        labelList faceMap_;
+        labelList pointToFace_;
 
 
     // Private Member Functions

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C

@@ -36,20 +36,12 @@ Foam::sampledTriSurfaceMesh::sampleField
 ) const
 {
     // One value per face
-    tmp<Field<Type> > tvalues(new Field<Type>(faceMap_.size()));
+    tmp<Field<Type> > tvalues(new Field<Type>(cellLabels_.size()));
     Field<Type>& values = tvalues();
 
-    forAll(faceMap_, i)
+    forAll(cellLabels_, triI)
     {
-        label cellI = cellLabels_[faceMap_[i]];
-        if (cellI != -1)
-        {
-            values[i] = vField[cellI];
-        }
-        else
-        {
-            values[i] = pTraits<Type>::zero;
-        }
+        values[triI] = vField[cellLabels_[triI]];
     }
 
     return tvalues;
@@ -64,24 +56,15 @@ Foam::sampledTriSurfaceMesh::interpolateField
 ) const
 {
     // One value per vertex
-    tmp<Field<Type> > tvalues(new Field<Type>(pointMap_.size()));
+    tmp<Field<Type> > tvalues(new Field<Type>(pointToFace_.size()));
     Field<Type>& values = tvalues();
 
-    forAll(pointMap_, i)
+    forAll(pointToFace_, pointI)
     {
-        label origPointI = pointMap_[i];
-        label origTriI = surface_.pointFaces()[origPointI][0];
+        label triI = pointToFace_[pointI];
+        label cellI = cellLabels_[triI];
 
-        label cellI = cellLabels_[origTriI];
-
-        if (cellI != -1)
-        {
-            values[i] = interpolator.interpolate(points()[i], cellI);
-        }
-        else
-        {
-            values[i] = pTraits<Type>::zero;
-        }
+        values[pointI] = interpolator.interpolate(points()[pointI], cellI);
     }
 
     return tvalues;

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/0/U

@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet           
+    {
+        type            fixedValue;
+        value           uniform (10 0 0);
+    }
+
+    outlet          
+    {
+        //type            zeroGradient;
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
+    }
+
+    upperWall       
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    lowerWall       
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    frontAndBack    
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/0/Ua

@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      Ua;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet           
+    {
+        type            fixedValue;
+        value           uniform (-1 0 0);
+    }
+
+    outlet          
+    {
+        //type            zeroGradient;
+        //type            inletOutlet;
+        //inletValue      uniform (0 0 0);
+        type            adjointOutletVelocity;
+        value           uniform (0 0 0);
+    }
+
+    upperWall       
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    lowerWall       
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    frontAndBack    
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/0/epsilon

@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 14.855;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 14.855;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    upperWall
+    {
+        type            epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 14.855;
+    }
+    lowerWall
+    {
+        type            epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 14.855;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/0/k

@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.375;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.375;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    upperWall
+    {
+        type            kqRWallFunction;
+        value           uniform 0.375;
+    }
+    lowerWall
+    {
+        type            kqRWallFunction;
+        value           uniform 0.375;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/0/nut

@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    upperWall
+    {
+        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+    lowerWall
+    {
+        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/0/p

@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet           
+    {
+        type            zeroGradient;
+    }
+
+    outlet          
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    upperWall       
+    {
+        type            zeroGradient;
+    }
+
+    lowerWall       
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack    
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/0/pa

@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      pa;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet           
+    {
+        type            zeroGradient;
+    }
+
+    outlet          
+    {
+        type            adjointOutletPressure;
+        value           uniform 0;
+    }
+
+    upperWall       
+    {
+        type            zeroGradient;
+    }
+
+    lowerWall       
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack    
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/constant/RASProperties

@@ -0,0 +1,24 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel        kEpsilon;
+
+turbulence      on;
+
+printCoeffs     on;
+
+// ************************************************************************* //

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/constant/polyMesh/blockMeshDict

@@ -0,0 +1,153 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.001;
+
+vertices        
+(
+    (-20.6 0 -0.5)
+    (-20.6 3 -0.5)
+    (-20.6 12.7 -0.5)
+    (-20.6 25.4 -0.5)
+    (0 -25.4 -0.5)
+    (0 -5 -0.5)
+    (0 0 -0.5)
+    (0 3 -0.5)
+    (0 12.7 -0.5)
+    (0 25.4 -0.5)
+    (206 -25.4 -0.5)
+    (206 -8.5 -0.5)
+    (206 0 -0.5)
+    (206 6.5 -0.5)
+    (206 17 -0.5)
+    (206 25.4 -0.5)
+    (290 -16.6 -0.5)
+    (290 -6.3 -0.5)
+    (290 0 -0.5)
+    (290 4.5 -0.5)
+    (290 11 -0.5)
+    (290 16.6 -0.5)
+    (-20.6 0 0.5)
+    (-20.6 3 0.5)
+    (-20.6 12.7 0.5)
+    (-20.6 25.4 0.5)
+    (0 -25.4 0.5)
+    (0 -5 0.5)
+    (0 0 0.5)
+    (0 3 0.5)
+    (0 12.7 0.5)
+    (0 25.4 0.5)
+    (206 -25.4 0.5)
+    (206 -8.5 0.5)
+    (206 0 0.5)
+    (206 6.5 0.5)
+    (206 17 0.5)
+    (206 25.4 0.5)
+    (290 -16.6 0.5)
+    (290 -6.3 0.5)
+    (290 0 0.5)
+    (290 4.5 0.5)
+    (290 11 0.5)
+    (290 16.6 0.5)
+);
+
+blocks          
+(
+    hex (0 6 7 1 22 28 29 23) (18 7 1) simpleGrading (0.5 1.8 1)
+    hex (1 7 8 2 23 29 30 24) (18 10 1) simpleGrading (0.5 4 1)
+    hex (2 8 9 3 24 30 31 25) (18 13 1) simpleGrading (0.5 0.25 1)
+    hex (4 10 11 5 26 32 33 27) (180 18 1) simpleGrading (4 1 1)
+    hex (5 11 12 6 27 33 34 28) (180 9 1) edgeGrading (4 4 4 4 0.5 1 1 0.5 1 1 1 1)
+    hex (6 12 13 7 28 34 35 29) (180 7 1) edgeGrading (4 4 4 4 1.8 1 1 1.8 1 1 1 1)
+    hex (7 13 14 8 29 35 36 30) (180 10 1) edgeGrading (4 4 4 4 4 1 1 4 1 1 1 1)
+    hex (8 14 15 9 30 36 37 31) (180 13 1) simpleGrading (4 0.25 1)
+    hex (10 16 17 11 32 38 39 33) (25 18 1) simpleGrading (2.5 1 1)
+    hex (11 17 18 12 33 39 40 34) (25 9 1) simpleGrading (2.5 1 1)
+    hex (12 18 19 13 34 40 41 35) (25 7 1) simpleGrading (2.5 1 1)
+    hex (13 19 20 14 35 41 42 36) (25 10 1) simpleGrading (2.5 1 1)
+    hex (14 20 21 15 36 42 43 37) (25 13 1) simpleGrading (2.5 0.25 1)
+);
+
+edges           
+(
+);
+
+patches         
+(
+    patch inlet 
+    (
+        (0 22 23 1)
+        (1 23 24 2)
+        (2 24 25 3)
+    )
+    patch outlet 
+    (
+        (16 17 39 38)
+        (17 18 40 39)
+        (18 19 41 40)
+        (19 20 42 41)
+        (20 21 43 42)
+    )
+    wall upperWall 
+    (
+        (3 25 31 9)
+        (9 31 37 15)
+        (15 37 43 21)
+    )
+    wall lowerWall 
+    (
+        (0 6 28 22)
+        (6 5 27 28)
+        (5 4 26 27)
+        (4 10 32 26)
+        (10 16 38 32)
+    )
+    empty frontAndBack 
+    (
+        (22 28 29 23)
+        (23 29 30 24)
+        (24 30 31 25)
+        (26 32 33 27)
+        (27 33 34 28)
+        (28 34 35 29)
+        (29 35 36 30)
+        (30 36 37 31)
+        (32 38 39 33)
+        (33 39 40 34)
+        (34 40 41 35)
+        (35 41 42 36)
+        (36 42 43 37)
+        (0 1 7 6)
+        (1 2 8 7)
+        (2 3 9 8)
+        (4 5 11 10)
+        (5 6 12 11)
+        (6 7 13 12)
+        (7 8 14 13)
+        (8 9 15 14)
+        (10 11 17 16)
+        (11 12 18 17)
+        (12 13 19 18)
+        (13 14 20 19)
+        (14 15 21 20)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/constant/polyMesh/boundary

@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+5
+(
+    inlet
+    {
+        type            patch;
+        physicalType    inlet;
+        nFaces          30;
+        startFace       24170;
+    }
+    outlet
+    {
+        type            patch;
+        physicalType    outlet;
+        nFaces          57;
+        startFace       24200;
+    }
+    upperWall
+    {
+        type            wall;
+        physicalType    wallFunctions;
+        nFaces          223;
+        startFace       24257;
+    }
+    lowerWall
+    {
+        type            wall;
+        physicalType    wallFunctions;
+        nFaces          250;
+        startFace       24480;
+    }
+    frontAndBack
+    {
+        type            empty;
+        physicalType    empty;
+        nFaces          24450;
+        startFace       24730;
+    }
+)
+
+// ************************************************************************* //

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/constant/transportProperties

@@ -0,0 +1,25 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              nu [0 2 -1 0 0 0 0] 1e-5;
+
+lambda          lambda   [0 -2 1 0 0 0 0] 1e5;
+alphaMax        alphaMax [0 0 -1 0 0 0 0] 200.0;
+
+// ************************************************************************* //

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/system/controlDict

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     adjointShapeOptimizationFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1000;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   100;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  12;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   12;
+
+runTimeModifiable yes;
+
+// ************************************************************************* //

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/system/fvSchemes

@@ -0,0 +1,76 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+
+    grad(p)         Gauss linear;
+    grad(U)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+
+    div(phi,U)      Gauss upwind;
+    div(phi,k)      Gauss upwind;
+    div(phi,epsilon) Gauss upwind;
+    div((nuEff*dev(grad(U).T()))) Gauss linear;
+
+    div(-phi,Ua)    Gauss upwind;
+    div((nuEff*dev(grad(Ua).T()))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+
+    laplacian(nuEff,U) Gauss linear corrected;
+    laplacian(1,p) Gauss linear corrected;
+    laplacian(1|A(U),p) Gauss linear corrected;
+    laplacian(DkEff,k) Gauss linear corrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+
+    laplacian(nuEff,Ua) Gauss linear corrected;
+    laplacian(1|A(Ua),pa) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+    interpolate(U)  linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p;
+    pa;
+}
+
+// ************************************************************************* //

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/system/fvSolution

@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "(p|pa)"
+    {
+        solver          GAMG;
+        tolerance       1e-08;
+        relTol          0.01;
+        smoother        GaussSeidel;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    "(U|Ua|k|epsilon)"
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        nSweeps         2;
+        tolerance       1e-08;
+        relTol          0.1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    pRefCell 0;
+    pRefValue 0;
+}
+
+relaxationFactors
+{
+    "(p|pa)"        0.3;
+    "(U|Ua)"        0.7;
+    "(k|epsilon)"   0.7;
+    alpha           0.1;
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/0/U

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/

tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/0/alpha1

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/

tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/polyMesh/boundary

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/

tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/pointDisplacement

@@ -37,7 +37,8 @@ boundaryField
         centreOfMass    (0.5 0.5 0.5);
         momentOfInertia (0.08622222 0.08622222 0.144);
         mass            9.6;
-        rhoInf          1;                      // for forces calculation
+        rhoInf          1;
+        report          on;
         value           uniform (0 0 0);
     }
 }

tutorials/multiphase/interDyMFoam/ras/floatingObject/constant/polyMesh/boundary

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6.x                                 |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/

tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict

@@ -17,7 +17,7 @@ FoamFile
 
 application     interDyMFoam;
 
-startFrom       startTime;
+startFrom       latestTime;
 
 startTime       0;
 
@@ -29,7 +29,7 @@ deltaT          0.01;
 
 writeControl    adjustableRunTime;
 
-writeInterval   0.025;
+writeInterval   0.1;
 
 purgeWrite      0;
 
@@ -47,9 +47,9 @@ runTimeModifiable yes;
 
 adjustTimeStep  yes;
 
-maxCo           0.2;
+maxCo           0.5;
 
-maxDeltaT       0.025;
+maxDeltaT       0.01;
 
 libs
 (

tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution

@@ -20,7 +20,7 @@ solvers
     cellDisplacement
     {
         solver          GAMG;
-        tolerance       1e-08;
+        tolerance       1e-5;
         relTol          0;
         smoother        GaussSeidel;
         cacheAgglomeration true;
@@ -35,7 +35,7 @@ solvers
         preconditioner
         {
             preconditioner  GAMG;
-            tolerance       1e-05;
+            tolerance       1e-5;
             relTol          0;
             smoother        DICGaussSeidel;
             nPreSweeps      0;
@@ -55,7 +55,7 @@ solvers
     p
     {
         solver          GAMG;
-        tolerance       1e-08;
+        tolerance       1e-8;
         relTol          0.01;
         smoother        DIC;
         nPreSweeps      0;
@@ -73,7 +73,7 @@ solvers
         preconditioner
         {
             preconditioner  GAMG;
-            tolerance       2e-09;
+            tolerance       1e-8;
             relTol          0;
             nVcycles        2;
             smoother        DICGaussSeidel;
@@ -86,7 +86,7 @@ solvers
             mergeLevels     1;
         }
 
-        tolerance       2e-09;
+        tolerance       1e-8;
         relTol          0;
         maxIter         20;
     }
@@ -95,18 +95,18 @@ solvers
     {
         solver          smoothSolver;
         smoother        GaussSeidel;
-        tolerance       1e-06;
+        tolerance       1e-6;
         relTol          0;
         nSweeps         1;
     }
 
-    "(k|epsilon|R|nuTilda)"
+    "(k|epsilon|omega|R|nuTilda)"
     {
-        $U;
-        tolerance       1e-08;
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-8;
         relTol          0;
     }
-
 }
 
 PISO
@@ -115,14 +115,14 @@ PISO
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
     nAlphaCorr      1;
-    nAlphaSubCycles 1;
-    cAlpha          1.5;
+    nAlphaSubCycles 3;
+    cAlpha          1;
     correctPhi      yes;
 }
 
 relaxationFactors
 {
-    U               1;
+    "(U|k|epsilon|omega|R|nuTilda)"               1;
 }