COMP: merge conflict

ReleaseNotes-dev

@@ -0,0 +1,130 @@
+#                            -*- mode: org; -*-
+#
+#+TITLE:  OpenFOAM release notes for version dev
+#+AUTHOR:                      OpenCFD Ltd.
+#+DATE:                            TBA
+#+LINK:                  http://www.openfoam.com
+#+OPTIONS: author:nil ^:{}
+# Copyright (c) 2010 OpenCFD Ltd.
+
+* Overview
+  OpenFOAM-dev is the latest major release of OpenFOAM including many new
+  developments a number of bug-fixes.  This release passes our standard tests
+  and the tutorials have been broadly checked.  Please report any bugs by
+  following the link: http://www.openfoam.com/bugs.
+
+* GNU/Linux version
+  This release of OpenFOAM is distributed primarily in 2 ways: (1) as a Debian
+  pack containing binaries and source; (2) from the SourceForge source code
+  repository (see [[./README.org][README]]).
+
+  The Ubuntu/Debian pack is available for 32 and 64 bit versions of the 10.04
+  LTS operating system using the system compiler and libraries that will be
+  installed automatically from standard Debian packs.
+
+  To use the source version from the SourceForge repository, we provide a source
+  pack of third-party packages that can be compiled on the user's system.  This
+  does not include =gcc=, since the system installed version is typically
+  sufficient, but includes =paraview-3.8.0=, =openmpi-1.4.1=, =scotch_5.1=,
+  =metis-5.0pre2=, =ParMetis-3.1= and =ParMGridGen-1.0=.
+
+* Library developments
+  There have been a number of developments to the libraries to support the
+  extension of functionality in solver and utility applications.
+*** Core library
+    + Large number of code refinements and consistency improvements to support
+      other developments.
+*** Turbulence modelling
+*** Thermo-physical Models
+*** DSMC
+*** Dynamic Mesh
+*** Numerics
+*** *Updated* command line help, e.g. `snappyHexMesh -help' now gives:
+
+    Usage: snappyHexMesh [OPTIONS]
+    options:
+      -case <dir>       specify alternate case directory, default is the cwd
+      -overwrite        overwrite existing mesh/results files
+      -parallel         run in parallel
+      -srcDoc           display source code in browser
+      -doc              display application documentation in browser
+      -help             print the usage
+
+*** *New* Surface film library
+    + Creation of films by particle addition, or initial film distribution
+    + Coupled with the lagrangian/intermediate cloud hierarchy library
+    + Hierarchical design, consisting of
+      + kinematic film: mass, momentum
+        + constant thermodynamic properties
+      + thermodynamic film: mass, momentum and enthalpy
+        + constant, or temperature dependant thermodynamic properties
+    + Sub-models:
+      + detachment/dripping whereby particles (re)enter the originating cloud
+        + particle sizes set according to PDF
+        + other properties set to ensure mass, momentum and energy conservation
+      + heat transfer to/from walls and film surface
+      + film evaporation and boiling
+    + Additional wall functions for primary region momentum and temperature
+      taking film into account
+    + Parallel aware
+*** *Updated* particle tracking algorithm
+
+* Solvers
+  A number of new solvers have been developed for a range of engineering
+  applications.  There has been a set of improvements to certain classes of
+  solver that are introduced in this release.
+*** *New* Solvers
+    + =reactingParcelFilmFoam=: Lagrangian cloud and film transport in a
+      reacting gas phase system
+*** Modifications to multiphase and buoyant solvers
+    + ...
+*** Modifications to solvers for sensible enthalpy
+    + ...
+*** Modifications to steady-state compressible solvers
+    + ...
+*** Other modifications
+    + ...
+
+* Boundary conditions
+  New boundary conditions have been introduced to support new applications in
+  OpenFOAM.
+  + *New* wall functions:
+    + kappatJayatillekeWallFunction: incompressible RAS thermal wall function
+
+* Utilities
+  There have been some utilities added and updated in this release.
+*** *New* utilities
+    + =extrudeToRegionMesh=: Extrude faceZones into separate mesh (as a
+      different region)
+      + used to e.g. extrude baffles (extrude internal faces) or create
+        liquid film regions
+      + if extruding internal faces:
+         + create baffles in original mesh with directMappedWall patches
+      + if extruding boundary faces:
+        + convert boundary faces to directMappedWall patches
+      + extrude edges of faceZone as a <zone>_sidePatch
+      + extrude edges inbetween different faceZones as a
+        (nonuniformTransform)cyclic <zoneA>_<zoneB>
+      + extrudes into master direction (i.e. away from the owner cell
+        if flipMap is false)
+*** Updated utilities
+    + ...
+
+* Post-processing
+  + =foamToEnsight=: new =-nodeValues= option to generate and output nodal
+    field data.
+  + Function objects:
+    + =residualControl=: new function object to allow users to terminate steady
+    state calculations when the defined residual levels are achieved
+    + =abortCalculation=: watches for presence of the named file in the
+    $FOAM_CASE directory and aborts the calculation if it is present
+    + =timeActivatedFileUpdate=: performs a file copy/replacement once a
+    specified time has been reached, e.g. to automagically change fvSchemes and
+    fvSolution during a calculation
+    + =streamLine=: generate streamlines; ouputs both trajectory and field data
+
+* New tutorials
+  There is a large number of new tutorials to support the new solvers in the
+  release.
+  + =reactingParcelFilmFoam= tutorials:
+    + multipleBoxes, hotBoxes, panel, evaporationTest

applications/solvers/combustion/fireFoam/UEqn.H

@@ -13,7 +13,9 @@ solve
   ==
     fvc::reconstruct
     (
-        fvc::interpolate(rho)*(g & mesh.Sf())
+        (
-      - fvc::snGrad(p)*mesh.magSf()
+          - ghf*fvc::snGrad(rho)
+          - fvc::snGrad(p_rgh)
+        )*mesh.magSf()
     )
 );

applications/solvers/combustion/fireFoam/combustionModels/Make/files

@@ -1,5 +1,5 @@
 combustionModel/combustionModel.C
-combustionModel/combustionModelNew.C
+combustionModel/newCombustionModel.C
 
 infinitelyFastChemistry/infinitelyFastChemistry.C
 

applications/solvers/combustion/fireFoam/combustionModels/combustionModel/newCombustionModel.C

@@ -29,19 +29,22 @@ License
 
 Foam::autoPtr<Foam::combustionModel> Foam::combustionModel::New
 (
-    const dictionary& propDict,
+    const dictionary& combustionProperties,
     const hsCombustionThermo& thermo,
     const compressible::turbulenceModel& turbulence,
     const surfaceScalarField& phi,
     const volScalarField& rho
 )
 {
-    const word modelType(propDict.lookup("combustionModel"));
+    word combustionModelTypeName = combustionProperties.lookup
+    (
+        "combustionModel"
+    );
 
-    Info<< "Selecting combustion model " << modelType << endl;
+    Info<< "Selecting combustion model " << combustionModelTypeName << endl;
 
     dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(modelType);
+        dictionaryConstructorTablePtr_->find(combustionModelTypeName);
 
     if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
@@ -49,14 +52,14 @@ Foam::autoPtr<Foam::combustionModel> Foam::combustionModel::New
         (
             "combustionModel::New"
         )   << "Unknown combustionModel type "
-            << modelType << nl << nl
+            << combustionModelTypeName << endl << endl
-            << "Valid combustionModels are : " << endl
+            << "Valid  combustionModels are : " << endl
-            << dictionaryConstructorTablePtr_->toc()
+            << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }
 
     return autoPtr<combustionModel>
-        (cstrIter()(propDict, thermo, turbulence, phi, rho));
+        (cstrIter()(combustionProperties, thermo, turbulence, phi, rho));
 }
 
 

applications/solvers/combustion/fireFoam/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.H

@@ -88,8 +88,9 @@ public:
         );
 
 
-    //- Destructor
+    // Destructor
-    virtual ~infinitelyFastChemistry();
+
+        virtual ~infinitelyFastChemistry();
 
 
     // Member Functions

applications/solvers/combustion/fireFoam/combustionModels/noCombustion/noCombustion.H

@@ -83,8 +83,9 @@ public:
         );
 
 
-    //- Destructor
+    // Destructor
-    virtual ~noCombustion();
+
+        virtual ~noCombustion();
 
     // Member Functions
 

applications/solvers/combustion/fireFoam/createFields.H

@@ -35,6 +35,7 @@ const volScalarField& psi = thermo.psi();
 volScalarField& ft = composition.Y("ft");
 volScalarField& fu = composition.Y("fu");
 
+
 Info<< "Reading field U\n" << endl;
 
 volVectorField U
@@ -73,7 +74,7 @@ IOdictionary combustionProperties
 Info<< "Creating combustion model\n" << endl;
 autoPtr<combustionModel> combustion
 (
-    combustionModel::New
+    combustionModel::combustionModel::New
     (
         combustionProperties,
         thermo,
@@ -83,6 +84,29 @@ autoPtr<combustionModel> combustion
     )
 );
 
+
+Info<< "Calculating field g.h\n" << endl;
+volScalarField gh("gh", g & mesh.C());
+surfaceScalarField ghf("gh", g & mesh.Cf());
+
+Info<< "Reading field p_rgh\n" << endl;
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Force p_rgh to be consistent with p
+p_rgh = p - rho*gh;
+
+
 volScalarField dQ
 (
     IOobject
@@ -103,15 +127,6 @@ volScalarField DpDt =
     fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
 
-Info<< "Calculating field g.h\n" << endl;
-volScalarField gh("gh", g & mesh.C());
-
-surfaceScalarField ghf("gh", g & mesh.Cf());
-
-p += rho*gh;
-
-thermo.correct();
-
 dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 
 

applications/solvers/combustion/fireFoam/pEqn.H

@@ -1,5 +1,3 @@
-bool closedVolume = false;
-
 rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
@@ -15,43 +13,41 @@ surfaceScalarField phiU
     )
 );
 
-phi = phiU + rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
+phi = phiU - rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
-       surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA);
+    surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA);
 
-        fvScalarMatrix pEqn
+    fvScalarMatrix p_rghEqn
+    (
+        fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh
+      + fvc::div(phi)
+      - fvm::laplacian(rhorUAf, p_rgh)
+    );
+
+    p_rghEqn.solve
+    (
+        mesh.solver
         (
-            fvm::ddt(psi,p)
+            p_rgh.select
-          + fvc::div(phi)
-          - fvm::laplacian(rhorUAf, p)
-        );
-
-        closedVolume = p.needReference();
-
-        pEqn.solve
-        (
-            mesh.solver
             (
-                p.select
                 (
-                    (
+                    finalIter
-                        finalIter
+                 && corr == nCorr-1
-                     && corr == nCorr-1
+                 && nonOrth == nNonOrthCorr
-                     && nonOrth == nNonOrthCorr
-                    )
                 )
             )
-        );
+        )
+    );
 
-        if (nonOrth == nNonOrthCorr)
+    if (nonOrth == nNonOrthCorr)
-        {
+    {
-            phi += pEqn.flux();
+        phi += p_rghEqn.flux();
-        }
+    }
 }
 
-DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+p = p_rgh + rho*gh;
 
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
@@ -59,12 +55,4 @@ DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
 U.correctBoundaryConditions();
 
-// For closed-volume cases adjust the pressure and density levels
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
-// to obey overall mass continuity
-if (closedVolume)
-{
-    p +=
-        (initialMass - fvc::domainIntegrate(thermo.psi()*p))
-       /fvc::domainIntegrate(thermo.psi());
-    rho = thermo.rho();
-}

applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H

@@ -39,9 +39,14 @@
 
     #include "compressibleCreatePhi.H"
 
-    dimensionedScalar pMin
+    dimensionedScalar rhoMax
     (
-        mesh.solutionDict().subDict("PIMPLE").lookup("pMin")
+        mesh.solutionDict().subDict("PIMPLE").lookup("rhoMax")
+    );
+
+    dimensionedScalar rhoMin
+    (
+        mesh.solutionDict().subDict("PIMPLE").lookup("rhoMin")
     );
 
     Info<< "Creating turbulence model\n" << endl;

applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H

@@ -102,6 +102,8 @@ else
     p.relax();
 
     rho = thermo.rho();
+    rho = max(rho, rhoMin);
+    rho = min(rho, rhoMax);
     rho.relax();
     Info<< "rho max/min : " << max(rho).value()
         << " " << min(rho).value() << endl;
@@ -112,8 +114,6 @@ 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
 /*

applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H

@@ -70,16 +70,29 @@
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 
     thermalPorousZones pZones(mesh);
+    Switch pressureImplicitPorosity(false);
 
     // nUCorrectors used for pressureImplicitPorosity
     int nUCorr = 0;
-    const bool pressureImplicitPorosity =
+    if (pZones.size())
-    (
+    {
-        pZones.size()
+        // nUCorrectors for pressureImplicitPorosity
-     && mesh.solutionDict().subDict("SIMPLE").readIfPresent
+        if (mesh.solutionDict().subDict("SIMPLE").found("nUCorrectors"))
-        (
+        {
-            "nUCorrectors",
+            nUCorr = readInt
-            nUCorr
+            (
-        )
+                mesh.solutionDict().subDict("SIMPLE").lookup("nUCorrectors")
-     && (nUCorr > 0)
+            );
-    );
+        }
+
+        if (nUCorr > 0)
+        {
+            pressureImplicitPorosity = true;
+            Info<< "Using pressure implicit porosity" << endl;
+        }
+        else
+        {
+            Info<< "Using pressure explicit porosity" << endl;
+        }
+    }
+

applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H

@@ -5,8 +5,8 @@
       - fvm::Sp(fvc::div(phi), h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)")
+        fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p, "div(U,p)"))
-      - (rho/psi)*fvc::div(phi/fvc::interpolate(rho))
+      - p*fvc::div(phi/fvc::interpolate(rho))
     );
 
     pZones.addEnthalpySource(thermo, rho, hEqn);

applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H

@@ -1,8 +1,3 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
 if (pressureImplicitPorosity)
 {
     U = trTU()&UEqn().H();

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -30,15 +30,7 @@ if (transonic)
 
         pEqn.setReference(pRefCell, pRefValue);
 
-        // Retain the residual from the first iteration
+        pEqn.solve();
-        if (nonOrth == 0)
-        {
-            pEqn.solve();
-        }
-        else
-        {
-            pEqn.solve();
-        }
 
         if (nonOrth == nNonOrthCorr)
         {
@@ -60,15 +52,7 @@ else
 
         pEqn.setReference(pRefCell, pRefValue);
 
-        // Retain the residual from the first iteration
+        pEqn.solve();
-        if (nonOrth == 0)
-        {
-            pEqn.solve();
-        }
-        else
-        {
-            pEqn.solve();
-        }
 
         if (nonOrth == nNonOrthCorr)
         {

applications/solvers/compressible/rhoSimplecFoam/createFields.H

@@ -43,9 +43,14 @@
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
 
-    dimensionedScalar pMin
+    dimensionedScalar rhoMax
     (
-        mesh.solutionDict().subDict("SIMPLE").lookup("pMin")
+        mesh.solutionDict().subDict("SIMPLE").lookup("rhoMax")
+    );
+
+    dimensionedScalar rhoMin
+    (
+        mesh.solutionDict().subDict("SIMPLE").lookup("rhoMin")
     );
 
     Info<< "Creating turbulence model\n" << endl;

applications/solvers/compressible/rhoSimplecFoam/pEqn.H

@@ -101,8 +101,6 @@ U -= (fvc::grad(p0)*(1.0/AU - 1.0/AtU) + fvc::grad(p)/AtU);
 
 U.correctBoundaryConditions();
 
-bound(p, pMin);
-
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
 if (closedVolume)
@@ -112,6 +110,8 @@ if (closedVolume)
 }
 
 rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
 
 if (!transonic)
 {

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H

@@ -12,7 +12,7 @@
     );
 
     TEqn.relax();
-    TEqn.solve();
+    TEqn.solve(mesh.solver(T.select(finalIter)));
 
     rhok = 1.0 - beta*(T - TRef);
 }

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H

@@ -18,9 +18,10 @@
             fvc::reconstruct
             (
                 (
-                    fvc::interpolate(rhok)*(g & mesh.Sf())
+                  - ghf*fvc::snGrad(rhok)
-                  - fvc::snGrad(p)*mesh.magSf()
+                  - fvc::snGrad(p_rgh)
-                )
+                )*mesh.magSf()
-            )
+            ),
+            mesh.solver(U.select(finalIter))
         );
     }

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C

@@ -87,7 +87,7 @@ int main(int argc, char *argv[])
 
             if (nOuterCorr != 1)
             {
-                p.storePrevIter();
+                p_rgh.storePrevIter();
             }
 
             #include "UEqn.H"

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H

@@ -14,12 +14,12 @@
         mesh
     );
 
-    Info<< "Reading field p\n" << endl;
+    Info<< "Reading field p_rgh\n" << endl;
-    volScalarField p
+    volScalarField p_rgh
     (
         IOobject
         (
-            "p",
+            "p_rgh",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -52,6 +52,18 @@
         incompressible::RASModel::New(U, phi, laminarTransport)
     );
 
+    // Kinematic density for buoyancy force
+    volScalarField rhok
+    (
+        IOobject
+        (
+            "rhok",
+            runTime.timeName(),
+            mesh
+        ),
+        1.0 - beta*(T - TRef)
+    );
+
     // kinematic turbulent thermal thermal conductivity m2/s
     Info<< "Reading field kappat\n" << endl;
     volScalarField kappat
@@ -67,25 +79,41 @@
         mesh
     );
 
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        p_rgh + rhok*gh
+    );
+
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell
     (
         p,
+        p_rgh,
         mesh.solutionDict().subDict("PIMPLE"),
         pRefCell,
         pRefValue
     );
 
-
+    if (p_rgh.needReference())
-    // Kinematic density for buoyancy force
+    {
-    volScalarField rhok
+        p += dimensionedScalar
-    (
-        IOobject
         (
-            "rhok",
+            "p",
-            runTime.timeName(),
+            p.dimensions(),
-            mesh
+            pRefValue - getRefCellValue(p, pRefCell)
-        ),
+        );
-        1.0 - beta*(T - TRef)
+    }
-    );
+

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H

@@ -7,22 +7,23 @@
     phi = (fvc::interpolate(U) & mesh.Sf())
         + fvc::ddtPhiCorr(rUA, U, phi);
 
-    surfaceScalarField buoyancyPhi =
+    surfaceScalarField buoyancyPhi = rUAf*ghf*fvc::snGrad(rhok)*mesh.magSf();
-        rUAf*fvc::interpolate(rhok)*(g & mesh.Sf());
+    phi -= buoyancyPhi;
-    phi += buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pEqn
+        fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rUAf, p) == fvc::div(phi)
+            fvm::laplacian(rUAf, p_rgh) == fvc::div(phi)
         );
 
-        pEqn.solve
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve
         (
             mesh.solver
             (
-                p.select
+                p_rgh.select
                 (
                     (
                         finalIter
@@ -36,17 +37,30 @@
         if (nonOrth == nNonOrthCorr)
         {
             // Calculate the conservative fluxes
-            phi -= pEqn.flux();
+            phi -= p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
-            p.relax();
+            p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rUAf);
+            U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rUAf);
             U.correctBoundaryConditions();
         }
     }
 
     #include "continuityErrs.H"
+
+    p = p_rgh + rhok*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rhok*gh;
+    }
 }

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H

@@ -96,29 +96,23 @@
         p_rgh + rhok*gh
     );
 
-    label p_rghRefCell = 0;
+    label pRefCell = 0;
-    scalar p_rghRefValue = 0.0;
+    scalar pRefValue = 0.0;
     setRefCell
     (
+        p,
         p_rgh,
         mesh.solutionDict().subDict("SIMPLE"),
-        p_rghRefCell,
+        pRefCell,
-        p_rghRefValue
+        pRefValue
     );
 
-    scalar pRefValue = 0.0;
-
     if (p_rgh.needReference())
     {
-        pRefValue = readScalar
-        (
-            mesh.solutionDict().subDict("SIMPLE").lookup("pRefValue")
-        );
-
         p += dimensionedScalar
         (
             "p",
             p.dimensions(),
-            pRefValue - getRefCellValue(p, p_rghRefCell)
+            pRefValue - getRefCellValue(p, pRefCell)
         );
     }

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H

@@ -18,17 +18,9 @@
             fvm::laplacian(rUAf, p_rgh) == fvc::div(phi)
         );
 
-        p_rghEqn.setReference(p_rghRefCell, p_rghRefValue);
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        // retain the residual from the first iteration
+        p_rghEqn.solve();
-        if (nonOrth == 0)
-        {
-            p_rghEqn.solve();
-        }
-        else
-        {
-            p_rghEqn.solve();
-        }
 
         if (nonOrth == nNonOrthCorr)
         {
@@ -55,7 +47,8 @@
         (
             "p",
             p.dimensions(),
-            pRefValue - getRefCellValue(p, p_rghRefCell)
+            pRefValue - getRefCellValue(p, pRefCell)
         );
+        p_rgh = p - rhok*gh;
     }
 }

applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H

@@ -17,8 +17,11 @@
          ==
             fvc::reconstruct
             (
-                fvc::interpolate(rho)*(g & mesh.Sf())
+                (
-              - fvc::snGrad(p)*mesh.magSf()
+                  - ghf*fvc::snGrad(rho)
-            )
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            ),
+            mesh.solver(U.select(finalIter))
         );
     }

applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C

@@ -80,7 +80,7 @@ int main(int argc, char *argv[])
 
             if (nOuterCorr != 1)
             {
-                p.storePrevIter();
+                p_rgh.storePrevIter();
             }
 
             #include "UEqn.H"

applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H

@@ -53,15 +53,30 @@
         )
     );
 
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+    Info<< "Reading field p_rgh\n" << endl;
+    volScalarField p_rgh
+    (
+        IOobject
+        (
+            "p_rgh",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // Force p_rgh to be consistent with p
+    p_rgh = p - rho*gh;
+
     Info<< "Creating field DpDt\n" << endl;
     volScalarField DpDt
     (
         "DpDt",
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
     );
-
-    thermo.correct();
-
-    dimensionedScalar initialMass = fvc::domainIntegrate(rho);
-
-    dimensionedScalar totalVolume = sum(mesh.V());

applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H

@@ -9,7 +9,7 @@
     );
 
     hEqn.relax();
-    hEqn.solve();
+    hEqn.solve(mesh.solver(h.select(finalIter)));
 
     thermo.correct();
 }

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -1,11 +1,9 @@
 {
-    bool closedVolume = p.needReference();
-
     rho = thermo.rho();
 
     // Thermodynamic density needs to be updated by psi*d(p) after the
     // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p;
+    thermo.rho() -= psi*p_rgh;
 
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
@@ -18,24 +16,23 @@
       + fvc::ddtPhiCorr(rUA, rho, U, phi)
     );
 
-    surfaceScalarField buoyancyPhi =
+    surfaceScalarField buoyancyPhi = -rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
-        rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
     phi += buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pEqn
+        fvScalarMatrix p_rghEqn
         (
-            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
           + fvc::div(phi)
-          - fvm::laplacian(rhorUAf, p)
+          - fvm::laplacian(rhorUAf, p_rgh)
         );
 
-        pEqn.solve
+        p_rghEqn.solve
         (
             mesh.solver
             (
-                p.select
+                p_rgh.select
                 (
                     (
                         finalIter
@@ -49,34 +46,25 @@
         if (nonOrth == nNonOrthCorr)
         {
             // Calculate the conservative fluxes
-            phi += pEqn.flux();
+            phi += p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
-            p.relax();
+            p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U += rUA*fvc::reconstruct((buoyancyPhi + pEqn.flux())/rhorUAf);
+            U += rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf);
             U.correctBoundaryConditions();
         }
     }
 
+    p = p_rgh + rho*gh;
+
     // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
+    thermo.rho() += psi*p_rgh;
 
     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
     #include "rhoEqn.H"
     #include "compressibleContinuityErrs.H"
-
-    // 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);
-        thermo.rho() = psi*p;
-        rho += (initialMass - fvc::domainIntegrate(rho))/totalVolume;
-    }
 }

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C

@@ -62,10 +62,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "hEqn.H"
-            for (int i=0; i<3; i++)
-            {
             #include "pEqn.H"
-            }
         }
 
         turbulence->correct();

applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H

@@ -23,20 +23,6 @@
     volScalarField& h = thermo.h();
     const volScalarField& psi = thermo.psi();
 
-    Info<< "Reading field p_rgh\n" << endl;
-    volScalarField p_rgh
-    (
-        IOobject
-        (
-            "p_rgh",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -53,7 +39,6 @@
 
     #include "compressibleCreatePhi.H"
 
-
     Info<< "Creating turbulence model\n" << endl;
     autoPtr<compressible::RASModel> turbulence
     (
@@ -66,40 +51,39 @@
         )
     );
 
+
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
     surfaceScalarField ghf("ghf", g & mesh.Cf());
 
-    p = p_rgh + rho*gh;
+    Info<< "Reading field p_rgh\n" << endl;
-    thermo.correct();
+    volScalarField p_rgh
-    rho = thermo.rho();
-    p_rgh = p - rho*gh;
-
-    label p_rghRefCell = 0;
-    scalar p_rghRefValue = 0.0;
-    setRefCell
     (
-        p_rgh,
+        IOobject
-        mesh.solutionDict().subDict("SIMPLE"),
+        (
-        p_rghRefCell,
+            "p_rgh",
-        p_rghRefValue
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
     );
 
+    // Force p_rgh to be consistent with p
+    p_rgh = p - rho*gh;
+
+
+    label pRefCell = 0;
     scalar pRefValue = 0.0;
-
+    setRefCell
-    if (p_rgh.needReference())
+    (
-    {
+        p,
-        pRefValue = readScalar
+        p_rgh,
-        (
+        mesh.solutionDict().subDict("SIMPLE"),
-            mesh.solutionDict().subDict("SIMPLE").lookup("pRefValue")
+        pRefCell,
-        );
+        pRefValue
-
+    );
-        p += dimensionedScalar
-        (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, p_rghRefCell)
-        );
-    }
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+    dimensionedScalar totalVolume = sum(mesh.V());

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -1,11 +1,12 @@
 {
     rho = thermo.rho();
+    rho.relax();
 
     volScalarField rUA = 1.0/UEqn().A();
     surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
 
     U = rUA*UEqn().H();
-    //UEqn.clear();
+    UEqn.clear();
 
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
     bool closedVolume = adjustPhi(phi, U, p_rgh);
@@ -20,7 +21,7 @@
             fvm::laplacian(rhorUAf, p_rgh) == fvc::div(phi)
         );
 
-        p_rghEqn.setReference(p_rghRefCell, p_rghRefValue);
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
         p_rghEqn.solve();
 
         if (nonOrth == nNonOrthCorr)
@@ -42,13 +43,13 @@
 
     p = p_rgh + rho*gh;
 
-    // For closed-volume cases adjust the pressure and density levels
+    // For closed-volume cases adjust the pressure level
     // to obey overall mass continuity
     if (closedVolume)
     {
         p += (initialMass - fvc::domainIntegrate(psi*p))
             /fvc::domainIntegrate(psi);
-        p_rgh == p - rho*gh;
+        p_rgh = p - rho*gh;
     }
 
     rho = thermo.rho();

applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C

@@ -58,7 +58,7 @@ int main(int argc, char *argv[])
 
         #include "readSIMPLEControls.H"
 
-        p.storePrevIter();
+        p_rgh.storePrevIter();
         rho.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -93,6 +93,8 @@ int main(int argc, char *argv[])
         // --- PIMPLE loop
         for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
         {
+            bool finalIter = oCorr == nOuterCorr-1;
+
             forAll(fluidRegions, i)
             {
                 Info<< "\nSolving for fluid region "

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H

@@ -13,7 +13,9 @@
      ==
         fvc::reconstruct
         (
-            fvc::interpolate(rho)*(g & mesh.Sf())
+            (
-          - fvc::snGrad(p)*mesh.magSf()
+              - ghf*fvc::snGrad(rho)
+              - fvc::snGrad(p_rgh)
+            )*mesh.magSf()
         )
     );

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H

@@ -6,12 +6,16 @@
     PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
     PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
     PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
-    PtrList<volScalarField> DpDtf(fluidRegions.size());
+    PtrList<volScalarField> p_rghFluid(fluidRegions.size());
+    PtrList<volScalarField> ghFluid(fluidRegions.size());
+    PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
 
     List<scalar> initialMassFluid(fluidRegions.size());
     List<label> pRefCellFluid(fluidRegions.size(),0);
     List<scalar> pRefValueFluid(fluidRegions.size(),0.0);
 
+    PtrList<dimensionedScalar> rhoMax(fluidRegions.size());
+    PtrList<dimensionedScalar> rhoMin(fluidRegions.size());
 
     // Populate fluid field pointer lists
     forAll(fluidRegions, i)
@@ -130,15 +134,74 @@
             ).ptr()
         );
 
+        Info<< "    Adding to ghFluid\n" << endl;
+        ghFluid.set
+        (
+            i,
+            new volScalarField("gh", gFluid[i] & fluidRegions[i].C())
+        );
+
+        Info<< "    Adding to ghfFluid\n" << endl;
+        ghfFluid.set
+        (
+            i,
+            new surfaceScalarField("ghf", gFluid[i] & fluidRegions[i].Cf())
+        );
+
+        p_rghFluid.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "p_rgh",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                fluidRegions[i]
+            )
+        );
+
+        // Force p_rgh to be consistent with p
+        p_rghFluid[i] = thermoFluid[i].p() - rhoFluid[i]*ghFluid[i];
+
         initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
 
         setRefCell
         (
             thermoFluid[i].p(),
+            p_rghFluid[i],
             fluidRegions[i].solutionDict().subDict("SIMPLE"),
             pRefCellFluid[i],
             pRefValueFluid[i]
         );
+
+        rhoMax.set
+        (
+            i,
+            new dimensionedScalar
+            (
+                fluidRegions[i].solutionDict().subDict("SIMPLE").lookup
+                (
+                    "rhoMax"
+                )
+            )
+        );
+
+        rhoMin.set
+        (
+            i,
+            new dimensionedScalar
+            (
+                fluidRegions[i].solutionDict().subDict("SIMPLE").lookup
+                (
+                    "rhoMin"
+                )
+            )
+        );
     }
 
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H

@@ -5,8 +5,8 @@
       - fvm::Sp(fvc::div(phi), h)
       - fvm::laplacian(turb.alphaEff(), h)
      ==
-        fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p))
+        fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)")
-      - p*fvc::div(phi/fvc::interpolate(rho))
+      - (rho/psi)*fvc::div(phi/fvc::interpolate(rho))
     );
 
     hEqn.relax();

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H

@@ -1,7 +1,9 @@
 {
     // From buoyantSimpleFoam
-
     rho = thermo.rho();
+    rho = max(rho, rhoMin[i]);
+    rho = min(rho, rhoMax[i]);
+    rho.relax();
 
     volScalarField rUA = 1.0/UEqn().A();
     surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
@@ -10,59 +12,54 @@
     UEqn.clear();
 
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
-    bool closedVolume = adjustPhi(phi, U, p);
+    bool closedVolume = adjustPhi(phi, U, p_rgh);
 
-    surfaceScalarField buoyancyPhi =
+    surfaceScalarField buoyancyPhi = rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
-        rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
+    phi -= buoyancyPhi;
-    phi += buoyancyPhi;
 
     // Solve pressure
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pEqn
+        fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rhorUAf, p) == fvc::div(phi)
+            fvm::laplacian(rhorUAf, p_rgh) == fvc::div(phi)
         );
 
-        pEqn.setReference(pRefCell, pRefValue);
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        // retain the residual from the first iteration
+        p_rghEqn.solve();
-        if (nonOrth == 0)
-        {
-            pEqn.solve();
-        }
-        else
-        {
-            pEqn.solve();
-        }
 
         if (nonOrth == nNonOrthCorr)
         {
-            // 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);
-            }
-
             // Calculate the conservative fluxes
-            phi -= pEqn.flux();
+            phi -= p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
-            p.relax();
+            p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rhorUAf);
+            U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf);
             U.correctBoundaryConditions();
         }
     }
 
+    p = p_rgh + rho*gh;
 
     #include "continuityErrs.H"
 
+    // For closed-volume cases adjust the pressure level
+    // to obey overall mass continuity
+    if (closedVolume)
+    {
+        p += (initialMass - fvc::domainIntegrate(psi*p))
+            /fvc::domainIntegrate(psi);
+        p_rgh = p - rho*gh;
+    }
+
     rho = thermo.rho();
+    rho = max(rho, rhoMin[i]);
+    rho = min(rho, rhoMax[i]);
     rho.relax();
 
     Info<< "Min/max rho:" << min(rho).value() << ' '

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H

@@ -5,7 +5,6 @@
     volScalarField& K = KFluid[i];
     volVectorField& U = UFluid[i];
     surfaceScalarField& phi = phiFluid[i];
-    const dimensionedVector& g = gFluid[i];
 
     compressible::turbulenceModel& turb = turbulence[i];
 
@@ -22,3 +21,7 @@
 
     const label pRefCell = pRefCellFluid[i];
     const scalar pRefValue = pRefValueFluid[i];
+
+    volScalarField& p_rgh = p_rghFluid[i];
+    const volScalarField& gh = ghFluid[i];
+    const surfaceScalarField& ghf = ghfFluid[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H

@@ -1,6 +1,6 @@
 //  Pressure-velocity SIMPLE corrector
 
-    p.storePrevIter();
+    p_rgh.storePrevIter();
     rho.storePrevIter();
     {
         #include "UEqn.H"

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H

@@ -16,8 +16,11 @@
          ==
             fvc::reconstruct
             (
-                fvc::interpolate(rho)*(g & mesh.Sf())
+                (
-              - fvc::snGrad(p)*mesh.magSf()
+                  - ghf*fvc::snGrad(rho)
-            )
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            ),
+            mesh.solver(U.select(finalIter))
         );
     }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H

@@ -6,6 +6,9 @@
     PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
     PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
     PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
+    PtrList<volScalarField> p_rghFluid(fluidRegions.size());
+    PtrList<volScalarField> ghFluid(fluidRegions.size());
+    PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
     PtrList<volScalarField> DpDtFluid(fluidRegions.size());
 
     List<scalar> initialMassFluid(fluidRegions.size());
@@ -129,6 +132,42 @@
             ).ptr()
         );
 
+        Info<< "    Adding to ghFluid\n" << endl;
+        ghFluid.set
+        (
+            i,
+            new volScalarField("gh", gFluid[i] & fluidRegions[i].C())
+        );
+
+        Info<< "    Adding to ghfFluid\n" << endl;
+        ghfFluid.set
+        (
+            i,
+            new surfaceScalarField("ghf", gFluid[i] & fluidRegions[i].Cf())
+        );
+
+        p_rghFluid.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "p_rgh",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                fluidRegions[i]
+            )
+        );
+
+        // Force p_rgh to be consistent with p
+        p_rghFluid[i] = thermoFluid[i].p() - rhoFluid[i]*ghFluid[i];
+
+        initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
+
         Info<< "    Adding to DpDtFluid\n" << endl;
         DpDtFluid.set
         (
@@ -147,6 +186,4 @@
                 )
             )
         );
-
-        initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
     }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H

@@ -7,16 +7,9 @@
      ==
         DpDt
     );
-    if (oCorr == nOuterCorr-1)
+
-    {
+    hEqn.relax();
-        hEqn.relax();
+    hEqn.solve(mesh.solver(h.select(finalIter)));
-        hEqn.solve(mesh.solver("hFinal"));
-    }
-    else
-    {
-        hEqn.relax();
-        hEqn.solve();
-    }
 
     thermo.correct();
 

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H

@@ -1,5 +1,5 @@
 {
-    bool closedVolume = p.needReference();
+    bool closedVolume = p_rgh.needReference();
 
     rho = thermo.rho();
 
@@ -17,34 +17,35 @@
         )
     );
 
-    phi = phiU + fvc::interpolate(rho)*(g & mesh.Sf())*rhorUAf;
+    phi = phiU - rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pEqn
+        fvScalarMatrix p_rghEqn
         (
-            fvm::ddt(psi, p)
+            fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh
           + fvc::div(phi)
-          - fvm::laplacian(rhorUAf, p)
+          - fvm::laplacian(rhorUAf, p_rgh)
         );
 
-        if
+        p_rghEqn.solve
         (
-            oCorr == nOuterCorr-1
+            mesh.solver
-         && corr == nCorr-1
+            (
-         && nonOrth == nNonOrthCorr
+                p_rgh.select
-        )
+                (
-        {
+                    (
-            pEqn.solve(mesh.solver(p.name() + "Final"));
+                        oCorr == nOuterCorr-1
-        }
+                     && corr == nCorr-1
-        else
+                     && nonOrth == nNonOrthCorr
-        {
+                    )
-            pEqn.solve(mesh.solver(p.name()));
+                )
-        }
+            )
+        );
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi += pEqn.flux();
+            phi += p_rghEqn.flux();
         }
     }
 
@@ -52,6 +53,8 @@
     U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
     U.correctBoundaryConditions();
 
+    p = p_rgh + rho*gh;
+
     // Update pressure substantive derivative
     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
@@ -65,9 +68,10 @@
     // to obey overall mass continuity
     if (closedVolume)
     {
-        p += (massIni - fvc::domainIntegrate(psi*p))
+        p += (initialMass - fvc::domainIntegrate(psi*p))
             /fvc::domainIntegrate(psi);
         rho = thermo.rho();
+        p_rgh = p - rho*gh;
     }
 
     // Update thermal conductivity

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPISOControls.H

@@ -1,17 +0,0 @@
-    const dictionary& piso = fluidRegions[i].solutionDict().subDict("PISO");
-
-    const int nOuterCorr =
-        piso.lookupOrDefault<int>("nOuterCorrectors", 1);
-
-    const int nCorr =
-        piso.lookupOrDefault<int>("nCorrectors", 1);
-
-    const int nNonOrthCorr =
-        piso.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-
-    const bool momentumPredictor =
-        piso.lookupOrDefault("momentumPredictor", true);
-
-    const bool transonic =
-        piso.lookupOrDefault("transonic", false);
-

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H

@@ -1,11 +1,10 @@
-    const fvMesh& mesh = fluidRegions[i];
+    fvMesh& mesh = fluidRegions[i];
 
     basicPsiThermo& thermo = thermoFluid[i];
     volScalarField& rho = rhoFluid[i];
     volScalarField& K = KFluid[i];
     volVectorField& U = UFluid[i];
     surfaceScalarField& phi = phiFluid[i];
-    const dimensionedVector& g = gFluid[i];
 
     compressible::turbulenceModel& turb = turbulence[i];
     volScalarField& DpDt = DpDtFluid[i];
@@ -14,4 +13,13 @@
     const volScalarField& psi = thermo.psi();
     volScalarField& h = thermo.h();
 
-    const dimensionedScalar massIni("massIni", dimMass, initialMassFluid[i]);
+    volScalarField& p_rgh = p_rghFluid[i];
+    const volScalarField& gh = ghFluid[i];
+    const surfaceScalarField& ghf = ghfFluid[i];
+
+    const dimensionedScalar initialMass
+    (
+        "initialMass",
+        dimMass,
+        initialMassFluid[i]
+    );

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H

@@ -1,3 +1,8 @@
+if (finalIter)
+{
+    mesh.data::add("finalIteration", true);
+}
+
 if (oCorr == 0)
 {
     #include "rhoEqn.H"
@@ -16,3 +21,8 @@ for (int corr=0; corr<nCorr; corr++)
 turb.correct();
 
 rho = thermo.rho();
+
+if (finalIter)
+{
+    mesh.data::remove("finalIteration");
+}

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/storeOldFluidFields.H

@@ -1,2 +1,2 @@
-    p.storePrevIter();
+    p_rgh.storePrevIter();
     rho.storePrevIter();

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H

@@ -8,9 +8,5 @@
             << solidRegions[i].name() << nl << endl;
 
         Info<< "    Adding to thermos\n" << endl;
-        thermos.set
+        thermos.set(i, basicSolidThermo::New(solidRegions[i]));
-        (
-            i,
-            basicSolidThermo::New(solidRegions[i])
-        );
     }

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H

@@ -1,3 +1,8 @@
+if (finalIter)
+{
+    mesh.data::add("finalIteration", true);
+}
+
 {
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
@@ -7,8 +12,13 @@
           - fvm::laplacian(K, T)
         );
         TEqn().relax();
-        TEqn().solve();
+        TEqn().solve(mesh.solver(T.select(finalIter)));
     }
 
     Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl;
 }
+
+if (finalIter)
+{
+    mesh.data::remove("finalIteration");
+}

applications/solvers/incompressible/pimpleFoam/UEqn.H

@@ -9,7 +9,7 @@ tmp<fvVectorMatrix> UEqn
 
 UEqn().relax();
 
-volScalarField rUA = 1.0/UEqn().A();
+volScalarField rAU = 1.0/UEqn().A();
 
 if (momentumPredictor)
 {
@@ -17,6 +17,6 @@ if (momentumPredictor)
 }
 else
 {
-    U = rUA*(UEqn().H() - fvc::grad(p));
+    U = rAU*(UEqn().H() - fvc::grad(p));
     U.correctBoundaryConditions();
 }

applications/solvers/incompressible/pimpleFoam/pEqn.H

@@ -1,4 +1,4 @@
-U = rUA*UEqn().H();
+U = rAU*UEqn().H();
 
 if (nCorr <= 1)
 {
@@ -6,7 +6,7 @@ if (nCorr <= 1)
 }
 
 phi = (fvc::interpolate(U) & mesh.Sf())
-    + fvc::ddtPhiCorr(rUA, U, phi);
+    + fvc::ddtPhiCorr(rAU, U, phi);
 
 adjustPhi(phi, U, p);
 
@@ -16,7 +16,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     // Pressure corrector
     fvScalarMatrix pEqn
     (
-        fvm::laplacian(rUA, p) == fvc::div(phi)
+        fvm::laplacian(rAU, p) == fvc::div(phi)
     );
 
     pEqn.setReference(pRefCell, pRefValue);
@@ -47,5 +47,5 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-U -= rUA*fvc::grad(p);
+U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createFields.H

@@ -27,7 +27,7 @@
         mesh
     );
 
-#   include "createPhi.H"
+    #include "createPhi.H"
 
 
     label pRefCell = 0;

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H

@@ -46,6 +46,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
 #include "continuityErrs.H"
 
+// Explicitly relax pressure for momentum corrector
 p.relax();
 
 // Make the fluxes relative to the mesh motion

applications/solvers/incompressible/porousSimpleFoam/createFields.H

@@ -43,16 +43,28 @@
 
 
     porousZones pZones(mesh);
+    Switch pressureImplicitPorosity(false);
 
     // nUCorrectors used for pressureImplicitPorosity
     int nUCorr = 0;
-    const bool pressureImplicitPorosity =
+    if (pZones.size())
-    (
+    {
-        pZones.size()
+        // nUCorrectors for pressureImplicitPorosity
-     && mesh.solutionDict().subDict("SIMPLE").readIfPresent
+        if (mesh.solutionDict().subDict("SIMPLE").found("nUCorrectors"))
-        (
+        {
-            "nUCorrectors",
+            nUCorr = readInt
-            nUCorr
+            (
-        )
+                mesh.solutionDict().subDict("SIMPLE").lookup("nUCorrectors")
-     && (nUCorr > 0)
+            );
-    );
+        }
+
+        if (nUCorr > 0)
+        {
+            pressureImplicitPorosity = true;
+            Info<< "Using pressure implicit porosity" << endl;
+        }
+        else
+        {
+            Info<< "Using pressure explicit porosity" << endl;
+        }
+    }

applications/solvers/incompressible/simpleFoam/pEqn.H

@@ -15,15 +15,8 @@
         );
 
         pEqn.setReference(pRefCell, pRefValue);
-        // retain the residual from the first iteration
+
-        if (nonOrth == 0)
+        pEqn.solve();
-        {
-            pEqn.solve();
-        }
-        else
-        {
-            pEqn.solve();
-        }
 
         if (nonOrth == nNonOrthCorr)
         {

applications/solvers/lagrangian/coalChemistryFoam/UEqn.H

@@ -13,5 +13,5 @@
 
     if (momentumPredictor)
     {
-        solve(UEqn == -fvc::grad(p));
+        solve(UEqn == -fvc::grad(p), mesh.solver(U.select(finalIter)));
     }

applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C

@@ -90,17 +90,28 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- PIMPLE loop
-        for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
         {
+            bool finalIter = oCorr == nOuterCorr - 1;
+            if (finalIter)
+            {
+                mesh.data::add("finalIteration", true);
+            }
+
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            for (int corr=0; corr<nCorr; corr++)
             {
                 #include "pEqn.H"
             }
+
+            if (finalIter)
+            {
+                mesh.data::remove("finalIteration");
+            }
         }
 
         turbulence->correct();

applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H

@@ -15,7 +15,7 @@
 
     hsEqn.relax();
 
-    hsEqn.solve();
+    hsEqn.solve(mesh.solver(hs.select(finalIter)));
 
     thermo.correct();
 

applications/solvers/lagrangian/coalChemistryFoam/pEqn.H

@@ -26,7 +26,20 @@ if (transonic)
             coalParcels.Srho()
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver
+            (
+                p.select
+                (
+                    (
+                        finalIter
+                     && corr == nCorr-1
+                     && nonOrth == nNonOrthCorr
+                    )
+                )
+            )
+        );
 
         if (nonOrth == nNonOrthCorr)
         {
@@ -54,7 +67,20 @@ else
             coalParcels.Srho()
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver
+            (
+                p.select
+                (
+                    (
+                        finalIter
+                     && corr == nCorr-1
+                     && nonOrth == nNonOrthCorr
+                    )
+                )
+            )
+        );
 
         if (nonOrth == nNonOrthCorr)
         {

applications/solvers/multiphase/interFoam/alphaCourantNo.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Global
-    CourantNo
+    alphaCourantNo
 
 Description
     Calculates and outputs the mean and maximum Courant Numbers.
@@ -39,17 +39,14 @@ scalar meanAlphaCoNum = 0.0;
 
 if (mesh.nInternalFaces())
 {
-    surfaceScalarField alphaf = fvc::interpolate(alpha1);
+    scalarField sumPhi =
+        pos(alpha1 - 0.01)*pos(0.99 - alpha1)
+       *fvc::surfaceSum(mag(phi))().internalField();
 
-    surfaceScalarField SfUfbyDelta =
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaT().value();
-        pos(alphaf - 0.01)*pos(0.99 - alphaf)
-       *mesh.surfaceInterpolation::deltaCoeffs()*mag(phi);
 
-    alphaCoNum = max(SfUfbyDelta/mesh.magSf())
+    meanAlphaCoNum =
-        .value()*runTime.deltaT().value();
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaT().value();
-
-    meanAlphaCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
-        .value()*runTime.deltaT().value();
 }
 
 Info<< "Interface Courant Number mean: " << meanAlphaCoNum

applications/test/CompactIOList/Make/files

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

applications/test/CompactIOList/Make/options

@@ -0,0 +1,3 @@
+EXE_INC =
+
+EXE_LIBS =

applications/test/CompactIOList/testCompactIOList.C

@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-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 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    testCompactIOList
+
+Description
+    Simple demonstration and test application for the CompactIOList container
+
+\*---------------------------------------------------------------------------*/
+
+#include "IOstreams.H"
+#include "argList.H"
+#include "Time.H"
+#include "polyMesh.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+//  Main program:
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+
+    IOstream::streamFormat format=IOstream::BINARY;
+    // IOstream::streamFormat format=IOstream::ASCII;
+
+    const label size = 20000000;
+
+    // Old format
+    // ~~~~~~~~~~
+
+    {
+        // Construct big faceList in old format
+        faceIOList faces2
+        (
+            IOobject
+            (
+                "faces2",
+                runTime.constant(),
+                polyMesh::meshSubDir,
+                runTime,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            size
+        );
+
+        const face f(identity(4));
+
+        forAll(faces2, i)
+        {
+            faces2[i] = f;
+        }
+
+        Info<< "Constructed faceList in = "
+            << runTime.cpuTimeIncrement() << " s" << nl << endl;
+
+
+        // Write binary
+        faces2.writeObject
+        (
+            format,
+            IOstream::currentVersion,
+            IOstream::UNCOMPRESSED
+        );
+
+        Info<< "Written old format faceList in = "
+            << runTime.cpuTimeIncrement() << " s" << nl << endl;
+
+        // Read
+        faceIOList faces3
+        (
+            IOobject
+            (
+                "faces2",
+                runTime.constant(),
+                polyMesh::meshSubDir,
+                runTime,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        );
+
+        Info<< "Read old format " << faces3.size() << " faceList in = "
+            << runTime.cpuTimeIncrement() << " s" << nl << endl;
+    }
+
+
+    // New format
+    // ~~~~~~~~~~
+
+    {
+        // Construct big faceList in new format
+        faceCompactIOList faces2
+        (
+            IOobject
+            (
+                "faces2",
+                runTime.constant(),
+                polyMesh::meshSubDir,
+                runTime,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            size
+        );
+
+        const face f(identity(4));
+
+        forAll(faces2, i)
+        {
+            faces2[i] = f;
+        }
+
+        Info<< "Constructed new format faceList in = "
+            << runTime.cpuTimeIncrement() << " s" << nl << endl;
+
+
+        // Write binary
+        faces2.writeObject
+        (
+            format,
+            IOstream::currentVersion,
+            IOstream::UNCOMPRESSED
+        );
+
+        Info<< "Written new format faceList in = "
+            << runTime.cpuTimeIncrement() << " s" << nl << endl;
+
+        // Read
+        faceCompactIOList faces3
+        (
+            IOobject
+            (
+                "faces2",
+                runTime.constant(),
+                polyMesh::meshSubDir,
+                runTime,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        );
+
+        Info<< "Read new format " << faces3.size() << " faceList in = "
+            << runTime.cpuTimeIncrement() << " s" << nl << endl;
+    }
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/IndirectList/IndirectListTest.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,6 +30,15 @@ Description
 
 using namespace Foam;
 
+template<class ListType>
+void printInfo(const ListType& lst)
+{
+    Info<< "addr: " << lst.addressing() << nl
+        << "list: " << lst << nl
+        << endl;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Main program:
 
@@ -42,8 +51,7 @@ int main(int argc, char *argv[])
         completeList[i] = 0.1*i;
     }
 
-    Info<< "raw : " << completeList << nl
+    Info<< "raw : " << completeList << nl << endl;
-        << endl;
 
 
     List<label> addresses(5);
@@ -53,11 +61,9 @@ int main(int argc, char *argv[])
     addresses[3] = 8;
     addresses[4] = 5;
 
-    IndirectList<double> idl(completeList, addresses);
+    IndirectList<double> idl1(completeList, addresses);
 
-    Info<< "addr: " << idl.addressing() << nl
+    printInfo(idl1);
-        << "list: " << idl() << nl
-        << endl;
 
     addresses[4] = 1;
     addresses[3] = 0;
@@ -65,11 +71,26 @@ int main(int argc, char *argv[])
     addresses[1] = 8;
     addresses[0] = 5;
 
-    idl.resetAddressing(addresses.xfer());
+    idl1.resetAddressing(addresses.xfer());
 
-    Info<< "addr: " << idl.addressing() << nl
+    printInfo(idl1);
-        << "list: " << idl() << nl
+
-        << endl;
+    // test copying
+    UIndirectList<double> uidl1(idl1);
+    IndirectList<double> idl2(uidl1);
+    IndirectList<double> idl3(idl2);
+
+    printInfo(uidl1);
+
+    idl1.resetAddressing(List<label>());
+//    idl2.resetAddressing(List<label>());
+
+    Info<<"after resetAddressing:" << nl << endl;
+
+    printInfo(uidl1);
+    printInfo(idl1);
+    printInfo(idl2);
+    printInfo(idl3);
 
     Info<< "End\n" << endl;
 

applications/test/IndirectList2/IndirectList2.H

@@ -1,164 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-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 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::IndirectList2
-
-Description
-    A List with indirect addressing.
-
-SourceFiles
-    IndirectListI.H
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef IndirectList2_H
-#define IndirectList2_H
-
-#include "UIndirectList.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                   Class IndirectListAddressing Declaration
-\*---------------------------------------------------------------------------*/
-
-//- A helper class for storing addresses.
-class IndirectListAddressing
-{
-    // Private data
-
-        //- Storage for the list addressing
-        List<label> addressing_;
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        IndirectListAddressing(const IndirectListAddressing&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const IndirectListAddressing&);
-
-
-protected:
-
-    // Constructors
-
-        //- Construct by copying the addressing array
-        explicit inline IndirectListAddressing(const UList<label>& addr);
-
-        //- Construct by transferring addressing array
-        explicit inline IndirectListAddressing(const Xfer<List<label> >& addr);
-
-
-    // Member Functions
-
-        // Access
-
-            //- Return the list addressing
-            inline const List<label>& addressing() const;
-
-        // Edit
-
-            //- Reset addressing
-            inline void resetAddressing(const UList<label>&);
-            inline void resetAddressing(const Xfer<List<label> >&);
-
-};
-
-
-/*---------------------------------------------------------------------------*\
-                        Class IndirectList2 Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class T>
-class IndirectList2
-:
-    private IndirectListAddressing,
-    public  UIndirectList<T>
-{
-    // Private Member Functions
-
-        //- Disable default assignment operator
-        void operator=(const IndirectList2<T>&);
-
-        //- Disable assignment from UIndirectList
-        void operator=(const UIndirectList<T>&);
-
-
-public:
-
-    // Constructors
-
-        //- Construct given the complete list and the addressing array
-        inline IndirectList2(const UList<T>&, const UList<label>&);
-
-        //- Construct given the complete list and by transferring addressing
-        inline IndirectList2(const UList<T>&, const Xfer<List<label> >&);
-
-        //- Copy constructor
-        inline IndirectList2(const IndirectList2<T>&);
-
-        //- Construct from UIndirectList
-        explicit inline IndirectList2(const UIndirectList<T>&);
-
-
-    // Member Functions
-
-
-        // Access
-
-            //- Return the list addressing
-            using UIndirectList<T>::addressing;
-
-
-        // Edit
-
-            //- Reset addressing
-            using IndirectListAddressing::resetAddressing;
-
-
-        // Member Operators
-
-            //- Assignment operator
-            using UIndirectList<T>::operator=;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "IndirectList2I.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/test/IndirectList2/IndirectList2I.H

@@ -1,136 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-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 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-
-inline Foam::IndirectListAddressing::IndirectListAddressing
-(
-    const UList<label>& addr
-)
-:
-    addressing_(addr)
-{}
-
-
-inline Foam::IndirectListAddressing::IndirectListAddressing
-(
-    const Xfer<List<label> >& addr
-)
-:
-    addressing_(addr)
-{}
-
-
-template<class T>
-inline Foam::IndirectList2<T>::IndirectList2
-(
-    const UList<T>& completeList,
-    const UList<label>& addr
-)
-:
-    IndirectListAddressing(addr),
-    UIndirectList<T>
-    (
-        completeList,
-        IndirectListAddressing::addressing()
-    )
-{}
-
-
-template<class T>
-inline Foam::IndirectList2<T>::IndirectList2
-(
-    const UList<T>& completeList,
-    const Xfer<List<label> >& addr
-)
-:
-    IndirectListAddressing(addr),
-    UIndirectList<T>
-    (
-        completeList,
-        IndirectListAddressing::addressing()
-    )
-{}
-
-
-template<class T>
-inline Foam::IndirectList2<T>::IndirectList2
-(
-    const IndirectList2<T>& lst
-)
-:
-    IndirectListAddressing(lst.addressing()),
-    UIndirectList<T>
-    (
-        lst.completeList(),
-        IndirectListAddressing::addressing()
-    )
-{}
-
-
-template<class T>
-inline Foam::IndirectList2<T>::IndirectList2
-(
-    const UIndirectList<T>& lst
-)
-:
-    IndirectListAddressing(lst.addressing()),
-    UIndirectList<T>
-    (
-        lst.completeList(),
-        IndirectListAddressing::addressing()
-    )
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-inline const Foam::List<Foam::label>&
-Foam::IndirectListAddressing::addressing() const
-{
-    return addressing_;
-}
-
-
-inline void Foam::IndirectListAddressing::resetAddressing
-(
-    const UList<label>& addr
-)
-{
-    addressing_ = addr;
-}
-
-
-inline void Foam::IndirectListAddressing::resetAddressing
-(
-    const Xfer<List<label> >& addr
-)
-{
-    addressing_.transfer(addr());
-}
-
-
-// ************************************************************************* //

applications/test/IndirectList2/Make/files

@@ -1,3 +0,0 @@
-IndirectListTest2.C
-
-EXE = $(FOAM_USER_APPBIN)/IndirectListTest2

applications/test/IndirectList2/Make/options

@@ -1,2 +0,0 @@
-/* EXE_INC = -I$(LIB_SRC)/cfdTools/include */
-/* EXE_LIBS = -lfiniteVolume */

applications/test/Polynomial/Make/options

@@ -1,3 +1,3 @@
-EXE_INC = \
+EXE_INC =
 
-EXE_LIBS = \
+EXE_LIBS =

applications/test/Polynomial/PolynomialTest.C

@@ -29,24 +29,30 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "IFstream.H"
+#include "IStringStream.H"
 #include "Polynomial.H"
 #include "Random.H"
+#include "cpuTime.H"
 
 using namespace Foam;
 
+const int PolySize = 8;
+const scalar coeff[] = { 0.11, 0.45, -0.94, 1.58, -2.58, 0.08, 3.15, -4.78 };
+const char* polyDef = "testPoly (0.11 0.45 -0.94 1.58 -2.58 0.08 3.15 -4.78)";
+
+
 scalar polyValue(const scalar x)
 {
     // Hard-coded polynomial 8 coeff (7th order)
     return
-        0.11
+        coeff[0]
-      + 0.45*x
+      + coeff[1]*x
-      - 0.94*sqr(x)
+      + coeff[2]*sqr(x)
-      + 1.58*pow3(x)
+      + coeff[3]*pow3(x)
-      - 2.58*pow4(x)
+      + coeff[4]*pow4(x)
-      + 0.08*pow5(x)
+      + coeff[5]*pow5(x)
-      + 3.15*pow6(x)
+      + coeff[6]*pow6(x)
-      - 4.78*x*pow6(x);
+      + coeff[7]*x*pow6(x);
 }
 
 
@@ -54,14 +60,44 @@ scalar intPolyValue(const scalar x)
 {
     // Hard-coded integrated form of above polynomial
     return
-        0.11*x
+        coeff[0]*x
-      + 0.45/2.0*sqr(x)
+      + coeff[1]/2.0*sqr(x)
-      - 0.94/3.0*pow3(x)
+      + coeff[2]/3.0*pow3(x)
-      + 1.58/4.0*pow4(x)
+      + coeff[3]/4.0*pow4(x)
-      - 2.58/5.0*pow5(x)
+      + coeff[4]/5.0*pow5(x)
-      + 0.08/6.0*pow6(x)
+      + coeff[5]/6.0*pow6(x)
-      + 3.15/7.0*x*pow6(x)
+      + coeff[6]/7.0*x*pow6(x)
-      - 4.78/8.0*x*x*pow6(x);
+      + coeff[7]/8.0*x*x*pow6(x);
+}
+
+
+scalar polyValue1(const scalar x)
+{
+    // "normal" evaluation using pow()
+    scalar value = coeff[0];
+
+    for (int i=1; i < PolySize; ++i)
+    {
+        value += coeff[i]*pow(x, i);
+    }
+
+    return value;
+}
+
+
+scalar polyValue2(const scalar x)
+{
+    // calculation avoiding pow()
+    scalar value = coeff[0];
+
+    scalar powX = x;
+    for (int i=1; i < PolySize; ++i)
+    {
+        value += coeff[i] * powX;
+        powX *= x;
+    }
+
+    return value;
 }
 
 
@@ -69,9 +105,14 @@ scalar intPolyValue(const scalar x)
 
 int main(int argc, char *argv[])
 {
-    IFstream is("polyTestInput");
+    const label n = 10000;
+    const label nIters = 1000;
+    scalar sum = 0.0;
 
-    Polynomial<8> poly("testPoly", is);
+    Info<< "null poly = " << (Polynomial<8>()) << nl << endl;
+
+    // Polynomial<8> poly("testPoly", IStringStream(polyDef)());
+    Polynomial<8> poly(coeff);
     Polynomial<9> intPoly(poly.integrate(0.0));
 
     Info<< "poly = " << poly << endl;
@@ -118,6 +159,62 @@ int main(int argc, char *argv[])
         Info<< endl;
     }
 
+
+    //
+    // test speed of Polynomial:
+    //
+    Info<< "start timing loops" << nl
+        << "~~~~~~~~~~~~~~~~~~" << endl;
+
+    cpuTime timer;
+
+    for (int loop = 0; loop < n; ++loop)
+    {
+        sum = 0.0;
+        for (label iter = 0; iter < nIters; ++iter)
+        {
+            sum += poly.evaluate(loop+iter);
+        }
+    }
+    Info<< "evaluate:     " << sum
+        << " in " << timer.cpuTimeIncrement() << " s\n";
+
+
+    for (int loop = 0; loop < n; ++loop)
+    {
+        sum = 0.0;
+        for (label iter = 0; iter < nIters; ++iter)
+        {
+            sum += polyValue(loop+iter);
+        }
+    }
+    Info<< "hard-coded 0: " << sum
+        << " in " << timer.cpuTimeIncrement() << " s\n";
+
+
+    for (int loop = 0; loop < n; ++loop)
+    {
+        sum = 0.0;
+        for (label iter = 0; iter < nIters; ++iter)
+        {
+            sum += polyValue1(loop+iter);
+        }
+    }
+    Info<< "hard-coded 1: " << sum
+        << " in " << timer.cpuTimeIncrement() << " s\n";
+
+    for (int loop = 0; loop < n; ++loop)
+    {
+        sum = 0.0;
+        for (label iter = 0; iter < nIters; ++iter)
+        {
+            sum += polyValue2(loop+iter);
+        }
+    }
+    Info<< "hard-coded 2: " << sum
+        << " in " << timer.cpuTimeIncrement() << " s\n";
+
+
     Info<< nl << "Done." << endl;
 
     return 0;

applications/test/Polynomial/polyTestInput

@@ -1,11 +0,0 @@
-    testPoly
-    (
-        0.11
-        0.45
-       -0.94
-        1.58
-       -2.58
-        0.08
-        3.15
-       -4.78
-    )

applications/utilities/mesh/conversion/ansysToFoam/ansysToFoam.L

@@ -475,7 +475,7 @@ int main(int argc, char *argv[])
         {-1, -1, -1, -1, -1, -1}, // 1
         {-1, -1, -1, -1, -1, -1}, // 2
         {-1, -1, -1, -1, -1, -1}, // 3
-        {-1,  2,  0,  3,  1, -1}, // tet (version 2.0)
+        { 3,  2,  0, -1,  1, -1}, // tet (version 2.0)
         { 0,  4,  3, -1,  2,  1}, // prism
         { 4,  2,  1,  3,  0,  5}, // hex
     };
@@ -500,12 +500,14 @@ int main(int argc, char *argv[])
             ++cellIter, ++faceIter
         )
         {
+            const cellShape& shape = cellShapes[cellMap[cellIter()]];
+
             patchFaces.append
             (
-                cellShapes[cellMap[cellIter()] ].faces()
+                shape.faces()
                 [
                     faceIndex
-                        [cellShapes[cellMap[cellIter()] ].nFaces()]
+                        [shape.nFaces()]
                         [faceIter()-1]
                 ]
             );

applications/utilities/postProcessing/dataConversion/foamToEnsight/Make/files

@@ -2,6 +2,5 @@ itoa.C
 ensightMesh.C
 ensightParticlePositions.C
 foamToEnsight.C
-ensightWriteBinary.C
 
 EXE = $(FOAM_APPBIN)/foamToEnsight

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

@@ -44,21 +44,12 @@ namespace Foam
 
 class cellSets
 {
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        cellSets(const cellSets&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const cellSets&);
-
-
 public:
 
-        label nHexesWedges;
-        label nPrisms;
-        label nPyrs;
         label nTets;
+        label nPyrs;
+        label nPrisms;
+        label nHexesWedges;
         label nPolys;
 
         labelList tets;
@@ -66,6 +57,7 @@ public:
         labelList prisms;
         labelList wedges;
         labelList hexes;
+        labelList hexesWedges;
         labelList polys;
 
 
@@ -74,10 +66,10 @@ public:
         //- Construct given the number ov cells
         cellSets(const label nCells)
         :
-            nHexesWedges(0),
-            nPrisms(0),
-            nPyrs(0),
             nTets(0),
+            nPyrs(0),
+            nPrisms(0),
+            nHexesWedges(0),
             nPolys(0),
 
             tets(nCells),
@@ -85,6 +77,7 @@ public:
             prisms(nCells),
             wedges(nCells),
             hexes(nCells),
+            hexesWedges(nCells),
             polys(nCells)
         {}
 };

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

@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-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 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::ensightAsciiStream
+
+Description
+
+SourceFiles
+    ensightAsciiStream.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ensightAsciiStream_H
+#define ensightAsciiStream_H
+
+#include "ensightStream.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class ensightAsciiStream Declaration
+\*---------------------------------------------------------------------------*/
+
+class ensightAsciiStream
+:
+    public ensightStream
+{
+    // Private data
+
+        //- Description of data_
+        OFstream str_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        ensightAsciiStream(const ensightAsciiStream&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const ensightAsciiStream&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        ensightAsciiStream(const fileName& f, const Time& runTime)
+        :
+            ensightStream(f),
+            str_
+            (
+                f,
+                runTime.writeFormat(),
+                runTime.writeVersion(),
+                IOstream::UNCOMPRESSED
+            )
+        {
+
+            str_.setf(ios_base::scientific, ios_base::floatfield);
+            str_.precision(5);
+        }
+
+
+    //- Destructor
+    virtual ~ensightAsciiStream()
+    {}
+
+
+    // Member Functions
+
+        virtual bool ascii() const
+        {
+            return true;
+        }
+
+        virtual void write(const char* c)
+        {
+            str_ << c << nl;
+        }
+
+        virtual void write(const int v)
+        {
+            str_ << setw(10) << v << nl;
+        }
+
+        virtual void write(const scalarField& sf)
+        {
+            forAll(sf, i)
+            {
+                if (mag(sf[i]) >= scalar(floatScalarVSMALL))
+                {
+                    str_ << setw(12) << sf[i] << nl;
+                }
+                else
+                {
+                    str_ << setw(12) << scalar(0) << nl;
+                }
+            }
+        }
+
+        virtual void write(const List<int>& sf)
+        {
+            forAll(sf, i)
+            {
+                str_ << setw(10) << sf[i];
+            }
+            str_<< nl;
+        }
+
+        virtual void writePartHeader(const label partI)
+        {
+            str_<< "part" << nl
+                << setw(10) << partI << nl;
+        }
+
+    // Member Operators
+
+    // Friend Functions
+
+    // Friend Operators
+
+    // IOstream Operators
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-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 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::ensightBinaryStream
+
+Description
+
+SourceFiles
+    ensightBinaryStream.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ensightBinaryStream_H
+#define ensightBinaryStream_H
+
+#include "ensightStream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class ensightBinaryStream Declaration
+\*---------------------------------------------------------------------------*/
+
+class ensightBinaryStream
+:
+    public ensightStream
+{
+    // Private data
+
+        //- Description of data_
+        autoPtr<std::ofstream> str_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        ensightBinaryStream(const ensightBinaryStream&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const ensightBinaryStream&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        ensightBinaryStream(const fileName& f, const Time&)
+        :
+            ensightStream(f),
+            str_
+            (
+                new std::ofstream
+                (
+                    f.c_str(),
+                    ios_base::out | ios_base::binary | ios_base::trunc
+                )
+            )
+        {}
+
+
+    //- Destructor
+    virtual ~ensightBinaryStream()
+    {}
+
+
+    // Member Functions
+
+        virtual bool ascii() const
+        {
+            return false;
+        }
+
+        virtual void write(const char* val)
+        {
+            char buffer[80] = {0};
+            strcpy(buffer, val);
+            str_().write(buffer, 80*sizeof(char));
+        }
+
+        virtual void write(const int val)
+        {
+            str_().write(reinterpret_cast<const char*>(&val), sizeof(int));
+        }
+
+        virtual void write(const scalarField& sf)
+        {
+            if (sf.size())
+            {
+                List<float> temp(sf.size());
+
+                forAll(sf, i)
+                {
+                    temp[i] = float(sf[i]);
+                }
+
+                str_().write
+                (
+                    reinterpret_cast<const char*>(temp.begin()),
+                    sf.size()*sizeof(float)
+                );
+            }
+        }
+
+        virtual void write(const List<int>& sf)
+        {
+            str_().write
+            (
+                reinterpret_cast<const char*>(sf.begin()),
+                sf.size()*sizeof(int)
+            );
+        }
+
+        virtual void writePartHeader(const label partI)
+        {
+            write("part");
+            write(partI);
+        }
+
+    // Member Operators
+
+    // Friend Functions
+
+    // Friend Operators
+
+    // IOstream Operators
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -29,99 +29,38 @@ License
 #include "OFstream.H"
 #include "IOmanip.H"
 #include "itoa.H"
-#include "ensightWriteBinary.H"
+#include "volPointInterpolation.H"
+#include "ensightBinaryStream.H"
+#include "ensightAsciiStream.H"
+#include "globalIndex.H"
 
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
-void writeData(const scalarField& sf, OFstream& ensightFile)
-{
-    forAll(sf, i)
-    {
-        if (mag( sf[i] ) >= scalar(floatScalarVSMALL))
-        {
-            ensightFile << setw(12) << sf[i] << nl;
-        }
-        else
-        {
-            ensightFile << setw(12) << scalar(0) << nl;
-        }
-    }
-}
-
-
 template<class Type>
-scalarField map
+void writeField
-(
-    const Field<Type>& vf,
-    const labelList& map,
-    const label cmpt
-)
-{
-    scalarField mf(map.size());
-
-    forAll(map, i)
-    {
-        mf[i] = component(vf[map[i]], cmpt);
-    }
-
-    return mf;
-}
-
-
-template<class Type>
-scalarField map
-(
-    const Field<Type>& vf,
-    const labelList& map1,
-    const labelList& map2,
-    const label cmpt
-)
-{
-    scalarField mf(map1.size() + map2.size());
-
-    forAll(map1, i)
-    {
-        mf[i] = component(vf[map1[i]], cmpt);
-    }
-
-    label offset = map1.size();
-
-    forAll(map2, i)
-    {
-        mf[i + offset] = component(vf[map2[i]], cmpt);
-    }
-
-    return mf;
-}
-
-
-template<class Type>
-void writeAllData
 (
     const char* key,
     const Field<Type>& vf,
-    const labelList& prims,
+    ensightStream& ensightFile
-    const label nPrims,
-    OFstream& ensightFile
 )
 {
-    if (nPrims)
+    if (returnReduce(vf.size(), sumOp<label>()) > 0)
     {
         if (Pstream::master())
         {
-            ensightFile << key << nl;
+            ensightFile.write(key);
 
             for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
             {
-                writeData(map(vf, prims, cmpt), ensightFile);
+                ensightFile.write(vf.component(cmpt));
 
                 for (int slave=1; slave<Pstream::nProcs(); slave++)
                 {
                     IPstream fromSlave(Pstream::scheduled, slave);
-                    scalarField data(fromSlave);
+                    scalarField slaveData(fromSlave);
-                    writeData(data, ensightFile);
+                    ensightFile.write(slaveData);
                 }
             }
         }
@@ -130,129 +69,7 @@ void writeAllData
             for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
             {
                 OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
-                toMaster<< map(vf, prims, cmpt);
+                toMaster<< vf.component(cmpt);
-            }
-        }
-    }
-}
-
-
-template<class Type>
-void writeAllDataBinary
-(
-    const char* key,
-    const Field<Type>& vf,
-    const labelList& prims,
-    const label nPrims,
-    std::ofstream& ensightFile
-)
-{
-    if (nPrims)
-    {
-        if (Pstream::master())
-        {
-            writeEnsDataBinary(key,ensightFile);
-
-            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
-            {
-                writeEnsDataBinary(map(vf, prims, cmpt), ensightFile);
-
-                for (int slave=1; slave<Pstream::nProcs(); 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);
-            }
-        }
-    }
-}
-
-
-template<class Type>
-void writeAllFaceData
-(
-    const char* key,
-    const labelList& prims,
-    const label nPrims,
-    const Field<Type>& pf,
-    OFstream& ensightFile
-)
-{
-    if (nPrims)
-    {
-        if (Pstream::master())
-        {
-            ensightFile << key << nl;
-
-            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
-            {
-                writeData(map(pf, prims, cmpt), ensightFile);
-
-                for (int slave=1; slave<Pstream::nProcs(); 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);
-            }
-        }
-    }
-}
-
-
-template<class Type>
-void writeAllFaceDataBinary
-(
-    const char* key,
-    const labelList& prims,
-    const label nPrims,
-    const Field<Type>& pf,
-    std::ofstream& ensightFile
-)
-{
-    if (nPrims)
-    {
-        if (Pstream::master())
-        {
-            writeEnsDataBinary(key,ensightFile);
-
-            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
-            {
-                writeEnsDataBinary(map(pf, prims, cmpt), ensightFile);
-
-                for (int slave=1; slave<Pstream::nProcs(); 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);
             }
         }
     }
@@ -267,97 +84,34 @@ bool writePatchField
     const Foam::label ensightPatchI,
     const Foam::faceSets& boundaryFaceSet,
     const Foam::ensightMesh::nFacePrimitives& nfp,
-    Foam::OFstream& ensightFile
+    ensightStream& ensightFile
 )
 {
     if (nfp.nTris || nfp.nQuads || nfp.nPolys)
     {
         if (Pstream::master())
         {
-            ensightFile
+            ensightFile.writePartHeader(ensightPatchI);
-                << "part" << nl
-                << setw(10) << ensightPatchI << nl;
         }
 
-        writeAllFaceData
+        writeField
         (
             "tria3",
-            boundaryFaceSet.tris,
+            Field<Type>(pf, boundaryFaceSet.tris),
-            nfp.nTris,
-            pf,
             ensightFile
         );
 
-        writeAllFaceData
+        writeField
         (
             "quad4",
-            boundaryFaceSet.quads,
+            Field<Type>(pf, boundaryFaceSet.quads),
-            nfp.nQuads,
-            pf,
             ensightFile
         );
 
-        writeAllFaceData
+        writeField
         (
             "nsided",
-            boundaryFaceSet.polys,
+            Field<Type>(pf, boundaryFaceSet.polys),
-            nfp.nPolys,
-            pf,
-            ensightFile
-        );
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-template<class Type>
-bool writePatchFieldBinary
-(
-    const Foam::Field<Type>& pf,
-    const Foam::label patchi,
-    const Foam::label ensightPatchI,
-    const Foam::faceSets& boundaryFaceSet,
-    const Foam::ensightMesh::nFacePrimitives& nfp,
-    std::ofstream& ensightFile
-)
-{
-    if (nfp.nTris || nfp.nQuads || nfp.nPolys)
-    {
-        if (Pstream::master())
-        {
-            writeEnsDataBinary("part",ensightFile);
-            writeEnsDataBinary(ensightPatchI,ensightFile);
-        }
-
-        writeAllFaceDataBinary
-        (
-            "tria3",
-            boundaryFaceSet.tris,
-            nfp.nTris,
-            pf,
-            ensightFile
-        );
-
-        writeAllFaceDataBinary
-        (
-            "quad4",
-            boundaryFaceSet.quads,
-            nfp.nQuads,
-            pf,
-            ensightFile
-        );
-
-        writeAllFaceDataBinary
-        (
-            "nsided",
-            boundaryFaceSet.polys,
-            nfp.nPolys,
-            pf,
             ensightFile
         );
 
@@ -380,6 +134,7 @@ void writePatchField
     const Foam::fileName& postProcPath,
     const Foam::word& prepend,
     const Foam::label timeIndex,
+    const bool binary,
     Foam::Ostream& ensightCaseFile
 )
 {
@@ -409,7 +164,7 @@ void writePatchField
 
     word timeFile = prepend + itoa(timeIndex);
 
-    OFstream *ensightFilePtr = NULL;
+    ensightStream* ensightFilePtr = NULL;
     if (Pstream::master())
     {
         if (timeIndex == 0)
@@ -426,20 +181,30 @@ void writePatchField
 
         // set the filename of the ensight file
         fileName ensightFileName(timeFile + "." + pfName);
-        ensightFilePtr = new OFstream
+
-        (
+        if (binary)
-            postProcPath/ensightFileName,
+        {
-            runTime.writeFormat(),
+            ensightFilePtr = new ensightBinaryStream
-            runTime.writeVersion(),
+            (
-            runTime.writeCompression()
+                postProcPath/ensightFileName,
-        );
+                runTime
+            );
+        }
+        else
+        {
+            ensightFilePtr = new ensightAsciiStream
+            (
+                postProcPath/ensightFileName,
+                runTime
+            );
+        }
     }
 
-    OFstream& ensightFile = *ensightFilePtr;
+    ensightStream& ensightFile = *ensightFilePtr;
 
     if (Pstream::master())
     {
-        ensightFile << pTraits<Type>::typeName << nl;
+        ensightFile.write(pTraits<Type>::typeName);
     }
 
     if (patchi >= 0)
@@ -477,17 +242,18 @@ void writePatchField
 
 
 template<class Type>
-void ensightFieldAscii
+void ensightField
 (
-    const Foam::IOobject& fieldObject,
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
     const Foam::ensightMesh& eMesh,
     const Foam::fileName& postProcPath,
     const Foam::word& prepend,
     const Foam::label timeIndex,
+    const bool binary,
     Foam::Ostream& ensightCaseFile
 )
 {
-    Info<< "Converting field " << fieldObject.name() << endl;
+    Info<< "Converting field " << vf.name() << endl;
 
     word timeFile = prepend + itoa(timeIndex);
 
@@ -508,27 +274,34 @@ void ensightFieldAscii
     const labelList& tets = meshCellSets.tets;
     const labelList& pyrs = meshCellSets.pyrs;
     const labelList& prisms = meshCellSets.prisms;
-    const labelList& wedges = meshCellSets.wedges;
+    const labelList& hexesWedges = meshCellSets.hexesWedges;
-    const labelList& hexes = meshCellSets.hexes;
     const labelList& polys = meshCellSets.polys;
 
-    OFstream *ensightFilePtr = NULL;
+    ensightStream* ensightFilePtr = NULL;
     if (Pstream::master())
     {
         // set the filename of the ensight file
-        fileName ensightFileName(timeFile + "." + fieldObject.name());
+        fileName ensightFileName(timeFile + "." + vf.name());
-        ensightFilePtr = new OFstream
+
-        (
+        if (binary)
-            postProcPath/ensightFileName,
+        {
-            runTime.writeFormat(),
+            ensightFilePtr = new ensightBinaryStream
-            runTime.writeVersion(),
+            (
-            IOstream::UNCOMPRESSED
+                postProcPath/ensightFileName,
-        );
+                runTime
+            );
+        }
+        else
+        {
+            ensightFilePtr = new ensightAsciiStream
+            (
+                postProcPath/ensightFileName,
+                runTime
+            );
+        }
     }
 
-    OFstream& ensightFile = *ensightFilePtr;
+    ensightStream& ensightFile = *ensightFilePtr;
-
-    GeometricField<Type, fvPatchField, volMesh> vf(fieldObject, mesh);
 
     if (patchNames.empty())
     {
@@ -548,80 +321,42 @@ void ensightFieldAscii
                     << nl;
             }
 
+            ensightFile.write(pTraits<Type>::typeName);
+            ensightFile.writePartHeader(1);
+        }
+
+        writeField
+        (
+            "hexa8",
+            Field<Type>(vf, hexesWedges),
             ensightFile
-                << pTraits<Type>::typeName << nl
+        );
-                << "part" << nl
-                << setw(10) << 1 << nl;
 
-            ensightFile.setf(ios_base::scientific, ios_base::floatfield);
+        writeField
-            ensightFile.precision(5);
-        }
-
-        if (meshCellSets.nHexesWedges)
-        {
-            if (Pstream::master())
-            {
-                ensightFile << "hexa8" << nl;
-
-                for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
-                {
-                    writeData
-                    (
-                        map(vf, hexes, wedges, cmpt),
-                        ensightFile
-                    );
-
-                    for (int slave=1; slave<Pstream::nProcs(); 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,
+            Field<Type>(vf, prisms),
-            prisms,
-            meshCellSets.nPrisms,
             ensightFile
         );
 
-        writeAllData
+        writeField
         (
             "pyramid5",
-            vf,
+            Field<Type>(vf, pyrs),
-            pyrs,
-            meshCellSets.nPyrs,
             ensightFile
         );
 
-        writeAllData
+        writeField
         (
             "tetra4",
-            vf,
+            Field<Type>(vf, tets),
-            tets,
-            meshCellSets.nTets,
             ensightFile
         );
 
-        writeAllData
+        writeField
         (
             "nfaced",
-            vf,
+            Field<Type>(vf, polys),
-            polys,
-            meshCellSets.nPolys,
             ensightFile
         );
     }
@@ -727,7 +462,6 @@ void ensightFieldAscii
             }
         }
     }
-
     if (Pstream::master())
     {
         delete ensightFilePtr;
@@ -736,59 +470,54 @@ void ensightFieldAscii
 
 
 template<class Type>
-void ensightFieldBinary
+void ensightPointField
 (
-    const Foam::IOobject& fieldObject,
+    const GeometricField<Type, pointPatchField, pointMesh>& pf,
     const Foam::ensightMesh& eMesh,
     const Foam::fileName& postProcPath,
     const Foam::word& prepend,
     const Foam::label timeIndex,
+    const bool binary,
     Foam::Ostream& ensightCaseFile
 )
 {
-    Info<< "Converting field (binary) " << fieldObject.name() << endl;
+    Info<< "Converting field " << pf.name() << endl;
 
     word timeFile = prepend + itoa(timeIndex);
 
     const fvMesh& mesh = eMesh.mesh();
-    //const Time& runTime = mesh.time();
-
-    const cellSets& meshCellSets = eMesh.meshCellSets();
-    const List<faceSets>& boundaryFaceSets = eMesh.boundaryFaceSets();
     const wordList& allPatchNames = eMesh.allPatchNames();
     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;
-    const labelList& prisms = meshCellSets.prisms;
-    const labelList& wedges = meshCellSets.wedges;
-    const labelList& hexes = meshCellSets.hexes;
-    const labelList& polys = meshCellSets.polys;
 
-    std::ofstream *ensightFilePtr = NULL;
+    ensightStream* ensightFilePtr = NULL;
     if (Pstream::master())
     {
         // set the filename of the ensight file
-        fileName ensightFileName(timeFile + "." + fieldObject.name());
+        fileName ensightFileName(timeFile + "." + pf.name());
-        ensightFilePtr = new std::ofstream
+
-        (
+        if (binary)
-            (postProcPath/ensightFileName).c_str(),
+        {
-            ios_base::out | ios_base::binary | ios_base::trunc
+            ensightFilePtr = new ensightBinaryStream
-        );
+            (
-        // Check on file opened?
+                postProcPath/ensightFileName,
+                mesh.time()
+            );
+        }
+        else
+        {
+            ensightFilePtr = new ensightAsciiStream
+            (
+                postProcPath/ensightFileName,
+                mesh.time()
+            );
+        }
     }
 
-    std::ofstream& ensightFile = *ensightFilePtr;
+    ensightStream& ensightFile = *ensightFilePtr;
 
-    GeometricField<Type, fvPatchField, volMesh> vf(fieldObject, mesh);
+    if (eMesh.patchNames().empty())
-
-    if (patchNames.empty())
     {
         eMesh.barrier();
 
@@ -800,192 +529,122 @@ void ensightFieldBinary
 
                 ensightCaseFile
                     << pTraits<Type>::typeName
-                    << " per element:            1       "
+                    << " per node:            1       "
-                    << setw(15) << vf.name()
+                    << setw(15) << pf.name()
-                    << (' ' + prepend + "***." + vf.name()).c_str()
+                    << (' ' + prepend + "***." + pf.name()).c_str()
                     << nl;
             }
 
-            writeEnsDataBinary(pTraits<Type>::typeName,ensightFile);
+            ensightFile.write(pTraits<Type>::typeName);
-            writeEnsDataBinary("part",ensightFile);
+            ensightFile.writePartHeader(1);
-            writeEnsDataBinary(1,ensightFile);
         }
 
-        if (meshCellSets.nHexesWedges)
+        writeField
-        {
-            if (Pstream::master())
-            {
-                writeEnsDataBinary("hexa8",ensightFile);
-
-                for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
-                {
-                    writeEnsDataBinary
-                    (
-                        map(vf, hexes, wedges, cmpt),
-                        ensightFile
-                    );
-
-                    for (int slave=1; slave<Pstream::nProcs(); 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
         (
-            "penta6",
+            "coordinates",
-            vf,
+            Field<Type>(pf.internalField(), eMesh.uniquePointMap()),
-            prisms,
-            meshCellSets.nPrisms,
-            ensightFile
-        );
-
-        writeAllDataBinary
-        (
-            "pyramid5",
-            vf,
-            pyrs,
-            meshCellSets.nPyrs,
-            ensightFile
-        );
-
-        writeAllDataBinary
-        (
-            "tetra4",
-            vf,
-            tets,
-            meshCellSets.nTets,
-            ensightFile
-        );
-
-        writeAllDataBinary
-        (
-            "nfaced",
-            vf,
-            polys,
-            meshCellSets.nPolys,
             ensightFile
         );
     }
 
-    label ensightPatchI = eMesh.patchPartOffset();
 
-    forAll(allPatchNames, patchi)
+     label ensightPatchI = eMesh.patchPartOffset();
-    {
+ 
-        const word& patchName = allPatchNames[patchi];
+     forAll(allPatchNames, patchi)
-
+     {
-        eMesh.barrier();
+         const word& patchName = allPatchNames[patchi];
-
+ 
-        if (patchNames.empty() || patchNames.found(patchName))
+         eMesh.barrier();
-        {
+ 
-            if
+         if (patchNames.empty() || patchNames.found(patchName))
-            (
+         {
-                writePatchFieldBinary
+             const fvPatch& p = mesh.boundary()[patchi];
-                (
+             if
-                    vf.boundaryField()[patchi],
+             (
-                    patchi,
+                 returnReduce(p.size(), sumOp<label>())
-                    ensightPatchI,
+               > 0
-                    boundaryFaceSets[patchi],
+             )
-                    nPatchPrims.find(patchName)(),
+             {
-                    ensightFile
+                 // Renumber the patch points/faces into unique points
-                )
+                 labelList pointToGlobal;
-            )
+                 labelList uniqueMeshPointLabels;
-            {
+                 autoPtr<globalIndex> globalPointsPtr =
-                ensightPatchI++;
+                 mesh.globalData().mergePoints
-            }
+                 (
-        }
+                     p.patch().meshPoints(),
-
+                     p.patch().meshPointMap(),
-    }
+                     pointToGlobal,
-
+                     uniqueMeshPointLabels
-    // write faceZones, if requested
+                 );
-    if (faceZoneNames.size())
+ 
-    {
+                 if (Pstream::master())
-        // Interpolates cell values to faces - needed only when exporting
+                 {
-        // faceZones...
+                     ensightFile.writePartHeader(ensightPatchI);
-        GeometricField<Type, fvsPatchField, surfaceMesh> sf
+                 }
-        (
+ 
-            linearInterpolate(vf)
+                 writeField
-        );
+                 (
-
+                     "coordinates",
-        forAllConstIter(wordHashSet, faceZoneNames, iter)
+                     Field<Type>(pf.internalField(), uniqueMeshPointLabels),
-        {
+                     ensightFile
-            const word& faceZoneName = iter.key();
+                 );
-
+ 
-            eMesh.barrier();
+                 ensightPatchI++;
-
+             }
-            label zoneID = mesh.faceZones().findZoneID(faceZoneName);
+         }
-
+     }
-            const faceZone& fz = mesh.faceZones()[zoneID];
+ 
-
+     // write faceZones, if requested
-            // Prepare data to write
+     if (faceZoneNames.size())
-            label nIncluded = 0;
+     {
-            forAll(fz, i)
+         forAllConstIter(wordHashSet, faceZoneNames, iter)
-            {
+         {
-                if (eMesh.faceToBeIncluded(fz[i]))
+             const word& faceZoneName = iter.key();
-                {
+ 
-                    ++nIncluded;
+             eMesh.barrier();
-                }
+ 
-            }
+             label zoneID = mesh.faceZones().findZoneID(faceZoneName);
-
+ 
-            Field<Type> values(nIncluded);
+             const faceZone& fz = mesh.faceZones()[zoneID];
-
+ 
-            // Loop on the faceZone and store the needed field values
+             if (returnReduce(fz().nPoints(), sumOp<label>()) > 0)
-            label j = 0;
+             {
-            forAll(fz, i)
+                 // Renumber the faceZone points/faces into unique points
-            {
+                 labelList pointToGlobal;
-                label faceI = fz[i];
+                 labelList uniqueMeshPointLabels;
-                if (mesh.isInternalFace(faceI))
+                 autoPtr<globalIndex> globalPointsPtr =
-                {
+                 mesh.globalData().mergePoints
-                    values[j] = sf[faceI];
+                 (
-                    ++j;
+                     fz().meshPoints(),
-                }
+                     fz().meshPointMap(),
-                else
+                     pointToGlobal,
-                {
+                     uniqueMeshPointLabels
-                    if (eMesh.faceToBeIncluded(faceI))
+                 );
-                    {
+ 
-                        label patchI = mesh.boundaryMesh().whichPatch(faceI);
+                 if (Pstream::master())
-                        const polyPatch& pp = mesh.boundaryMesh()[patchI];
+                 {
-                        label patchFaceI = pp.whichFace(faceI);
+                     ensightFile.writePartHeader(ensightPatchI);
-                        Type value = sf.boundaryField()[patchI][patchFaceI];
+                 }
-                        values[j] = value;
+ 
-                        ++j;
+                 writeField
-                    }
+                 (
-                }
+                     "coordinates",
-            }
+                     Field<Type>
-
+                     (
-            if
+                         pf.internalField(),
-            (
+                         uniqueMeshPointLabels
-                writePatchFieldBinary
+                     ),
-                (
+                     ensightFile
-                    values,
+                 );
-                    zoneID,
+ 
-                    ensightPatchI,
+                 ensightPatchI++;
-                    faceZoneFaceSets[zoneID],
+             }
-                    nFaceZonePrims.find(faceZoneName)(),
+         }
-                    ensightFile
+     }
-                )
-            )
-            {
-                ensightPatchI++;
-            }
-        }
-    }
 
     if (Pstream::master())
     {
-        ensightFile.close();
+        delete ensightFilePtr;
     }
 }
 
@@ -999,30 +658,42 @@ void ensightField
     const Foam::word& prepend,
     const Foam::label timeIndex,
     const bool binary,
+    const bool nodeValues,
     Foam::Ostream& ensightCaseFile
 )
 {
-    if (binary)
+    // Read field
+    GeometricField<Type, fvPatchField, volMesh> vf(fieldObject, eMesh.mesh());
+
+    if (nodeValues)
     {
-        ensightFieldBinary<Type>
+        tmp<GeometricField<Type, pointPatchField, pointMesh> > pfld
         (
-            fieldObject,
+            volPointInterpolation::New(eMesh.mesh()).interpolate(vf)
+        );
+        pfld().rename(vf.name());
+
+        ensightPointField<Type>
+        (
+            pfld,
             eMesh,
             postProcPath,
             prepend,
             timeIndex,
+            binary,
             ensightCaseFile
         );
     }
     else
     {
-        ensightFieldAscii<Type>
+        ensightField<Type>
         (
-            fieldObject,
+            vf,
             eMesh,
             postProcPath,
             prepend,
             timeIndex,
+            binary,
             ensightCaseFile
         );
     }

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

@@ -35,7 +35,6 @@ SourceFiles
 #define ensightField_H
 
 #include "ensightMesh.H"
-#include "HashSet.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,6 +47,7 @@ void ensightField
     const Foam::word& prepend,
     const Foam::label timeIndex,
     const bool binary,
+    const bool nodeValues,
     Foam::Ostream& ensightCaseFile
 );
 

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

@@ -38,10 +38,13 @@ SourceFiles
 #include "faceSets.H"
 #include "HashTable.H"
 #include "HashSet.H"
-#include "fvMesh.H"
-#include "OFstream.H"
-#include <fstream>
 #include "PackedBoolList.H"
+#include "wordReList.H"
+#include "scalarField.H"
+#include "cellShapeList.H"
+#include "cellList.H"
+
+#include <fstream>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,6 +54,7 @@ namespace Foam
 class fvMesh;
 class argList;
 class globalIndex;
+class ensightStream;
 
 /*---------------------------------------------------------------------------*\
                            Class ensightMesh Declaration
@@ -82,8 +86,19 @@ private:
         //- Reference to the OpenFOAM mesh
         const fvMesh& mesh_;
 
+        //- Suppress patches
+        const bool noPatches_;
+
+        //- Output selected patches only
+        const bool patches_;
+        const wordReList patchPatterns_;
+
+        //- Output selected faceZones
+        const bool faceZones_;
+        const wordReList faceZonePatterns_;
+
         //- Set binary file output
-        bool binary_;
+        const bool binary_;
 
         //- The ensight part id for the first patch
         label patchPartOffset_;
@@ -109,6 +124,19 @@ private:
         PackedBoolList boundaryFaceToBeIncluded_;
 
 
+        // Parallel merged points
+
+            //- Global numbering for merged points
+            autoPtr<globalIndex> globalPointsPtr_;
+
+            //- From mesh point to global merged point
+            labelList pointToGlobal_;
+
+            //- Local points that are unique
+            labelList uniquePointMap_;
+
+
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -120,7 +148,7 @@ private:
         void writePoints
         (
             const scalarField& pointsComponent,
-            OFstream& ensightGeometryFile
+            ensightStream& ensightGeometryFile
         ) const;
 
         cellShapeList map
@@ -141,14 +169,14 @@ private:
         void writePrims
         (
             const cellShapeList& cellShapes,
-            OFstream& ensightGeometryFile
+            ensightStream& ensightGeometryFile
         ) const;
 
         void writePolysNFaces
         (
             const labelList& polys,
             const cellList& cellFaces,
-            OFstream& ensightGeometryFile
+            ensightStream& ensightGeometryFile
         ) const;
 
         void writePolysNPointsPerFace
@@ -156,7 +184,7 @@ private:
             const labelList& polys,
             const cellList& cellFaces,
             const faceList& faces,
-            OFstream& ensightGeometryFile
+            ensightStream& ensightGeometryFile
         ) const;
 
         void writePolysPoints
@@ -164,13 +192,13 @@ private:
             const labelList& polys,
             const cellList& cellFaces,
             const faceList& faces,
-            OFstream& ensightGeometryFile
+            ensightStream& ensightGeometryFile
         ) const;
 
         void writeAllPolys
         (
             const labelList& pointToGlobal,
-            OFstream& ensightGeometryFile
+            ensightStream& ensightGeometryFile
         ) const;
 
         void writeAllPrims
@@ -178,13 +206,13 @@ private:
             const char* key,
             const label nPrims,
             const cellShapeList& cellShapes,
-            OFstream& ensightGeometryFile
+            ensightStream& ensightGeometryFile
         ) const;
 
         void writeFacePrims
         (
             const faceList& patchFaces,
-            OFstream& ensightGeometryFile
+            ensightStream& ensightGeometryFile
         ) const;
 
         void writeAllFacePrims
@@ -193,19 +221,19 @@ private:
             const labelList& prims,
             const label nPrims,
             const faceList& patchFaces,
-            OFstream& ensightGeometryFile
+            ensightStream& ensightGeometryFile
         ) const;
 
         void writeNSidedNPointsPerFace
         (
             const faceList& patchFaces,
-            OFstream& ensightGeometryFile
+            ensightStream& ensightGeometryFile
         ) const;
 
         void writeNSidedPoints
         (
             const faceList& patchFaces,
-            OFstream& ensightGeometryFile
+            ensightStream& ensightGeometryFile
         ) const;
 
         void writeAllNSided
@@ -213,14 +241,14 @@ private:
             const labelList& prims,
             const label nPrims,
             const faceList& patchFaces,
-            OFstream& ensightGeometryFile
+            ensightStream& ensightGeometryFile
         ) const;
 
         void writeAllInternalPoints
         (
             const pointField& uniquePoints,
             const label nPoints,
-            OFstream& ensightGeometryFile
+            ensightStream& ensightGeometryFile
         ) const;
 
         void writeAllPatchPoints
@@ -229,123 +257,7 @@ private:
             const word& patchName,
             const pointField& uniquePoints,
             const label nPoints,
-            OFstream& ensightGeometryFile
+            ensightStream& 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,
-            const word& prepend,
-            const label timeIndex,
-            Ostream& ensightCaseFile,
-            const labelList& pointToGlobal,
-            const pointField& uniquePoints,
-            const globalIndex& globalPoints
-        ) const;
-
-        void writeBinary
-        (
-            const fileName& postProcPath,
-            const word& prepend,
-            const label timeIndex,
-            Ostream& ensightCaseFile,
-            const labelList& pointToGlobal,
-            const pointField& uniquePoints,
-            const globalIndex& globalPoints
-        ) const;
-
-        void writePrimsBinary
-        (
-            const cellShapeList& cellShapes,
-            std::ofstream& ensightGeometryFile
-        ) const;
-
-        void writeAllPrimsBinary
-        (
-            const char* key,
-            const label nPrims,
-            const cellShapeList& cellShapes,
-            std::ofstream& ensightGeometryFile
-        ) const;
-
-        void writePolysNFacesBinary
-        (
-            const labelList& polys,
-            const cellList& cellFaces,
-            std::ofstream& ensightGeometryFile
-        ) const;
-
-        void writePolysNPointsPerFaceBinary
-        (
-            const labelList& polys,
-            const cellList& cellFaces,
-            const faceList& faces,
-            std::ofstream& ensightGeometryFile
-        ) const;
-
-        void writePolysPointsBinary
-        (
-            const labelList& polys,
-            const cellList& cellFaces,
-            const faceList& faces,
-            std::ofstream& ensightGeometryFile
-        ) const;
-
-        void writeAllPolysBinary
-        (
-            const labelList& pointToGlobal,
-            std::ofstream& ensightGeometryFile
-        ) const;
-
-        void writeAllFacePrimsBinary
-        (
-            const char* key,
-            const labelList& prims,
-            const label nPrims,
-            const faceList& patchFaces,
-            std::ofstream& ensightGeometryFile
-        ) const;
-
-        void writeFacePrimsBinary
-        (
-            const faceList& patchFaces,
-            std::ofstream& ensightGeometryFile
-        ) const;
-
-        void writeNSidedPointsBinary
-        (
-            const faceList& patchFaces,
-            std::ofstream& ensightGeometryFile
-        ) const;
-
-        void writeNSidedNPointsPerFaceBinary
-        (
-            const faceList& patchFaces,
-            std::ofstream& ensightGeometryFile
-        ) const;
-
-        void writeAllNSidedBinary
-        (
-            const labelList& prims,
-            const label nPrims,
-            const faceList& patchFaces,
-            std::ofstream& ensightGeometryFile
         ) const;
 
 public:
@@ -355,8 +267,12 @@ public:
         //- Construct from fvMesh
         ensightMesh
         (
-            const fvMesh&,
+            const fvMesh& mesh,
-            const argList& args,
+            const bool noPatches,
+            const bool patches,
+            const wordReList& patchPatterns,
+            const bool faceZones,
+            const wordReList& faceZonePatterns,
             const bool binary
         );
 
@@ -420,8 +336,35 @@ public:
                 return patchPartOffset_;
             }
 
+
+        // Parallel point merging
+
+            //- Global numbering for merged points
+            const globalIndex& globalPoints() const
+            {
+                return globalPointsPtr_();
+            }
+
+            //- From mesh point to global merged point
+            const labelList& pointToGlobal() const
+            {
+                return pointToGlobal_;
+            }
+
+            //- Local points that are unique
+            const labelList& uniquePointMap() const
+            {
+                return uniquePointMap_;
+            }
+
+
+
+
     // Other
 
+        //- Update for new mesh
+        void correct();
+
         //- 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;

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

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-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 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::ensightStream
+
+Description
+    Abstract base class for writing Ensight data
+
+SourceFiles
+    ensightStream.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ensightStream_H
+#define ensightStream_H
+
+#include "fileName.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+
+/*---------------------------------------------------------------------------*\
+                         Class ensightStream Declaration
+\*---------------------------------------------------------------------------*/
+
+class ensightStream
+{
+    // Private data
+
+        const fileName name_;
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        ensightStream(const ensightStream&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const ensightStream&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        ensightStream(const fileName& f)
+        :
+            name_(f)
+        {}
+
+
+    //- Destructor
+    virtual ~ensightStream()
+    {}
+
+
+    // Member Functions
+
+        const fileName& name() const
+        {
+            return name_;
+        }
+
+        virtual bool ascii() const = 0;
+
+        virtual void write(const char*) = 0;
+
+        virtual void write(const int) = 0;
+
+        virtual void write(const scalarField&) = 0;
+
+        virtual void write(const List<int>&) = 0;
+
+        virtual void writePartHeader(const label) = 0;
+
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -44,15 +44,6 @@ namespace Foam
 
 class faceSets
 {
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        faceSets(const faceSets&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const faceSets&);
-
-
 public:
 
         label nTris;

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

@@ -104,6 +104,11 @@ int main(int argc, char *argv[])
         "write in ASCII format instead of 'C Binary'"
     );
     argList::addBoolOption
+    (
+        "nodeValues",
+        "write values in nodes"
+    );
+    argList::addBoolOption
     (
         "noPatches",
         "suppress writing any patches"
@@ -126,6 +131,7 @@ int main(int argc, char *argv[])
 
     // Check options
     const bool binary = !args.optionFound("ascii");
+    const bool nodeValues = args.optionFound("nodeValues");
 
 #   include "createTime.H"
 
@@ -191,7 +197,29 @@ int main(int argc, char *argv[])
     OFstream& ensightCaseFile = *ensightCaseFilePtr;
 
     // Construct the EnSight mesh
-    ensightMesh eMesh(mesh, args, binary);
+    const bool selectedPatches = args.optionFound("patches");
+    wordReList patchPatterns;
+    if (selectedPatches)
+    {
+        patchPatterns = wordReList(args.optionLookup("patches")());
+    }
+    const bool selectedZones = args.optionFound("faceZones");
+    wordReList zonePatterns;
+    if (selectedZones)
+    {
+        zonePatterns = wordReList(args.optionLookup("faceZones")());
+    }
+
+    ensightMesh eMesh
+    (
+        mesh,
+        args.optionFound("noPatches"),
+        selectedPatches,
+        patchPatterns,
+        selectedZones,
+        zonePatterns,
+        binary
+    );
 
     // Set Time to the last time before looking for the lagrangian objects
     runTime.setTime(Times.last(), Times.size()-1);
@@ -313,6 +341,11 @@ int main(int argc, char *argv[])
 
         polyMesh::readUpdateState meshState = mesh.readUpdate();
 
+        if (meshState != polyMesh::UNCHANGED)
+        {
+            eMesh.correct();
+        }
+
         if (timeIndex == 0 || (meshState != polyMesh::UNCHANGED))
         {
             eMesh.write
@@ -371,6 +404,7 @@ int main(int argc, char *argv[])
                         prepend,
                         timeIndex,
                         binary,
+                        nodeValues,
                         ensightCaseFile
                     );
                 }
@@ -384,6 +418,7 @@ int main(int argc, char *argv[])
                         prepend,
                         timeIndex,
                         binary,
+                        nodeValues,
                         ensightCaseFile
                     );
                 }
@@ -397,6 +432,7 @@ int main(int argc, char *argv[])
                         prepend,
                         timeIndex,
                         binary,
+                        nodeValues,
                         ensightCaseFile
                     );
                 }
@@ -410,6 +446,7 @@ int main(int argc, char *argv[])
                         prepend,
                         timeIndex,
                         binary,
+                        nodeValues,
                         ensightCaseFile
                     );
                 }
@@ -423,6 +460,7 @@ int main(int argc, char *argv[])
                         prepend,
                         timeIndex,
                         binary,
+                        nodeValues,
                         ensightCaseFile
                     );
                 }

applications/utilities/postProcessing/sampling/sample/sample.C

@@ -116,7 +116,7 @@ int main(int argc, char *argv[])
         sampledSets::typeName,
         mesh,
         sampleDict,
-        IOobject::MUST_READ,
+        IOobject::MUST_READ_IF_MODIFIED,
         true
     );
 
@@ -125,7 +125,7 @@ int main(int argc, char *argv[])
         sampledSurfaces::typeName,
         mesh,
         sampleDict,
-        IOobject::MUST_READ,
+        IOobject::MUST_READ_IF_MODIFIED,
         true
     );
 

applications/utilities/postProcessing/velocityField/Co/Co.C

@@ -25,7 +25,7 @@ Application
     Co
 
 Description
-    Calculates and writes the Co number as a surfaceScalarField obtained
+    Calculates and writes the Co number as a volScalarField obtained
     from field phi.
 
     The -noWrite option just outputs the max values without writing the
@@ -38,52 +38,6 @@ Description
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-    tmp<volScalarField> Co(const surfaceScalarField& Cof)
-    {
-        const fvMesh& mesh = Cof.mesh();
-
-        tmp<volScalarField> tCo
-        (
-            new volScalarField
-            (
-                IOobject
-                (
-                    "Co",
-                    mesh.time().timeName(),
-                    mesh
-                ),
-                mesh,
-                dimensionedScalar("0", Cof.dimensions(), 0)
-            )
-        );
-
-        volScalarField& Co = tCo();
-
-        // Set local references to mesh data
-        const unallocLabelList& owner = mesh.owner();
-        const unallocLabelList& neighbour = mesh.neighbour();
-
-        forAll(owner, facei)
-        {
-            label own = owner[facei];
-            label nei = neighbour[facei];
-
-            Co[own] = max(Co[own], Cof[facei]);
-            Co[nei] = max(Co[nei], Cof[facei]);
-        }
-
-        forAll(Co.boundaryField(), patchi)
-        {
-            Co.boundaryField()[patchi] = Cof.boundaryField()[patchi];
-        }
-
-        return tCo;
-    }
-}
-
-
 void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
 {
     bool writeResults = !args.optionFound("noWrite");
@@ -98,15 +52,28 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
 
     if (phiHeader.headerOk())
     {
-        autoPtr<surfaceScalarField> CoPtr;
+        volScalarField Co
+        (
+            IOobject
+            (
+                "Co",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ
+            ),
+            mesh,
+            dimensionedScalar("0", dimless, 0),
+            zeroGradientFvPatchScalarField::typeName
+        );
 
         Info<< "    Reading phi" << endl;
         surfaceScalarField phi(phiHeader, mesh);
-        Info<< "    Calculating Co" << endl;
 
         if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
         {
-            // compressible
+            Info<< "    Calculating compressible Co" << endl;
+
+            Info<< "    Reading rho" << endl;
             volScalarField rho
             (
                 IOobject
@@ -119,60 +86,34 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
                 mesh
             );
 
-            CoPtr.set
+            Co.dimensionedInternalField() =
-            (
+                (0.5*runTime.deltaT())
-                new surfaceScalarField
+               *fvc::surfaceSum(mag(phi))().dimensionedInternalField()
-                (
+               /(rho*mesh.V());
-                    IOobject
+            Co.correctBoundaryConditions();
-                    (
-                        "Cof",
-                        runTime.timeName(),
-                        mesh,
-                        IOobject::NO_READ
-                    ),
-                    (
-                        mesh.surfaceInterpolation::deltaCoeffs()
-                      * (mag(phi)/(fvc::interpolate(rho)*mesh.magSf()))
-                      * runTime.deltaT()
-                    )
-                )
-            );
         }
         else if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
         {
-            // incompressible
+            Info<< "    Calculating incompressible Co" << endl;
-            CoPtr.set
+
-            (
+            Co.dimensionedInternalField() =
-                new surfaceScalarField
+                (0.5*runTime.deltaT())
-                (
+               *fvc::surfaceSum(mag(phi))().dimensionedInternalField()
-                    IOobject
+               /mesh.V();
-                    (
+            Co.correctBoundaryConditions();
-                        "Cof",
-                        runTime.timeName(),
-                        mesh,
-                        IOobject::NO_READ
-                    ),
-                    (
-                        mesh.surfaceInterpolation::deltaCoeffs()
-                      * (mag(phi)/mesh.magSf())
-                      * runTime.deltaT()
-                    )
-                )
-            );
         }
         else
         {
             FatalErrorIn(args.executable())
                 << "Incorrect dimensions of phi: " << phi.dimensions()
-                    << abort(FatalError);
+                << abort(FatalError);
         }
 
-        Info<< "Co max : " << max(CoPtr()).value() << endl;
+        Info<< "Co max : " << max(Co).value() << endl;
 
         if (writeResults)
         {
-            CoPtr().write();
+            Co.write();
-            Co(CoPtr())().write();
         }
     }
     else
@@ -183,4 +124,5 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
     Info<< "\nEnd\n" << endl;
 }
 
+
 // ************************************************************************* //

bin/tools/CleanFunctions

@@ -76,6 +76,11 @@ cleanCase()
            sets/streamLines \
            > /dev/null 2>&1
 
+    if [ -e constant/polyMesh/blockMeshDict.m4 ]
+    then
+        rm -f constant/polyMesh/blockMeshDict > /dev/null 2>&1
+    fi
+
     for f in `find . -name "*Dict"`
     do
         sed -e /arguments/d $f > temp.$$

doc/changes/splitCyclic.txt

@@ -20,7 +20,8 @@ The disadvantages:
 - a patch-wise loop now might need to store data to go to the neighbour half
 since it is no longer handled in a single patch.
 - decomposed cyclics now require overlapping communications so will
-only work in non-blocking mode. Hence the underlying message passing library
+only work in 'nonBlocking' mode or 'blocking' (=buffered) mode but not
+in 'scheduled' mode. The underlying message passing library
 will require overlapping communications with message tags.
 - it is quite a code-change and there might be some oversights.
 - once converted (see foamUpgradeCyclics below) cases are not backwards
@@ -103,19 +104,14 @@ type 'processorCyclic'.
 
 
 - processor patches use overlapping communication using a different message
-tag. This maps straight through into the MPI message tag.
+tag. This maps straight through into the MPI message tag. Each processor
-See processorCyclicPolyPatch::tag(). This needs to be calculated the
+'interface' (processorPolyPatch, processorFvPatch, etc.) has a 'tag()'
-same on both sides so is calculated as
+to use for communication.
-        Pstream::nProcs()*max(myProcNo, neighbProcNo)
-      + min(myProcNo, neighbProcNo)
-which is
-- unique
-- commutative
-- does not interfere with the default tag (= 1)
 
 
 - when constructing a GeometricField from a dictionary it will explicitly
 check for non-existing entries for cyclic patches and exit with an error message
-warning to run foamUpgradeCyclics.
+warning to run foamUpgradeCyclics. (1.7.x will check if you are trying
+to run a case which has split cyclics)
 
 

src/OpenFOAM/Make/files

@@ -85,6 +85,7 @@ $(primitiveLists)/vectorListIOList.C
 $(primitiveLists)/sphericalTensorList.C
 $(primitiveLists)/symmTensorList.C
 $(primitiveLists)/tensorList.C
+$(primitiveLists)/hashedWordList.C
 
 Streams = db/IOstreams
 $(Streams)/token/tokenIO.C
@@ -296,6 +297,7 @@ primitiveShapes = meshes/primitiveShapes
 $(primitiveShapes)/line/line.C
 $(primitiveShapes)/plane/plane.C
 $(primitiveShapes)/triangle/intersection.C
+$(primitiveShapes)/objectHit/pointIndexHitIOList.C
 
 meshShapes = meshes/meshShapes
 $(meshShapes)/edge/edge.C

src/OpenFOAM/containers/Lists/IndirectList/IndirectList.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,11 @@ Class
     Foam::IndirectList
 
 Description
-    A List with indirect addressing
+    A List with indirect addressing.
+
+SeeAlso
+    Foam::UIndirectList for a version without any allocation for the
+    addressing.
 
 SourceFiles
     IndirectListI.H
@@ -35,24 +39,79 @@ SourceFiles
 #ifndef IndirectList_H
 #define IndirectList_H
 
-#include "List.H"
+#include "UIndirectList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+/*---------------------------------------------------------------------------*\
+                   Class IndirectListAddressing Declaration
+\*---------------------------------------------------------------------------*/
+
+//- A helper class for storing addresses.
+class IndirectListAddressing
+{
+    // Private data
+
+        //- Storage for the list addressing
+        List<label> addressing_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        IndirectListAddressing(const IndirectListAddressing&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const IndirectListAddressing&);
+
+
+protected:
+
+    // Constructors
+
+        //- Construct by copying the addressing array
+        explicit inline IndirectListAddressing(const UList<label>& addr);
+
+        //- Construct by transferring addressing array
+        explicit inline IndirectListAddressing(const Xfer<List<label> >& addr);
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the list addressing
+            inline const List<label>& addressing() const;
+
+        // Edit
+
+            //- Reset addressing
+            inline void resetAddressing(const UList<label>&);
+            inline void resetAddressing(const Xfer<List<label> >&);
+
+};
+
+
 /*---------------------------------------------------------------------------*\
                         Class IndirectList Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class T>
 class IndirectList
+:
+    private IndirectListAddressing,
+    public  UIndirectList<T>
 {
-    // Private data
+    // Private Member Functions
 
-        UList<T>& completeList_;
+        //- Disable default assignment operator
-        List<label> addressing_;
+        void operator=(const IndirectList<T>&);
+
+        //- Disable assignment from UIndirectList
+        void operator=(const UIndirectList<T>&);
 
 
 public:
@@ -65,59 +124,32 @@ public:
         //- Construct given the complete list and by transferring addressing
         inline IndirectList(const UList<T>&, const Xfer<List<label> >&);
 
+        //- Copy constructor
+        inline IndirectList(const IndirectList<T>&);
+
+        //- Construct from UIndirectList
+        explicit inline IndirectList(const UIndirectList<T>&);
+
 
     // Member Functions
 
+
         // Access
 
-            //- Return the number of elements in the list
-            inline label size() const;
-
-            //- Return true if the list is empty (ie, size() is zero).
-            inline bool empty() const;
-
-            //- Return the first element of the list.
-            inline T& first();
-
-            //- Return first element of the list.
-            inline const T& first() const;
-
-            //- Return the last element of the list.
-            inline T& last();
-
-            //- Return the last element of the list.
-            inline const T& last() const;
-
-            //- Return the complete list
-            inline const UList<T>& completeList() const;
-
             //- Return the list addressing
-            inline const List<label>& addressing() const;
+            using UIndirectList<T>::addressing;
 
 
         // Edit
 
             //- Reset addressing
-            inline void resetAddressing(const UList<label>&);
+            using IndirectListAddressing::resetAddressing;
-            inline void resetAddressing(const Xfer<List<label> >&);
 
 
         // Member Operators
 
-            //- Return the addressed elements as a List
+            //- Assignment operator
-            inline List<T> operator()() const;
+            using UIndirectList<T>::operator=;
-
-            //- Return non-const access to an element
-            inline T& operator[](const label);
-
-            //- Return const access to an element
-            inline const T& operator[](const label) const;
-
-            //- Assignment from UList of addressed elements
-            inline void operator=(const UList<T>&);
-
-            //- Assignment of all entries to the given value
-            inline void operator=(const T&);
 };
 
 

src/OpenFOAM/containers/Lists/IndirectList/IndirectListI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,25 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+
+inline Foam::IndirectListAddressing::IndirectListAddressing
+(
+    const UList<label>& addr
+)
+:
+    addressing_(addr)
+{}
+
+
+inline Foam::IndirectListAddressing::IndirectListAddressing
+(
+    const Xfer<List<label> >& addr
+)
+:
+    addressing_(addr)
+{}
+
+
 template<class T>
 inline Foam::IndirectList<T>::IndirectList
 (
@@ -32,8 +51,12 @@ inline Foam::IndirectList<T>::IndirectList
     const UList<label>& addr
 )
 :
-    completeList_(const_cast<UList<T>&>(completeList)),
+    IndirectListAddressing(addr),
-    addressing_(addr)
+    UIndirectList<T>
+    (
+        completeList,
+        IndirectListAddressing::addressing()
+    )
 {}
 
 
@@ -44,71 +67,55 @@ inline Foam::IndirectList<T>::IndirectList
     const Xfer<List<label> >& addr
 )
 :
-    completeList_(const_cast<UList<T>&>(completeList)),
+    IndirectListAddressing(addr),
-    addressing_(addr)
+    UIndirectList<T>
+    (
+        completeList,
+        IndirectListAddressing::addressing()
+    )
+{}
+
+
+template<class T>
+inline Foam::IndirectList<T>::IndirectList
+(
+    const IndirectList<T>& lst
+)
+:
+    IndirectListAddressing(lst.addressing()),
+    UIndirectList<T>
+    (
+        lst.completeList(),
+        IndirectListAddressing::addressing()
+    )
+{}
+
+
+template<class T>
+inline Foam::IndirectList<T>::IndirectList
+(
+    const UIndirectList<T>& lst
+)
+:
+    IndirectListAddressing(lst.addressing()),
+    UIndirectList<T>
+    (
+        lst.completeList(),
+        IndirectListAddressing::addressing()
+    )
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class T>
+inline const Foam::List<Foam::label>&
-inline Foam::label Foam::IndirectList<T>::size() const
+Foam::IndirectListAddressing::addressing() const
-{
-    return addressing_.size();
-}
-
-
-template<class T>
-inline bool Foam::IndirectList<T>::empty() const
-{
-    return addressing_.empty();
-}
-
-
-template<class T>
-inline T& Foam::IndirectList<T>::first()
-{
-    return completeList_[addressing_.first()];
-}
-
-
-template<class T>
-inline const T& Foam::IndirectList<T>::first() const
-{
-    return completeList_[addressing_.first()];
-}
-
-
-template<class T>
-inline T& Foam::IndirectList<T>::last()
-{
-    return completeList_[addressing_.last()];
-}
-
-
-template<class T>
-inline const T& Foam::IndirectList<T>::last() const
-{
-    return completeList_[addressing_.last()];
-}
-
-
-template<class T>
-inline const Foam::UList<T>& Foam::IndirectList<T>::completeList() const
-{
-    return completeList_;
-}
-
-
-template<class T>
-inline const Foam::List<Foam::label>& Foam::IndirectList<T>::addressing() const
 {
     return addressing_;
 }
 
 
-template<class T>
+inline void Foam::IndirectListAddressing::resetAddressing
-inline void Foam::IndirectList<T>::resetAddressing
 (
     const UList<label>& addr
 )
@@ -117,8 +124,7 @@ inline void Foam::IndirectList<T>::resetAddressing
 }
 
 
-template<class T>
+inline void Foam::IndirectListAddressing::resetAddressing
-inline void Foam::IndirectList<T>::resetAddressing
 (
     const Xfer<List<label> >& addr
 )
@@ -127,63 +133,4 @@ inline void Foam::IndirectList<T>::resetAddressing
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-template<class T>
-inline Foam::List<T> Foam::IndirectList<T>::operator()() const
-{
-    List<T> result(size());
-
-    forAll(*this, i)
-    {
-        result[i] = operator[](i);
-    }
-
-    return result;
-}
-
-
-template<class T>
-inline T& Foam::IndirectList<T>::operator[](const label i)
-{
-    return completeList_[addressing_[i]];
-}
-
-
-template<class T>
-inline const T& Foam::IndirectList<T>::operator[](const label i) const
-{
-    return completeList_[addressing_[i]];
-}
-
-
-template<class T>
-inline void Foam::IndirectList<T>::operator=(const UList<T>& ae)
-{
-    if (addressing_.size() != ae.size())
-    {
-        FatalErrorIn("IndirectList<T>::operator=(const UList<T>&)")
-            << "Addressing and list of addressed elements "
-               "have different sizes: "
-            << addressing_.size() << " " << ae.size()
-            << abort(FatalError);
-    }
-
-    forAll(addressing_, i)
-    {
-        completeList_[addressing_[i]] = ae[i];
-    }
-}
-
-
-template<class T>
-inline void Foam::IndirectList<T>::operator=(const T& t)
-{
-    forAll(addressing_, i)
-    {
-        completeList_[addressing_[i]] = t;
-    }
-}
-
-
 // ************************************************************************* //

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

@@ -266,24 +266,6 @@ Foam::List<T>::List(const SLList<T>& lst)
 }
 
 
-// Construct as copy of IndirectList<T>
-template<class T>
-Foam::List<T>::List(const IndirectList<T>& lst)
-:
-    UList<T>(NULL, lst.size())
-{
-    if (this->size_)
-    {
-        this->v_ = new T[this->size_];
-
-        forAll(*this, i)
-        {
-            this->operator[](i) = lst[i];
-        }
-    }
-}
-
-
 // Construct as copy of UIndirectList<T>
 template<class T>
 Foam::List<T>::List(const UIndirectList<T>& lst)
@@ -517,25 +499,6 @@ void Foam::List<T>::operator=(const SLList<T>& lst)
 }
 
 
-// Assignment operator. Takes linear time.
-template<class T>
-void Foam::List<T>::operator=(const IndirectList<T>& lst)
-{
-    if (lst.size() != this->size_)
-    {
-        if (this->v_) delete[] this->v_;
-        this->v_ = 0;
-        this->size_ = lst.size();
-        if (this->size_) this->v_ = new T[this->size_];
-    }
-
-    forAll(*this, i)
-    {
-        this->operator[](i) = lst[i];
-    }
-}
-
-
 // Assignment operator. Takes linear time.
 template<class T>
 void Foam::List<T>::operator=(const UIndirectList<T>& lst)

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

@@ -131,9 +131,6 @@ public:
         //- Construct as copy of SLList<T>
         explicit List(const SLList<T>&);
 
-        //- Construct as copy of IndirectList<T>
-        explicit List(const IndirectList<T>&);
-
         //- Construct as copy of UIndirectList<T>
         explicit List(const UIndirectList<T>&);
 
@@ -219,9 +216,6 @@ public:
         //- Assignment from SLList operator. Takes linear time.
         void operator=(const SLList<T>&);
 
-        //- Assignment from IndirectList operator. Takes linear time.
-        void operator=(const IndirectList<T>&);
-
         //- Assignment from UIndirectList operator. Takes linear time.
         void operator=(const UIndirectList<T>&);
 

src/OpenFOAM/db/dlLibraryTable/dlLibraryTable.C

@@ -82,9 +82,8 @@ bool Foam::dlLibraryTable::open(const fileName& functionLibName)
         }
         else
         {
-            if (!loadedLibraries.found(functionLibPtr))
+            if (loadedLibraries.insert(functionLibPtr, functionLibName))
             {
-                loadedLibraries.insert(functionLibPtr, functionLibName);
                 return true;
             }
             else

src/OpenFOAM/meshes/meshShapes/face/face.H

@@ -126,6 +126,14 @@ class face
 
 public:
 
+    //- Return types for classify
+    enum proxType
+    {
+        NONE,
+        POINT,  // Close to point
+        EDGE    // Close to edge
+    };
+
     // Static data members
 
         static const char* const typeName;
@@ -249,6 +257,20 @@ public:
             const pointField& meshPoints
         ) const;
 
+        //- Return nearest point to face and classify it:
+        //  + near point (nearType=POINT, nearLabel=0, 1, 2)
+        //  + near edge (nearType=EDGE, nearLabel=0, 1, 2)
+        //    Note: edges are counted from starting vertex so
+        //    e.g. edge n is from f[n] to f[0], where the face has n + 1
+        //    points
+        pointHit nearestPointClassify
+        (
+            const point& p,
+            const pointField& meshPoints,
+            label& nearType,
+            label& nearLabel
+        ) const;
+
         //- Return contact sphere diameter
         scalar contactSphereDiameter
         (

src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C

@@ -181,6 +181,22 @@ Foam::pointHit Foam::face::nearestPoint
     const point& p,
     const pointField& meshPoints
 ) const
+{
+    // Dummy labels
+    label nearType = -1;
+    label nearLabel = -1;
+
+    return nearestPointClassify(p, meshPoints, nearType, nearLabel);
+}
+
+
+Foam::pointHit Foam::face::nearestPointClassify
+(
+    const point& p,
+    const pointField& meshPoints,
+    label& nearType,
+    label& nearLabel
+) const
 {
     const face& f = *this;
     point ctr = centre(meshPoints);
@@ -188,6 +204,9 @@ Foam::pointHit Foam::face::nearestPoint
     // Initialize to miss, distance=GREAT
     pointHit nearest(p);
 
+    nearType = -1;
+    nearLabel = -1;
+
     label nPoints = f.size();
 
     point nextPoint = ctr;
@@ -196,8 +215,10 @@ Foam::pointHit Foam::face::nearestPoint
     {
         nextPoint = meshPoints[f[fcIndex(pI)]];
 
+        label tmpNearType = -1;
+        label tmpNearLabel = -1;
+
         // Note: for best accuracy, centre point always comes last
-        //
         triPointRef tri
         (
             meshPoints[f[pI]],
@@ -205,12 +226,42 @@ Foam::pointHit Foam::face::nearestPoint
             ctr
         );
 
-        pointHit curHit = tri.nearestPoint(p);
+        pointHit curHit = tri.nearestPointClassify
+        (
+            p,
+            tmpNearType,
+            tmpNearLabel
+        );
 
         if (Foam::mag(curHit.distance()) < Foam::mag(nearest.distance()))
         {
             nearest.setDistance(curHit.distance());
 
+            // Assume at first that the near type is NONE on the
+            // triangle (i.e. on the face of the triangle) then it is
+            // therefore also for the face.
+
+            nearType = NONE;
+
+            if (tmpNearType == triPointRef::EDGE && tmpNearLabel == 0)
+            {
+                // If the triangle edge label is 0, then this is also
+                // an edge of the face, if not, it is on the face
+
+                nearType = EDGE;
+
+                nearLabel = pI;
+            }
+            else if (tmpNearType == triPointRef::POINT && tmpNearLabel < 2)
+            {
+                // If the triangle point label is 0 or 1, then this is
+                // also a point of the face, if not, it is on the face
+
+                nearType = POINT;
+
+                nearLabel = pI + tmpNearLabel;
+            }
+
             if (curHit.hit())
             {
                 nearest.setHit();

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C

@@ -27,7 +27,10 @@ License
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-const Foam::scalar Foam::polyMeshTetDecomposition::minTetQuality = 1e-9;
+// Note: the use of this tolerance is ad-hoc, there may be extreme
+// cases where the resulting tetrahedra still have particle tracking
+// problems.
+const Foam::scalar Foam::polyMeshTetDecomposition::minTetQuality = SMALL;
 
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H

@@ -686,8 +686,14 @@ public:
 
         // Useful derived info
 
-            //- Is the point in the cell bounding box
+            //- Is the point in the cell bounding box, option relative
-            bool pointInCellBB(const point& p, label celli) const;
+            //  tolerance to increase the effective size of the boundBox
+            bool pointInCellBB
+            (
+                const point& p,
+                label celli,
+                scalar tol = 0
+            ) const;
 
             //- Is the point in the cell
             bool pointInCell(const point& p, label celli) const;

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C

@@ -30,16 +30,33 @@ License
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 // Is the point in the cell bounding box
-bool Foam::primitiveMesh::pointInCellBB(const point& p, label celli) const
+bool Foam::primitiveMesh::pointInCellBB
+(
+    const point& p,
+    label celli,
+    scalar tol
+) const
 {
-    return boundBox
+    boundBox bb
     (
         cells()[celli].points
         (
             faces(),
             points()
-        )
+        ),
-    ).contains(p);
+        false
+    );
+
+    if (tol > SMALL)
+    {
+        bb = boundBox
+        (
+            bb.min() - tol*bb.span(),
+            bb.max() + tol*bb.span()
+        );
+    }
+
+    return bb.contains(p);
 }