Merge branch 'master' into splitCyclic

Conflicts:
	src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
	src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
	src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
	src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
	src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
	src/parallel/decompose/scotchDecomp/scotchDecomp.C
	src/parallel/parMetisDecomp/parMetisDecomp.C
	src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C

.gitignore

@@ -55,6 +55,10 @@ doc/[Dd]oxygen/man
 /*.html
 /doc/*.html
 
+# untracked configuration files
+/etc/prefs.csh
+/etc/prefs.sh
+
 # source packages - anywhere
 *.tar.bz2
 *.tar.gz

README

@@ -2,7 +2,7 @@
 #
 #+TITLE:             OpenFOAM README for version 1.6
 #+AUTHOR:                      OpenCFD Ltd.
-#+DATE:                       November 2009
+#+DATE:                         March 2010
 #+LINK:                  http://www.opencfd.co.uk
 #+OPTIONS: author:nil ^:{}
 
@@ -23,10 +23,10 @@
   section "Running OpenFOAM in 32-bit mode".
 
 *** Qt (from http://trolltech.com/products/qt)
-    The ParaView 3.6.1 visualisation package requires Qt to be installed on the
+    The ParaView 3.7.0 visualisation package requires Qt to be installed on the
     system.  ParaView's producers state that ParaView is only officially
-    supported on Qt version 4.3.x.  However, we have found in limited tests that
-    ParaView works satisfactorily with newer versions of Qt than 4.3.x.  To
+    supported on Qt version 4.6.x.  However, we have found in limited tests that
+    ParaView works satisfactorily with Qt than 4.5.x.  To
     check whether Qt4 is installed, and the version, type:
     + qmake --version
 
@@ -44,13 +44,14 @@
     + openSUSE-10.3: Version 4.3.1
     + openSUSE-11.0: Version 4.4.0
     + openSUSE-11.1: Version 4.4.3
+    + openSUSE-11.2: Version 4.5.3
 
     Compilation and running of ParaView has been successful using the libraries
     downloaded in the "libqt4-dev" package on ubuntu.
 
     If you don't have an appropriate version of Qt installed you can download
-    the sources from TrollTech e.g.:
-    ftp://ftp.trolltech.com/qt/source/qt-x11-opensource-src-4.3.5.tar.bz2
+    the sources e.g.:
+    http://get.qt.nokia.com/qt/source/qt-everywhere-opensource-src-4.6.2.tar.gz
     and compile and install in /usr/local or some other location that does not
     conflict with the pre-installed version.
 
@@ -108,11 +109,11 @@
 * Building from Sources (Optional)
   If you cannot find an appropriate binary pack for your platform, you can build
   the complete OpenFOAM from the source-pack.  You will first need to compile or
-  obtain a recent version of gcc (we recomend gcc-4.3.?) for your platform,
+  obtain a recent version of gcc (we recommend gcc-4.4.?) for your platform,
   which may be obtained from http://gcc.gnu.org/.
 
   Install the compiler in
-  $WM_THIRD_PARTY_DIR/gcc-<GCC_VERSION>/platforms/$WM_ARCH$WM_COMPILER_ARCH/
+  $WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER_ARCH/gcc-<GCC_VERSION>
   and change the gcc version number in $WM_PROJECT_DIR/etc/settings.sh and
   $WM_PROJECT_DIR/etc/settings.csh appropriately and finally update the
   environment variables as in section 3.
@@ -157,15 +158,16 @@
   Refer to the OpenFOAM User Guide at http://www.OpenFOAM.org/doc/user.html for
   more information.
 
-* Compiling Paraview 3.6.1 and the PV3FoamReader module
+* Compiling Paraview 3.7.0 and the PV3FoamReader module
   If there are problems encountered with ParaView, then it may be necessary to
   compile ParaView from sources.  The compilation
   is a fairly simple process using the makeParaView script
   (found in ThirdParty directory), which has worked in our tests with other
-  packages supplied in the ThirdParty directory, namely cmake-2.6.4 and
-  gcc-4.3.3.  Execute the following:
+  packages supplied in the ThirdParty directory, namely cmake-2.8.0 and
+  gcc-4.4.3.  Execute the following:
   + cd $WM_THIRD_PARTY_DIR
-  + rm -rf paraview-3.6.1/platforms
+  + rm -rf paraview-3.7.0/platforms
+  + rm -rf platforms/*/paraview-3.7.0
   + ./makeParaView
 
   The PV3blockMeshReader and the PV3FoamReader ParaView plugins are compiled
@@ -177,7 +179,7 @@
 *** Compiling Paraview with a local version of Qt
     If the user still encounters problems with ParaView, it may relate to the
     version of Qt, in which case, it is recommended that the user first
-    downloads a supported version of Qt /e.g./ 4.3.5 as described in the section
+    downloads a supported version of Qt /e.g./ 4.5.3 as described in the section
     on "Qt".  The user should unpack the source pack in the $WM_THIRD_PARTY_DIR.
     Then the user can build Qt by executing from within $WM_THIRD_PARTY_DIR:
     + ./makeQt

TODO

@@ -24,6 +24,8 @@ OK - amg.
 Tested on unitTestCases/singleCyclic/
 
 - initTransfer in GAMGprocessorInterfaces using nonblocking+tags
+untested.
+
 - test createPatch pointSync
 - pointFields on cyclics. volPointInterpolation.
 - jumpCyclics

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/compressible/rhoPorousMRFPimpleFoam/Make/files

@@ -0,0 +1,4 @@
+rhoPorousMRFPimpleFoam.C
+
+EXE = $(FOAM_APPBIN)/rhoPorousMRFPimpleFoam
+

applications/solvers/compressible/rhoPorousMRFPimpleFoam/Make/options

@@ -0,0 +1,15 @@
+EXE_INC = \
+    -I../rhoPimpleFoam \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lfiniteVolume \
+    -lmeshTools

applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H

@@ -0,0 +1,39 @@
+// Solve the Momentum equation
+
+tmp<fvVectorMatrix> UEqn
+(
+    pZones.ddt(rho, U)
+  + fvm::div(phi, U)
+  + turbulence->divDevRhoReff(U)
+);
+
+if (oCorr == nOuterCorr-1)
+{
+    UEqn().relax(1);
+}
+else
+{
+    UEqn().relax();
+}
+
+mrfZones.addCoriolis(rho, UEqn());
+pZones.addResistance(UEqn());
+
+volScalarField rUA = 1.0/UEqn().A();
+
+if (momentumPredictor)
+{
+    if (oCorr == nOuterCorr-1)
+    {
+        solve(UEqn() == -fvc::grad(p), mesh.solver("UFinal"));
+    }
+    else
+    {
+        solve(UEqn() == -fvc::grad(p));
+    }
+}
+else
+{
+    U = rUA*(UEqn().H() - fvc::grad(p));
+    U.correctBoundaryConditions();
+}

applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H

@@ -0,0 +1,70 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<basicPsiThermo> pThermo
+    (
+        basicPsiThermo::New(mesh)
+    );
+    basicPsiThermo& thermo = pThermo();
+
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "compressibleCreatePhi.H"
+
+    dimensionedScalar pMin
+    (
+        mesh.solutionDict().subDict("PIMPLE").lookup("pMin")
+    );
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    //dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+
+    Info<< "Creating field DpDt\n" << endl;
+    volScalarField DpDt =
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+    MRFZones mrfZones(mesh);
+    mrfZones.correctBoundaryVelocity(U);
+
+    porousZones pZones(mesh);
+    Switch pressureImplicitPorosity(false);

applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H

@@ -0,0 +1,123 @@
+rho = thermo.rho();
+
+volScalarField rUA = 1.0/UEqn().A();
+U = rUA*UEqn().H();
+
+if (nCorr <= 1)
+{
+    UEqn.clear();
+}
+
+if (transonic)
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        )
+    );
+    mrfZones.relativeFlux(fvc::interpolate(psi), phid);
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        if
+        (
+            oCorr == nOuterCorr-1
+            && corr == nCorr-1
+            && nonOrth == nNonOrthCorr
+        )
+        {
+            pEqn.solve(mesh.solver("pFinal"));
+        }
+        else
+        {
+            pEqn.solve();
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi == pEqn.flux();
+        }
+    }
+}
+else
+{
+    phi =
+        fvc::interpolate(rho)*
+        (
+            (fvc::interpolate(U) & mesh.Sf())
+        //+ fvc::ddtPhiCorr(rUA, rho, U, phi)
+        );
+    mrfZones.relativeFlux(fvc::interpolate(rho), phi);
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        // Pressure corrector
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phi)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        if
+        (
+            oCorr == nOuterCorr-1
+         && corr == nCorr-1
+         && nonOrth == nNonOrthCorr
+        )
+        {
+            pEqn.solve(mesh.solver("pFinal"));
+        }
+        else
+        {
+            pEqn.solve();
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+//if (oCorr != nOuterCorr-1)
+{
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    rho = thermo.rho();
+    rho.relax();
+    Info<< "rho max/min : " << max(rho).value()
+        << " " << min(rho).value() << endl;
+}
+
+U -= rUA*fvc::grad(p);
+U.correctBoundaryConditions();
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+bound(p, pMin);
+
+// For closed-volume cases adjust the pressure and density levels
+// to obey overall mass continuity
+/*
+if (closedVolume)
+{
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /fvc::domainIntegrate(psi);
+}
+*/

applications/solvers/compressible/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C

@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    rhoPorousMRFPimpleFoam
+
+Description
+    Transient solver for laminar or turbulent flow of compressible fluids
+    with support for porous media and MRF for HVAC and similar applications.
+
+    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
+    pseudo-transient simulations.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "basicPsiThermo.H"
+#include "turbulenceModel.H"
+#include "bound.H"
+#include "MRFZones.H"
+#include "porousZones.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+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;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "readPIMPLEControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        if (nOuterCorr != 1)
+        {
+            p.storePrevIter();
+            rho.storePrevIter();
+        }
+
+        #include "rhoEqn.H"
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        {
+            #include "UEqn.H"
+            #include "hEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<nCorr; corr++)
+            {
+                #include "pEqn.H"
+            }
+
+            turbulence->correct();
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    createPhia
+
+Description
+    Creates and initialises the face-flux field phia.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef createPhia_H
+#define createPhia_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Info<< "Reading/calculating face flux field phia\n" << endl;
+
+surfaceScalarField phia
+(
+    IOobject
+    (
+        "phia",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    linearInterpolate(Ua) & mesh.Sf()
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqns.H

@@ -59,7 +59,17 @@
                 alpharScheme
             );
 
-        MULES::explicitSolve(oneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0);
+        MULES::explicitSolve
+        (
+            geometricOneField(),
+            alpha1,
+            phi,
+            phiAlpha1,
+            Sp,
+            Su,
+            1,
+            0
+        );
 
         surfaceScalarField rho1f = fvc::interpolate(rho1);
         surfaceScalarField rho2f = fvc::interpolate(rho2);

applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H

@@ -59,7 +59,7 @@
                 alpharScheme
             );
 
-        MULES::explicitSolve(oneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0);
+        MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0);
 
         surfaceScalarField rho1f = fvc::interpolate(rho1);
         surfaceScalarField rho2f = fvc::interpolate(rho2);

applications/solvers/multiphase/interMixingFoam/alphaEqns.H

@@ -70,12 +70,12 @@
         MULES::limiter
         (
             allLambda,
-            oneField(),
+            geometricOneField(),
             alpha1,
             phiAlpha1BD,
             phiAlpha1,
-            zeroField(),
-            zeroField(),
+            zero(),
+            zero(),
             1,
             0,
             3
@@ -107,12 +107,12 @@
         MULES::limiter
         (
             allLambda,
-            oneField(),
+            geometricOneField(),
             alpha2,
             phiAlpha2BD,
             phiAlpha2,
-            zeroField(),
-            zeroField(),
+            zero(),
+            zero(),
             1,
             0,
             3

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H

@@ -51,8 +51,8 @@
         );
 
         //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
-        //MULES::explicitSolve(oneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
-        MULES::implicitSolve(oneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
+        //MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
+        MULES::implicitSolve(geometricOneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
 
         rhoPhi +=
             (runTime.deltaT()/totalDeltaT)

applications/test/nearWallDist-wave/testWallDistData.C

@@ -58,7 +58,7 @@ int main(int argc, char *argv[])
             mesh
         ),
         mesh,
-        dimensionedVector("n", dimLength, wallPoint::greatPoint)
+        dimensionedVector("n", dimLength, point::max)
     );
 
     // Fill wall patches with unit normal

applications/test/syncTools/syncToolsTest.C

@@ -187,8 +187,6 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
     );
     const pointField& localPoints = allBoundary.localPoints();
 
-    const point greatPoint(GREAT, GREAT, GREAT);
-
 
     // Point data
     // ~~~~~~~~~~
@@ -196,7 +194,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
     {
         // Create some data. Use slightly perturbed positions.
         Map<vector> sparseData;
-        pointField fullData(mesh.nPoints(), greatPoint);
+        pointField fullData(mesh.nPoints(), point::max);
 
         forAll(localPoints, i)
         {
@@ -222,7 +220,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
             mesh,
             fullData,
             minEqOp<vector>(),
-            greatPoint,
+            point::max,
             true                    // apply separation
         );
 
@@ -232,7 +230,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
         {
             const point& fullPt = fullData[meshPointI];
 
-            if (fullPt != greatPoint)
+            if (fullPt != point::max)
             {
                 const point& sparsePt = sparseData[meshPointI];
 
@@ -272,7 +270,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
     {
         // Create some data. Use slightly perturbed positions.
         EdgeMap<vector> sparseData;
-        pointField fullData(mesh.nEdges(), greatPoint);
+        pointField fullData(mesh.nEdges(), point::max);
 
         const edgeList& edges = allBoundary.edges();
         const labelList meshEdges = allBoundary.meshEdges
@@ -307,7 +305,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
             mesh,
             fullData,
             minEqOp<vector>(),
-            greatPoint,
+            point::max,
             true
         );
 
@@ -317,7 +315,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
         {
             const point& fullPt = fullData[meshEdgeI];
 
-            if (fullPt != greatPoint)
+            if (fullPt != point::max)
             {
                 const point& sparsePt = sparseData[mesh.edges()[meshEdgeI]];
 
@@ -364,8 +362,6 @@ void testPointSync(const polyMesh& mesh, Random& rndGen)
 {
     Info<< nl << "Testing point-wise data synchronisation." << endl;
 
-    const point greatPoint(GREAT, GREAT, GREAT);
-
     // Test position.
 
     {
@@ -379,7 +375,7 @@ void testPointSync(const polyMesh& mesh, Random& rndGen)
             mesh,
             syncedPoints,
             minEqOp<point>(),
-            greatPoint,
+            point::max,
             true
         );
 
@@ -444,8 +440,6 @@ void testEdgeSync(const polyMesh& mesh, Random& rndGen)
 
     const edgeList& edges = mesh.edges();
 
-    const point greatPoint(GREAT, GREAT, GREAT);
-
     // Test position.
 
     {
@@ -463,7 +457,7 @@ void testEdgeSync(const polyMesh& mesh, Random& rndGen)
             mesh,
             syncedMids,
             minEqOp<point>(),
-            greatPoint,
+            point::max,
             true
         );
 

applications/utilities/mesh/generation/extrudeToRegionMesh/Make/files

@@ -0,0 +1,7 @@
+patchPointEdgeCirculator.C
+createShellMesh.C
+extrudeToRegionMesh.C
+
+EXE = $(FOAM_APPBIN)/extrudeToRegionMesh
+
+

applications/utilities/mesh/generation/extrudeToRegionMesh/Make/options

@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    /* -I$(LIB_SRC)/surfMesh/lnInclude */ \
+    /* -I$(LIB_SRC)/lagrangian/basic/lnInclude */ \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    /* -lsurfMesh */ \
+    /* -llagrangian */ \
+    -lmeshTools \
+    -ldynamicMesh
+

applications/utilities/mesh/generation/extrudeToRegionMesh/createShellMesh.C

@@ -0,0 +1,588 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "createShellMesh.H"
+#include "polyTopoChange.H"
+#include "meshTools.H"
+#include "mapPolyMesh.H"
+#include "polyAddPoint.H"
+#include "polyAddFace.H"
+#include "polyModifyFace.H"
+#include "polyAddCell.H"
+#include "patchPointEdgeCirculator.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::createShellMesh, 0);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::createShellMesh::calcPointRegions
+(
+    const primitiveFacePatch& patch,
+    const PackedBoolList& nonManifoldEdge,
+    faceList& pointRegions,
+    labelList& regionPoints
+)
+{
+    pointRegions.setSize(patch.size());
+    forAll(pointRegions, faceI)
+    {
+        const face& f = patch.localFaces()[faceI];
+        pointRegions[faceI].setSize(f.size(), -1);
+    }
+
+    label nRegions = 0;
+
+    forAll(pointRegions, faceI)
+    {
+        const face& f = patch.localFaces()[faceI];
+
+        forAll(f, fp)
+        {
+            if (pointRegions[faceI][fp] == -1)
+            {
+                // Found unassigned point. Distribute current region.
+                label pointI = f[fp];
+                label edgeI = patch.faceEdges()[faceI][fp];
+
+                patchPointEdgeCirculator circ
+                (
+                    patch,
+                    nonManifoldEdge,
+                    edgeI,
+                    findIndex(patch.edgeFaces()[edgeI], faceI),
+                    pointI
+                );
+
+                for
+                (
+                    patchPointEdgeCirculator iter = circ.begin();
+                    iter != circ.end();
+                    ++iter
+                )
+                {
+                    label face2 = iter.faceID();
+
+                    if (face2 != -1)
+                    {
+                        const face& f2 = patch.localFaces()[face2];
+                        label fp2 = findIndex(f2, pointI);
+                        label& region = pointRegions[face2][fp2];
+                        if (region != -1)
+                        {   
+                            FatalErrorIn
+                            (
+                                "createShellMesh::calcPointRegions(..)"
+                            )   << "On point " << pointI
+                                << " at:" << patch.localPoints()[pointI]
+                                << " found region:" << region
+                                << abort(FatalError);
+                        }
+                        region = nRegions;
+                    }
+                }
+
+                nRegions++;
+            }
+        }
+    }
+
+
+    // From region back to originating point (many to one, a point might
+    // have multiple regions though)
+    regionPoints.setSize(nRegions);
+    forAll(pointRegions, faceI)
+    {
+        const face& f = patch.localFaces()[faceI];
+
+        forAll(f, fp)
+        {
+            regionPoints[pointRegions[faceI][fp]] = f[fp];
+        }
+    }
+
+
+    if (debug)
+    {
+        const labelListList& pointFaces = patch.pointFaces();
+        forAll(pointFaces, pointI)
+        {
+            label region = -1;
+            const labelList& pFaces = pointFaces[pointI];
+            forAll(pFaces, i)
+            {
+                label faceI = pFaces[i];
+                const face& f = patch.localFaces()[faceI];
+                label fp = findIndex(f, pointI);
+
+                if (region == -1)
+                {
+                    region = pointRegions[faceI][fp];
+                }
+                else if (region != pointRegions[faceI][fp])
+                {
+                    Pout<< "Non-manifold point:" << pointI
+                        << " at " << patch.localPoints()[pointI]
+                        << " region:" << region
+                        << " otherRegion:" << pointRegions[faceI][fp]
+                        << endl;
+
+                }
+            }
+        }
+    }
+}
+
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::createShellMesh::createShellMesh
+(
+    const primitiveFacePatch& patch,
+    const faceList& pointRegions,
+    const labelList& regionPoints
+)
+:
+    patch_(patch),
+    pointRegions_(pointRegions),
+    regionPoints_(regionPoints)
+{
+    if (pointRegions_.size() != patch_.size())
+    {
+        FatalErrorIn("createShellMesh::createShellMesh(..)")
+            << "nFaces:" << patch_.size()
+            << " pointRegions:" << pointRegions.size()
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+void Foam::createShellMesh::setRefinement
+(
+    const pointField& thickness,
+    const labelList& topPatchID,
+    const labelList& bottomPatchID,
+    const labelListList& extrudeEdgePatches,
+    polyTopoChange& meshMod
+)
+{
+    if (thickness.size() != regionPoints_.size())
+    {
+        FatalErrorIn("createShellMesh::setRefinement(..)")
+            << "nRegions:" << regionPoints_.size()
+            << " thickness:" << thickness.size()
+            << exit(FatalError);
+    }
+
+    if
+    (
+        topPatchID.size() != patch_.size()
+     && bottomPatchID.size() != patch_.size()
+    )
+    {
+        FatalErrorIn("createShellMesh::setRefinement(..)")
+            << "nFaces:" << patch_.size()
+            << " topPatchID:" << topPatchID.size()
+            << " bottomPatchID:" << bottomPatchID.size()
+            << exit(FatalError);
+    }
+
+    if (extrudeEdgePatches.size() != patch_.nEdges())
+    {
+        FatalErrorIn("createShellMesh::setRefinement(..)")
+            << "nEdges:" << patch_.nEdges()
+            << " extrudeEdgePatches:" << extrudeEdgePatches.size()
+            << exit(FatalError);
+    }
+
+
+
+    // From cell to patch (trivial)
+    DynamicList<label> cellToFaceMap(patch_.size());
+    // From face to patch+turning index
+    DynamicList<label> faceToFaceMap(2*patch_.size()+patch_.nEdges());
+    // From face to patch edge index
+    DynamicList<label> faceToEdgeMap(patch_.nEdges()+patch_.nEdges());
+    // From point to patch point index
+    DynamicList<label> pointToPointMap(2*patch_.nPoints());
+
+
+    // Introduce new cell for every face
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelList addedCells(patch_.size());
+    forAll(patch_, faceI)
+    {
+        addedCells[faceI] = meshMod.addCell
+        (
+            -1,                     // masterPointID
+            -1,                     // masterEdgeID
+            -1,                     // masterFaceID
+            cellToFaceMap.size(),   // masterCellID
+            -1                      // zoneID
+        );
+        cellToFaceMap.append(faceI);
+    }
+
+
+    // Introduce original points
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Original point numbers in local point ordering so no need to store.
+    forAll(patch_.localPoints(), pointI)
+    {
+        //label addedPointI =
+        meshMod.addPoint
+        (
+            patch_.localPoints()[pointI],   // point
+            pointToPointMap.size(),         // masterPointID 
+            -1,                             // zoneID
+            true                            // inCell
+        );
+        pointToPointMap.append(pointI);
+
+        //Pout<< "Added bottom point " << addedPointI
+        //    << " at " << patch_.localPoints()[pointI]
+        //    << "  from point " << pointI
+        //    << endl;
+    }
+
+
+    // Introduce new points (one for every region)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelList addedPoints(regionPoints_.size());
+    forAll(regionPoints_, regionI)
+    {
+        label pointI = regionPoints_[regionI];
+        point extrudedPt = patch_.localPoints()[pointI] + thickness[regionI];
+
+        addedPoints[regionI] = meshMod.addPoint
+        (
+            extrudedPt,             // point
+            pointToPointMap.size(), // masterPointID - used only addressing
+            -1,                     // zoneID
+            true                    // inCell
+        );
+        pointToPointMap.append(pointI);
+
+        //Pout<< "Added top point " << addedPoints[regionI]
+        //    << " at " << extrudedPt
+        //    << "  from point " << pointI
+        //    << endl;
+    }
+
+
+    // Add face on bottom side
+    forAll(patch_.localFaces(), faceI)
+    {
+        meshMod.addFace
+        (
+            patch_.localFaces()[faceI].reverseFace(),// vertices
+            addedCells[faceI],          // own
+            -1,                         // nei
+            -1,                         // masterPointID
+            -1,                         // masterEdgeID
+            faceToFaceMap.size(),       // masterFaceID : current faceI
+            true,                       // flipFaceFlux
+            bottomPatchID[faceI],// patchID
+            -1,                         // zoneID
+            false                       // zoneFlip
+        );
+        faceToFaceMap.append(-faceI-1); // points to flipped original face
+        faceToEdgeMap.append(-1);
+
+        //Pout<< "Added bottom face "
+        //    << patch_.localFaces()[faceI].reverseFace()
+        //    << " own " << addedCells[faceI]
+        //    << "  at " << patch_.faceCentres()[faceI]
+        //    << endl;
+
+    }
+
+    // Add face on top
+    forAll(patch_.localFaces(), faceI)
+    {
+        // Get face in original ordering
+        const face& f = patch_.localFaces()[faceI];
+
+        // Pick up point based on region
+        face newF(f.size());
+        forAll(f, fp)
+        {
+            label region = pointRegions_[faceI][fp];
+            newF[fp] = addedPoints[region];
+        }
+
+        meshMod.addFace
+        (
+            newF,                       // vertices
+            addedCells[faceI],          // own
+            -1,                         // nei
+            -1,                         // masterPointID
+            -1,                         // masterEdgeID
+            faceToFaceMap.size(),       // masterFaceID : current faceI
+            false,                      // flipFaceFlux
+            topPatchID[faceI],          // patchID
+            -1,                         // zoneID
+            false                       // zoneFlip
+        );
+        faceToFaceMap.append(faceI+1);  // unflipped
+        faceToEdgeMap.append(-1);
+
+        //Pout<< "Added top face " << newF
+        //    << " own " << addedCells[faceI]
+        //    << "  at " << patch_.faceCentres()[faceI]
+        //    << endl;
+    }
+
+
+    // Add side faces
+    // ~~~~~~~~~~~~~~
+    // Note that we loop over edges multiple times so for edges with
+    // two cyclic faces they get added in two passes (for correct ordering)
+
+    // Pass1. Internal edges and first face of other edges
+    forAll(extrudeEdgePatches, edgeI)
+    {
+        const labelList& eFaces = patch_.edgeFaces()[edgeI];
+        const labelList& ePatches = extrudeEdgePatches[edgeI];
+
+        if (ePatches.size() == 0)
+        {
+            // internal face.
+
+            if (eFaces.size() != 2)
+            {
+                FatalErrorIn("createShellMesh::setRefinement(..)")
+                    << "edge:" << edgeI
+                    << " not internal but does not have side-patches defined."
+                    << exit(FatalError);
+            }
+
+            // Extrude
+
+            // Make face pointing in to eFaces[0] so out of new master face
+            const face& f = patch_.localFaces()[eFaces[0]];
+            const edge& e = patch_.edges()[edgeI];
+
+            label fp0 = findIndex(f, e[0]);
+            label fp1 = f.fcIndex(fp0);
+
+            if (f[fp1] != e[1])
+            {
+                fp1 = fp0;
+                fp0 = f.rcIndex(fp1);
+            }
+
+            face newF(4);
+            newF[0] = f[fp0];
+            newF[1] = f[fp1];
+            newF[2] = addedPoints[pointRegions_[eFaces[0]][fp1]];
+            newF[3] = addedPoints[pointRegions_[eFaces[0]][fp0]];
+
+            label minCellI = addedCells[eFaces[0]];
+            label maxCellI = addedCells[eFaces[1]];
+
+            if (minCellI > maxCellI)
+            {
+                // Swap
+                Swap(minCellI, maxCellI);
+                newF = newF.reverseFace();
+            }
+
+            //Pout<< "for internal edge:" << e
+            //    << " at:" << patch_.localPoints()[e[0]]
+            //    << patch_.localPoints()[e[1]]
+            //    << " adding face:" << newF
+            //    << " from f:" << f
+            //    << " inbetween " << minCellI << " and " << maxCellI << endl;
+
+            // newF already outwards pointing.
+            meshMod.addFace
+            (
+                newF,       // vertices
+                minCellI,   // own
+                maxCellI,   // nei
+                -1,         // masterPointID
+                -1,         // masterEdgeID
+                faceToFaceMap.size(),   // masterFaceID
+                false,      // flipFaceFlux
+                -1,         // patchID
+                -1,         // zoneID
+                false       // zoneFlip
+            );
+            faceToFaceMap.append(0);
+            faceToEdgeMap.append(edgeI);
+        }
+        else
+        {
+            if (eFaces.size() != ePatches.size())
+            {
+                FatalErrorIn("createShellMesh::setRefinement(..)")
+                    << "external/feature edge:" << edgeI
+                    << " has " << eFaces.size() << " connected extruded faces "
+                    << " but only " << ePatches.size()
+                    << " boundary faces defined." << exit(FatalError);
+            }
+
+
+            // Extrude eFaces[0]
+            label minFaceI = eFaces[0];
+
+            // Make face pointing in to eFaces[0] so out of new master face
+            const face& f = patch_.localFaces()[minFaceI];
+
+            const edge& e = patch_.edges()[edgeI];
+            label fp0 = findIndex(f, e[0]);
+            label fp1 = f.fcIndex(fp0);
+
+            if (f[fp1] != e[1])
+            {
+                fp1 = fp0;
+                fp0 = f.rcIndex(fp1);
+            }
+
+            face newF(4);
+            newF[0] = f[fp0];
+            newF[1] = f[fp1];
+            newF[2] = addedPoints[pointRegions_[minFaceI][fp1]];
+            newF[3] = addedPoints[pointRegions_[minFaceI][fp0]];
+
+
+            //Pout<< "for external edge:" << e
+            //    << " at:" << patch_.localPoints()[e[0]]
+            //    << patch_.localPoints()[e[1]]
+            //    << " adding first patch face:" << newF
+            //    << " from:" << f
+            //    << " into patch:" << ePatches[0]
+            //    << " own:" << addedCells[minFaceI]
+            //    << endl;
+
+
+            // newF already outwards pointing.
+            meshMod.addFace
+            (
+                newF,                   // vertices
+                addedCells[minFaceI],   // own
+                -1,                     // nei
+                -1,                     // masterPointID
+                -1,                     // masterEdgeID
+                faceToFaceMap.size(),   // masterFaceID
+                false,                  // flipFaceFlux
+                ePatches[0],            // patchID
+                -1,                     // zoneID
+                false                   // zoneFlip
+            );
+            faceToFaceMap.append(0);
+            faceToEdgeMap.append(edgeI);
+        }
+    }
+
+    // Pass2. Other faces of boundary edges
+    forAll(extrudeEdgePatches, edgeI)
+    {
+        const labelList& eFaces = patch_.edgeFaces()[edgeI];
+        const labelList& ePatches = extrudeEdgePatches[edgeI];
+
+        if (ePatches.size() >= 2)
+        {
+            for (label i = 1; i < ePatches.size(); i++)
+            {
+                // Extrude eFaces[i]
+                label minFaceI = eFaces[i];
+
+                // Make face pointing in to eFaces[0] so out of new master face
+                const face& f = patch_.localFaces()[minFaceI];
+
+                const edge& e = patch_.edges()[edgeI];
+                label fp0 = findIndex(f, e[0]);
+                label fp1 = f.fcIndex(fp0);
+
+                if (f[fp1] != e[1])
+                {
+                    fp1 = fp0;
+                    fp0 = f.rcIndex(fp1);
+                }
+
+                face newF(4);
+                newF[0] = f[fp0];
+                newF[1] = f[fp1];
+                newF[2] = addedPoints[pointRegions_[minFaceI][fp1]];
+                newF[3] = addedPoints[pointRegions_[minFaceI][fp0]];
+
+                //Pout<< "for external edge:" << e
+                //    << " at:" << patch_.localPoints()[e[0]]
+                //    << patch_.localPoints()[e[1]]
+                //    << " adding patch face:" << newF
+                //    << " from:" << f
+                //    << " into patch:" << ePatches[i]
+                //    << endl;
+
+                // newF already outwards pointing.
+                meshMod.addFace
+                (
+                    newF,                   // vertices
+                    addedCells[minFaceI],   // own
+                    -1,                     // nei
+                    -1,                     // masterPointID
+                    -1,                     // masterEdgeID
+                    faceToFaceMap.size(),   // masterFaceID
+                    false,                  // flipFaceFlux
+                    ePatches[i],            // patchID
+                    -1,                     // zoneID
+                    false                   // zoneFlip
+                );
+                faceToFaceMap.append(0);
+                faceToEdgeMap.append(edgeI);
+            }
+        }
+    }
+
+
+    cellToFaceMap_.transfer(cellToFaceMap);
+    faceToFaceMap_.transfer(faceToFaceMap);
+    faceToEdgeMap_.transfer(faceToEdgeMap);
+    pointToPointMap_.transfer(pointToPointMap);
+}
+
+
+void Foam::createShellMesh::updateMesh(const mapPolyMesh& map)
+{
+    inplaceReorder(map.reverseCellMap(), cellToFaceMap_);
+    inplaceReorder(map.reverseFaceMap(), faceToFaceMap_);
+    inplaceReorder(map.reverseFaceMap(), faceToEdgeMap_);
+    inplaceReorder(map.reversePointMap(), pointToPointMap_);
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/extrudeToRegionMesh/createShellMesh.H

@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::createShellMesh
+
+Description
+    Creates mesh by extruding a patch.
+
+SourceFiles
+    createShellMesh.C
+
+    Extrudes into thickness direction.
+    - bottom faces originate from reversed original faces (have turning index)
+    - top faces originate from original faces (no turning index)
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef createShellMesh_H
+#define createShellMesh_H
+
+#include "primitiveFacePatch.H"
+#include "PackedBoolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class mapPolyMesh;
+class polyTopoChange;
+
+/*---------------------------------------------------------------------------*\
+                           Class createShellMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class createShellMesh
+{
+    // Private data
+
+        //- Reference to patch to extrude
+        const primitiveFacePatch& patch_;
+
+        //- Per point on face the region
+        const faceList& pointRegions_;
+
+        //- From region to original patch point
+        const labelList& regionPoints_;
+
+
+        labelList cellToFaceMap_;
+
+        labelList faceToFaceMap_;
+
+        labelList faceToEdgeMap_;
+
+        labelList pointToPointMap_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        createShellMesh(const createShellMesh&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const createShellMesh&);
+
+
+public:
+
+    //- Runtime type information
+    ClassName("createShellMesh");
+
+
+    // Constructors
+
+        //- Construct from mesh.
+        createShellMesh
+        (
+            const primitiveFacePatch& patch,
+            const faceList& pointRegions,
+            const labelList& regionPoints
+        );
+
+
+    // Member Functions
+
+
+        // Access
+
+            //- From region cell to patch face
+            const labelList& cellToFaceMap() const
+            {
+                return cellToFaceMap_;
+            }
+
+            //- From region face to patch face. Contains turning index:
+            //      > 0 : face in same orientation as patch face. face will
+            //            be in top patch
+            //      < 0 : face in opposite orientation as patch face. face will
+            //            be in bottom patch
+            //      = 0 : for all side faces
+            const labelList& faceToFaceMap() const
+            {
+                return faceToFaceMap_;
+            }
+
+            //- From region side-face to patch edge. -1 for non-edge faces.
+            const labelList& faceToEdgeMap() const
+            {
+                return faceToEdgeMap_;
+            }
+
+            //- From region point to patch point.
+            const labelList& pointToPointMap() const
+            {
+                return pointToPointMap_;
+            }
+
+
+        // Other
+
+            //- Helper: calculate point regions
+            static void calcPointRegions
+            (
+                const primitiveFacePatch& patch,
+                const PackedBoolList& nonManifoldEdge,
+                faceList& pointRegions,
+                labelList& regionPoints
+            );
+
+            //- Play commands into polyTopoChange to create layer mesh.
+            void setRefinement
+            (
+                const pointField& thickness,
+                const labelList& topPatchID,
+                const labelList& bottomPatchID,
+                const labelListList& extrudeEdgePatches,
+                polyTopoChange& meshMod
+            );
+
+            //- Update any locally stored mesh information.
+            void updateMesh(const mapPolyMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/extrudeToRegionMesh/patchPointEdgeCirculator.C

@@ -0,0 +1,68 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "patchPointEdgeCirculator.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::patchPointEdgeCirculator
+Foam::patchPointEdgeCirculator::endConstIter
+(
+    *reinterpret_cast<primitiveFacePatch*>(0),  // primitiveFacePatch
+    *reinterpret_cast<PackedBoolList*>(0),      // PackedBoolList
+    -1,                                         // edgeID
+    -1,                                         // index
+    -1                                          // pointID
+);
+
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const InfoProxy<patchPointEdgeCirculator>& ip
+)
+{
+    const patchPointEdgeCirculator& c = ip.t_;
+
+    const pointField& pts = c.patch_.localPoints();
+    const edge& e = c.patch_.edges()[c.edgeID_];
+    label faceI = c.faceID();
+
+    os  << "around point:" << c.pointID_
+        << " coord:" << pts[c.pointID_]
+        << " at edge:" << c.edgeID_
+        << " coord:" << pts[e.otherVertex(c.pointID_)]
+        << " in direction of face:" << faceI;
+
+    if (faceI != -1)
+    {
+        os<< " fc:" << c.patch_.localFaces()[faceI].centre(pts);
+    }
+    return os;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/extrudeToRegionMesh/patchPointEdgeCirculator.H

@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::patchPointEdgeCirculator
+
+Description
+    Walks from starting edge/face around point on patch.
+
+    -# Use in-place: \n
+        @code
+            patchPointEdgeCirculator circ(..);
+
+            // Walk
+            do
+            {
+                Info<< "edge:" << circ.edgeID() << endl;
+                ++circ;
+            }
+            while (circ != circ.end());
+        @endcode
+
+    -# Use like STL iterator: \n
+        @code
+            patchPointEdgeCirculator circ(..);
+            for
+            (
+                patchPointEdgeCirculator iter = circ.begin();
+                iter != circ.end();
+                ++iter
+            )
+            {
+                Info<< "edge:" << iter.edgeID() << endl;
+            }
+        @endcode
+
+
+SourceFiles
+    patchPointEdgeCirculator.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef patchPointEdgeCirculator_H
+#define patchPointEdgeCirculator_H
+
+#include "face.H"
+#include "ListOps.H"
+#include "primitiveFacePatch.H"
+#include "PackedBoolList.H"
+#include "InfoProxy.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+/*---------------------------------------------------------------------------*\
+                           Class patchPointEdgeCirculator Declaration
+\*---------------------------------------------------------------------------*/
+
+class patchPointEdgeCirculator
+{
+    // Static data members
+
+        //- end iterator
+        static const patchPointEdgeCirculator endConstIter;
+
+
+    // Private data
+
+        //- Patch
+        const primitiveFacePatch& patch_;
+
+        const PackedBoolList& nonManifoldEdge_;
+
+        //- Current edge
+        label edgeID_;
+
+        //- Current direction (face, expressed as index into edgeFaces()[edgeID]
+        label index_;
+
+        //- Point
+        label pointID_;
+
+        //- Starting edge so we know when to stop.
+        label startEdgeID_;
+
+
+    // Private Member Functions
+
+        //- Set to end() iterator
+        inline void setEnd();
+
+        //- Set edgeID_ to be the other edge on the face that uses the
+        //  point.
+        inline void otherEdge(const label cellI);
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        inline patchPointEdgeCirculator
+        (
+            const primitiveFacePatch&,
+            const PackedBoolList& nonManifoldEdge,
+            const label edgeID,
+            const label index,
+            const label pointID
+        );
+
+        //- Construct as copy
+        inline patchPointEdgeCirculator(const patchPointEdgeCirculator&);
+
+
+    // Member Functions
+
+        inline label edgeID() const;
+
+        inline label index() const;
+
+        inline label pointID() const;
+
+        //- Helper: get face according to the index.
+        //  Returns -1 if on end.
+        inline label faceID() const;
+
+        //- Walk back until non-manifold edge (if any) or minimum edge.
+        inline void setCanonical();
+
+
+    // Member Operators
+
+        inline void operator=(const patchPointEdgeCirculator& iter);
+
+        inline bool operator==(const patchPointEdgeCirculator& iter) const;
+
+        inline bool operator!=(const patchPointEdgeCirculator& iter) const;
+
+        //- Step to next face.
+        inline patchPointEdgeCirculator& operator++();
+
+        //- iterator set to the beginning face. For internal edges this is
+        //  the current face. For boundary edges this is the first boundary face
+        //  reached from walking back (i.e. in opposite direction to ++)
+        inline patchPointEdgeCirculator begin() const;
+        inline patchPointEdgeCirculator cbegin() const;
+
+        //- iterator set to beyond the end of the walk.
+        inline const patchPointEdgeCirculator& end() const;
+        inline const patchPointEdgeCirculator& cend() const;
+
+    // Info
+
+        //- Return info proxy.
+        //  Used to print token information to a stream
+        InfoProxy<patchPointEdgeCirculator> info() const
+        {
+            return *this;
+        }
+
+        friend Ostream& operator<<
+        (
+            Ostream&,
+            const InfoProxy<patchPointEdgeCirculator>&
+        );
+};
+
+Ostream& operator<<(Ostream&, const InfoProxy<patchPointEdgeCirculator>&);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "patchPointEdgeCirculatorI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/extrudeToRegionMesh/patchPointEdgeCirculatorI.H

@@ -0,0 +1,309 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::patchPointEdgeCirculator::setEnd()
+{
+    edgeID_ = -1;
+    pointID_ = -1;
+}
+
+
+// Cross face to other edge on point
+void Foam::patchPointEdgeCirculator::otherEdge(const label faceI)
+{
+    const labelList& fEdges = patch_.faceEdges()[faceI];
+    const face& f = patch_.localFaces()[faceI];
+    label fp = findIndex(f, pointID_);
+
+    if (fEdges[fp] == edgeID_)
+    {
+        edgeID_ = fEdges[f.rcIndex(fp)];
+    }
+    else
+    {
+        // Check for now
+        if (fEdges[f.rcIndex(fp)] != edgeID_)
+        {
+            FatalErrorIn("patchPointEdgeCirculator::otherEdge(const label)")
+                << "face:" << faceI
+                << " verts:" << f
+                << " edges:" << fEdges
+                << " looking for other edge than " << edgeID_
+                << abort(FatalError);
+        }
+        edgeID_ = fEdges[fp];
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+//- Construct from components
+Foam::patchPointEdgeCirculator::patchPointEdgeCirculator
+(
+    const primitiveFacePatch& patch,
+    const PackedBoolList& nonManifoldEdge,
+    const label edgeID,
+    const label index,
+
+    const label pointID
+)
+:
+    patch_(patch),
+    nonManifoldEdge_(nonManifoldEdge),
+    edgeID_(edgeID),
+    index_(index),
+    pointID_(pointID),
+    startEdgeID_(edgeID_)
+{
+    if (edgeID_ != -1)
+    {
+        const edge& e = patch_.edges()[edgeID_];
+
+        if (pointID_ != e[0] && pointID_ != e[1])
+        {
+            FatalErrorIn
+            (
+                "patchPointEdgeCirculator::patchPointEdgeCirculator(..)"
+            )   << "edge " << edgeID_ << " verts " << e
+                << " does not contain point " << pointID_ << abort(FatalError);
+        }
+    }
+}
+
+
+//- Construct copy
+Foam::patchPointEdgeCirculator::patchPointEdgeCirculator
+(
+    const patchPointEdgeCirculator& circ
+)
+:
+    patch_(circ.patch_),
+    nonManifoldEdge_(circ.nonManifoldEdge_),
+    edgeID_(circ.edgeID_),
+    index_(circ.index_),
+    pointID_(circ.pointID_),
+    startEdgeID_(circ.startEdgeID_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::label Foam::patchPointEdgeCirculator::edgeID() const
+{
+    return edgeID_;
+}
+
+
+Foam::label Foam::patchPointEdgeCirculator::index() const
+{
+    return index_;
+}
+
+
+Foam::label Foam::patchPointEdgeCirculator::pointID() const
+{
+    return pointID_;
+}
+
+
+Foam::label Foam::patchPointEdgeCirculator::faceID() const
+{
+    if (edgeID_ != -1 && index_ != -1)
+    {
+        return patch_.edgeFaces()[edgeID_][index_];
+    }
+    else
+    {
+        return -1;
+    }
+}
+
+
+void Foam::patchPointEdgeCirculator::operator=
+(
+    const patchPointEdgeCirculator& circ
+)
+{
+    edgeID_ = circ.edgeID_;
+    index_ = circ.index_;
+    pointID_ = circ.pointID_;
+    startEdgeID_ = circ.startEdgeID_;
+}
+
+
+bool Foam::patchPointEdgeCirculator::operator==
+(
+    const patchPointEdgeCirculator& circ
+) const
+{
+    // Do just enough to have setEnd() produce an iterator equal to end().
+    // Could include the direction(i.e. index_) to make sure that two
+    // circulators around same point but in different direction are not equal.
+    return edgeID_ == circ.edgeID_ && pointID_ == circ.pointID_;
+}
+
+
+bool Foam::patchPointEdgeCirculator::operator!=
+(
+    const patchPointEdgeCirculator& circ
+) const
+{
+    return !(*this == circ);
+}
+
+
+void Foam::patchPointEdgeCirculator::setCanonical()
+{
+    if (edgeID_ == -1)
+    {
+        FatalErrorIn("patchPointEdgeCirculator::setCanonical()")
+            << "Already reached end(). Cannot walk any further."
+            << abort(FatalError);
+    }
+
+    label minEdgeID = edgeID_;
+    label minIndex = index_;
+
+    while (true)
+    {
+        if (nonManifoldEdge_[edgeID_])
+        {
+            break;
+        }
+
+        // Step back
+        const labelList& eFaces = patch_.edgeFaces()[edgeID_];
+
+        if (eFaces.size() != 2)
+        {
+            FatalErrorIn("patchPointEdgeCirculator::setCanonical()")
+                << "problem eFaces:" << eFaces << abort(FatalError);
+        }
+
+        label faceI = (index_ == 0 ? eFaces[1] : eFaces[0]);
+
+        // Step to other edge on face
+        otherEdge(faceI);
+
+        // Update index
+        index_ = findIndex(patch_.edgeFaces()[edgeID_], faceI);
+
+        if (edgeID_ < minEdgeID)
+        {
+            minEdgeID = edgeID_;
+            minIndex = index_;
+        }
+
+        if (edgeID_ == startEdgeID_)
+        {
+            edgeID_ = minEdgeID;
+            index_ = minIndex;
+            break;
+        }
+    }
+
+    startEdgeID_ = edgeID_;
+}
+
+
+//- Step to next edge.
+Foam::patchPointEdgeCirculator::patchPointEdgeCirculator&
+Foam::patchPointEdgeCirculator::operator++()
+{
+    if (index_ == -1)
+    {
+        setEnd();
+    }
+    else
+    {
+        // Step to other edge on face
+        label faceI = patch_.edgeFaces()[edgeID_][index_];
+        otherEdge(faceI);
+
+        if (edgeID_ == startEdgeID_)
+        {
+            setEnd();
+        }
+        else if (nonManifoldEdge_[edgeID_])
+        {
+            // Reached non-manifold edge. Cannot walk further.
+            // Mark so it gets set to end next time.
+            index_ = -1;
+        }
+        else
+        {
+            const labelList& eFaces = patch_.edgeFaces()[edgeID_];
+
+            if (eFaces.size() != 2)
+            {
+                FatalErrorIn("patchPointEdgeCirculator:::operator++()")
+                    << "problem eFaces:" << eFaces << abort(FatalError);
+            }
+            // Point to the face that is not faceI
+            index_ = (eFaces[0] != faceI ? 0 : 1);
+        }
+    }
+
+    return *this;
+}
+
+
+Foam::patchPointEdgeCirculator Foam::patchPointEdgeCirculator::begin() const
+{
+    patchPointEdgeCirculator iter(*this);
+    iter.setCanonical();
+
+    return iter;
+}
+
+
+Foam::patchPointEdgeCirculator Foam::patchPointEdgeCirculator::cbegin() const
+{
+    patchPointEdgeCirculator iter(*this);
+    iter.setCanonical();
+
+    return iter;
+}
+
+
+const Foam::patchPointEdgeCirculator& Foam::patchPointEdgeCirculator::end()
+ const
+{
+    return endConstIter;
+}
+
+const Foam::patchPointEdgeCirculator& Foam::patchPointEdgeCirculator::cend()
+ const
+{
+    return endConstIter;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/snappyHexMesh/Make/options

@@ -1,6 +1,6 @@
 EXE_INC = \
     /* -g -DFULLDEBUG -O0 */ \
-    -I$(LIB_SRC)/parallel/decompositionMethods/lnInclude \
+    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/mesh/autoMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
@@ -11,6 +11,7 @@ EXE_INC = \
 EXE_LIBS = \
     -lfiniteVolume \
     -ldecompositionMethods \
+    -L$(FOAM_MPI_LIBBIN) -lptscotchDecomp \
     -lmeshTools \
     -ldynamicMesh \
     -lautoMesh

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

@@ -347,7 +347,7 @@ int main(int argc, char *argv[])
             << "You have selected decomposition method "
             << decomposer.typeName
             << " which is not parallel aware." << endl
-            << "Please select one that is (hierarchical, parMetis)"
+            << "Please select one that is (hierarchical, ptscotch)"
             << exit(FatalError);
     }
 

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/mesh/manipulation/createPatch/createPatch.C

@@ -259,13 +259,34 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
 
 void separateList
 (
-    const vector& separation,
+    const vectorField& separation,
     UList<vector>& field
 )
 {
-    forAll(field, i)
+    if (separation.size() == 1)
     {
-        field[i] += separation;
+        // Single value for all.
+
+        forAll(field, i)
+        {
+            field[i] += separation[0];
+        }
+    }
+    else if (separation.size() == field.size())
+    {
+        forAll(field, i)
+        {
+            field[i] += separation[i];
+        }
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "separateList(const vectorField&, UList<vector>&)"
+        )   << "Sizes of field and transformation not equal. field:"
+            << field.size() << " transformation:" << separation.size()
+            << abort(FatalError);
     }
 }
 
@@ -423,8 +444,7 @@ void syncPoints
             else if (cycPatch.separated())
             {
                 hasTransformation = true;
-                const vector& v = cycPatch.separation();
-                separateList(v, half0Values);
+                separateList(cycPatch.separation(), half0Values);
             }
 
             forAll(coupledPoints, i)
@@ -540,7 +560,7 @@ int main(int argc, char *argv[])
     dumpCyclicMatch("initial_", mesh);
 
     // Read patch construct info from dictionary
-    PtrList<dictionary> patchSources(dict.lookup("patchInfo"));
+    PtrList<dictionary> patchSources(dict.lookup("patches"));
 
 
 
@@ -585,7 +605,7 @@ int main(int argc, char *argv[])
 
             if (destPatchI == -1)
             {
-                dictionary patchDict(dict.subDict("dictionary"));
+                dictionary patchDict(dict.subDict("patchInfo"));
 
                 destPatchI = allPatches.size();
 
@@ -772,7 +792,7 @@ int main(int argc, char *argv[])
                         << " separation[0] was "
                         << cpp.separation()[0] << endl;
 
-                    if (isA<cyclicPolyPatch>(pp))
+                    if (isA<cyclicPolyPatch>(pp) && pp.size())
                     {
                         const cyclicPolyPatch& cycpp =
                             refCast<const cyclicPolyPatch>(pp);
@@ -783,14 +803,44 @@ int main(int argc, char *argv[])
                             Info<< "On cyclic translation patch " << pp.name()
                                 << " forcing uniform separation of "
                                 << cycpp.separationVector() << endl;
-                            const_cast<vector&>(cpp.separation()) =
-                                cycpp.separationVector();
+                            const_cast<vectorField&>(cpp.separation()) =
+                                pointField(1, cycpp.separationVector());
+                        }
+                        else
+                        {
+                            const cyclicPolyPatch& nbr = cycpp.neighbPatch();
+                            const_cast<vectorField&>(cpp.separation()) =
+                                pointField
+                                (
+                                    1,
+                                    nbr[0].centre(mesh.points())
+                                  - cycpp[0].centre(mesh.points())
+                                );
                         }
                     }
                     Info<< "On coupled patch " << pp.name()
                         << " forcing uniform separation of "
                         << cpp.separation() << endl;
                 }
+                else if (!cpp.parallel())
+                {
+                    Info<< "On coupled patch " << pp.name()
+                        << " forcing uniform rotation of "
+                        << cpp.forwardT()[0] << endl;
+
+                    const_cast<tensorField&>
+                    (
+                        cpp.forwardT()
+                    ).setSize(1);
+                    const_cast<tensorField&>
+                    (
+                        cpp.reverseT()
+                    ).setSize(1);
+
+                    Info<< "On coupled patch " << pp.name()
+                        << " forcing uniform rotation of "
+                        << cpp.forwardT() << endl;
+                }
             }
         }
 

applications/utilities/mesh/manipulation/createPatch/createPatchDict

@@ -48,14 +48,14 @@ matchTolerance 1E-3;
 pointSync true;
 
 // Patches to create.
-patchInfo
+patches
 (
     {
         // Name of new patch
         name sidePatches;
 
         // Type of new patch
-        dictionary
+        patchInfo
         {
             type cyclic;
 

applications/utilities/mesh/manipulation/mirrorMesh/mirrorFvMesh.C

@@ -28,20 +28,6 @@ License
 #include "Time.H"
 #include "plane.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-const Foam::label Foam::mirrorFvMesh::cellRenumber[8][8] =
-{
-    {-1, -1, -1, -1, -1, -1, -1, -1},    // unknown
-    {-1, -1, -1, -1, -1, -1, -1, -1},    //
-    {-1, -1, -1, -1, -1, -1, -1, -1},    //
-    { 0,  3,  2,  1,  4,  7,  6,  5},    // hex
-    { 2,  1,  0,  5,  4,  3,  6, -1},    // wedge
-    { 0,  2,  1,  3,  5,  4, -1, -1},    // prism
-    { 0,  3,  2,  1,  4, -1, -1, -1},    // pyramid
-    { 2,  1,  0,  3, -1, -1, -1, -1},    // tet
-};
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::mirrorFvMesh::mirrorFvMesh(const IOobject& io)

applications/utilities/mesh/manipulation/mirrorMesh/mirrorFvMesh.H

@@ -71,12 +71,6 @@ class mirrorFvMesh
 public:
 
 
-    // Static data
-
-        //- Cell renumnering table
-        static const label cellRenumber[8][8];
-
-
     // Constructors
 
         //- Construct from IOobject

applications/utilities/mesh/manipulation/renumberMesh/Make/options

@@ -2,7 +2,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/parallel/decompositionMethods/lnInclude
+    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude
 
 EXE_LIBS = \
     -lmeshTools \

applications/utilities/mesh/manipulation/setSet/setSet.C

@@ -960,7 +960,14 @@ int main(int argc, char *argv[])
 #               endif
             }
 
-            if (rawLine.empty() || rawLine[0] == '#')
+            // Strip off anything after #
+            string::size_type i = rawLine.find_first_of("#");
+            if (i != string::npos)
+            {
+                rawLine = rawLine(0, i);
+            }
+
+            if (rawLine.empty())
             {
                 continue;
             }

applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C

@@ -68,6 +68,104 @@ Description
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+label addPointZone(const polyMesh& mesh, const word& name)
+{
+    label zoneID = mesh.pointZones().findZoneID(name);
+
+    if (zoneID != -1)
+    {
+        Info<< "Reusing existing pointZone "
+            << mesh.pointZones()[zoneID].name()
+            << " at index " << zoneID << endl;
+    }
+    else
+    {
+        pointZoneMesh& pointZones = const_cast<polyMesh&>(mesh).pointZones();
+        zoneID = pointZones.size();
+        Info<< "Adding pointZone " << name << " at index " << zoneID << endl;
+
+        pointZones.setSize(zoneID+1);
+        pointZones.set
+        (
+            zoneID,
+            new pointZone
+            (
+                name,
+                labelList(0),
+                zoneID,
+                pointZones
+            )
+        );
+    }
+    return zoneID;
+}
+
+
+label addFaceZone(const polyMesh& mesh, const word& name)
+{
+    label zoneID = mesh.faceZones().findZoneID(name);
+
+    if (zoneID != -1)
+    {
+        Info<< "Reusing existing faceZone " << mesh.faceZones()[zoneID].name()
+            << " at index " << zoneID << endl;
+    }
+    else
+    {
+        faceZoneMesh& faceZones = const_cast<polyMesh&>(mesh).faceZones();
+        zoneID = faceZones.size();
+        Info<< "Adding faceZone " << name << " at index " << zoneID << endl;
+
+        faceZones.setSize(zoneID+1);
+        faceZones.set
+        (
+            zoneID,
+            new faceZone
+            (
+                name,
+                labelList(0),
+                boolList(),
+                zoneID,
+                faceZones
+            )
+        );
+    }
+    return zoneID;
+}
+
+
+label addCellZone(const polyMesh& mesh, const word& name)
+{
+    label zoneID = mesh.cellZones().findZoneID(name);
+
+    if (zoneID != -1)
+    {
+        Info<< "Reusing existing cellZone " << mesh.cellZones()[zoneID].name()
+            << " at index " << zoneID << endl;
+    }
+    else
+    {
+        cellZoneMesh& cellZones = const_cast<polyMesh&>(mesh).cellZones();
+        zoneID = cellZones.size();
+        Info<< "Adding cellZone " << name << " at index " << zoneID << endl;
+
+        cellZones.setSize(zoneID+1);
+        cellZones.set
+        (
+            zoneID,
+            new cellZone
+            (
+                name,
+                labelList(0),
+                zoneID,
+                cellZones
+            )
+        );
+    }
+    return zoneID;
+}
+
+
 // Checks whether patch present
 void checkPatch(const polyBoundaryMesh& bMesh, const word& name)
 {
@@ -211,29 +309,20 @@ int main(int argc, char *argv[])
     polyTopoChanger stitcher(mesh);
     stitcher.setSize(1);
 
-    DynamicList<pointZone*> pz;
-    DynamicList<faceZone*> fz;
-    DynamicList<cellZone*> cz;
+    mesh.pointZones().clearAddressing();
+    mesh.faceZones().clearAddressing();
+    mesh.cellZones().clearAddressing();
 
     if (perfectCover)
     {
         // Add empty zone for resulting internal faces
-        fz.append
-        (
-            new faceZone
-            (
-                cutZoneName,
-                isf,
-                boolList(masterPatch.size(), false),
-                0,
-                mesh.faceZones()
-            )
-        );
+        label cutZoneID = addFaceZone(mesh, cutZoneName);
 
-        // Note: make sure to add the zones BEFORE constructing polyMeshModifier
-        // (since looks up various zones at construction time)
-        Info<< "Adding point and face zones" << endl;
-        mesh.addZones(pz.shrink(), fz.shrink(), cz.shrink());
+        mesh.faceZones()[cutZoneID].resetAddressing
+        (
+            isf,
+            boolList(masterPatch.size(), false)
+        );
 
         // Add the perfect interface mesh modifier
         stitcher.set
@@ -252,27 +341,15 @@ int main(int argc, char *argv[])
     }
     else
     {
-        pz.append
-        (
-            new pointZone
-            (
-                mergePatchName + "CutPointZone",
-                labelList(0),
-                0,
-                mesh.pointZones()
-            )
-        );
+        label pointZoneID = addPointZone(mesh, mergePatchName + "CutPointZone");
+        mesh.pointZones()[pointZoneID] = labelList(0);
 
-        fz.append
+        label masterZoneID = addFaceZone(mesh, mergePatchName + "MasterZone");
+
+        mesh.faceZones()[masterZoneID].resetAddressing
         (
-            new faceZone
-            (
-                mergePatchName + "MasterZone",
-                isf,
-                boolList(masterPatch.size(), false),
-                0,
-                mesh.faceZones()
-            )
+            isf,
+            boolList(masterPatch.size(), false)
         );
 
         // Slave patch
@@ -284,42 +361,27 @@ int main(int argc, char *argv[])
 
         labelList osf(slavePatch.size());
 
-        forAll(osf, i)
+        forAll (osf, i)
         {
             osf[i] = slavePatch.start() + i;
         }
 
-        fz.append
+        label slaveZoneID = addFaceZone(mesh, mergePatchName + "SlaveZone");
+        mesh.faceZones()[slaveZoneID].resetAddressing
         (
-            new faceZone
-            (
-                mergePatchName + "SlaveZone",
-                osf,
-                boolList(slavePatch.size(), false),
-                1,
-                mesh.faceZones()
-            )
+            osf,
+            boolList(slavePatch.size(), false)
         );
 
         // Add empty zone for cut faces
-        fz.append
+        label cutZoneID = addFaceZone(mesh, cutZoneName);
+        mesh.faceZones()[cutZoneID].resetAddressing
         (
-            new faceZone
-            (
-                cutZoneName,
-                labelList(0),
-                boolList(0, false),
-                2,
-                mesh.faceZones()
-            )
+            labelList(0),
+            boolList(0, false)
         );
 
 
-        // Note: make sure to add the zones BEFORE constructing polyMeshModifier
-        // (since looks up various zones at construction time)
-        Info<< "Adding point and face zones" << endl;
-        mesh.addZones(pz.shrink(), fz.shrink(), cz.shrink());
-
         // Add the sliding interface mesh modifier
         stitcher.set
         (

applications/utilities/parallelProcessing/decomposePar/Make/options

@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I$(LIB_SRC)/parallel/decompositionMethods/lnInclude \
+    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
@@ -7,6 +7,6 @@ EXE_INC = \
 EXE_LIBS = \
     -lfiniteVolume \
     -lgenericPatchFields \
-    -ldecompositionMethods \
+    -ldecompositionMethods -lmetisDecomp -lscotchDecomp \
     -llagrangian \
     -lmeshTools

applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C

@@ -161,15 +161,13 @@ void Foam::domainDecomposition::distributeCells()
         // somewhere in the middle of the domain which might not be anywhere
         // near any of the cells.
 
-        const point greatPoint(GREAT, GREAT, GREAT);
-
-        pointField regionCentres(globalRegion.nRegions(), greatPoint);
+        pointField regionCentres(globalRegion.nRegions(), point::max);
 
         forAll(globalRegion, cellI)
         {
             label regionI = globalRegion[cellI];
 
-            if (regionCentres[regionI] == greatPoint)
+            if (regionCentres[regionI] == point::max)
             {
                 regionCentres[regionI] = cellCentres()[cellI];
             }

applications/utilities/parallelProcessing/decomposePar/domainDecompositionMesh.C

@@ -40,7 +40,6 @@ Description
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-<<<<<<< HEAD
 void Foam::domainDecomposition::append(labelList& lst, const label elem)
 {
     label sz = lst.size();
@@ -413,38 +412,6 @@ void Foam::domainDecomposition::decomposeMesh()
     }
 
 
-    // Set the patch map. No filterPatches allowed.
-    forAll(procBoundaryAddressing_, procI)
-    {
-        label nNormal = procPatchSize_[procI].size();
-
-        const labelListList& curSubPatchIDs = 
-            procProcessorPatchSubPatchIDs_[procI];
-
-//        label nInterProc = procProcessorPatchSize_[procI].size();
-        label nInterProc = 0;
-        forAll(curSubPatchIDs, procPatchI)
-        {
-            nInterProc += curSubPatchIDs[procPatchI].size();
-        }
-
-        procBoundaryAddressing_[procI].setSize(nNormal + nInterProc);
-
-        for (label patchI = 0; patchI < nNormal; patchI++)
-        {
-            procBoundaryAddressing_[procI][patchI] = patchI;
-        }
-        label patchI = nNormal;
-        forAll(curSubPatchIDs, procPatchI)
-        {
-            forAll(curSubPatchIDs[procPatchI], i)
-            {
-                procBoundaryAddressing_[procI][patchI++] =
-                    curSubPatchIDs[procPatchI][i];
-            }
-        }
-    }
-
 //XXXXXXX
 // Print a bit
     forAll(procPatchStartIndex_, procI)
@@ -465,13 +432,6 @@ void Foam::domainDecomposition::decomposeMesh()
     }
     Info<< endl;
 
-    forAll(procBoundaryAddressing_, procI)
-    {
-        Info<< "Processor:" << procI << endl;
-        Info<< "    patchMap:" << procBoundaryAddressing_[procI] << endl;
-    }
-    Info<< endl;
-
     forAll(procNeighbourProcessors_, procI)
     {
         Info<< "Processor " << procI << endl;

applications/utilities/parallelProcessing/reconstructPar/Make/options

@@ -1,7 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I$(LIB_SRC)/parallel/reconstruct/lnInclude
+    -I$(LIB_SRC)/parallel/reconstruct/reconstruct/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/utilities/parallelProcessing/redistributeMeshPar/Make/options

@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I$(LIB_SRC)/parallel/decompositionMethods/lnInclude \
+    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude
@@ -7,5 +7,6 @@ EXE_INC = \
 EXE_LIBS = \
     -lfiniteVolume \
     -ldecompositionMethods \
+    -L$(FOAM_MPI_LIBBIN) -lptscotchDecomp \
     -lmeshTools \
     -ldynamicMesh

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

@@ -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
@@ -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
@@ -109,15 +112,21 @@ int main(int argc, char *argv[])
     argList::addOption
     (
         "patches",
-        "wordList",
-        "specify particular patches to write - eg '(inlet outlet)'. "
+        "wordReList",
+        "specify particular patches to write - eg '(outlet \"inlet.*\")'. "
         "An empty list suppresses writing the internalMesh."
     );
+    argList::addOption
+    (
+        "faceZones",
+        "wordReList",
+        "specify faceZones to write - eg '( slice \"mfp-.*\" )'."
+    );
 
 #   include "setRootCase.H"
 
     // Check options
-    bool binary = !args.optionFound("ascii");
+    const bool binary = !args.optionFound("ascii");
 
 #   include "createTime.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;
         }
 
 

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

@@ -159,13 +159,6 @@ Note
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-
-static const label VTK_TETRA      = 10;
-static const label VTK_PYRAMID    = 14;
-static const label VTK_WEDGE      = 13;
-static const label VTK_HEXAHEDRON = 12;
-
-
 template<class GeoField>
 void print(const char* msg, Ostream& os, const PtrList<GeoField>& flds)
 {
@@ -229,10 +222,7 @@ labelList getSelectedPatches
 }
 
 
-
-
-
-
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Main program:
 
 int main(int argc, char *argv[])
@@ -285,11 +275,11 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 
-    bool doWriteInternal = !args.optionFound("noInternal");
-    bool doFaceZones     = !args.optionFound("noFaceZones");
-    bool doLinks         = !args.optionFound("noLinks");
-    bool binary          = !args.optionFound("ascii");
-    bool useTimeName     = args.optionFound("useTimeName");
+    const bool doWriteInternal = !args.optionFound("noInternal");
+    const bool doFaceZones     = !args.optionFound("noFaceZones");
+    const bool doLinks         = !args.optionFound("noLinks");
+    const bool binary          = !args.optionFound("ascii");
+    const bool useTimeName     = args.optionFound("useTimeName");
 
     if (binary && (sizeof(floatScalar) != 4 || sizeof(label) != 4))
     {
@@ -299,7 +289,7 @@ int main(int argc, char *argv[])
             << exit(FatalError);
     }
 
-    bool nearCellValue = args.optionFound("nearCellValue");
+    const bool nearCellValue = args.optionFound("nearCellValue");
 
     if (nearCellValue)
     {
@@ -308,7 +298,7 @@ int main(int argc, char *argv[])
             << nl << endl;
     }
 
-    bool noPointValues = args.optionFound("noPointValues");
+    const bool noPointValues = args.optionFound("noPointValues");
 
     if (noPointValues)
     {
@@ -316,7 +306,7 @@ int main(int argc, char *argv[])
             << "Outputting cell values only" << nl << endl;
     }
 
-    bool allPatches = args.optionFound("allPatches");
+    const bool allPatches = args.optionFound("allPatches");
 
     List<wordRe> excludePatches;
     if (args.optionFound("excludePatches"))
@@ -392,15 +382,8 @@ int main(int argc, char *argv[])
 
         Info<< "Time: " << runTime.timeName() << endl;
 
-        word timeDesc = "";
-        if (useTimeName)
-        {
-            timeDesc = runTime.timeName();
-        }
-        else
-        {
-            timeDesc = name(runTime.timeIndex());
-        }
+        word timeDesc =
+            useTimeName ? runTime.timeName() : Foam::name(runTime.timeIndex());
 
         // Check for new polyMesh/ and update mesh, fvMeshSubset and cell
         // decomposition.
@@ -470,10 +453,7 @@ int main(int argc, char *argv[])
         IOobjectList objects(mesh, runTime.timeName());
 
         HashSet<word> selectedFields;
-        if (args.optionFound("fields"))
-        {
-            args.optionLookup("fields")() >> selectedFields;
-        }
+        args.optionReadIfPresent("fields", selectedFields);
 
         // Construct the vol fields (on the original mesh if subsetted)
 

applications/utilities/postProcessing/dataConversion/foamToVTK/internalWriter.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

applications/utilities/postProcessing/dataConversion/foamToVTK/internalWriter.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

applications/utilities/postProcessing/dataConversion/foamToVTK/internalWriterTemplates.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

applications/utilities/postProcessing/dataConversion/foamToVTK/lagrangianWriter.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

applications/utilities/postProcessing/dataConversion/foamToVTK/lagrangianWriter.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

applications/utilities/postProcessing/dataConversion/foamToVTK/lagrangianWriterTemplates.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

applications/utilities/postProcessing/dataConversion/foamToVTK/patchWriter.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

applications/utilities/postProcessing/dataConversion/foamToVTK/patchWriter.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

applications/utilities/postProcessing/dataConversion/foamToVTK/patchWriterTemplates.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

applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.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

applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.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

applications/utilities/postProcessing/dataConversion/foamToVTK/vtkMesh.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

applications/utilities/postProcessing/dataConversion/foamToVTK/vtkMesh.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

applications/utilities/postProcessing/dataConversion/foamToVTK/vtkTopo.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
@@ -33,10 +33,10 @@ Description
 #include "polyMesh.H"
 #include "cellShape.H"
 #include "cellModeller.H"
+#include "Swap.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
 :
     mesh_(mesh),
@@ -61,6 +61,9 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
     // Number of additional cells generated by the decomposition of polyhedra
     label nAddCells = 0;
 
+    // face owner is needed to determine the face orientation
+    const labelList& owner = mesh.faceOwner();
+
     // Scan for cells which need to be decomposed and count additional points
     // and cells
 
@@ -71,7 +74,7 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
         if
         (
             model != hex
-//         && model != wedge            // See above.
+            // && model != wedge    // See above.
          && model != prism
          && model != pyr
          && model != tet
@@ -111,7 +114,7 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
     cellTypes_.setSize(cellShapes.size() + nAddCells);
 
     // Set counters for additional points and additional cells
-    label api = 0, aci = 0;
+    label addPointI = 0, addCellI = 0;
 
     forAll(cellShapes, cellI)
     {
@@ -134,15 +137,13 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
         }
         else if (cellModel == prism)
         {
-            vtkVerts.setSize(6);
-            vtkVerts[0] = cellShape[0];
-            vtkVerts[1] = cellShape[2];
-            vtkVerts[2] = cellShape[1];
-            vtkVerts[3] = cellShape[3];
-            vtkVerts[4] = cellShape[5];
-            vtkVerts[5] = cellShape[4];
+            // VTK has a different node order for VTK_WEDGE
+            // their triangles point outwards!
+            vtkVerts = cellShape;
+
+            Foam::Swap(vtkVerts[1], vtkVerts[2]);
+            Foam::Swap(vtkVerts[4], vtkVerts[5]);
 
-            // VTK calls this a wedge.
             cellTypes_[cellI] = VTK_WEDGE;
         }
         else if (cellModel == tetWedge)
@@ -175,34 +176,28 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
 //        }
         else if (cellModel == hex)
         {
-            vtkVerts.setSize(8);
-            vtkVerts[0] = cellShape[0];
-            vtkVerts[1] = cellShape[1];
-            vtkVerts[2] = cellShape[2];
-            vtkVerts[3] = cellShape[3];
-            vtkVerts[4] = cellShape[4];
-            vtkVerts[5] = cellShape[5];
-            vtkVerts[6] = cellShape[6];
-            vtkVerts[7] = cellShape[7];
+            vtkVerts = cellShape;
 
             cellTypes_[cellI] = VTK_HEXAHEDRON;
         }
         else
         {
             // Polyhedral cell. Decompose into tets + prisms.
-            // (see dxFoamExec/createDxConnections.C)
 
             // Mapping from additional point to cell
-            addPointCellLabels_[api] = cellI;
+            addPointCellLabels_[addPointI] = cellI;
+
+            // The new vertex from the cell-centre
+            const label newVertexLabel = mesh_.nPoints() + addPointI;
 
             // Whether to insert cell in place of original or not.
             bool substituteCell = true;
 
             const labelList& cFaces = mesh_.cells()[cellI];
-
             forAll(cFaces, cFaceI)
             {
                 const face& f = mesh_.faces()[cFaces[cFaceI]];
+                const bool isOwner = (owner[cFaces[cFaceI]] == cellI);
 
                 // Number of triangles and quads in decomposition
                 label nTris = 0;
@@ -216,75 +211,95 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
                 label quadi = 0;
                 f.trianglesQuads(mesh_.points(), trii, quadi, triFcs, quadFcs);
 
-                forAll(quadFcs, quadi)
+                forAll(quadFcs, quadI)
                 {
-                    label thisCellI = -1;
+                    label thisCellI;
 
                     if (substituteCell)
                     {
                         thisCellI = cellI;
-
                         substituteCell = false;
                     }
                     else
                     {
-                        thisCellI = mesh_.nCells() + aci;
-
-                        superCells_[aci] = cellI;
-
-                        aci++;
+                        thisCellI = mesh_.nCells() + addCellI;
+                        superCells_[addCellI++] = cellI;
                     }
 
                     labelList& addVtkVerts = vertLabels_[thisCellI];
 
                     addVtkVerts.setSize(5);
 
-                    const face& quad = quadFcs[quadi];
+                    const face& quad = quadFcs[quadI];
 
-                    addVtkVerts[0] = quad[0];
-                    addVtkVerts[1] = quad[1];
-                    addVtkVerts[2] = quad[2];
-                    addVtkVerts[3] = quad[3];
-                    addVtkVerts[4] = mesh_.nPoints() + api;
+                    // Ensure we have the correct orientation for the
+                    // base of the primitive cell shape.
+                    // If the cell is face owner, the orientation needs to be
+                    // flipped.
+                    // At the moment, VTK doesn't actually seem to care if
+                    // negative cells are defined, but we'll do it anyhow
+                    // (for safety).
+                    if (isOwner)
+                    {
+                        addVtkVerts[0] = quad[3];
+                        addVtkVerts[1] = quad[2];
+                        addVtkVerts[2] = quad[1];
+                        addVtkVerts[3] = quad[0];
+                    }
+                    else
+                    {
+                        addVtkVerts[0] = quad[0];
+                        addVtkVerts[1] = quad[1];
+                        addVtkVerts[2] = quad[2];
+                        addVtkVerts[3] = quad[3];
+                    }
+                    addVtkVerts[4] = newVertexLabel;
 
                     cellTypes_[thisCellI] = VTK_PYRAMID;
                 }
 
-                forAll(triFcs, trii)
+                forAll(triFcs, triI)
                 {
-                    label thisCellI = -1;
+                    label thisCellI;
 
                     if (substituteCell)
                     {
                         thisCellI = cellI;
-
                         substituteCell = false;
                     }
                     else
                     {
-                        thisCellI = mesh_.nCells() + aci;
-
-                        superCells_[aci] = cellI;
-
-                        aci++;
+                        thisCellI = mesh_.nCells() + addCellI;
+                        superCells_[addCellI++] = cellI;
                     }
 
 
                     labelList& addVtkVerts = vertLabels_[thisCellI];
 
-                    const face& tri = triFcs[trii];
+                    const face& tri = triFcs[triI];
 
                     addVtkVerts.setSize(4);
-                    addVtkVerts[0] = tri[0];
-                    addVtkVerts[1] = tri[1];
-                    addVtkVerts[2] = tri[2];
-                    addVtkVerts[3] = mesh_.nPoints() + api;
+
+                    // See note above about the orientation.
+                    if (isOwner)
+                    {
+                        addVtkVerts[0] = tri[2];
+                        addVtkVerts[1] = tri[1];
+                        addVtkVerts[2] = tri[0];
+                    }
+                    else
+                    {
+                        addVtkVerts[0] = tri[0];
+                        addVtkVerts[1] = tri[1];
+                        addVtkVerts[2] = tri[2];
+                    }
+                    addVtkVerts[3] = newVertexLabel;
 
                     cellTypes_[thisCellI] = VTK_TETRA;
                 }
             }
 
-            api++;
+            addPointI++;
         }
     }
 

applications/utilities/postProcessing/dataConversion/foamToVTK/vtkTopo.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
@@ -80,7 +80,7 @@ public:
 
    // Public static data
 
-        // this must be consistent with the enumeration in "vtkCell.H"
+        //- equivalent to enumeration in "vtkCellType.h"
         enum vtkTypes
         {
             VTK_TRIANGLE   = 5,
@@ -88,9 +88,10 @@ public:
             VTK_QUAD       = 9,
 
             VTK_TETRA      = 10,
-            VTK_PYRAMID    = 14,
-            VTK_WEDGE      = 13,
             VTK_HEXAHEDRON = 12,
+            VTK_WEDGE      = 13,
+            VTK_PYRAMID    = 14,
+            VTK_POLYHEDRON = 42
         };
 
     // Constructors

applications/utilities/postProcessing/dataConversion/foamToVTK/writeFaceSet.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

applications/utilities/postProcessing/dataConversion/foamToVTK/writeFaceSet.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

applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.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

applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.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

applications/utilities/postProcessing/dataConversion/foamToVTK/writeFunsTemplates.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

applications/utilities/postProcessing/dataConversion/foamToVTK/writePatchGeom.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

applications/utilities/postProcessing/dataConversion/foamToVTK/writePatchGeom.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

applications/utilities/postProcessing/dataConversion/foamToVTK/writePointSet.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

applications/utilities/postProcessing/dataConversion/foamToVTK/writePointSet.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

applications/utilities/postProcessing/dataConversion/foamToVTK/writeSurfFields.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

applications/utilities/postProcessing/dataConversion/foamToVTK/writeSurfFields.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

applications/utilities/postProcessing/graphics/PV3Readers/Allwclean

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wclean libso vtkPV3Readers
+PV3blockMeshReader/Allwclean
+PV3FoamReader/Allwclean
+
+# ----------------------------------------------------------------- end-of-file

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/PV3FoamReader/PV3FoamReader_SM.xml

@@ -91,7 +91,8 @@
       animateable="0">
       <BooleanDomain name="bool"/>
       <Documentation>
-        Use vtkPolyhedron instead of decomposing polyhedra
+        Use vtkPolyhedron instead of decomposing polyhedra.
+        !!Actually uses vtkConvexPointSet until this is properly supported in VTK!!
       </Documentation>
     </IntVectorProperty>
 

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/PV3FoamReader/pqPV3FoamReaderPanel.cxx

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/PV3FoamReader/pqPV3FoamReaderPanel.h

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/Make/options

@@ -2,12 +2,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I$(ParaView_DIR)/VTK \
-    -I$(ParaView_INST_DIR) \
-    -I$(ParaView_INST_DIR)/VTK \
-    -I$(ParaView_INST_DIR)/VTK/Common \
-    -I$(ParaView_INST_DIR)/VTK/Filtering \
-    -I$(ParaView_INST_DIR)/VTK/Rendering \
+    -I$(ParaView_DIR)/include/paraview-$(ParaView_MAJOR) \
     -I../../vtkPV3Readers/lnInclude \
     -I../PV3FoamReader
 

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H

@@ -356,8 +356,15 @@ class vtkPV3Foam
             //- Zone info
             void updateInfoZones(vtkDataArraySelection*);
 
-            //- Read zone names for zoneType from file
-            wordList readZoneNames(const word& zoneType);
+            //- Get non-empty zone names for zoneType from file
+            wordList getZoneNames(const word& zoneType) const;
+
+            //- Get non-empty zone names from mesh info
+            template<class ZoneType>
+            wordList getZoneNames
+            (
+                const ZoneMesh<ZoneType, polyMesh>&
+            ) const;
 
             //- Add objects of Type to paraview array selection
             template<class Type>

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshVolume.C

@@ -31,6 +31,7 @@ License
 #include "fvMesh.H"
 #include "cellModeller.H"
 #include "vtkOpenFOAMPoints.H"
+#include "Swap.H"
 
 // VTK includes
 #include "vtkCellArray.h"
@@ -45,6 +46,13 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
     polyDecomp& decompInfo
 )
 {
+    const cellModel& tet = *(cellModeller::lookup("tet"));
+    const cellModel& pyr = *(cellModeller::lookup("pyr"));
+    const cellModel& prism = *(cellModeller::lookup("prism"));
+    const cellModel& wedge = *(cellModeller::lookup("wedge"));
+    const cellModel& tetWedge = *(cellModeller::lookup("tetWedge"));
+    const cellModel& hex = *(cellModeller::lookup("hex"));
+
     vtkUnstructuredGrid* vtkmesh = vtkUnstructuredGrid::New();
 
     if (debug)
@@ -53,36 +61,27 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
         printMemory();
     }
 
+    const cellShapeList& cellShapes = mesh.cellShapes();
+
     // Number of additional points needed by the decomposition of polyhedra
     label nAddPoints = 0;
 
     // Number of additional cells generated by the decomposition of polyhedra
     label nAddCells = 0;
 
+    // face owner is needed to determine the face orientation
+    const labelList& owner = mesh.faceOwner();
+
     labelList& superCells = decompInfo.superCells();
     labelList& addPointCellLabels = decompInfo.addPointCellLabels();
 
-    const cellModel& tet = *(cellModeller::lookup("tet"));
-    const cellModel& pyr = *(cellModeller::lookup("pyr"));
-    const cellModel& prism = *(cellModeller::lookup("prism"));
-    const cellModel& wedge = *(cellModeller::lookup("wedge"));
-    const cellModel& tetWedge = *(cellModeller::lookup("tetWedge"));
-    const cellModel& hex = *(cellModeller::lookup("hex"));
-
-    // Scan for cells which need to be decomposed and count additional points
-    // and cells
-    if (debug)
-    {
-        Info<< "... building cell-shapes" << endl;
-    }
-    const cellShapeList& cellShapes = mesh.cellShapes();
-
     if (debug)
     {
         Info<< "... scanning" << endl;
     }
 
-    // count number of cells to decompose
+    // Scan for cells which need to be decomposed and count additional points
+    // and cells
     if (!reader_->GetUseVTKPolyhedron())
     {
         forAll(cellShapes, cellI)
@@ -105,10 +104,10 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
                 {
                     const face& f = mesh.faces()[cFaces[cFaceI]];
 
-                    label nFacePoints = f.size();
+                    label nQuads = 0;
+                    label nTris = 0;
+                    f.nTrianglesQuads(mesh.points(), nTris, nQuads);
 
-                    label nQuads = (nFacePoints - 2)/2;
-                    label nTris = (nFacePoints - 2)%2;
                     nAddCells += nQuads + nTris;
                 }
 
@@ -201,8 +200,8 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
         }
         else if (cellModel == prism)
         {
-            // VTK has a different node order - their triangles point outwards!
-
+            // VTK has a different node order for VTK_WEDGE
+            // their triangles point outwards!
             nodeIds[0] = cellShape[0];
             nodeIds[1] = cellShape[2];
             nodeIds[2] = cellShape[1];
@@ -349,29 +348,34 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
             // Mapping from additional point to cell
             addPointCellLabels[addPointI] = cellI;
 
-            // Insert the new vertex from the cell-centre
-            label newVertexLabel = mesh.nPoints() + addPointI;
+            // The new vertex from the cell-centre
+            const label newVertexLabel = mesh.nPoints() + addPointI;
             vtkInsertNextOpenFOAMPoint(vtkpoints, mesh.C()[cellI]);
 
             // Whether to insert cell in place of original or not.
             bool substituteCell = true;
 
             const labelList& cFaces = mesh.cells()[cellI];
-
             forAll(cFaces, cFaceI)
             {
                 const face& f = mesh.faces()[cFaces[cFaceI]];
+                const bool isOwner = (owner[cFaces[cFaceI]] == cellI);
 
-                label nFacePoints = f.size();
+                // Number of triangles and quads in decomposition
+                label nTris = 0;
+                label nQuads = 0;
+                f.nTrianglesQuads(mesh.points(), nTris, nQuads);
 
-                label nQuads = (nFacePoints - 2)/2;
-                label nTris = (nFacePoints - 2)%2;
+                // Do actual decomposition into triFcs and quadFcs.
+                faceList triFcs(nTris);
+                faceList quadFcs(nQuads);
+                label trii = 0;
+                label quadi = 0;
+                f.trianglesQuads(mesh.points(), trii, quadi, triFcs, quadFcs);
 
-                label qpi = 0;
-
-                for (label quadi=0; quadi<nQuads; quadi++)
+                forAll(quadFcs, quadI)
                 {
-                    label thisCellI = -1;
+                    label thisCellI;
 
                     if (substituteCell)
                     {
@@ -384,10 +388,29 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
                         superCells[addCellI++] = cellI;
                     }
 
-                    nodeIds[0] = f[0];
-                    nodeIds[1] = f[qpi + 1];
-                    nodeIds[2] = f[qpi + 2];
-                    nodeIds[3] = f[qpi + 3];
+                    const face& quad = quadFcs[quadI];
+
+                    // Ensure we have the correct orientation for the
+                    // base of the primitive cell shape.
+                    // If the cell is face owner, the orientation needs to be
+                    // flipped.
+                    // At the moment, VTK doesn't actually seem to care if
+                    // negative cells are defined, but we'll do it anyhow
+                    // (for safety).
+                    if (isOwner)
+                    {
+                        nodeIds[0] = quad[3];
+                        nodeIds[1] = quad[2];
+                        nodeIds[2] = quad[1];
+                        nodeIds[3] = quad[0];
+                    }
+                    else
+                    {
+                        nodeIds[0] = quad[0];
+                        nodeIds[1] = quad[1];
+                        nodeIds[2] = quad[2];
+                        nodeIds[3] = quad[3];
+                    }
                     nodeIds[4] = newVertexLabel;
                     vtkmesh->InsertNextCell
                     (
@@ -395,13 +418,11 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
                         5,
                         nodeIds
                     );
-
-                    qpi += 2;
                 }
 
-                if (nTris)
+                forAll(triFcs, triI)
                 {
-                    label thisCellI = -1;
+                    label thisCellI;
 
                     if (substituteCell)
                     {
@@ -414,10 +435,23 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
                         superCells[addCellI++] = cellI;
                     }
 
-                    nodeIds[0] = f[0];
-                    nodeIds[1] = f[qpi + 1];
-                    nodeIds[2] = f[qpi + 2];
+                    const face& tri = triFcs[triI];
+
+                    // See note above about the orientation.
+                    if (isOwner)
+                    {
+                        nodeIds[0] = tri[2];
+                        nodeIds[1] = tri[1];
+                        nodeIds[2] = tri[0];
+                    }
+                    else
+                    {
+                        nodeIds[0] = tri[0];
+                        nodeIds[1] = tri[1];
+                        nodeIds[2] = tri[2];
+                    }
                     nodeIds[3] = newVertexLabel;
+
                     vtkmesh->InsertNextCell
                     (
                         VTK_TETRA,

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C

@@ -82,9 +82,31 @@ public:
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::wordList Foam::vtkPV3Foam::readZoneNames(const word& zoneType)
+template<class ZoneType>
+Foam::wordList Foam::vtkPV3Foam::getZoneNames
+(
+    const ZoneMesh<ZoneType, polyMesh>& zmesh
+) const
 {
-    wordList zoneNames;
+    wordList names(zmesh.size());
+    label nZone = 0;
+
+    forAll(zmesh, zoneI)
+    {
+        if (zmesh[zoneI].size())
+        {
+            names[nZone++] = zmesh[zoneI].name();
+        }
+    }
+    names.setSize(nZone);
+
+    return names;
+}
+
+
+Foam::wordList Foam::vtkPV3Foam::getZoneNames(const word& zoneType) const
+{
+    wordList names;
 
     // mesh not loaded - read from file
     IOobject ioObj
@@ -107,14 +129,14 @@ Foam::wordList Foam::vtkPV3Foam::readZoneNames(const word& zoneType)
     {
         zonesEntries zones(ioObj);
 
-        zoneNames.setSize(zones.size());
+        names.setSize(zones.size());
         forAll(zones, zoneI)
         {
-            zoneNames[zoneI] = zones[zoneI].keyword();
+            names[zoneI] = zones[zoneI].keyword();
         }
     }
 
-    return zoneNames;
+    return names;
 }
 
 
@@ -312,11 +334,11 @@ void Foam::vtkPV3Foam::updateInfoZones
     // ~~~~~~~~~~~~~~~~~~~~~
     if (meshPtr_)
     {
-        namesLst = meshPtr_->cellZones().names();
+        namesLst = getZoneNames(meshPtr_->cellZones());
     }
     else
     {
-        namesLst = readZoneNames("cellZones");
+        namesLst = getZoneNames("cellZones");
     }
 
     arrayRangeCellZones_.reset(arraySelection->GetNumberOfArrays());
@@ -335,11 +357,11 @@ void Foam::vtkPV3Foam::updateInfoZones
     // ~~~~~~~~~~~~~~~~~~~~~
     if (meshPtr_)
     {
-        namesLst = meshPtr_->faceZones().names();
+        namesLst = getZoneNames(meshPtr_->faceZones());
     }
     else
     {
-        namesLst = readZoneNames("faceZones");
+        namesLst = getZoneNames("faceZones");
     }
 
     arrayRangeFaceZones_.reset(arraySelection->GetNumberOfArrays());
@@ -358,11 +380,11 @@ void Foam::vtkPV3Foam::updateInfoZones
     // ~~~~~~~~~~~~~~~~~~~~~~
     if (meshPtr_)
     {
-        namesLst = meshPtr_->pointZones().names();
+        namesLst = getZoneNames(meshPtr_->pointZones());
     }
     else
     {
-        namesLst = readZoneNames("pointZones");
+        namesLst = getZoneNames("pointZones");
     }
 
     arrayRangePointZones_.reset(arraySelection->GetNumberOfArrays());

applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/Make/options

@@ -1,12 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/mesh/blockMesh/lnInclude \
-    -I$(ParaView_DIR)/VTK \
-    -I$(ParaView_INST_DIR) \
-    -I$(ParaView_INST_DIR)/VTK \
-    -I$(ParaView_INST_DIR)/VTK/Common \
-    -I$(ParaView_INST_DIR)/VTK/Filtering \
-    -I$(ParaView_INST_DIR)/VTK/Rendering \
+    -I$(ParaView_DIR)/include/paraview-$(ParaView_MAJOR) \
     -I../../vtkPV3Readers/lnInclude \
     -I../PV3blockMeshReader
 

applications/utilities/postProcessing/graphics/PV3Readers/vtkPV3Readers/Make/options

@@ -1,10 +1,5 @@
 EXE_INC = \
-    -I$(ParaView_DIR)/VTK \
-    -I$(ParaView_INST_DIR) \
-    -I$(ParaView_INST_DIR)/VTK \
-    -I$(ParaView_INST_DIR)/VTK/Common \
-    -I$(ParaView_INST_DIR)/VTK/Filtering \
-    -I$(ParaView_INST_DIR)/VTK/Rendering
+    -I$(ParaView_DIR)/include/paraview-$(ParaView_MAJOR)
 
 LIB_LIBS = \
     $(GLIBS)

bin/doxyFilt

@@ -3,7 +3,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

bin/engridFoam

@@ -3,7 +3,7 @@
 # =========                 |
 # \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
 #  \\    /   O peration     |
-#   \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+#   \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
 #    \\/     M anipulation  |
 #-------------------------------------------------------------------------------
 # License
@@ -38,7 +38,7 @@ Usage: ${0##*/} [OPTION]
 options:
   -help
 
-* start engrid using the paraview-$ParaView_VERSION libraries
+* start engrid using the paraview-$ParaView_MAJOR libraries
   passes through engrid options unmodified
 
 USAGE
@@ -48,13 +48,8 @@ USAGE
 # report usage
 [ "$1" = "-h" -o "$1" = "-help" ] && usage
 
-# set the major version "<digits>.<digits>"
-ParaView_MAJOR_VERSION=$(echo $ParaView_VERSION | \
-    sed -e 's/^\([0-9][0-9]*\.[0-9][0-9]*\).*$/\1/')
-
-
-bindir=$WM_THIRD_PARTY_DIR/engrid/platforms/$WM_ARCH
-libdir="$ParaView_DIR/lib/paraview-$ParaView_MAJOR_VERSION"
+bindir=$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/engrid
+libdir="$ParaView_DIR/lib/paraview-${ParaView_MAJOR:-unknown}"
 
 [ -x $bindir/engrid ] || usage "engrid executable not found in $bindir"
 [ -d $libdir ] || usage "paraview libraries not found"

bin/foamAllHC

@@ -3,7 +3,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

bin/foamCheckJobs

@@ -3,7 +3,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