Merge branch 'master' into cvm

ReleaseNotes-dev

@@ -105,7 +105,15 @@
     + Additional wall functions for primary region momentum and temperature
       taking film into account
     + Parallel aware
-*** *New* ptscotch decomposition method
+*** *New* ptscotch decomposition method.
+*** *New* multiLevel decomposition method.
+    Decomposes in levels, e.g. first decompose onto number of nodes and
+    then onto number of cores per node. This will minimise off-node
+    communication. Each level can use any of the other decomposition methods
+*** *New* structured decomposition method.
+    Does a 2D decomposition of a mesh. Valid only for an 'extruded' mesh, i.e.
+    columns of cells originating from a patch. Bases decomposition on this
+    patch and assigns the cells according to the patch decomposition.
 *** *Updated* scotch decomposition method to run in parallel by doing
     decomposition on the master. Unfortunately scotch and ptscotch cannot
     be linked in to the same executable.
@@ -212,7 +220,8 @@
     }
     #+END_SRC
     See also [[./doc/changes/dynamicCode.org]]
-
+  + cyclicSlip: cyclic with point motion constrained to tangential plane.
+    Can be used with any mesh movement, e.g. snapping in snappyHexMesh.
 * Utilities
   There have been some utilities added and updated in this release.
 *** *New* utilities
@@ -229,7 +238,8 @@
         (nonuniformTransform)cyclic <zoneA>_<zoneB>
       + extrudes into master direction (i.e. away from the owner cell
         if flipMap is false)
-    + =topoSet=: replacement of cellSet,faceSet,pointSet utilities.
+    + =topoSet=: replacement of cellSet,faceSet,pointSet utilities. Multiple
+      commands operating on different sets.
       Comparable to a dictionary driven =setSet= utility.
 *** Updated utilities
     + =setFields=: optionally use faceSets to set patch values (see
@@ -244,6 +254,8 @@
       - works in parallel
     + =snappyHexMesh=:
       + extrude across multi-processor boundaries
+      + specify type of patches created during meshing
+      + handle cyclics in initial mesh (non parallel meshing only)
       + preserve faceZones shape during layering
       + combining coincident patch faces is now default after snapping
       + *Warning*:

applications/solvers/combustion/fireFoam/Make/files

@@ -1,4 +1,3 @@
 fireFoam.C
 
 EXE = $(FOAM_APPBIN)/fireFoam
-

applications/solvers/combustion/fireFoam/Make/options

@@ -1,20 +1,52 @@
 EXE_INC = \
-    -I./combustionModels/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solid/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidChemistryModel/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/regionModels/pyrolysisModels/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
 
 
 EXE_LIBS = \
-    -lcombustionModels \
-    -lspecie \
-    -lreactionThermophysicalModels \
-    -lbasicThermophysicalModels \
     -lfiniteVolume \
-    -lcompressibleTurbulenceModel \
+    -lmeshTools \
-    -lcompressibleLESModels \
     -lcompressibleRASModels \
-    -lradiationModels
+    -lcompressibleLESModels \
+    -lspecie \
+    -lbasicThermophysicalModels \
+    -lsolidProperties \
+    -lsolidMixtureProperties \
+    -lthermophysicalFunctions \
+    -lreactionThermophysicalModels \
+    -lSLGThermo \
+    -lchemistryModel \
+    -lsolidChemistryModel \
+    -lcombustionModels \
+    -lregionModels \
+    -lradiationModels \
+    -lsurfaceFilmModels \
+    -lpyrolysisModels \
+    -llagrangianIntermediate \
+    -lODE \
+    -lsampling

applications/solvers/combustion/fireFoam/UEqn.H

@@ -1,21 +1,26 @@
-fvVectorMatrix UEqn
+    fvVectorMatrix UEqn
-(
-    fvm::ddt(rho, U)
-  + fvm::div(phi, U)
-  + turbulence->divDevRhoReff(U)
-);
-
-UEqn.relax();
-
-solve
-(
-    UEqn
-  ==
-    fvc::reconstruct
     (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        parcels.SU(U)
+    );
+
+    UEqn.relax();
+
+    if (pimple.momentumPredictor())
+    {
+        solve
         (
-          - ghf*fvc::snGrad(rho)
+            UEqn
-          - fvc::snGrad(p_rgh)
+          ==
-        )*mesh.magSf()
+            fvc::reconstruct
-    )
+            (
-);
+                (
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            )
+        );
+    }

applications/solvers/combustion/fireFoam/YhsEqn.H

@@ -0,0 +1,70 @@
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_hs)")
+    )
+);
+{
+    radiation->correct();
+    combustion->correct();
+    dQ = combustion->dQ();
+    label inertIndex = -1;
+    volScalarField Yt(0.0*Y[0]);
+
+    forAll(Y, i)
+    {
+        if (Y[i].name() != inertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+            fvScalarMatrix R(combustion->R(Yi));
+
+            fvScalarMatrix YiEqn
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->alphaEff(), Yi)
+              ==
+                parcels.SYi(i, Yi)
+              + surfaceFilm.Srho(i)
+              + R
+            );
+
+            YiEqn.relax();
+            YiEqn.solve(mesh.solver("Yi"));
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+        else
+        {
+            inertIndex = i;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+      + dQ
+      + radiation->Shs(thermo)
+      + parcels.Sh(hs)
+      + surfaceFilm.Sh()
+    );
+
+    hsEqn.relax();
+    hsEqn.solve();
+
+    thermo.correct();
+
+    Info<< "min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
+}

applications/solvers/combustion/fireFoam/createClouds.H

@@ -0,0 +1,9 @@
+Info<< "\nConstructing reacting cloud" << endl;
+basicReactingCloud parcels
+(
+    "reactingCloud1",
+    rho,
+    U,
+    g,
+    slgThermo
+);

applications/solvers/combustion/fireFoam/createFields.H

@@ -1,147 +1,154 @@
-Info<< "Reading thermophysical properties\n" << endl;
+    Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<hsCombustionThermo> pThermo
+    autoPtr<hsCombustionThermo> pThermo
-(
-    hsCombustionThermo::New(mesh)
-);
-
-hsCombustionThermo& thermo = pThermo();
-
-basicMultiComponentMixture& composition = thermo.composition();
-
-volScalarField rho
-(
-    IOobject
     (
-        "rho",
+        hsCombustionThermo::New(mesh)
-        runTime.timeName(),
+    );
+    hsCombustionThermo& thermo = pThermo();
+
+    SLGThermo slgThermo(mesh, thermo);
+
+    basicMultiComponentMixture& composition = thermo.composition();
+    PtrList<volScalarField>& Y = composition.Y();
+
+    const word inertSpecie(thermo.lookup("inertSpecie"));
+
+    Info<< "Creating field rho\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    volScalarField& p = thermo.p();
+    volScalarField& hs = thermo.hs();
+    const volScalarField& T = thermo.T();
+    const volScalarField& psi = thermo.psi();
+
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "compressibleCreatePhi.H"
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    IOdictionary combustionProperties
+    (
+        IOobject
+        (
+            "combustionProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    );
+
+    Info<< "Creating combustion model\n" << endl;
+    autoPtr<combustionModel> combustion
+    (
+        combustionModel::combustionModel::New
+        (
+            combustionProperties,
+            thermo,
+            turbulence(),
+            phi,
+            rho
+        )
+    );
+
+    volScalarField dQ
+    (
+        IOobject
+        (
+            "dQ",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
         mesh,
-        IOobject::NO_READ,
+        dimensionedScalar("dQ", dimMass/pow3(dimTime)/dimLength, 0.0)
-        IOobject::AUTO_WRITE
+    );
-    ),
-    thermo.rho()
-);
 
-dimensionedScalar stoicRatio
+    Info<< "Creating field DpDt\n" << endl;
-(
+    volScalarField DpDt
-    thermo.lookup("stoichiometricAirFuelMassRatio")
-);
-
-volScalarField& p = thermo.p();
-volScalarField& hs = thermo.hs();
-
-const volScalarField& psi = thermo.psi();
-
-volScalarField& ft = composition.Y("ft");
-volScalarField& fu = composition.Y("fu");
-
-
-Info<< "Reading field U\n" << endl;
-
-volVectorField U
-(
-    IOobject
     (
-        "U",
+        "DpDt",
-        runTime.timeName(),
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
-        mesh,
+    );
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
 
-#include "compressibleCreatePhi.H"
 
-Info<< "Creating turbulence model\n" << endl;
+    Info<< "Calculating field g.h\n" << endl;
-autoPtr<compressible::turbulenceModel> turbulence
+    volScalarField gh("gh", g & mesh.C());
-(
-    compressible::turbulenceModel::New(rho, U, phi, thermo)
-);
 
-IOdictionary combustionProperties
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
-(
+
-    IOobject
+    volScalarField p_rgh
     (
-        "combustionProperties",
+        IOobject
-        runTime.constant(),
+        (
-        mesh,
+            "p_rgh",
-        IOobject::MUST_READ_IF_MODIFIED,
+            runTime.timeName(),
-        IOobject::NO_WRITE
+            mesh,
-    )
+            IOobject::MUST_READ,
-);
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
 
-Info<< "Creating combustion model\n" << endl;
+    // Force p_rgh to be consistent with p
-autoPtr<combustionModel> combustion
+    p_rgh = p - rho*gh;
-(
+
-    combustionModel::combustionModel::New
+    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+    forAll(Y, i)
+    {
+        fields.add(Y[i]);
+    }
+    fields.add(hs);
+
+    IOdictionary additionalControlsDict
     (
-        combustionProperties,
+        IOobject
-        thermo,
+        (
-        turbulence(),
+            "additionalControls",
-        phi,
+            runTime.constant(),
-        rho
+            mesh,
-    )
+            IOobject::MUST_READ_IF_MODIFIED,
-);
+            IOobject::NO_WRITE
+        )
+    );
 
-
+    Switch solvePrimaryRegion
-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",
+        additionalControlsDict.lookup("solvePrimaryRegion")
-        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
-    (
-        "dQ",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh,
-    dimensionedScalar("dQ", dimMass/pow3(dimTime)/dimLength, 0.0)
-);
-
-
-Info<< "Creating field DpDt\n" << endl;
-volScalarField DpDt
-(
-    fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
-);
-
-
-dimensionedScalar initialMass = fvc::domainIntegrate(rho);
-
-
-multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
-
-if (composition.contains("ft"))
-{
-    fields.add(composition.Y("ft"));
-}
-
-if (composition.contains("fu"))
-{
-    fields.add(composition.Y("fu"));
-}
-
-fields.add(hs);

applications/solvers/combustion/fireFoam/createPyrolysisModel.H

@@ -0,0 +1,6 @@
+Info<< "Creating pyrolysis model" << endl;
+
+autoPtr<regionModels::pyrolysisModels::pyrolysisModel> pyrolysis
+(
+    regionModels::pyrolysisModels::pyrolysisModel::New(mesh)
+);

applications/solvers/combustion/fireFoam/createSurfaceFilmModel.H

@@ -0,0 +1,6 @@
+Info<< "\nConstructing surface film model" << endl;
+
+typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;
+
+autoPtr<filmModelType> tsurfaceFilm(filmModelType::New(mesh, g));
+filmModelType& surfaceFilm = tsurfaceFilm();

applications/solvers/combustion/fireFoam/fireFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,15 +25,22 @@ Application
     fireFoam
 
 Description
-    Transient Solver for Fires and turbulent diffusion flames
+    Transient PIMPLE solver for Fires and turbulent diffusion flames with
+    reacting Lagrangian parcels, surface film and pyrolysis modelling.
 
 \*---------------------------------------------------------------------------*/
 
+#include "mapDistribute.H"
 #include "fvCFD.H"
-#include "hsCombustionThermo.H"
 #include "turbulenceModel.H"
-#include "combustionModel.H"
+#include "basicReactingCloud.H"
+#include "surfaceFilmModel.H"
+#include "pyrolysisModel.H"
 #include "radiationModel.H"
+#include "SLGThermo.H"
+#include "hsCombustionThermo.H"
+#include "solidChemistryModel.H"
+#include "combustionModel.H"
 #include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -41,19 +48,25 @@ Description
 int main(int argc, char *argv[])
 {
     #include "setRootCase.H"
+
     #include "createTime.H"
     #include "createMesh.H"
+    #include "readChemistryProperties.H"
     #include "readGravitationalAcceleration.H"
-    #include "initContinuityErrs.H"
     #include "createFields.H"
+    #include "createClouds.H"
+    #include "createSurfaceFilmModel.H"
+    #include "createPyrolysisModel.H"
     #include "createRadiationModel.H"
+    #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
+    #include "readPyrolysisTimeControls.H"
 
     pimpleControl pimple(mesh);
 
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -61,46 +74,56 @@ int main(int argc, char *argv[])
     {
         #include "readTimeControls.H"
         #include "compressibleCourantNo.H"
+        #include "solidRegionDiffusionNo.H"
+        #include "setMultiRegionDeltaT.H"
         #include "setDeltaT.H"
 
         runTime++;
+
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "rhoEqn.H"
+        parcels.evolve();
 
-        // --- Pressure-velocity PIMPLE corrector loop
+        surfaceFilm.evolve();
-        for (pimple.start(); pimple.loop(); pimple++)
+
+        pyrolysis->evolve();
+
+        if (solvePrimaryRegion)
         {
-            #include "UEqn.H"
+            #include "rhoEqn.H"
 
-            #include "ftEqn.H"
+            // --- PIMPLE loop
-            #include "fuhsEqn.H"
+            for (pimple.start(); pimple.loop(); pimple++)
-
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
-                #include "pEqn.H"
+                #include "UEqn.H"
+                #include "YhsEqn.H"
+
+                // --- PISO loop
+                for (int corr=1; corr<=pimple.nCorr(); corr++)
+                {
+                    #include "pEqn.H"
+                }
+
+                if (pimple.turbCorr())
+                {
+                    turbulence->correct();
+                }
             }
 
-            if (pimple.turbCorr())
+            rho = thermo.rho();
-            {
-                turbulence->correct();
-            }
         }
 
-        rho = thermo.rho();
-
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
-
     }
 
-    Info<< "End\n" << endl;
+    Info<< "End" << endl;
 
-    return 0;
+    return(0);
 }
 
+
 // ************************************************************************* //

applications/solvers/combustion/fireFoam/ftEqn.H

@@ -1,25 +0,0 @@
-tmp<fv::convectionScheme<scalar> > mvConvection
-(
-    fv::convectionScheme<scalar>::New
-    (
-        mesh,
-        fields,
-        phi,
-        mesh.divScheme("div(phi,ft_b_h)")
-    )
-);
-
-{
-    fvScalarMatrix ftEqn
-    (
-        fvm::ddt(rho, ft)
-      + mvConvection->fvmDiv(phi, ft)
-      - fvm::laplacian(turbulence->alphaEff(), ft)
-    );
-
-    ftEqn.relax();
-    ftEqn.solve();
-}
-
-Info<< "max(ft) = " << max(ft).value() << endl;
-Info<< "min(ft) = " << min(ft).value() << endl;

applications/solvers/combustion/fireFoam/fuhsEqn.H

@@ -1,47 +0,0 @@
-{
-    // Solve fuel equation
-    // ~~~~~~~~~~~~~~~~~~~
-    fvScalarMatrix R(combustion->R(fu));
-
-    {
-        fvScalarMatrix fuEqn
-        (
-            fvm::ddt(rho, fu)
-          + mvConvection->fvmDiv(phi, fu)
-          - fvm::laplacian(turbulence->alphaEff(), fu)
-         ==
-            R
-        );
-
-        fuEqn.relax();
-        fuEqn.solve();
-    }
-
-    Info<< "max(fu) = " << max(fu).value() << endl;
-    Info<< "min(fu) = " << min(fu).value() << endl;
-
-
-    // Solve sensible enthalpy equation
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    radiation->correct();
-    dQ = combustion->dQ(R);
-
-    {
-        fvScalarMatrix hsEqn
-        (
-            fvm::ddt(rho, hs)
-          + mvConvection->fvmDiv(phi,hs)
-          - fvm::laplacian(turbulence->alphaEff(), hs)
-         ==
-            DpDt
-          + dQ
-          + radiation->Shs(thermo)
-        );
-
-        hsEqn.relax();
-        hsEqn.solve();
-    }
-
-    thermo.correct();
-    combustion->correct();
-}

applications/solvers/combustion/fireFoam/pEqn.H

@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf(rAU.name() + 'f', fvc::interpolate(rho*rAU));
 U = rAU*UEqn.H();
 
 surfaceScalarField phiU
@@ -17,13 +17,15 @@ phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
 for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
 {
-    surfaceScalarField rhorAUf(fvc::interpolate(rho*rAU));
-
     fvScalarMatrix p_rghEqn
     (
-        fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh
+        fvc::ddt(psi, rho)*gh
       + fvc::div(phi)
+      + fvm::ddt(psi, p_rgh)
       - fvm::laplacian(rhorAUf, p_rgh)
+     ==
+        parcels.Srho()
+      + surfaceFilm.Srho()
     );
 
     p_rghEqn.solve

applications/solvers/combustion/fireFoam/readChemistryProperties.H

@@ -0,0 +1,23 @@
+Info<< "Reading chemistry properties\n" << endl;
+
+IOdictionary chemistryProperties
+(
+    IOobject
+    (
+        "chemistryProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE,
+        false
+    )
+);
+
+Switch turbulentReaction(chemistryProperties.lookup("turbulentReaction"));
+
+dimensionedScalar Cmix("Cmix", dimless, 1.0);
+
+if (turbulentReaction)
+{
+    chemistryProperties.lookup("Cmix") >> Cmix;
+}

applications/solvers/combustion/fireFoam/readPyrolysisTimeControls.H

@@ -0,0 +1,34 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2011 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/>.
+
+Global
+    readPyrolysisTimeControls
+
+Description
+
+
+\*---------------------------------------------------------------------------*/
+
+scalar maxDi = pyrolysis->maxDiff();
+
+// ************************************************************************* //

applications/solvers/combustion/fireFoam/rhoEqn.H

@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2011 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/>.
+
+Global
+    rhoEqn
+
+Description
+    Solve the continuity for density.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    solve
+    (
+        fvm::ddt(rho)
+      + fvc::div(phi)
+      ==
+        parcels.Srho(rho)
+      + surfaceFilm.Srho()
+    );
+}
+
+// ************************************************************************* //

applications/solvers/combustion/fireFoam/setMultiRegionDeltaT.H

@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2011 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/>.
+
+Global
+    setMultiRegionDeltaT
+
+Description
+    Reset the timestep to maintain a constant maximum Courant numbers.
+    Reduction of time-step is immediate, but increase is damped to avoid
+    unstable oscillations.
+
+\*---------------------------------------------------------------------------*/
+
+if (adjustTimeStep)
+{
+    if (CoNum == -GREAT)
+    {
+        CoNum = SMALL;
+    }
+
+    if (DiNum == -GREAT)
+    {
+        DiNum = SMALL;
+    }
+
+
+    const scalar TFactorFluid = maxCo/(CoNum + SMALL);
+    const scalar TFactorSolid = maxDi/(DiNum + SMALL);
+    const scalar TFactorFilm = maxCo/(surfaceFilm.CourantNumber() + SMALL);
+
+    const scalar dt0 = runTime.deltaTValue();
+
+    runTime.setDeltaT
+    (
+        min
+        (
+            dt0*min(min(TFactorFluid, min(TFactorFilm, TFactorSolid)), 1.2),
+            maxDeltaT
+        )
+    );
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/fireFoam/solidRegionDiffusionNo.H

@@ -0,0 +1 @@
+scalar DiNum =  pyrolysis->solidRegionDiffNo();

applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/Make/files

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

applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/Make/options

@@ -0,0 +1,20 @@
+EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/thermoBaffleModels/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lthermoBaffleModels \
+    -lregionModels

applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/UEqn.H

@@ -0,0 +1,25 @@
+    // Solve the Momentum equation
+
+    tmp<fvVectorMatrix> UEqn
+    (
+        fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+    );
+
+    UEqn().relax();
+
+    if (simple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn()
+         ==
+            fvc::reconstruct
+            (
+                (
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            )
+        );
+    }

applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/buoyantBaffleSimpleFoam.C

@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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
+    buoyantBaffleSimpleFoam
+
+Description
+    Steady-state solver for buoyant, turbulent flow of compressible fluids
+    using thermal baffles
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "basicPsiThermo.H"
+#include "RASModel.H"
+#include "fixedGradientFvPatchFields.H"
+#include "simpleControl.H"
+#include "thermoBaffleModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+
+    simpleControl simple(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        p_rgh.storePrevIter();
+        rho.storePrevIter();
+
+        // Pressure-velocity SIMPLE corrector
+        {
+            #include "UEqn.H"
+            #include "hEqn.H"
+            #include "pEqn.H"
+        }
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/createFields.H

@@ -0,0 +1,94 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<basicPsiThermo> pThermo
+    (
+        basicPsiThermo::New(mesh)
+    );
+    basicPsiThermo& thermo = pThermo();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "compressibleCreatePhi.H"
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::RASModel> turbulence
+    (
+        compressible::RASModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+
+    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;
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell
+    (
+        p,
+        p_rgh,
+        mesh.solutionDict().subDict("SIMPLE"),
+        pRefCell,
+        pRefValue
+    );
+
+    autoPtr<regionModels::thermoBaffleModels::thermoBaffleModel> baffles
+    (
+        regionModels::thermoBaffleModels::thermoBaffleModel::New(mesh)
+    );
+
+    dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+    dimensionedScalar totalVolume = sum(mesh.V());

applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/hEqn.H

@@ -0,0 +1,17 @@
+{
+    fvScalarMatrix hEqn
+    (
+        fvm::div(phi, h)
+      - fvm::Sp(fvc::div(phi), h)
+      - fvm::laplacian(turbulence->alphaEff(), h)
+     ==
+        fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p))
+      - p*fvc::div(phi/fvc::interpolate(rho))
+    );
+
+    hEqn.relax();
+    hEqn.solve();
+
+    baffles->evolve();
+    thermo.correct();
+}

applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/pEqn.H

@@ -0,0 +1,59 @@
+{
+    rho = thermo.rho();
+    rho.relax();
+
+    volScalarField rAU(1.0/UEqn().A());
+    surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
+
+    U = rAU*UEqn().H();
+    UEqn.clear();
+
+    phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
+    bool closedVolume = adjustPhi(phi, U, p_rgh);
+
+    surfaceScalarField buoyancyPhi(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+    phi -= buoyancyPhi;
+
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    {
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        p_rghEqn.solve();
+
+        if (nonOrth == simple.nNonOrthCorr())
+        {
+            // Calculate the conservative fluxes
+            phi -= p_rghEqn.flux();
+
+            // Explicitly relax pressure for momentum corrector
+            p_rgh.relax();
+
+            // Correct the momentum source with the pressure gradient flux
+            // calculated from the relaxed pressure
+            U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
+            U.correctBoundaryConditions();
+        }
+    }
+
+    #include "continuityErrs.H"
+
+    p = p_rgh + rho*gh;
+
+    // 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.relax();
+    Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
+        << endl;
+}

applications/solvers/heatTransfer/chtMultiRegionFoam/coupledDerivedFvPatchFields/Make/files

@@ -1,4 +1,3 @@
-turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C
 externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C
 
 LIB = $(FOAM_LIBBIN)/libcoupledDerivedFvPatchFields

applications/solvers/heatTransfer/chtMultiRegionFoam/coupledDerivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C

@@ -202,7 +202,7 @@ void Foam::externalWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    if(oldMode_ == fixedHeatFlux)
+    if (oldMode_ == fixedHeatFlux)
     {
         this->refGrad() = q_/K(*this);
         this->refValue() = 0.0;

applications/solvers/heatTransfer/chtMultiRegionFoam/coupledDerivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.H

@@ -54,8 +54,6 @@ SourceFiles
 #ifndef solidWallHeatFluxTemperatureFvPatchScalarField_H
 #define solidWallHeatFluxTemperatureFvPatchScalarField_H
 
-//#include "fixedGradientFvPatchFields.H"
-
 #include "mixedFvPatchFields.H"
 #include "temperatureCoupledBase.H"
 

applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H

@@ -88,7 +88,7 @@
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
 
-    surfaceScalarField ghf("gh", g & mesh.Cf());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
 
     volScalarField p_rgh
     (

applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H

@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf(rAU.name(), fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf(rAU.name() + 'f', fvc::interpolate(rho*rAU));
 U = rAU*UEqn.H();
 
 surfaceScalarField phiU

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

@@ -132,6 +132,7 @@ Usage
 #include "surfaceFields.H"
 #include "syncTools.H"
 #include "cyclicPolyPatch.H"
+#include "wedgePolyPatch.H"
 #include "nonuniformTransformCyclicPolyPatch.H"
 #include "extrudeModel.H"
 
@@ -951,6 +952,12 @@ int main(int argc, char *argv[])
     // Region
     const word shellRegionName(dict.lookup("region"));
     const wordList zoneNames(dict.lookup("faceZones"));
+    wordList zoneShadowNames(0);
+    if (dict.found("faceZonesShadow"))
+    {
+        dict.lookup("faceZonesShadow") >> zoneShadowNames;
+    }
+
     const Switch oneD(dict.lookup("oneD"));
     const Switch adaptMesh(dict.lookup("adaptMesh"));
 
@@ -1007,6 +1014,47 @@ int main(int argc, char *argv[])
         }
     }
 
+    labelList zoneShadowIDs;
+    if (zoneShadowNames.size())
+    {
+        zoneShadowIDs.setSize(zoneShadowNames.size());
+        forAll(zoneShadowNames, i)
+        {
+            zoneShadowIDs[i] = faceZones.findZoneID(zoneShadowNames[i]);
+            if (zoneShadowIDs[i] == -1)
+            {
+                FatalErrorIn(args.executable())
+                    << "Cannot find zone " << zoneShadowNames[i] << endl
+                    << "Valid zones are " << faceZones.names()
+                    << exit(FatalError);
+            }
+        }
+    }
+
+    label nShadowFaces = 0;
+    forAll(zoneShadowIDs, i)
+    {
+        nShadowFaces += faceZones[zoneShadowIDs[i]].size();
+    }
+
+    labelList extrudeMeshShadowFaces(nShadowFaces);
+    boolList zoneShadowFlipMap(nShadowFaces);
+    labelList zoneShadowID(nShadowFaces);
+
+    nShadowFaces = 0;
+    forAll(zoneShadowIDs, i)
+    {
+        const faceZone& fz = faceZones[zoneShadowIDs[i]];
+        forAll(fz, j)
+        {
+            extrudeMeshShadowFaces[nShadowFaces] = fz[j];
+            zoneShadowFlipMap[nShadowFaces] = fz.flipMap()[j];
+            zoneShadowID[nShadowFaces] = zoneShadowIDs[i];
+            nShadowFaces++;
+        }
+    }
+
+
 
     // Collect faces to extrude and per-face information
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1047,6 +1095,19 @@ int main(int argc, char *argv[])
         << nl
         << endl;
 
+     // Check nExtrudeFaces = nShadowFaces
+    if (zoneShadowNames.size())
+    {
+        if (nExtrudeFaces != nShadowFaces)
+        {
+            FatalErrorIn(args.executable())
+                << "Extruded faces " << nExtrudeFaces << endl
+                << "is different from shadow faces. " << nShadowFaces
+                << "This is not permitted " << endl
+                << exit(FatalError);
+        }
+    }
+
 
     // Determine corresponding mesh edges
     const labelList extrudeMeshEdges
@@ -1135,7 +1196,7 @@ int main(int argc, char *argv[])
                 << '\t' << patches[interRegionBottomPatch[i]].type()
                 << nl;
         }
-        else
+        else if (zoneShadowNames.size() == 0)
         {
             interRegionTopPatch[i] = addPatch<polyPatch>
             (
@@ -1159,6 +1220,31 @@ int main(int argc, char *argv[])
                 << '\t' << patches[interRegionBottomPatch[i]].type()
                 << nl;
         }
+        else if (zoneShadowNames.size() > 0) //patch using shadow face zones.
+        {
+            interRegionTopPatch[i] = addPatch<directMappedWallPolyPatch>
+            (
+                mesh,
+                zoneShadowNames[i] + "_top"
+            );
+            nCoupled++;
+            Info<< interRegionTopPatch[i]
+                << '\t' << patches[interRegionTopPatch[i]].name()
+                << '\t' << patches[interRegionTopPatch[i]].type()
+                << nl;
+
+            interRegionBottomPatch[i] = addPatch<directMappedWallPolyPatch>
+            (
+                mesh,
+                interName
+            );
+            nCoupled++;
+            Info<< interRegionBottomPatch[i]
+                << '\t' << patches[interRegionBottomPatch[i]].name()
+                << '\t' << patches[interRegionBottomPatch[i]].type()
+                << nl;
+        }
+
     }
     Info<< "Added " << nCoupled << " inter-region patches." << nl
         << endl;
@@ -1216,13 +1302,32 @@ int main(int argc, char *argv[])
         if (oneD)
         {
             // Reuse single empty patch.
-            word patchName = "oneDEmptPatch";
+            word patchType = dict.lookup("oneDPolyPatchType");
-
+            word patchName;
-            zoneSidePatch[zoneI] = addPatch<emptyPolyPatch>
+            if (patchType == "emptyPolyPatch")
-            (
+            {
-                mesh,
+                patchName = "oneDEmptyPatch";
-                patchName
+                zoneSidePatch[zoneI] = addPatch<emptyPolyPatch>
-            );
+                (
+                    mesh,
+                    patchName
+                );
+            }
+            else if (patchType == "wedgePolyPatch")
+            {
+                patchName = "oneDWedgePatch";
+                zoneSidePatch[zoneI] = addPatch<wedgePolyPatch>
+                (
+                    mesh,
+                    patchName
+                );
+            }
+            else
+            {
+                FatalErrorIn(args.executable())
+                << "Type " << patchType << " does not exist "
+                << exit(FatalError);
+            }
 
             Info<< zoneSidePatch[zoneI] << '\t' << patchName << nl;
 
@@ -1335,12 +1440,16 @@ int main(int argc, char *argv[])
 
         if (oneD)
         {
-            nonManifoldEdge[edgeI] = 1;
+            //nonManifoldEdge[edgeI] = 1; //To fill the space
             ePatches.setSize(eFaces.size());
             forAll(eFaces, i)
             {
                 ePatches[i] = zoneSidePatch[zoneID[eFaces[i]]];
             }
+            if (eFaces.size() != 2)
+            {
+                nonManifoldEdge[edgeI] = 1;
+            }
         }
         else if (eFaces.size() == 2)
         {
@@ -1754,48 +1863,91 @@ int main(int argc, char *argv[])
             }
         }
 
-        // Add faces (using same points) to be in top patch
+        if (zoneShadowNames.size() > 0) //if there is a top faceZone specified
-        forAll(extrudeMeshFaces, zoneFaceI)
         {
-            label meshFaceI = extrudeMeshFaces[zoneFaceI];
+            forAll(extrudeMeshFaces, zoneFaceI)
-            bool flip = zoneFlipMap[zoneFaceI];
-            const face& f = mesh.faces()[meshFaceI];
-
-            if (!flip)
             {
-                if (mesh.isInternalFace(meshFaceI))
+                label meshFaceI = extrudeMeshShadowFaces[zoneFaceI];
+                label zoneI = zoneShadowID[zoneFaceI];
+                bool flip = zoneShadowFlipMap[zoneFaceI];
+                const face& f = mesh.faces()[meshFaceI];
+
+                if (!flip)
+                {
+                    meshMod.modifyFace
+                    (
+                        f,                          // modified face
+                        meshFaceI,                  // face being modified
+                        mesh.faceOwner()[meshFaceI],// owner
+                        -1,                         // neighbour
+                        false,                      // face flip
+                        extrudeTopPatchID[zoneFaceI],// patch for face
+                        zoneI,                      // zone for face
+                        flip                        // face flip in zone
+                    );
+                }
+                else if (mesh.isInternalFace(meshFaceI))
+                {
+                    meshMod.modifyFace
+                    (
+                        f.reverseFace(),                // modified face
+                        meshFaceI,                      // label modified face
+                        mesh.faceNeighbour()[meshFaceI],// owner
+                        -1,                             // neighbour
+                        true,                           // face flip
+                        extrudeTopPatchID[zoneFaceI],   // patch for face
+                        zoneI,                          // zone for face
+                        !flip                           // face flip in zone
+                    );
+                }
+            }
+
+        }
+        else
+        {
+            // Add faces (using same points) to be in top patch
+            forAll(extrudeMeshFaces, zoneFaceI)
+            {
+                label meshFaceI = extrudeMeshFaces[zoneFaceI];
+                bool flip = zoneFlipMap[zoneFaceI];
+                const face& f = mesh.faces()[meshFaceI];
+
+                if (!flip)
+                {
+                    if (mesh.isInternalFace(meshFaceI))
+                    {
+                        meshMod.addFace
+                        (
+                            f.reverseFace(),                // modified face
+                            mesh.faceNeighbour()[meshFaceI],// owner
+                            -1,                             // neighbour
+                            -1,                             // master point
+                            -1,                             // master edge
+                            meshFaceI,                      // master face
+                            true,                           // flip flux
+                            extrudeTopPatchID[zoneFaceI],   // patch for face
+                            -1,                             // zone for face
+                            false                           //face flip in zone
+                        );
+                    }
+                }
+                else
                 {
                     meshMod.addFace
                     (
-                        f.reverseFace(),                // modified face
+                        f,                              // face
-                        mesh.faceNeighbour()[meshFaceI],// owner
+                        mesh.faceOwner()[meshFaceI],    // owner
                         -1,                             // neighbour
                         -1,                             // master point
                         -1,                             // master edge
                         meshFaceI,                      // master face
-                        true,                           // flip flux
+                        false,                          // flip flux
                         extrudeTopPatchID[zoneFaceI],   // patch for face
                         -1,                             // zone for face
-                        false                           // face flip in zone
+                        false                           // zone flip
                     );
                 }
             }
-            else
-            {
-                meshMod.addFace
-                (
-                    f,                              // face
-                    mesh.faceOwner()[meshFaceI],    // owner
-                    -1,                             // neighbour
-                    -1,                             // master point
-                    -1,                             // master edge
-                    meshFaceI,                      // master face
-                    false,                          // flip flux
-                    extrudeTopPatchID[zoneFaceI],   // patch for face
-                    -1,                             // zone for face
-                    false                           // zone flip
-                );
-            }
         }
 
         // Change the mesh. Change points directly (no inflation).

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -144,6 +144,13 @@ castellatedMeshControls
                 }
             }
 
+            // Optional specification of patch type (default is wall). No
+            // constraint types (cyclic, symmetry) etc. are allowed.
+            patchInfo
+            {
+                type patch;
+            }
+
             //- Optional angle to detect small-large cell situation
             //  perpendicular to the surface. Is the angle of face w.r.t.
             //  the local surface normal. Use on flat(ish) surfaces only.

applications/utilities/miscellaneous/foamDebugSwitches/Make/options

@@ -43,10 +43,10 @@ EXE_LIBS = \
     -lmeshTools \
     -lmolecularMeasurements \
     -lmolecule \
-    -lmultiphaseInterFoam \
+/*    -lmultiphaseInterFoam  */ \
     -lODE \
     -lOpenFOAM \
-    -lphaseModel \
+/*    -lphaseModel  */ \
     -lpotential \
     -lradiationModels \
     -lrandomProcesses \
@@ -62,7 +62,6 @@ EXE_LIBS = \
     -lsurfaceFilmModels \
     -lsurfMesh \
     -lsystemCall \
-    -ltabulatedWallFunctions \
     -lthermalPorousZone \
     -lthermophysicalFunctions \
     -ltopoChangerFvMesh \

applications/utilities/parallelProcessing/decomposePar/Make/files

@@ -2,6 +2,7 @@ decomposePar.C
 domainDecomposition.C
 domainDecompositionMesh.C
 domainDecompositionDistribute.C
+dimFieldDecomposer.C
 fvFieldDecomposer.C
 pointFieldDecomposer.C
 lagrangianFieldDecomposer.C

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -81,6 +81,7 @@ Usage
 #include "pointFields.H"
 
 #include "readFields.H"
+#include "dimFieldDecomposer.H"
 #include "fvFieldDecomposer.H"
 #include "pointFieldDecomposer.H"
 #include "lagrangianFieldDecomposer.H"
@@ -354,6 +355,25 @@ int main(int argc, char *argv[])
     readFields(mesh, objects, volTensorFields);
 
 
+    // Construct the dimensioned fields
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    PtrList<DimensionedField<scalar, volMesh> > dimScalarFields;
+    readFields(mesh, objects, dimScalarFields);
+
+    PtrList<DimensionedField<vector, volMesh> > dimVectorFields;
+    readFields(mesh, objects, dimVectorFields);
+
+    PtrList<DimensionedField<sphericalTensor, volMesh> >
+        dimSphericalTensorFields;
+    readFields(mesh, objects, dimSphericalTensorFields);
+
+    PtrList<DimensionedField<symmTensor, volMesh> > dimSymmTensorFields;
+    readFields(mesh, objects, dimSymmTensorFields);
+
+    PtrList<DimensionedField<tensor, volMesh> > dimTensorFields;
+    readFields(mesh, objects, dimTensorFields);
+
+
     // Construct the surface fields
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     PtrList<surfaceScalarField> surfaceScalarFields;
@@ -727,19 +747,6 @@ int main(int argc, char *argv[])
         );
 
         // FV fields
-        if
-        (
-            volScalarFields.size()
-         || volVectorFields.size()
-         || volSphericalTensorFields.size()
-         || volSymmTensorFields.size()
-         || volTensorFields.size()
-         || surfaceScalarFields.size()
-         || surfaceVectorFields.size()
-         || surfaceSphericalTensorFields.size()
-         || surfaceSymmTensorFields.size()
-         || surfaceTensorFields.size()
-        )
         {
             fvFieldDecomposer fieldDecomposer
             (
@@ -763,16 +770,25 @@ int main(int argc, char *argv[])
             fieldDecomposer.decomposeFields(surfaceTensorFields);
         }
 
+        // Dimensioned fields
+        {
+            dimFieldDecomposer fieldDecomposer
+            (
+                mesh,
+                procMesh,
+                faceProcAddressing,
+                cellProcAddressing
+            );
+
+            fieldDecomposer.decomposeFields(dimScalarFields);
+            fieldDecomposer.decomposeFields(dimVectorFields);
+            fieldDecomposer.decomposeFields(dimSphericalTensorFields);
+            fieldDecomposer.decomposeFields(dimSymmTensorFields);
+            fieldDecomposer.decomposeFields(dimTensorFields);
+        }
+
 
         // Point fields
-        if
-        (
-            pointScalarFields.size()
-         || pointVectorFields.size()
-         || pointSphericalTensorFields.size()
-         || pointSymmTensorFields.size()
-         || pointTensorFields.size()
-        )
         {
             labelIOList pointProcAddressing
             (
@@ -822,21 +838,6 @@ int main(int argc, char *argv[])
                 );
 
                 // Lagrangian fields
-                if
-                (
-                    lagrangianLabelFields[cloudI].size()
-                 || lagrangianLabelFieldFields[cloudI].size()
-                 || lagrangianScalarFields[cloudI].size()
-                 || lagrangianScalarFieldFields[cloudI].size()
-                 || lagrangianVectorFields[cloudI].size()
-                 || lagrangianVectorFieldFields[cloudI].size()
-                 || lagrangianSphericalTensorFields[cloudI].size()
-                 || lagrangianSphericalTensorFieldFields[cloudI].size()
-                 || lagrangianSymmTensorFields[cloudI].size()
-                 || lagrangianSymmTensorFieldFields[cloudI].size()
-                 || lagrangianTensorFields[cloudI].size()
-                 || lagrangianTensorFieldFields[cloudI].size()
-                )
                 {
                     fieldDecomposer.decomposeFields
                     (

applications/utilities/parallelProcessing/decomposePar/decomposeParDict

@@ -37,6 +37,7 @@ method          scotch;
 // method          metis;
 // method          manual;
 // method          multiLevel;
+// method          structured;  // does 2D decomposition of structured mesh
 
 multiLevelCoeffs
 {
@@ -108,6 +109,15 @@ manualCoeffs
 }
 
 
+structuredCoeffs
+{
+    // Patches to do 2D decomposition on. Structured mesh only; cells have
+    // to be in 'columns' on top of patches.
+    patches     (bottomPatch);
+}
+
+
+
 //// Is the case distributed
 //distributed     yes;
 //// Per slave (so nProcs-1 entries) the directory above the case.

applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.C

@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "dimFieldDecomposer.H"
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::dimFieldDecomposer::dimFieldDecomposer
+(
+    const fvMesh& completeMesh,
+    const fvMesh& procMesh,
+    const labelList& faceAddressing,
+    const labelList& cellAddressing
+)
+:
+    completeMesh_(completeMesh),
+    procMesh_(procMesh),
+    faceAddressing_(faceAddressing),
+    cellAddressing_(cellAddressing)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::dimFieldDecomposer::~dimFieldDecomposer()
+{}
+
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.H

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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::dimFieldDecomposer
+
+Description
+    Dimensioned field decomposer.
+
+SourceFiles
+    dimFieldDecomposer.C
+    dimFieldDecomposerDecomposeFields.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dimFieldDecomposer_H
+#define dimFieldDecomposer_H
+
+#include "fvMesh.H"
+#include "fvPatchFieldMapper.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class IOobjectList;
+
+/*---------------------------------------------------------------------------*\
+                    Class fvFieldDecomposer Declaration
+\*---------------------------------------------------------------------------*/
+
+class dimFieldDecomposer
+{
+private:
+
+    // Private data
+
+        //- Reference to complete mesh
+        const fvMesh& completeMesh_;
+
+        //- Reference to processor mesh
+        const fvMesh& procMesh_;
+
+        //- Reference to face addressing
+        const labelList& faceAddressing_;
+
+        //- Reference to cell addressing
+        const labelList& cellAddressing_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        dimFieldDecomposer(const dimFieldDecomposer&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const dimFieldDecomposer&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        dimFieldDecomposer
+        (
+            const fvMesh& completeMesh,
+            const fvMesh& procMesh,
+            const labelList& faceAddressing,
+            const labelList& cellAddressing
+        );
+
+
+    //- Destructor
+    ~dimFieldDecomposer();
+
+
+    // Member Functions
+
+        //- Decompose field
+        template<class Type>
+        tmp<DimensionedField<Type, volMesh> > decomposeField
+        (
+            const DimensionedField<Type, volMesh>& field
+        ) const;
+
+
+        //- Decompose llist of fields
+        template<class GeoField>
+        void decomposeFields(const PtrList<GeoField>& fields) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "dimFieldDecomposerDecomposeFields.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposerDecomposeFields.C

@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "dimFieldDecomposer.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::DimensionedField<Type, Foam::volMesh> >
+Foam::dimFieldDecomposer::decomposeField
+(
+    const DimensionedField<Type, volMesh>& field
+) const
+{
+    // Create and map the internal field values
+    Field<Type> mappedField(field, cellAddressing_);
+
+    // Create the field for the processor
+    return tmp<DimensionedField<Type, volMesh> >
+    (
+        new DimensionedField<Type, volMesh>
+        (
+            IOobject
+            (
+                field.name(),
+                procMesh_.time().timeName(),
+                procMesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            procMesh_,
+            field.dimensions(),
+            mappedField
+        )
+    );
+}
+
+
+template<class GeoField>
+void Foam::dimFieldDecomposer::decomposeFields
+(
+    const PtrList<GeoField>& fields
+) const
+{
+    forAll(fields, fieldI)
+    {
+        decomposeField(fields[fieldI])().write();
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -860,6 +860,16 @@ bool Foam::domainDecomposition::writeDecomposition()
     scalar avgProcPatches = scalar(totProcPatches)/nProcs_;
     scalar avgProcFaces = scalar(totProcFaces)/nProcs_;
 
+    // In case of all faces on one processor. Just to avoid division by 0.
+    if (totProcPatches == 0)
+    {
+        avgProcPatches = 1;
+    }
+    if (totProcFaces == 0)
+    {
+        avgProcFaces = 1;
+    }
+
     Info<< nl
         << "Number of processor faces = " << totProcFaces/2 << nl
         << "Max number of cells = " << maxProcCells

applications/utilities/parallelProcessing/decomposePar/readFields.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,7 @@ void Foam::readFields
     PtrList<GeoField>& fields
 )
 {
-    // Search list of objects for volScalarFields
+    // Search list of objects for fields of type GeomField
     IOobjectList fieldObjects(objects.lookupClass(GeoField::typeName));
 
     // Remove the cellDist field
@@ -45,7 +45,7 @@ void Foam::readFields
         fieldObjects.erase(celDistIter);
     }
 
-    // Construct the vol scalar fields
+    // Construct the fields
     fields.setSize(fieldObjects.size());
 
     label fieldI = 0;

applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C

@@ -475,6 +475,16 @@ void printMeshData(const polyMesh& mesh)
     scalar avgProcPatches = scalar(totProcPatches)/Pstream::nProcs();
     scalar avgProcFaces = scalar(totProcFaces)/Pstream::nProcs();
 
+    // In case of all faces on one processor. Just to avoid division by 0.
+    if (totProcPatches == 0)
+    {
+        avgProcPatches = 1;
+    }
+    if (totProcFaces == 0)
+    {
+        avgProcFaces = 1;
+    }
+
     Info<< nl
         << "Number of processor faces = " << totProcFaces/2 << nl
         << "Max number of cells = " << maxProcCells

applications/utilities/surface/surfaceSubset/surfaceSubset.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -249,7 +249,11 @@ int main(int argc, char *argv[])
         // search engine
         indexedOctree<treeDataTriSurface> selectTree
         (
-            treeDataTriSurface(selectSurf),
+            treeDataTriSurface
+            (
+                selectSurf,
+                indexedOctree<treeDataTriSurface>::perturbTol()
+            ),
             bb.extend(rndGen, 1E-4),    // slightly randomize bb
             8,      // maxLevel
             10,     // leafsize

bin/paraFoam

@@ -229,7 +229,8 @@ else
         echo "created temporary '$caseFile'"
     }
 
-    paraview --data="$caseFile" "$@"
+    # For now filter out any ld.so errors. Caused by non-system compiler?
+    paraview --data="$caseFile" "$@" 2>&1 | fgrep -v 'Inconsistency detected by ld.so'
 fi
 
 

bin/tools/MakefileDirs

@@ -36,7 +36,7 @@
 application: $(FOAM_TARGETS)
 
 $(FOAM_TARGETS):
-	@(cd $@ && $(FOAM_APP))
+	+@(cd $@ && $(FOAM_APP))
 
 
 #------------------------------------------------------------------------------

bin/tools/RunFunctions

@@ -61,7 +61,14 @@ runParallel()
         nProcs=$1
         shift
         echo "Running $APP_RUN in parallel on $PWD using $nProcs processes"
-        ( mpirun -np $nProcs $APP_RUN -parallel $* < /dev/null > log.$APP_NAME 2>&1 )
+
+        #if [ "$WM_SCHEDULER" ]
+        #then
+        #    echo "$PWD: $WM_SCHEDULER -np $nProcs" 1>&2
+        #    $WM_SCHEDULER -np $nProcs "( mpirun -np $nProcs $APP_RUN -parallel $* < /dev/null > log.$APP_NAME 2>&1 )"
+        #else
+            ( mpirun -np $nProcs $APP_RUN -parallel $* < /dev/null > log.$APP_NAME 2>&1 )
+        #fi
     fi
 }
 

src/Allwmake

@@ -60,6 +60,7 @@ wmake $makeType randomProcesses
 thermophysicalModels/Allwmake $*
 transportModels/Allwmake $*
 turbulenceModels/Allwmake $*
+wmake $makeType combustionModels
 regionModels/Allwmake $*
 lagrangian/Allwmake $*
 postProcessing/Allwmake $*

src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C

@@ -714,7 +714,7 @@ Foam::point Foam::indexedOctree<Type>::pushPoint
 )
 {
     // Get local length scale.
-    const vector perturbVec = perturbTol_*(bb.span());
+    const vector perturbVec = perturbTol_*bb.span();
 
     point perturbedPt(pt);
 

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

@@ -525,7 +525,6 @@ void Foam::globalMeshData::calcGlobalPointSlaves() const
 
     // Calculate connected points for master points.
     globalPoints globalData(mesh_, coupledPatch(), true, true);
-    globalPointNumberingPtr_.reset(new globalIndex(globalData.globalIndices()));
 
     globalPointSlavesPtr_.reset
     (
@@ -1564,6 +1563,42 @@ void Foam::globalMeshData::calcGlobalPointBoundaryCells() const
 }
 
 
+void Foam::globalMeshData::calcGlobalCoPointSlaves() const
+{
+    if (debug)
+    {
+        Pout<< "globalMeshData::calcGlobalCoPointSlaves() :"
+            << " calculating coupled master to collocated"
+            << " slave point addressing." << endl;
+    }
+
+    // Calculate connected points for master points.
+    globalPoints globalData(mesh_, coupledPatch(), true, false);
+
+    globalCoPointSlavesPtr_.reset
+    (
+        new labelListList
+        (
+            globalData.pointPoints().xfer()
+        )
+    );
+    globalCoPointSlavesMapPtr_.reset
+    (
+        new mapDistribute
+        (
+            globalData.map().xfer()
+        )
+    );
+
+    if (debug)
+    {
+        Pout<< "globalMeshData::calcGlobalCoPointSlaves() :"
+            << " finished calculating coupled master to collocated"
+            << " slave point addressing." << endl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from polyMesh
@@ -1639,6 +1674,10 @@ void Foam::globalMeshData::clearOut()
     globalPointBoundaryCellsPtr_.clear();
     globalPointTransformedBoundaryCellsPtr_.clear();
     globalPointBoundaryCellsMapPtr_.clear();
+
+    // Other: collocated points
+    globalCoPointSlavesPtr_.clear();
+    globalCoPointSlavesMapPtr_.clear();
 }
 
 
@@ -1971,7 +2010,10 @@ const Foam::globalIndex& Foam::globalMeshData::globalPointNumbering() const
 {
     if (!globalPointNumberingPtr_.valid())
     {
-        calcGlobalPointSlaves();
+        globalPointNumberingPtr_.reset
+        (
+            new globalIndex(coupledPatch().nPoints())
+        );
     }
     return globalPointNumberingPtr_();
 }
@@ -2161,6 +2203,26 @@ const
 }
 
 
+const Foam::labelListList& Foam::globalMeshData::globalCoPointSlaves() const
+{
+    if (!globalCoPointSlavesPtr_.valid())
+    {
+        calcGlobalCoPointSlaves();
+    }
+    return globalCoPointSlavesPtr_();
+}
+
+
+const Foam::mapDistribute& Foam::globalMeshData::globalCoPointSlavesMap() const
+{
+    if (!globalCoPointSlavesMapPtr_.valid())
+    {
+        calcGlobalCoPointSlaves();
+    }
+    return globalCoPointSlavesMapPtr_();
+}
+
+
 Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
 (
     labelList& pointToGlobal,
@@ -2168,33 +2230,70 @@ Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
 ) const
 {
     const indirectPrimitivePatch& cpp = coupledPatch();
-    const labelListList& pointSlaves = globalPointSlaves();
+    const globalIndex& globalCoupledPoints = globalPointNumbering();
-    const mapDistribute& pointSlavesMap = globalPointSlavesMap();
+    // Use collocated only
+    const labelListList& pointSlaves = globalCoPointSlaves();
+    const mapDistribute& pointSlavesMap = globalCoPointSlavesMap();
+
+
+    // Points are either
+    // - master with slaves
+    // - slave with a master
+    // - other (since e.g. non-collocated cyclics not connected)
+
+    labelList masterGlobalPoint(cpp.nPoints(), -1);
+    forAll(masterGlobalPoint, pointI)
+    {
+        const labelList& slavePoints = pointSlaves[pointI];
+        if (slavePoints.size() > 0)
+        {
+            masterGlobalPoint[pointI] = globalCoupledPoints.toGlobal(pointI);
+        }
+    }
+
+    // Sync by taking max
+    syncData
+    (
+        masterGlobalPoint,
+        pointSlaves,
+        labelListList(cpp.nPoints()),   // no transforms
+        pointSlavesMap,
+        maxEqOp<label>()
+    );
 
 
     // 1. Count number of masters on my processor.
-    label nCoupledMaster = 0;
+    label nMaster = 0;
     PackedBoolList isMaster(mesh_.nPoints(), 1);
     forAll(pointSlaves, pointI)
     {
-        const labelList& slavePoints = pointSlaves[pointI];
+        if (masterGlobalPoint[pointI] == -1)
-
-        if (slavePoints.size() > 0)
         {
-            nCoupledMaster++;
+            // unconnected point (e.g. from separated cyclic)
+            nMaster++;
+        }
+        else if
+        (
+            masterGlobalPoint[pointI]
+         == globalCoupledPoints.toGlobal(pointI)
+        )
+        {
+            // connected master
+            nMaster++;
         }
         else
         {
+            // connected slave point
             isMaster[cpp.meshPoints()[pointI]] = 0;
         }
     }
 
-    label myUniquePoints = mesh_.nPoints() - cpp.nPoints() + nCoupledMaster;
+    label myUniquePoints = mesh_.nPoints() - cpp.nPoints() + nMaster;
 
     //Pout<< "Points :" << nl
     //    << "    mesh             : " << mesh_.nPoints() << nl
     //    << "    of which coupled : " << cpp.nPoints() << nl
-    //    << "    of which master  : " << nCoupledMaster << nl
+    //    << "    of which master  : " << nMaster << nl
     //    << endl;
 
 
@@ -2206,7 +2305,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
     pointToGlobal.setSize(mesh_.nPoints());
     pointToGlobal = -1;
     uniquePoints.setSize(myUniquePoints);
-    label nMaster = 0;
+    nMaster = 0;
 
     forAll(isMaster, meshPointI)
     {
@@ -2245,11 +2344,10 @@ Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
         // On slave copy master index into overal map.
         forAll(pointSlaves, pointI)
         {
-            const labelList& slaves = pointSlaves[pointI];
+            label meshPointI = cpp.meshPoints()[pointI];
 
-            if (slaves.size() == 0)
+            if (!isMaster[meshPointI])
             {
-                label meshPointI = cpp.meshPoints()[pointI];
                 pointToGlobal[meshPointI] = masterToGlobal[pointI];
             }
         }
@@ -2268,8 +2366,8 @@ Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
 ) const
 {
     const indirectPrimitivePatch& cpp = coupledPatch();
-    const labelListList& pointSlaves = globalPointSlaves();
+    const labelListList& pointSlaves = globalCoPointSlaves();
-    const mapDistribute& pointSlavesMap = globalPointSlavesMap();
+    const mapDistribute& pointSlavesMap = globalCoPointSlavesMap();
 
 
     // The patch points come in two variants:
@@ -2280,19 +2378,18 @@ Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
     // coupled points to be the master but this master point is not
     // necessarily on the patch itself! (it might just be connected to the
     // patch point via coupled patches).
-    // So this means that all master point loops should be over the
-    // master-slave structure, not over the patch points and that the unique
-    // point returned is a mesh point.
-    // (unless we want to do the whole master-slave analysis again for the
-    //  current patch only).
 
 
-    // Determine mapping from coupled point to patch point and vice versa
+    // Determine mapping:
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // - from patch point to coupled point (or -1)
+    // - from coupled point to global patch point
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    globalIndex globalPPoints(meshPoints.size());
 
     labelList patchToCoupled(meshPoints.size(), -1);
     label nCoupled = 0;
-    labelList coupledToPatch(pointSlavesMap.constructSize(), -1);
+    labelList coupledToGlobalPatch(pointSlavesMap.constructSize(), -1);
 
     // Note: loop over patch since usually smaller
     forAll(meshPoints, patchPointI)
@@ -2304,7 +2401,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
         if (iter != cpp.meshPointMap().end())
         {
             patchToCoupled[patchPointI] = iter();
-            coupledToPatch[iter()] = patchPointI;
+            coupledToGlobalPatch[iter()] = globalPPoints.toGlobal(patchPointI);
             nCoupled++;
         }
     }
@@ -2314,132 +2411,162 @@ Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
     //    << "    of which on coupled patch:" << nCoupled << endl;
 
 
-    // Pull coupled-to-patch information to master
+    // Determine master of connected points
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
+    // Problem is that the coupled master might not be on the patch. So
-    pointSlavesMap.distribute(coupledToPatch);
+    // work out the best patch-point master from all connected points.
-
+    // - if the coupled master is on the patch it becomes the patch-point master
-
+    // - else the slave with the lowest numbered patch point label
-    // Check on master whether point is anywhere on patch
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // List of master points that are on the patch
-    DynamicList<label> masterPoints(pointSlaves.size());
 
+    // Get all data on master
+    pointSlavesMap.distribute(coupledToGlobalPatch);
     forAll(pointSlaves, coupledPointI)
     {
         const labelList& slaves = pointSlaves[coupledPointI];
 
         if (slaves.size() > 0)
         {
-            // I am master. Is this point on the patch on myself or on any
+            // I am master. What is the best candidate for patch-point master
-            // any slave?
+            label masterI = labelMax;
-            if (coupledToPatch[coupledPointI] != -1)
+            if (coupledToGlobalPatch[coupledPointI] != -1)
             {
-                masterPoints.append(coupledPointI);
+                // I am master and on the coupled patch. Use me.
+                masterI = coupledToGlobalPatch[coupledPointI];
             }
             else
             {
+                // Get min of slaves as master.
                 forAll(slaves, i)
                 {
-                    if (coupledToPatch[slaves[i]] != -1)
+                    label slavePp = coupledToGlobalPatch[slaves[i]];
+                    if (slavePp != -1 && slavePp < masterI)
                     {
-                        masterPoints.append(coupledPointI);
+                        masterI = slavePp;
-                        break;
                     }
                 }
             }
+
+            if (masterI != labelMax)
+            {
+                // Push back
+                coupledToGlobalPatch[coupledPointI] = masterI;
+                forAll(slaves, i)
+                {
+                    coupledToGlobalPatch[slaves[i]] = masterI;
+                }
+            }
         }
     }
+    pointSlavesMap.reverseDistribute(cpp.nPoints(), coupledToGlobalPatch);
 
 
-    // Create global indexing
+    // Generate compact numbering for master points
-    // ~~~~~~~~~~~~~~~~~~~~~~
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 1. patch points that are not on coupled patch:
+    // Now coupledToGlobalPatch is the globalIndex of the master point.
-    //      meshPoints.size()-nCoupled
+    // Now every processor can check whether they hold it and generate a
-    // 2. master points that are on patch:
+    // compact numbering.
-    //      masterPoints.size()
-    label myUniquePoints = meshPoints.size()-nCoupled+masterPoints.size();
-    autoPtr<globalIndex> globalPointsPtr(new globalIndex(myUniquePoints));
 
-    //Pout<< "CoupledPatch:" << nl
+    label nMasters = 0;
-    //    << "    points:" << cpp.nPoints() << nl
+    forAll(meshPoints, patchPointI)
-    //    << "    of which on patch:" << masterPoints.size() << endl;
-
-
-    // Allocate unique points
-    // ~~~~~~~~~~~~~~~~~~~~~~
-
-    pointToGlobal.setSize(meshPoints.size());
-    pointToGlobal = -1;
-    uniqueMeshPoints.setSize(myUniquePoints);
-
-    // Allocate globals for uncoupled patch points
-    label nMaster = 0;
-    forAll(patchToCoupled, patchPointI)
     {
         if (patchToCoupled[patchPointI] == -1)
         {
-            // Allocate global point
+            nMasters++;
-            label globalI = globalPointsPtr().toGlobal(nMaster);
-            pointToGlobal[patchPointI] = globalI;
-            uniqueMeshPoints[nMaster] = meshPoints[patchPointI];
-            nMaster++;
         }
-    }
+        else
-
-    // Allocate globals for master
-    labelList masterToGlobal(pointSlavesMap.constructSize(), -456);
-
-    forAll(masterPoints, i)
-    {
-        label coupledPointI = masterPoints[i];
-
-        // Allocate global point
-        label globalI = globalPointsPtr().toGlobal(nMaster);
-        if (coupledToPatch[coupledPointI] != -1)
         {
-            pointToGlobal[coupledToPatch[coupledPointI]] = globalI;
+            label coupledPointI = patchToCoupled[patchPointI];
-        }
+            if
-        uniqueMeshPoints[nMaster] = cpp.meshPoints()[coupledPointI];
+            (
-        nMaster++;
+                globalPPoints.toGlobal(patchPointI)
-
+             == coupledToGlobalPatch[coupledPointI]
-        // Put global into slave slots
+            )
-        const labelList& slaves = pointSlaves[coupledPointI];
-        masterToGlobal[coupledPointI] = globalI;    // not really necessary
-        forAll(slaves, i)
-        {
-            masterToGlobal[slaves[i]] = globalI;
-        }
-    }
-
-
-    if (nMaster != myUniquePoints)
-    {
-        FatalErrorIn("globalMeshData::mergePoints(..)")
-            << "problem." << abort(FatalError);
-    }
-
-
-    // Send back (from slave slots) to originating processor
-    pointSlavesMap.reverseDistribute(cpp.nPoints(), masterToGlobal);
-
-    // On slaves take over global number
-    forAll(patchToCoupled, patchPointI)
-    {
-        label coupledPointI = patchToCoupled[patchPointI];
-
-        if (coupledPointI != -1)
-        {
-            const labelList& slaves = pointSlaves[coupledPointI];
-
-            if (slaves.size() == 0)
             {
-                pointToGlobal[patchPointI] = masterToGlobal[coupledPointI];
+                // I am the master
+                nMasters++;
             }
         }
     }
 
+    autoPtr<globalIndex> globalPointsPtr(new globalIndex(nMasters));
+
+    //Pout<< "Patch:" << nl
+    //    << "    points:" << meshPoints.size() << nl
+    //    << "    of which on coupled patch:" << nCoupled << nl
+    //    << "    of which master:" << nMasters << endl;
+
+
+
+    // Push back compact numbering for master point onto slaves
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    pointToGlobal.setSize(meshPoints.size());
+    pointToGlobal = -1;
+    uniqueMeshPoints.setSize(nMasters);
+
+    // Sync master in global point numbering so all know the master point.
+    // Initialise globalMaster to be -1 except at a globalMaster.
+    labelList globalMaster(cpp.nPoints(), -1);
+
+    nMasters = 0;
+    forAll(meshPoints, patchPointI)
+    {
+        if (patchToCoupled[patchPointI] == -1)
+        {
+            uniqueMeshPoints[nMasters++] = meshPoints[patchPointI];
+        }
+        else
+        {
+            label coupledPointI = patchToCoupled[patchPointI];
+            if
+            (
+                globalPPoints.toGlobal(patchPointI)
+             == coupledToGlobalPatch[coupledPointI]
+            )
+            {
+                globalMaster[coupledPointI] =
+                    globalPointsPtr().toGlobal(nMasters);
+                uniqueMeshPoints[nMasters++] = meshPoints[patchPointI];
+            }
+        }
+    }
+
+
+    // Sync by taking max
+    syncData
+    (
+        globalMaster,
+        pointSlaves,
+        labelListList(cpp.nPoints()),   // no transforms
+        pointSlavesMap,
+        maxEqOp<label>()
+    );
+
+
+    // Now everyone has the master point in globalPointsPtr numbering. Fill
+    // in the pointToGlobal map.
+    nMasters = 0;
+    forAll(meshPoints, patchPointI)
+    {
+        if (patchToCoupled[patchPointI] == -1)
+        {
+            pointToGlobal[patchPointI] = globalPointsPtr().toGlobal(nMasters++);
+        }
+        else
+        {
+            label coupledPointI = patchToCoupled[patchPointI];
+            pointToGlobal[patchPointI] = globalMaster[coupledPointI];
+
+            if
+            (
+                globalPPoints.toGlobal(patchPointI)
+             == coupledToGlobalPatch[coupledPointI]
+            )
+            {
+                nMasters++;
+            }
+        }
+    }
 
     return globalPointsPtr;
 }

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H

@@ -202,7 +202,7 @@ class globalMeshData
                 globalPointTransformedBoundaryFacesPtr_;
             mutable autoPtr<mapDistribute> globalPointBoundaryFacesMapPtr_;
 
-            // Coupled point to collocated boundary cells
+            // Coupled point to boundary cells
 
             mutable autoPtr<labelList> boundaryCellsPtr_;
             mutable autoPtr<globalIndex> globalBoundaryCellNumberingPtr_;
@@ -212,6 +212,12 @@ class globalMeshData
             mutable autoPtr<mapDistribute> globalPointBoundaryCellsMapPtr_;
 
 
+            // Other: coupled point to coupled COLLOCATED points
+            mutable autoPtr<labelListList> globalCoPointSlavesPtr_;
+            mutable autoPtr<mapDistribute> globalCoPointSlavesMapPtr_;
+
+
+
         // Globally shared point addressing
 
             //- Total number of global points
@@ -303,6 +309,18 @@ class globalMeshData
             //- Calculate global point to global boundary cell addressing.
             void calcGlobalPointBoundaryCells() const;
 
+        // Other
+
+            // Point to collocated points. Note that not all points on
+            // coupled patches now have a master! (since points on either
+            // side of a cyclic are not connected). So check whether the map
+            // reaches all points and decide who is master, slave and who is
+            // its own master. Maybe store as well?
+
+            void calcGlobalCoPointSlaves() const;
+            const labelListList& globalCoPointSlaves() const;
+            const mapDistribute& globalCoPointSlavesMap() const;
+
 
         //- Disallow default bitwise copy construct
         globalMeshData(const globalMeshData&);
@@ -547,7 +565,7 @@ public:
 
             // Other
 
-                //- Helper for merging mesh point data.
+                //- Helper for merging (collocated!) mesh point data.
                 //  Determines:
                 //  - my unique indices
                 //  - global numbering over all unique indices
@@ -559,11 +577,11 @@ public:
                     labelList& uniquePoints
                 ) const;
 
-                //- Helper for merging patch point data.
+                //- Helper for merging (collocated!) patch point data.
                 //  Takes maps from:
                 //  local points to/from mesh. Determines
-                //  - my unique points. These are mesh points, not patch points
+                //  - my unique points. These are mesh point indices, not patch
-                //    since the master might not be on the patch.
+                //    point indices.
                 //  - global numbering over all unique indices.
                 //  - the global number for all local points.
                 autoPtr<globalIndex> mergePoints

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

@@ -533,7 +533,7 @@ void Foam::globalPoints::sendPatchPoints
             // Send to neighbour
             if (debug)
             {
-                Pout<< " Sending to "
+                Pout<< " Sending from " << pp.name() << " to "
                     << procPatch.neighbProcNo() << "   point information:"
                     << patchFaces.size() << endl;
             }
@@ -589,7 +589,8 @@ void Foam::globalPoints::receivePatchPoints
 
             if (debug)
             {
-                Pout<< " Received from "
+                Pout<< " On " << pp.name()
+                    << " Received from "
                     << procPatch.neighbProcNo() << "   point information:"
                     << patchFaces.size() << endl;
             }

src/combustionModels/Make/files

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

src/combustionModels/Make/options

@@ -1,9 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(FOAM_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-    -I$(FOAM_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume

src/combustionModels/combustionModel/combustionModel.C

@@ -1,108 +1,153 @@
-/*---------------------------------------------------------------------------*\
+/*---------------------------------------------------------------------------*\
-  =========                 |
+  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
+   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
+     \\/     M anipulation  |
--------------------------------------------------------------------------------
+-------------------------------------------------------------------------------
-License
+License
-    This file is part of OpenFOAM.
+
-
+    OpenFOAM is free software: you can redistribute it and/or modify it
-    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
-    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
-    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
-    (at your option) any later version.
+
-
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
-    for more details.
+
-
+    You should have received a copy of the GNU General Public License
-    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
-
+\*---------------------------------------------------------------------------*/
-\*---------------------------------------------------------------------------*/
+
-
+#include "combustionModel.H"
-#include "combustionModel.H"
+#include "surfaceFields.H"
-#include "fvm.H"
+#include "fvScalarMatrix.H"
-
+
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
+
-namespace Foam
+namespace Foam
-{
+{
-    defineTypeNameAndDebug(combustionModel, 0);
+    defineTypeNameAndDebug(combustionModel, 0);
-    defineRunTimeSelectionTable(combustionModel, dictionary);
+    defineRunTimeSelectionTable(combustionModel, dictionary);
-}
+};
-
+
-
+
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
+
-Foam::combustionModel::combustionModel
+Foam::combustionModel::combustionModel
-(
+(
-    const dictionary& combustionProperties,
+    const dictionary& combustionProps,
-    const hsCombustionThermo& thermo,
+    hsCombustionThermo& thermo,
-    const compressible::turbulenceModel& turbulence,
+    const compressible::turbulenceModel& turbulence,
-    const surfaceScalarField& phi,
+    const surfaceScalarField& phi,
-    const volScalarField& rho
+    const volScalarField& rho
-)
+)
-:
+:
-    combustionModelCoeffs_
+    coeffs_(dictionary::null),
-    (
+    thermo_(thermo),
-        combustionProperties.subDict
+    turbulence_(turbulence),
-        (
+    mesh_(phi.mesh()),
-            word(combustionProperties.lookup("combustionModel")) + "Coeffs"
+    phi_(phi),
-        )
+    rho_(rho)
-    ),
+{}
-    thermo_(thermo),
+
-    turbulence_(turbulence),
+
-    mesh_(phi.mesh()),
+Foam::combustionModel::combustionModel
-    phi_(phi),
+(
-    rho_(rho),
+    const word& modelType,
-    stoicRatio_(thermo.lookup("stoichiometricAirFuelMassRatio")),
+    const dictionary& combustionProps,
-    s_(thermo.lookup("stoichiometricOxygenFuelMassRatio")),
+    hsCombustionThermo& thermo,
-    qFuel_(thermo_.lookup("qFuel")),
+    const compressible::turbulenceModel& turbulence,
-    composition_(thermo.composition())
+    const surfaceScalarField& phi,
-{}
+    const volScalarField& rho
-
+)
-
+:
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+    coeffs_(combustionProps.subDict(modelType + "Coeffs")),
-
+    thermo_(thermo),
-Foam::combustionModel::~combustionModel()
+    turbulence_(turbulence),
-{}
+    mesh_(phi.mesh()),
-
+    phi_(phi),
-
+    rho_(rho)
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+{}
-
+
-Foam::tmp<Foam::fvScalarMatrix>
+
-Foam::combustionModel::combustionModel::R(volScalarField& fu) const
+// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
-{
+
-    const basicMultiComponentMixture& composition = thermo_.composition();
+Foam::combustionModel::~combustionModel()
-    const volScalarField& ft = composition.Y("ft");
+{}
-    volScalarField fres(composition.fres(ft, stoicRatio_.value()));
+
-    volScalarField wFuelNorm(this->wFuelNorm()*pos(fu - fres));
+
-
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-    return wFuelNorm*fres - fvm::Sp(wFuelNorm, fu);
+
-}
+void Foam::combustionModel::correct()
-
+{
-
+    // do nothing
-Foam::tmp<Foam::volScalarField> Foam::combustionModel::combustionModel::dQ
+}
-(
+
-    const fvScalarMatrix& Rfu
+
-) const
+Foam::tmp<Foam::fvScalarMatrix> Foam::combustionModel::R
-{
+(
-    const basicMultiComponentMixture& composition = thermo_.composition();
+    volScalarField& Y
-    const volScalarField& fu = composition.Y("fu");
+) const
-
+{
-    return (-qFuel_)*(Rfu & fu);
+    return tmp<fvScalarMatrix>
-}
+    (
-
+        new fvScalarMatrix(Y, dimMass/dimTime*Y.dimensions())
-
+    );
-bool Foam::combustionModel::read(const dictionary& combustionProperties)
+}
-{
+
-    combustionModelCoeffs_ = combustionProperties.subDict(type() + "Coeffs");
+
-
+Foam::tmp<Foam::volScalarField> Foam::combustionModel::dQ() const
-    return true;
+{
-}
+    return tmp<Foam::volScalarField>
-
+    (
-
+        new volScalarField
-// ************************************************************************* //
+        (
+            IOobject
+            (
+                "dQ",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0.0)
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::combustionModel::wFuelNorm() const
+{
+    return tmp<Foam::volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "wFuelNorm",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimMass/dimTime/pow3(dimLength), 0.0)
+        )
+    );
+}
+
+
+bool Foam::combustionModel::read(const dictionary& combustionProps)
+{
+    coeffs_ = combustionProps.subDict(type() + "Coeffs");
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/combustionModels/combustionModel/combustionModel.H

@@ -2,11 +2,10 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 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
@@ -25,7 +24,8 @@ Class
     Foam::combustionModel
 
 Description
-    Base class for all non-premixed combustion models.
+    Base class for all non-premixed combustion models based on single step
+    chemistry
 
 SourceFiles
     combustionModel.C
@@ -38,7 +38,6 @@ SourceFiles
 #include "IOdictionary.H"
 #include "hsCombustionThermo.H"
 #include "turbulenceModel.H"
-#include "multivariateSurfaceInterpolationScheme.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -47,7 +46,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                        Class combustionModel Declaration
+                     Class combustionModel Declaration
 \*---------------------------------------------------------------------------*/
 
 class combustionModel
@@ -58,10 +57,10 @@ protected:
     // Protected data
 
         //- Dictionary of coefficients for the particular model
-        dictionary combustionModelCoeffs_;
+        dictionary coeffs_;
 
-        //- Reference to the thermodynamic
+        //- Reference to the thermodynamics
-        const hsCombustionThermo& thermo_;
+        hsCombustionThermo& thermo_;
 
         //- Reference to the turbulence model
         const compressible::turbulenceModel& turbulence_;
@@ -75,15 +74,6 @@ protected:
         //- Reference to the density field
         const volScalarField& rho_;
 
-        //- Stoichiometric air-fuel mass ratio
-        dimensionedScalar stoicRatio_;
-
-        //- Stoichiometric oxygen-fuel mass ratio
-        dimensionedScalar s_;
-
-        //- Heat of combustion (J/Kg)
-        dimensionedScalar qFuel_;
-
 
 private:
 
@@ -95,8 +85,6 @@ private:
         //- Disallow default bitwise assignment
         void operator=(const combustionModel&);
 
-        const basicMultiComponentMixture& composition_;
-
 
 public:
 
@@ -113,7 +101,7 @@ public:
             dictionary,
             (
                 const dictionary& combustionProperties,
-                const hsCombustionThermo& thermo,
+                hsCombustionThermo& thermo,
                 const compressible::turbulenceModel& turbulence,
                 const surfaceScalarField& phi,
                 const volScalarField& rho
@@ -134,7 +122,7 @@ public:
         static autoPtr<combustionModel> New
         (
             const dictionary& combustionProperties,
-            const hsCombustionThermo& thermo,
+            hsCombustionThermo& thermo,
             const compressible::turbulenceModel& turbulence,
             const surfaceScalarField& phi,
             const volScalarField& rho
@@ -143,11 +131,22 @@ public:
 
     // Constructors
 
+        //- Construct null from components
+        combustionModel
+        (
+            const dictionary& combustionProps,
+            hsCombustionThermo& thermo,
+            const compressible::turbulenceModel& turbulence,
+            const surfaceScalarField& phi,
+            const volScalarField& rho
+        );
+
         //- Construct from components
         combustionModel
         (
+            const word& modelType,
             const dictionary& combustionProperties,
-            const hsCombustionThermo& thermo,
+            hsCombustionThermo& thermo,
             const compressible::turbulenceModel& turbulence,
             const surfaceScalarField& phi,
             const volScalarField& rho
@@ -162,39 +161,32 @@ public:
 
         // Access functions
 
-            //- Access composition
-            const basicMultiComponentMixture& composition() const
-            {
-                return composition_;
-            }
-
             //- Access combustion dictionary
-            const dictionary combustionModelCoeffs() const
+            const dictionary coeffs() const
             {
-                return combustionModelCoeffs_;
+                return coeffs_;
             }
 
-            //- Access heat of combustion
-            const dimensionedScalar qFuel() const
-            {
-                return qFuel_;
-            }
 
-        //- Return normalised consumption rate of (fu - fres)
+    // Evolution
-        virtual tmp<volScalarField> wFuelNorm() const = 0;
-
-        //- Fuel consumption rate matrix i.e. source-term for the fuel equation
-        virtual tmp<fvScalarMatrix> R(volScalarField& fu) const;
-
-        //- Heat-release rate calculated from the given
-        //  fuel consumption rate matrix
-        virtual tmp<volScalarField> dQ(const fvScalarMatrix& Rfu) const;
 
         //- Correct combustion rate
-        virtual void correct() = 0;
+        virtual void correct();
+
+        //- Fuel consumption rate matrix, i.e. source term for fuel equation
+        virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
+
+        //- Heat release rate calculated from fuel consumption rate matrix
+        virtual tmp<volScalarField> dQ() const;
+
+        //- Return normalised consumption rate of (fu - fres)
+        virtual tmp<Foam::volScalarField> wFuelNorm() const;
+
+
+    // I-O
 
         //- Update properties from given dictionary
-        virtual bool read(const dictionary& combustionProperties) = 0;
+        virtual bool read(const dictionary& combustionProps);
 };
 
 

src/combustionModels/combustionModel/combustionModelNew.C

@@ -2,11 +2,10 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 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
@@ -30,7 +29,7 @@ License
 Foam::autoPtr<Foam::combustionModel> Foam::combustionModel::New
 (
     const dictionary& combustionProperties,
-    const hsCombustionThermo& thermo,
+    hsCombustionThermo& thermo,
     const compressible::turbulenceModel& turbulence,
     const surfaceScalarField& phi,
     const volScalarField& rho
@@ -54,7 +53,7 @@ Foam::autoPtr<Foam::combustionModel> Foam::combustionModel::New
         )   << "Unknown combustionModel type "
             << combustionModelTypeName << endl << endl
             << "Valid  combustionModels are : " << endl
-            << dictionaryConstructorTablePtr_->sortedToc()
+            << dictionaryConstructorTablePtr_->toc()
             << exit(FatalError);
     }
 

src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.C

@@ -2,11 +2,10 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 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
@@ -25,6 +24,7 @@ License
 
 #include "infinitelyFastChemistry.H"
 #include "addToRunTimeSelectionTable.H"
+#include "fvmSup.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -39,23 +39,40 @@ namespace combustionModels
         infinitelyFastChemistry,
         dictionary
     );
-}
+};
-}
+};
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::combustionModels::infinitelyFastChemistry::infinitelyFastChemistry
 (
-    const dictionary& combustionProperties,
+    const dictionary& combustionProps,
-    const hsCombustionThermo& thermo,
+    hsCombustionThermo& thermo,
     const compressible::turbulenceModel& turbulence,
     const surfaceScalarField& phi,
     const volScalarField& rho
 )
 :
-    combustionModel(combustionProperties, thermo, turbulence, phi, rho),
+    combustionModel(typeName, combustionProps, thermo, turbulence, phi, rho),
-    C_(readScalar(combustionModelCoeffs_.lookup("C")))
+    C_(readScalar(coeffs_.lookup("C"))),
+    singleMixture_
+    (
+        dynamic_cast<singleStepReactingMixture<gasThermoPhysics>&>(thermo)
+    ),
+    wFuelNorm_
+    (
+        IOobject
+        (
+            "wFuelNorm",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero", dimMass/pow3(dimLength)/dimTime, 0.0)
+    )
 {}
 
 
@@ -68,23 +85,66 @@ Foam::combustionModels::infinitelyFastChemistry::~infinitelyFastChemistry()
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void Foam::combustionModels::infinitelyFastChemistry::correct()
-{}
+{
+    singleMixture_.fresCorrect();
+
+    const label fuelI = singleMixture_.fuelIndex();
+
+    const volScalarField& YFuel = thermo_.composition().Y()[fuelI];
+
+    const dimensionedScalar s = singleMixture_.s();
+
+    if (thermo_.composition().contains("O2"))
+    {
+        const volScalarField& YO2 = thermo_.composition().Y("O2");
+        wFuelNorm_ == rho_/(mesh_.time().deltaT()*C_)*min(YFuel, YO2/s.value());
+    }
+}
+
+
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::combustionModels::infinitelyFastChemistry::R(volScalarField& Y) const
+{
+    const label specieI = thermo_.composition().species()[Y.name()];
+
+    const label fNorm = singleMixture_.specieProd()[specieI];
+
+    const volScalarField fres(singleMixture_.fres(specieI));
+
+    const volScalarField wSpecie
+    (
+        wFuelNorm_*singleMixture_.specieStoichCoeffs()[specieI]
+      / max(fNorm*(Y - fres), scalar(0.001))
+    );
+
+    return -fNorm*wSpecie*fres + fNorm*fvm::Sp(wSpecie, Y);
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::infinitelyFastChemistry::dQ() const
+{
+    const label fuelI = singleMixture_.fuelIndex();
+    volScalarField& YFuel = thermo_.composition().Y(fuelI);
+
+    return -singleMixture_.qFuel()*(R(YFuel) & YFuel);
+}
 
 
 Foam::tmp<Foam::volScalarField>
 Foam::combustionModels::infinitelyFastChemistry::wFuelNorm() const
 {
-    return rho_/(mesh_.time().deltaT()*C_);
+    return wFuelNorm_;
 }
 
 
 bool Foam::combustionModels::infinitelyFastChemistry::read
 (
-    const dictionary& combustionProperties
+    const dictionary& combustionProps
 )
 {
-    combustionModel::read(combustionProperties);
+    combustionModel::read(combustionProps);
-    combustionModelCoeffs_.lookup("C") >> C_ ;
+    coeffs_.lookup("C") >> C_ ;
 
     return true;
 }

src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.H

@@ -2,11 +2,10 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 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
@@ -37,8 +36,9 @@ SourceFiles
 #ifndef infinitelyFastChemistry_H
 #define infinitelyFastChemistry_H
 
-#include "fvc.H"
 #include "combustionModel.H"
+#include "singleStepReactingMixture.H"
+#include "thermoPhysicsTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,7 +48,7 @@ namespace combustionModels
 {
 
 /*---------------------------------------------------------------------------*\
-                Class infinitelyFastChemistry Declaration
+                  Class infinitelyFastChemistry Declaration
 \*---------------------------------------------------------------------------*/
 
 class infinitelyFastChemistry
@@ -60,6 +60,13 @@ class infinitelyFastChemistry
         //- Model constant
         scalar C_;
 
+        //- Reference to singleStepReactingMixture mixture
+        singleStepReactingMixture<gasThermoPhysics>& singleMixture_;
+
+        //- Normalised consumption rate of (fu - fres)
+        volScalarField wFuelNorm_;
+
+
     // Private Member Functions
 
         //- Disallow copy construct
@@ -80,29 +87,39 @@ public:
         //- Construct from components
         infinitelyFastChemistry
         (
-            const dictionary& combustionProperties,
+            const dictionary& combustionProps,
-            const hsCombustionThermo& thermo,
+            hsCombustionThermo& thermo,
             const compressible::turbulenceModel& turbulence,
             const surfaceScalarField& phi,
             const volScalarField& rho
         );
 
 
-    // Destructor
+    //- Destructor
-
+    virtual ~infinitelyFastChemistry();
-        virtual ~infinitelyFastChemistry();
 
 
     // Member Functions
 
-        //- Update properties from given dictionary
+        // Evolution
-        virtual bool read(const dictionary& combustionProperties);
 
-        //- Correct combustion rate
+            //- Correct combustion rate
-        virtual void correct();
+            virtual void correct();
 
-        //- Return normalised consumption rate of (fu - fres)
+            //- Fuel consumption rate matrix, i.e. source term for fuel equation
-        virtual tmp<volScalarField> wFuelNorm() const;
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
+
+            //- Heat release rate calculated from fuel consumption rate matrix
+            virtual tmp<volScalarField> dQ() const;
+
+            //- Return normalised consumption rate of (fu - fres)
+            virtual tmp<volScalarField> wFuelNorm() const;
+
+
+    // I-O
+
+            //- Update properties from given dictionary
+            virtual bool read(const dictionary& combustionProperties);
 };
 
 

src/combustionModels/noCombustion/noCombustion.C

@@ -2,11 +2,10 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 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
@@ -39,22 +38,22 @@ namespace combustionModels
         noCombustion,
         dictionary
     );
-}
+};
-}
+};
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::combustionModels::noCombustion::noCombustion
 (
-    const dictionary& combustionProperties,
+    const dictionary& combustionProps,
-    const hsCombustionThermo& thermo,
+    hsCombustionThermo& thermo,
     const compressible::turbulenceModel& turbulence,
     const surfaceScalarField& phi,
     const volScalarField& rho
 )
 :
-    combustionModel(combustionProperties, thermo, turbulence, phi, rho)
+    combustionModel(combustionProps, thermo, turbulence, phi, rho)
 {}
 
 
@@ -64,39 +63,4 @@ Foam::combustionModels::noCombustion::~noCombustion()
 {}
 
 
-void Foam::combustionModels::noCombustion::correct()
-{}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::noCombustion::wFuelNorm() const
-{
-    return tmp<Foam::volScalarField>
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "wFuelNorm",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedScalar("wFuelNorm", dimMass/dimTime/pow3(dimLength), 0.0)
-        )
-    );
-}
-
-
-bool Foam::combustionModels::noCombustion::read
-(
-    const dictionary& combustionProperties
-)
-{
-    return combustionModel::read(combustionProperties);
-}
-
-
 // ************************************************************************* //

src/combustionModels/noCombustion/noCombustion.H

@@ -2,11 +2,10 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 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
@@ -25,7 +24,7 @@ Class
     Foam::combustionModel::noCombustion
 
 Description
-    No combustion
+    Dummy combustion model for 'none' option
 
 SourceFiles
     noCombustion.C
@@ -45,16 +44,13 @@ namespace combustionModels
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class noCombustion Declaration
+                      Class noCombustion Declaration
 \*---------------------------------------------------------------------------*/
 
 class noCombustion
 :
     public combustionModel
 {
-    // Private data
-
-
     // Private Member Functions
 
         //- Disallow copy construct
@@ -67,7 +63,7 @@ class noCombustion
 public:
 
     //- Runtime type information
-    TypeName("noCombustion");
+    TypeName("none");
 
 
     // Constructors
@@ -76,27 +72,15 @@ public:
         noCombustion
         (
             const dictionary& combustionProperties,
-            const hsCombustionThermo& thermo,
+            hsCombustionThermo& thermo,
             const compressible::turbulenceModel& turbulence,
             const surfaceScalarField& phi,
             const volScalarField& rho
         );
 
 
-    // Destructor
+    //- Destructor
-
+    virtual ~noCombustion();
-        virtual ~noCombustion();
-
-    // Member Functions
-
-        //- Update properties from given dictionary
-        virtual bool read(const dictionary& combustionProperties);
-
-        //- Correct combustion rate
-        virtual void correct();
-
-        //- Return normalised consumption rate of (fu - fres)
-        virtual tmp<volScalarField> wFuelNorm() const;
 };
 
 

src/finiteVolume/Make/files

@@ -120,6 +120,7 @@ $(derivedFvPatchFields)/directMappedFixedInternalValue/directMappedFixedInternal
 $(derivedFvPatchFields)/directMappedFixedPushedInternalValue/directMappedFixedPushedInternalValueFvPatchFields.C
 $(derivedFvPatchFields)/directMappedFixedValue/directMappedFixedValueFvPatchFields.C
 $(derivedFvPatchFields)/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.C
+$(derivedFvPatchFields)/directMappedFlowRate/directMappedFlowRateFvPatchVectorField.C
 $(derivedFvPatchFields)/fan/fanFvPatchFields.C
 $(derivedFvPatchFields)/buoyantPressure/buoyantPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C

src/finiteVolume/fields/fvPatchFields/derived/directMappedFlowRate/directMappedFlowRateFvPatchVectorField.C

@@ -0,0 +1,208 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2006-2011 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "directMappedFlowRateFvPatchVectorField.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "directMappedPatchBase.H"
+#include "mapDistribute.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::directMappedFlowRateFvPatchVectorField::
+directMappedFlowRateFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<vector>(p, iF),
+    nbrPhiName_("none"),
+    phiName_("phi"),
+    rhoName_("rho")
+{}
+
+
+Foam::directMappedFlowRateFvPatchVectorField::
+directMappedFlowRateFvPatchVectorField
+(
+    const directMappedFlowRateFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchField<vector>(ptf, p, iF, mapper),
+    nbrPhiName_(ptf.nbrPhiName_),
+    phiName_(ptf.phiName_),
+    rhoName_(ptf.rhoName_)
+{}
+
+
+Foam::directMappedFlowRateFvPatchVectorField::
+directMappedFlowRateFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchField<vector>(p, iF, dict),
+    nbrPhiName_(dict.lookupOrDefault<word>("nbrPhi", "phi")),
+    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
+{}
+
+
+Foam::directMappedFlowRateFvPatchVectorField::
+directMappedFlowRateFvPatchVectorField
+(
+    const directMappedFlowRateFvPatchVectorField& ptf
+)
+:
+    fixedValueFvPatchField<vector>(ptf),
+    nbrPhiName_(ptf.nbrPhiName_),
+    phiName_(ptf.phiName_),
+    rhoName_(ptf.rhoName_)
+{}
+
+
+Foam::directMappedFlowRateFvPatchVectorField::
+directMappedFlowRateFvPatchVectorField
+(
+    const directMappedFlowRateFvPatchVectorField& ptf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<vector>(ptf, iF),
+    nbrPhiName_(ptf.nbrPhiName_),
+    phiName_(ptf.phiName_),
+    rhoName_(ptf.rhoName_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::directMappedFlowRateFvPatchVectorField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    // Get the coupling information from the directMappedPatchBase
+    const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
+    (
+        patch().patch()
+    );
+    const polyMesh& nbrMesh = mpp.sampleMesh();
+    const fvPatch& nbrPatch = refCast<const fvMesh>
+    (
+        nbrMesh
+    ).boundary()[mpp.samplePolyPatch().index()];
+
+    scalarList phi =
+        nbrPatch.lookupPatchField<surfaceScalarField, scalar>(nbrPhiName_);
+
+    mpp.map().distribute(phi);
+
+
+    const surfaceScalarField& phiName =
+        db().lookupObject<surfaceScalarField>(phiName_);
+
+    scalarField U(-phi/patch().magSf());
+
+    vectorField n(patch().nf());
+
+    if (phiName.dimensions() == dimVelocity*dimArea)
+    {
+        // volumetric flow-rate
+        operator==(n*U);
+    }
+    else if (phiName.dimensions() == dimDensity*dimVelocity*dimArea)
+    {
+        const fvPatchField<scalar>& rhop =
+            patch().lookupPatchField<volScalarField, scalar>(rhoName_);
+
+        // mass flow-rate
+        operator==(n*U/rhop);
+
+        if (debug)
+        {
+            scalar phi = gSum(rhop*(*this) & patch().Sf());
+            Info<< patch().boundaryMesh().mesh().name() << ':'
+                << patch().name() << ':'
+                << this->dimensionedInternalField().name() << " <- "
+                << nbrMesh.name() << ':'
+                << nbrPatch.name() << ':'
+                << this->dimensionedInternalField().name() << " :"
+                << " mass flux[Kg/s]:" << -phi
+                << endl;
+        }
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "directMappedFlowRateFvPatchVectorField::updateCoeffs()"
+        )   << "dimensions of " << phiName_ << " are incorrect" << nl
+            << "    on patch " << this->patch().name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << nl << exit(FatalError);
+    }
+
+    fixedValueFvPatchField<vector>::updateCoeffs();
+}
+
+
+void Foam::directMappedFlowRateFvPatchVectorField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchField<vector>::write(os);
+    writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
+    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    os.writeKeyword("nbrPhi") << nbrPhiName_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+   makePatchTypeField
+   (
+       fvPatchVectorField,
+       directMappedFlowRateFvPatchVectorField
+   );
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/directMappedFlowRate/directMappedFlowRateFvPatchVectorField.H

@@ -0,0 +1,176 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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::directMappedFlowRateFvPatchVectorField
+
+Description
+    Describes a volumetric/mass flow normal vector boundary condition by its
+    magnitude as an integral over its area.
+
+    The inlet mass flux is taken from the neighbor region.
+
+    phi is used to determine if the flow is compressible or incompressible.
+
+    The basis of the patch (volumetric or mass) is determined by the
+    dimensions of the flux, phi.
+    The current density is used to correct the velocity when applying the
+    mass basis.
+
+    Example of the boundary condition specification:
+    @verbatim
+    inlet
+    {
+        type            directMappedFlowRate;
+        phi             phi;
+        rho             rho;
+        neigPhi         neigPhiName_;    // Volumetric/mass flow rate
+                                         // [m3/s or kg/s]
+        value           uniform (0 0 0); // placeholder
+    }
+    @endverbatim
+
+SourceFiles
+    directMappedFlowRateFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directMappedFlowRateFvPatchVectorField_H
+#define directMappedFlowRateFvPatchVectorField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+               Class flowRateInletVelocityFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class directMappedFlowRateFvPatchVectorField
+:
+    public fixedValueFvPatchVectorField
+{
+    // Private data
+
+        //- Name of the neighbor flux setting the inlet mass flux
+        word nbrPhiName_;
+
+        //- Name of the local mass flux
+        word phiName_;
+
+        //- Name of the density field used to normalize the mass flux
+        word rhoName_;
+
+
+public:
+
+   //- Runtime type information
+   TypeName("directMappedFlowRate");
+
+
+   // Constructors
+
+        //- Construct from patch and internal field
+        directMappedFlowRateFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        directMappedFlowRateFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  directMappedFlowRateFvPatchVectorField
+        //  onto a new patch
+        directMappedFlowRateFvPatchVectorField
+        (
+            const directMappedFlowRateFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        directMappedFlowRateFvPatchVectorField
+        (
+            const directMappedFlowRateFvPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new directMappedFlowRateFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        directMappedFlowRateFvPatchVectorField
+        (
+            const directMappedFlowRateFvPatchVectorField&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchVectorField> clone
+        (
+            const DimensionedField<vector, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new directMappedFlowRateFvPatchVectorField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/fvm/fvmSup.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -54,6 +54,7 @@ Foam::fvm::Su
     return tfvm;
 }
 
+
 template<class Type>
 Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::Su
@@ -67,6 +68,7 @@ Foam::fvm::Su
     return tfvm;
 }
 
+
 template<class Type>
 Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::Su
@@ -80,6 +82,7 @@ Foam::fvm::Su
     return tfvm;
 }
 
+
 template<class Type>
 Foam::zeroField
 Foam::fvm::Su
@@ -117,6 +120,7 @@ Foam::fvm::Sp
     return tfvm;
 }
 
+
 template<class Type>
 Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::Sp
@@ -130,6 +134,7 @@ Foam::fvm::Sp
     return tfvm;
 }
 
+
 template<class Type>
 Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::Sp
@@ -169,6 +174,7 @@ Foam::fvm::Sp
     return tfvm;
 }
 
+
 template<class Type>
 Foam::zeroField
 Foam::fvm::Sp
@@ -209,6 +215,7 @@ Foam::fvm::SuSp
     return tfvm;
 }
 
+
 template<class Type>
 Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::SuSp
@@ -222,6 +229,7 @@ Foam::fvm::SuSp
     return tfvm;
 }
 
+
 template<class Type>
 Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::SuSp
@@ -235,6 +243,7 @@ Foam::fvm::SuSp
     return tfvm;
 }
 
+
 template<class Type>
 Foam::zeroField
 Foam::fvm::SuSp

src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -219,7 +219,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve()
 template<class Type>
 Foam::tmp<Foam::Field<Type> > Foam::fvMatrix<Type>::residual() const
 {
-    tmp<Field<Type> > tres(source_);
+    tmp<Field<Type> > tres(new Field<Type>(source_));
     Field<Type>& res = tres();
 
     addBoundarySource(res);

src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.C

@@ -44,13 +44,10 @@ void Foam::cyclicFvPatch::makeWeights(scalarField& w) const
 {
     const cyclicFvPatch& nbrPatch = neighbFvPatch();
 
-    const scalarField& magFa = magSf();
-    const scalarField& nbrMagFa = nbrPatch.magSf();
-
     const scalarField deltas(nf() & fvPatch::delta());
     const scalarField nbrDeltas(nbrPatch.nf() & nbrPatch.fvPatch::delta());
 
-    forAll(magFa, facei)
+    forAll(deltas, facei)
     {
         scalar di = deltas[facei];
         scalar dni = nbrDeltas[facei];

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C

@@ -130,7 +130,6 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
             dimensionedScalar("zero", dimEnergy/dimTemperature, 0.0)
         )
     )
-
 {
     if (this->solution().active())
     {

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H

@@ -167,7 +167,7 @@ public:
         // Private data
 
             //- Local copy of carrier specific heat field
-            //  Cp not stored on acrrier thermo, but returned as tmp<...>
+            //  Cp not stored on carrier thermo, but returned as tmp<...>
             const volScalarField Cp_;
 
 

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

@@ -52,7 +52,49 @@ void Foam::FacePostProcessing<CloudType>::applyToFace
                 faceI = j;
                 return;
             }
-        }        
+        }
+    }
+}
+
+
+template<class CloudType>
+void Foam::FacePostProcessing<CloudType>::makeLogFile
+(
+    const word& zoneName,
+    const label zoneI,
+    const label nFaces,
+    const scalar totArea
+)
+{
+    // Create the output file if not already created
+    if (log_)
+    {
+        if (debug)
+        {
+            Info<< "Creating output file." << endl;
+        }
+
+        if (Pstream::master())
+        {
+            const fileName logDir = outputDir_/this->owner().time().timeName();
+
+            // Create directory if does not exist
+            mkDir(logDir);
+
+            // Open new file at start up
+            outputFilePtr_.set
+            (
+                zoneI,
+                new OFstream(logDir/(type() + '_' + zoneName + ".dat"))
+            );
+
+            outputFilePtr_[zoneI]
+                << "# Source    : " << type() << nl
+                << "# Face zone : " << zoneName << nl
+                << "# Faces     : " << nFaces << nl
+                << "# Area      : " << totArea << nl
+                << "# Time" << tab << "mass" << tab << "massFlux" << endl;
+        }
     }
 }
 
@@ -61,8 +103,9 @@ template<class CloudType>
 void Foam::FacePostProcessing<CloudType>::write()
 {
     const fvMesh& mesh = this->owner().mesh();
-    const faceZoneMesh& fzm = this->owner().mesh().faceZones();
+    const Time& time = mesh.time();
-    const scalar dt = this->owner().time().deltaTValue();
+    const faceZoneMesh& fzm = mesh.faceZones();
+    const scalar dt = time.deltaTValue();
 
     totalTime_ += dt;
 
@@ -80,8 +123,11 @@ void Foam::FacePostProcessing<CloudType>::write()
     Info<< "particleFaceFlux output:" << nl;
 
     List<scalarField> zoneMassTotal(mass_.size());
-    forAll(zoneMassTotal, zoneI)
+    List<scalarField> zoneMassFlux(massFlux_.size());
+    forAll(faceZoneIDs_, zoneI)
     {
+        const word& zoneName = fzm[faceZoneIDs_[zoneI]].name();
+
         scalarListList allProcMass(Pstream::nProcs());
         allProcMass[procI] = massTotal_[zoneI];
         Pstream::gatherList(allProcMass);
@@ -90,15 +136,8 @@ void Foam::FacePostProcessing<CloudType>::write()
             (
                 allProcMass, accessOp<scalarList>()
             );
+        const scalar sumMassTotal = sum(zoneMassTotal[zoneI]);
 
-        const word& zoneName = fzm[faceZoneIDs_[zoneI]].name();
-        Info<< "    " << zoneName << " total mass      = "
-            << sum(zoneMassTotal[zoneI]) << nl;
-    }
-
-    List<scalarField> zoneMassFlux(massFlux_.size());
-    forAll(zoneMassFlux, zoneI)
-    {
         scalarListList allProcMassFlux(Pstream::nProcs());
         allProcMassFlux[procI] = massFlux_[zoneI];
         Pstream::gatherList(allProcMassFlux);
@@ -107,10 +146,19 @@ void Foam::FacePostProcessing<CloudType>::write()
             (
                 allProcMassFlux, accessOp<scalarList>()
             );
+        const scalar sumMassFlux = sum(zoneMassFlux[zoneI]);
 
-        const word& zoneName = fzm[faceZoneIDs_[zoneI]].name();
+        Info<< "    " << zoneName
-        Info<< "    " << zoneName << " average mass flux = "
+            << ": total mass = " << sumMassTotal
-            << sum(zoneMassFlux[zoneI]) << nl;
+            << "; average mass flux = " << sumMassFlux
+            << nl;
+
+        if (outputFilePtr_.set(zoneI))
+        {
+            OFstream& os = outputFilePtr_[zoneI];
+            os  << time.timeName() << token::TAB << sumMassTotal << token::TAB
+                <<  sumMassFlux<< endl;
+        }
     }
 
     Info<< endl;
@@ -118,23 +166,6 @@ void Foam::FacePostProcessing<CloudType>::write()
 
     if (surfaceFormat_ != "none")
     {
-        fileName outputDir = mesh.time().path();
-
-        if (Pstream::parRun())
-        {
-            // Put in undecomposed case (Note: gives problems for
-            // distributed data running)
-            outputDir =
-                outputDir/".."/"postProcessing"/cloud::prefix/
-                this->owner().name()/mesh.time().timeName();
-        }
-        else
-        {
-            outputDir =
-                outputDir/"postProcessing"/cloud::prefix/
-                this->owner().name()/mesh.time().timeName();
-        }
-
         forAll(faceZoneIDs_, zoneI)
         {
             const faceZone& fZone = fzm[faceZoneIDs_[zoneI]];
@@ -189,7 +220,7 @@ void Foam::FacePostProcessing<CloudType>::write()
 
                 writer->write
                 (
-                    outputDir,
+                    outputDir_/time.timeName(),
                     fZone.name(),
                     allPoints,
                     allFaces,
@@ -200,7 +231,7 @@ void Foam::FacePostProcessing<CloudType>::write()
 
                 writer->write
                 (
-                    outputDir,
+                    outputDir_/time.timeName(),
                     fZone.name(),
                     allPoints,
                     allFaces,
@@ -247,15 +278,34 @@ Foam::FacePostProcessing<CloudType>::FacePostProcessing
     totalTime_(0.0),
     mass_(),
     massTotal_(),
-    massFlux_()
+    massFlux_(),
+    log_(this->coeffDict().lookup("log")),
+    outputFilePtr_(),
+    outputDir_(owner.mesh().time().path())
 {
     wordList faceZoneNames(this->coeffDict().lookup("faceZones"));
     mass_.setSize(faceZoneNames.size());
     massTotal_.setSize(faceZoneNames.size());
     massFlux_.setSize(faceZoneNames.size());
 
+    outputFilePtr_.setSize(faceZoneNames.size());
+
+    if (Pstream::parRun())
+    {
+        // Put in undecomposed case (Note: gives problems for
+        // distributed data running)
+        outputDir_ =
+            outputDir_/".."/"postProcessing"/cloud::prefix/owner.name();
+    }
+    else
+    {
+        outputDir_ =
+            outputDir_/"postProcessing"/cloud::prefix/owner.name();
+    }
+
     DynamicList<label> zoneIDs;
     const faceZoneMesh& fzm = owner.mesh().faceZones();
+    const surfaceScalarField& magSf = owner.mesh().magSf();
     forAll(faceZoneNames, i)
     {
         const word& zoneName = faceZoneNames[i];
@@ -268,7 +318,16 @@ Foam::FacePostProcessing<CloudType>::FacePostProcessing
             mass_[i].setSize(nFaces, 0.0);
             massTotal_[i].setSize(nFaces, 0.0);
             massFlux_[i].setSize(nFaces, 0.0);
-            Info<< "        " << zoneName << " faces: " << nFaces;
+            Info<< "        " << zoneName << " faces: " << nFaces << nl;
+
+            scalar totArea = 0.0;
+            forAll(fz, j)
+            {
+                totArea += magSf[fz[j]];
+            }
+            totArea = returnReduce(totArea, sumOp<scalar>());
+
+            makeLogFile(zoneName, i, nFaces, totArea);
         }
     }
 
@@ -291,7 +350,10 @@ Foam::FacePostProcessing<CloudType>::FacePostProcessing
     totalTime_(pff.totalTime_),
     mass_(pff.mass_),
     massTotal_(pff.massTotal_),
-    massFlux_(pff.massFlux_)
+    massFlux_(pff.massFlux_),
+    log_(pff.log_),
+    outputFilePtr_(),
+    outputDir_(pff.outputDir_)
 {}
 
 

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

@@ -85,9 +85,27 @@ class FacePostProcessing
         //- Mass flux storage
         List<scalarField> massFlux_;
 
+        //- Flag to indicate whether data should be written to file
+        Switch log_;
+
+        //- Output file pointer per zone
+        PtrList<OFstream> outputFilePtr_;
+
+        //- Output directory
+        fileName outputDir_;
+
 
     // Private Member Functions
 
+        //- Helper function to create log files
+        void makeLogFile
+        (
+            const word& zoneName,
+            const label zoneI,
+            const label nFaces,
+            const scalar totArea
+        );
+
         //- Return index into storage lists if valid zone and face
         void applyToFace
         (

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

@@ -30,6 +30,82 @@ License
 
 using namespace Foam::constant;
 
+// * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::ConeNozzleInjection<CloudType>::setInjectionMethod()
+{
+    word injectionMethodType = this->coeffDict().lookup("injectionMethod");
+    if (injectionMethodType == "disc")
+    {
+        injectionMethod_ = imDisc;
+    }
+    else if (injectionMethodType == "point")
+    {
+        injectionMethod_ = imPoint;
+
+        // Set/cache the injector cell
+        this->findCellAtPosition
+        (
+            injectorCell_,
+            tetFaceI_,
+            tetPtI_,
+            position_,
+            false
+        );
+    }
+    else
+    {
+        FatalErrorIn("Foam::InjectionModel<CloudType>::setInjectionMethod()")
+            << "injectionMethod must be either 'point' or 'disc'"
+            << exit(FatalError);
+    }
+}
+
+
+template<class CloudType>
+void Foam::ConeNozzleInjection<CloudType>::setFlowType()
+{
+    word flowType = this->coeffDict().lookup("flowType");
+    if (flowType == "constantVelocity")
+    {
+        this->coeffDict().lookup("UMag") >> UMag_;
+        flowType_ = ftConstantVelocity;
+    }
+    else if (flowType == "pressureDrivenVelocity")
+    {
+        Pinj_.reset
+        (
+            DataEntry<scalar>::New
+            (
+                "Pinj",
+                this->coeffDict()
+            ).ptr()
+        );
+        flowType_ = ftPressureDrivenVelocity;
+    }
+    else if (flowType == "flowRateAndDischarge")
+    {
+        Cd_.reset
+        (
+            DataEntry<scalar>::New
+            (
+                "Cd",
+                this->coeffDict()
+            ).ptr()
+        );
+        flowType_ = ftFlowRateAndDischarge;
+    }
+    else
+    {
+        FatalErrorIn("Foam::InjectionModel<CloudType>::setFlowType()")
+            << "flowType must be either 'constantVelocity', "
+            <<"'pressureDrivenVelocity' or 'flowRateAndDischarge'"
+            << exit(FatalError);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -41,14 +117,9 @@ Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
 :
     InjectionModel<CloudType>(dict, owner, typeName),
     injectionMethod_(imPoint),
-    outerNozzleDiameter_
+    flowType_(ftConstantVelocity),
-    (
+    outerDiameter_(readScalar(this->coeffDict().lookup("outerDiameter"))),
-        readScalar(this->coeffDict().lookup("outerNozzleDiameter"))
+    innerDiameter_(readScalar(this->coeffDict().lookup("innerDiameter"))),
-    ),
-    innerNozzleDiameter_
-    (
-        readScalar(this->coeffDict().lookup("innerNozzleDiameter"))
-    ),
     duration_(readScalar(this->coeffDict().lookup("duration"))),
     position_(this->coeffDict().lookup("position")),
     injectorCell_(-1),
@@ -67,14 +138,6 @@ Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
             this->coeffDict()
         )
     ),
-    Cd_
-    (
-        DataEntry<scalar>::New
-        (
-            "Cd",
-            this->coeffDict()
-        )
-    ),
     thetaInner_
     (
         DataEntry<scalar>::New
@@ -101,9 +164,13 @@ Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
     ),
     tanVec1_(vector::zero),
     tanVec2_(vector::zero),
-    normal_(vector::zero)
+    normal_(vector::zero),
+
+    UMag_(0.0),
+    Cd_(NULL),
+    Pinj_(NULL)
 {
-    if (innerNozzleDiameter_ >= outerNozzleDiameter_)
+    if (innerDiameter_ >= outerDiameter_)
     {
         FatalErrorIn
         (
@@ -116,38 +183,9 @@ Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
          << exit(FatalError);
     }
 
-    word injectionMethodType = this->coeffDict().lookup("injectionMethod");
+    setInjectionMethod();
 
-    if (injectionMethodType == "disc")
+    setFlowType();
-    {
-        injectionMethod_ = imDisc;
-    }
-    else if (injectionMethodType == "point")
-    {
-        injectionMethod_ = imPoint;
-
-        // Set/cache the injector cell
-        this->findCellAtPosition
-        (
-            injectorCell_,
-            tetFaceI_,
-            tetPtI_,
-            position_,
-            false
-        );
-    }
-    else
-    {
-        FatalErrorIn
-        (
-            "Foam::InjectionModel<CloudType>::InjectionModel"
-            "("
-                "const dictionary&, "
-                "CloudType&"
-            ")"
-        )<< "injectionMethod must be either 'point' or 'disc'" << nl
-         << exit(FatalError);
-    }
 
     cachedRandom& rndGen = this->owner().rndGen();
 
@@ -182,21 +220,23 @@ Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
 :
     InjectionModel<CloudType>(im),
     injectionMethod_(im.injectionMethod_),
-    outerNozzleDiameter_(im.outerNozzleDiameter_),
+    outerDiameter_(im.outerDiameter_),
-    innerNozzleDiameter_(im.innerNozzleDiameter_),
+    innerDiameter_(im.innerDiameter_),
     duration_(im.duration_),
     position_(im.position_),
     injectorCell_(im.injectorCell_),
     direction_(im.direction_),
     parcelsPerSecond_(im.parcelsPerSecond_),
     volumeFlowRate_(im.volumeFlowRate_().clone().ptr()),
-    Cd_(im.Cd_().clone().ptr()),
     thetaInner_(im.thetaInner_().clone().ptr()),
     thetaOuter_(im.thetaOuter_().clone().ptr()),
     sizeDistribution_(im.sizeDistribution_().clone().ptr()),
     tanVec1_(im.tanVec1_),
     tanVec2_(im.tanVec1_),
-    normal_(im.normal_)
+    normal_(im.normal_),
+    UMag_(im.UMag_),
+    Cd_(im.Cd_().clone().ptr()),
+    Pinj_(im.Pinj_().clone().ptr())
 {}
 
 
@@ -283,8 +323,8 @@ void Foam::ConeNozzleInjection<CloudType>::setPositionAndCell
         case imDisc:
         {
             scalar frac = rndGen.sample01<scalar>();
-            scalar dr = outerNozzleDiameter_ - innerNozzleDiameter_;
+            scalar dr = outerDiameter_ - innerDiameter_;
-            scalar r = 0.5*(innerNozzleDiameter_ + frac*dr);
+            scalar r = 0.5*(innerDiameter_ + frac*dr);
             position = position_ + r*normal_;
 
             this->findCellAtPosition
@@ -344,13 +384,38 @@ void Foam::ConeNozzleInjection<CloudType>::setProperties
     dirVec += normal;
     dirVec /= mag(dirVec);
 
-    scalar Ao = 0.25*mathematical::pi*outerNozzleDiameter_*outerNozzleDiameter_;
+    switch (flowType_)
-    scalar Ai = 0.25*mathematical::pi*innerNozzleDiameter_*innerNozzleDiameter_;
+    {
-    scalar massFlowRate =
+        case ftConstantVelocity:
-        this->massTotal()*volumeFlowRate_().value(t)/this->volumeTotal();
+        {
+            parcel.U() = UMag_*dirVec;
+            break;
+        }
+        case ftPressureDrivenVelocity:
+        {
+            scalar pAmbient = this->owner().pAmbient();
+            scalar rho = parcel.rho();
+            scalar UMag = ::sqrt(2.0*(Pinj_().value(t) - pAmbient)/rho);
+            parcel.U() = UMag*dirVec;
+            break;
+        }
+        case ftFlowRateAndDischarge:
+        {
+            scalar Ao = 0.25*mathematical::pi*outerDiameter_*outerDiameter_;
+            scalar Ai = 0.25*mathematical::pi*innerDiameter_*innerDiameter_;
+            scalar massFlowRate =
+                this->massTotal()
+               *volumeFlowRate_().value(t)
+               /this->volumeTotal();
 
-    scalar Umag = massFlowRate/(parcel.rho()*Cd_().value(t)*(Ao - Ai));
+            scalar Umag = massFlowRate/(parcel.rho()*Cd_().value(t)*(Ao - Ai));
-    parcel.U() = Umag*dirVec;
+            parcel.U() = Umag*dirVec;
+            break;
+        }
+        default:
+        {
+        }
+    }
 
     // set particle diameter
     parcel.d() = sizeDistribution_->sample();

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

@@ -32,16 +32,18 @@ Description
       - injector position
       - direction (along injection axis)
       - parcel flow rate
-      - discharge coefficient, Cd
       - inner and outer cone angles
 
     - Parcel diameters obtained by size distribution model
 
     - Parcel velocity is calculated as:
 
-        U = V_dot/(A * Cd), where V_dot is the volume flow rate
+        - Constant velocity
-
+            U = <specified by user>
-    Based on the old 'unitInjection' model
+        - Pressure driven velocity
+            U = sqrt(2*(Pinj - Pamb)/rho)
+        - Flow rate and discharge
+            U = V_dot/(A*Cd)
 
 SourceFiles
     ConeNozzleInjection.C
@@ -83,19 +85,30 @@ public:
         imDisc
     };
 
+    //- Flow type enumeration
+    enum flowType
+    {
+        ftConstantVelocity,
+        ftPressureDrivenVelocity,
+        ftFlowRateAndDischarge
+    };
+
 
 private:
 
     // Private data
 
-        //- point/disc injection method
+        //- Point/disc injection method
         injectionMethod injectionMethod_;
 
+        //- Flow type
+        flowType flowType_;
+
         //- Outer nozzle diameter [m]
-        const scalar outerNozzleDiameter_;
+        const scalar outerDiameter_;
 
         //- Inner nozzle diameter [m]
-        const scalar innerNozzleDiameter_;
+        const scalar innerDiameter_;
 
         //- Injection duration [s]
         const scalar duration_;
@@ -121,9 +134,6 @@ private:
         //- Volume flow rate of parcels to introduce relative to SOI [m^3/s]
         const autoPtr<DataEntry<scalar> > volumeFlowRate_;
 
-        //- Discharge coefficient, relative to SOI [m/s]
-        const autoPtr<DataEntry<scalar> > Cd_;
-
         //- Inner cone angle relative to SOI [deg]
         const autoPtr<DataEntry<scalar> > thetaInner_;
 
@@ -146,6 +156,27 @@ private:
             vector normal_;
 
 
+        // Velocity model coefficients
+
+            //- Constant velocity [m/s]
+            scalar UMag_;
+
+            //- Discharge coefficient, relative to SOI [m/s]
+            autoPtr<DataEntry<scalar> > Cd_;
+
+            //- Injection pressure [Pa]
+            autoPtr<DataEntry<scalar> > Pinj_;
+
+
+    // Private Member Functions
+
+        //- Set the injection type
+        void setInjectionMethod();
+
+        //- Set the injection flow type
+        void setFlowType();
+
+
 public:
 
     //- Runtime type information

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

@@ -205,6 +205,9 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
                         tetPtI
                     );
 
+                // Check/set new parcel thermo properties
+                td.cloud().setParcelThermoProperties(*pPtr, 0.0);
+
                 setParcelProperties(*pPtr, j);
 
                 // Check new parcel properties

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

@@ -345,7 +345,7 @@ void Foam::SprayCloud<CloudType>::motion(TrackData& td)
             parcelType& p = iter();
             if (p.mass() < VSMALL)
             {
-                deleteParticle(p);
+                this->deleteParticle(p);
             }
         }
     }

src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloudI.H

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

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelI.H

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

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelIO.C

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

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

@@ -648,7 +648,7 @@ Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
         -GREAT              // null value
     );
 
-    return snapParams.snapTol()*maxEdgeLen;
+    return scalarField(snapParams.snapTol()*maxEdgeLen);
 }
 
 

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

@@ -773,7 +773,7 @@ void Foam::autoSnapDriver::determineAllFeatures
     labelList nearPointFeat;
     labelList nearPointIndex;
     {
-        scalarField snapDistSqr = sqr(snapDist);
+        scalarField snapDistSqr(sqr(snapDist));
         features.findNearestEdge
         (
             pp.localPoints(),

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -45,6 +45,7 @@ License
 #include "slipPointPatchFields.H"
 #include "fixedValuePointPatchFields.H"
 #include "calculatedPointPatchFields.H"
+#include "cyclicSlipPointPatchFields.H"
 #include "processorPointPatch.H"
 #include "globalIndex.H"
 #include "meshTools.H"
@@ -1458,6 +1459,10 @@ Foam::tmp<Foam::pointVectorField> Foam::meshRefinement::makeDisplacementField
         {
             patchFieldTypes[patchI] = calculatedPointPatchVectorField::typeName;
         }
+        else if (isA<cyclicPointPatch>(pointPatches[patchI]))
+        {
+            patchFieldTypes[patchI] = cyclicSlipPointPatchVectorField::typeName;
+        }
     }
 
     // Note: time().timeName() instead of meshRefinement::timeName() since

src/meshTools/Make/files

@@ -159,7 +159,10 @@ twoDPointCorrector/twoDPointCorrector.C
 directMapped/directMappedPolyPatch/directMappedPatchBase.C
 directMapped/directMappedPolyPatch/directMappedPolyPatch.C
 directMapped/directMappedPolyPatch/directMappedWallPolyPatch.C
+directMapped/directMappedPolyPatch/directMappedVariableThicknessWallPolyPatch.C
 directMapped/directMappedPointPatch/directMappedPointPatch.C
 directMapped/directMappedPointPatch/directMappedWallPointPatch.C
 
+
+
 LIB = $(FOAM_LIBBIN)/libmeshTools

src/meshTools/directMapped/directMappedPolyPatch/directMappedVariableThicknessWallPolyPatch.C

@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "directMappedVariableThicknessWallPolyPatch.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(directMappedVariableThicknessWallPolyPatch, 0);
+
+    addToRunTimeSelectionTable
+    (
+        polyPatch,
+        directMappedVariableThicknessWallPolyPatch,
+        word
+    );
+
+    addToRunTimeSelectionTable
+    (
+        polyPatch,
+        directMappedVariableThicknessWallPolyPatch,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //
+
+Foam::directMappedVariableThicknessWallPolyPatch::
+directMappedVariableThicknessWallPolyPatch
+(
+    const word& name,
+    const label size,
+    const label start,
+    const label index,
+    const polyBoundaryMesh& bm
+)
+:
+    directMappedWallPolyPatch(name, size, start, index, bm),
+    thickness_(size)
+{}
+
+
+Foam::directMappedVariableThicknessWallPolyPatch::
+directMappedVariableThicknessWallPolyPatch
+(
+    const word& name,
+    const label size,
+    const label start,
+    const label index,
+    const word& sampleRegion,
+    const directMappedPatchBase::sampleMode mode,
+    const word& samplePatch,
+    const vectorField& offset,
+    const polyBoundaryMesh& bm
+)
+:
+    directMappedWallPolyPatch(name, size, start, index, bm),
+    thickness_(size)
+{}
+
+
+Foam::directMappedVariableThicknessWallPolyPatch::
+directMappedVariableThicknessWallPolyPatch
+(
+    const word& name,
+    const label size,
+    const label start,
+    const label index,
+    const word& sampleRegion,
+    const directMappedPatchBase::sampleMode mode,
+    const word& samplePatch,
+    const vector& offset,
+    const polyBoundaryMesh& bm
+)
+:
+    directMappedWallPolyPatch(name, size, start, index, bm),
+    thickness_(size)
+{}
+
+
+Foam::directMappedVariableThicknessWallPolyPatch::
+directMappedVariableThicknessWallPolyPatch
+(
+    const word& name,
+    const dictionary& dict,
+    const label index,
+    const polyBoundaryMesh& bm
+)
+:
+    directMappedWallPolyPatch(name, dict, index, bm),
+    thickness_(scalarField("thickness", dict, this->size()))
+{}
+
+
+Foam::directMappedVariableThicknessWallPolyPatch::
+directMappedVariableThicknessWallPolyPatch
+(
+    const directMappedVariableThicknessWallPolyPatch& pp,
+    const polyBoundaryMesh& bm
+)
+:
+    directMappedWallPolyPatch(pp, bm),
+    thickness_(pp.thickness_)
+{}
+
+
+Foam::directMappedVariableThicknessWallPolyPatch::
+directMappedVariableThicknessWallPolyPatch
+(
+    const directMappedVariableThicknessWallPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const label newSize,
+    const label newStart
+)
+:
+    directMappedWallPolyPatch(pp, bm, index, newSize, newStart),
+    thickness_(newSize)
+{}
+
+
+Foam::directMappedVariableThicknessWallPolyPatch::
+directMappedVariableThicknessWallPolyPatch
+(
+    const directMappedVariableThicknessWallPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const labelUList& mapAddressing,
+    const label newStart
+)
+:
+    directMappedWallPolyPatch(pp, bm, index, mapAddressing, newStart),
+    thickness_(pp.size())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::directMappedVariableThicknessWallPolyPatch::
+~directMappedVariableThicknessWallPolyPatch()
+{
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::directMappedVariableThicknessWallPolyPatch::
+write(Foam::Ostream& os) const
+{
+    os.writeKeyword("thickness") << thickness_ << token::END_STATEMENT << nl;
+}
+
+// ************************************************************************* //

src/meshTools/directMapped/directMappedPolyPatch/directMappedVariableThicknessWallPolyPatch.H

@@ -0,0 +1,236 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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::directMappedVariableThicknessWallPolyPatch
+
+Description
+    Foam::directMappedVariableThicknessWallPolyPatch
+
+SourceFiles
+    directMappedVariableThicknessWallPolyPatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directMappedVariableThicknessWallPolyPatch_H
+#define directMappedVariableThicknessWallPolyPatch_H
+
+#include "wallPolyPatch.H"
+#include "directMappedWallPolyPatch.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+        Class directMappedVariableThicknessWallPolyPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class directMappedVariableThicknessWallPolyPatch
+:
+    public directMappedWallPolyPatch
+{
+
+    // Private data
+
+        //- Thickness
+        scalarList thickness_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("directMappedWallVariableThickness");
+
+
+    // Constructors
+
+        //- Construct from components
+        directMappedVariableThicknessWallPolyPatch
+        (
+            const word& name,
+            const label size,
+            const label start,
+            const label index,
+            const polyBoundaryMesh& bm
+        );
+
+        //- Construct from components
+        directMappedVariableThicknessWallPolyPatch
+        (
+            const word& name,
+            const label size,
+            const label start,
+            const label index,
+            const word& sampleRegion,
+            const directMappedPatchBase::sampleMode mode,
+            const word& samplePatch,
+            const vectorField& offset,
+            const polyBoundaryMesh& bm
+        );
+
+        //- Construct from components. Uniform offset.
+        directMappedVariableThicknessWallPolyPatch
+        (
+            const word& name,
+            const label size,
+            const label start,
+            const label index,
+            const word& sampleRegion,
+            const directMappedPatchBase::sampleMode mode,
+            const word& samplePatch,
+            const vector& offset,
+            const polyBoundaryMesh& bm
+        );
+
+        //- Construct from dictionary
+        directMappedVariableThicknessWallPolyPatch
+        (
+            const word& name,
+            const dictionary& dict,
+            const label index,
+            const polyBoundaryMesh& bm
+        );
+
+        //- Construct as copy, resetting the boundary mesh
+        directMappedVariableThicknessWallPolyPatch
+        (
+            const directMappedVariableThicknessWallPolyPatch&,
+            const polyBoundaryMesh&
+        );
+
+        //- Construct given the original patch and resetting the
+        //  face list and boundary mesh information
+        directMappedVariableThicknessWallPolyPatch
+        (
+            const directMappedVariableThicknessWallPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const label newSize,
+            const label newStart
+        );
+
+        //- Construct given the original patch and a map
+        directMappedVariableThicknessWallPolyPatch
+        (
+            const directMappedVariableThicknessWallPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const label newStart
+        );
+
+        //- Construct and return a clone, resetting the boundary mesh
+        virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new directMappedVariableThicknessWallPolyPatch(*this, bm)
+            );
+        }
+
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const label newSize,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new directMappedVariableThicknessWallPolyPatch
+                (
+                    *this,
+                    bm,
+                    index,
+                    newSize,
+                    newStart
+                )
+            );
+        }
+
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new directMappedVariableThicknessWallPolyPatch
+                (
+                    *this,
+                    bm,
+                    index,
+                    mapAddressing,
+                    newStart
+                )
+            );
+        }
+
+
+    //- Destructor
+    virtual ~directMappedVariableThicknessWallPolyPatch();
+
+
+    // Member functions
+
+        //- Return non const thickness
+        scalarList& thickness()
+        {
+            return thickness_;
+        }
+
+
+        //- Return const thickness
+        const scalarList& thickness() const
+        {
+            return thickness_;
+        }
+
+
+        //- Write the polyPatch data as a dictionary
+        void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/indexedOctree/treeDataTriSurface.C

@@ -179,9 +179,14 @@ defineTypeNameAndDebug(Foam::treeDataTriSurface, 0);
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
-Foam::treeDataTriSurface::treeDataTriSurface(const triSurface& surface)
+Foam::treeDataTriSurface::treeDataTriSurface
+(
+    const triSurface& surface,
+    const scalar planarTol
+)
 :
-    surface_(surface)
+    surface_(surface),
+    planarTol_(planarTol)
 {}
 
 
@@ -437,7 +442,7 @@ bool Foam::treeDataTriSurface::intersects
         dir,
         points,
         intersection::HALF_RAY,
-        indexedOctree<treeDataTriSurface>::perturbTol()
+        planarTol_
     );
 
     if (inter.hit() && inter.distance() <= 1)

src/meshTools/indexedOctree/treeDataTriSurface.H

@@ -55,8 +55,11 @@ class treeDataTriSurface
 {
     // Private data
 
+        //- Reference to triSurface
         const triSurface& surface_;
 
+        //- Tolerance to use for intersection tests
+        const scalar planarTol_;
 
     // Private Member Functions
 
@@ -83,8 +86,9 @@ public:
 
     // Constructors
 
-        //- Construct from triSurface. Holds reference.
+        //- Construct from triSurface and tolerance for intersection
-        treeDataTriSurface(const triSurface&);
+        //  tests. Holds reference.
+        treeDataTriSurface(const triSurface&, const scalar planarTol);
 
 
     // Member Functions

src/meshTools/searchableSurface/triSurfaceMesh.C

@@ -507,7 +507,7 @@ void Foam::triSurfaceMesh::movePoints(const pointField& newPoints)
 
 
 const Foam::indexedOctree<Foam::treeDataTriSurface>&
-    Foam::triSurfaceMesh::tree() const
+Foam::triSurfaceMesh::tree() const
 {
     if (tree_.empty())
     {
@@ -543,7 +543,7 @@ const Foam::indexedOctree<Foam::treeDataTriSurface>&
         (
             new indexedOctree<treeDataTriSurface>
             (
-                treeDataTriSurface(*this),
+                treeDataTriSurface(*this, tolerance_),
                 bb,
                 maxTreeDepth_,  // maxLevel
                 10,             // leafsize
@@ -559,7 +559,7 @@ const Foam::indexedOctree<Foam::treeDataTriSurface>&
 
 
 const Foam::indexedOctree<Foam::treeDataEdge>&
- Foam::triSurfaceMesh::edgeTree() const
+Foam::triSurfaceMesh::edgeTree() const
 {
     if (edgeTree_.empty())
     {

src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -64,7 +64,11 @@ Foam::triSurfaceSearch::triSurfaceSearch(const triSurface& surface)
     (
         new indexedOctree<treeDataTriSurface>
         (
-            treeDataTriSurface(surface_),
+            treeDataTriSurface
+            (
+                surface_,
+                indexedOctree<treeDataTriSurface>::perturbTol()
+            ),
             treeBb,
             8,      // maxLevel
             10,     // leafsize

src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C

@@ -448,8 +448,8 @@ void Foam::decompositionMethod::calcCellCells
     // Count number of faces (internal + coupled)
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Number of faces per cell
+    // Number of faces per coarse cell
-    labelList nFacesPerCell(mesh.nCells(), 0);
+    labelList nFacesPerCell(nCoarse, 0);
 
     for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
     {
@@ -481,7 +481,11 @@ void Foam::decompositionMethod::calcCellCells
             {
                 label own = agglom[faceOwner[faceI]];
                 label globalNei = globalNeighbour[bFaceI];
-                if (cellPair.insert(labelPair(own, globalNei)))
+                if
+                (
+                    globalAgglom.toGlobal(own) != globalNei
+                 && cellPair.insert(labelPair(own, globalNei))
+                )
                 {
                     nFacesPerCell[own]++;
                 }

src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.H

@@ -111,7 +111,7 @@ public:
         }
 
         //- Return for every coordinate the wanted processor number. Use the
-        //  mesh connectivity (if needed)
+        //  mesh connectivity (if needed). See note on weights in scotchDecomp.H
         virtual labelList decompose
         (
             const polyMesh& mesh,
@@ -122,7 +122,7 @@ public:
         //- Return for every coordinate the wanted processor number. Gets
         //  passed agglomeration map (from fine to coarse cells) and coarse cell
         //  location. Can be overridden by decomposers that provide this
-        //  functionality natively.
+        //  functionality natively. See note on weights in scotchDecomp.H
         virtual labelList decompose
         (
             const polyMesh& mesh,
@@ -138,6 +138,7 @@ public:
         //    from 0 at processor0 and then incrementing all through the
         //    processors)
         //  - the connections are across coupled patches
+        //  See note on weights in scotchDecomp.H
         virtual labelList decompose
         (
             const labelListList& globalCellCells,

src/parallel/decompose/scotchDecomp/scotchDecomp.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,9 +107,9 @@ public:
 
         //- Return for every coordinate the wanted processor number. Use the
         //  mesh connectivity (if needed)
-        //  Weights get truncated to convert into integer
+        //  Weights get normalised with minimum weight and truncated to
-        //  so e.g. 3.5 is seen as 3. The overall sum of weights
+        //  convert into integer so e.g. 3.5 is seen as 3. The overall sum
-        //  might otherwise overflow.
+        //  of weights might otherwise overflow.
         virtual labelList decompose
         (
             const polyMesh& mesh,

src/regionModels/Allwmake

@@ -4,9 +4,10 @@ makeType=${1:-libso}
 set -x
 
 wmake $makeType regionModel
-#wmake $makeType pyrolysisModels
+wmake $makeType pyrolysisModels
 wmake $makeType surfaceFilmModels
-#wmake $makeType regionCoupling
+wmake $makeType thermoBaffleModels
+wmake $makeType regionCoupling
 
 
 # ----------------------------------------------------------------- end-of-file

src/regionModels/pyrolysisModels/Make/files

@@ -0,0 +1,7 @@
+/* Pyrolysis models */
+pyrolysisModel/pyrolysisModel.C
+pyrolysisModel/pyrolysisModelNew.C
+reactingOneDim/reactingOneDim.C
+noPyrolysis/noPyrolysis.C
+
+LIB = $(FOAM_LIBBIN)/libpyrolysisModels

src/regionModels/pyrolysisModels/Make/options

@@ -0,0 +1,25 @@
+
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solid/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidChemistryModel/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/LES/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude
+
+
+EXE_LIBS = \
+    -lregionModels \
+    -lsolidChemistryModel \
+    -lsolidThermo \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lcompressibleLESModels

src/regionModels/pyrolysisModels/noPyrolysis/noPyrolysis.C

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "noPyrolysis.H"
+#include "addToRunTimeSelectionTable.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace pyrolysisModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(noPyrolysis, 0);
+addToRunTimeSelectionTable(pyrolysisModel, noPyrolysis, mesh);
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+bool noPyrolysis::read()
+{
+    if (pyrolysisModel::read())
+    {
+        // no additional info to read
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+noPyrolysis::noPyrolysis(const word& modelType, const fvMesh& mesh)
+:
+    pyrolysisModel(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+noPyrolysis::~noPyrolysis()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+const volScalarField& noPyrolysis::rho() const
+{
+    FatalErrorIn("const volScalarField& noPyrolysis::rho() const")
+        << "rho field not available for " << type() << abort(FatalError);
+    return volScalarField::null();
+}
+
+
+const volScalarField& noPyrolysis::T() const
+{
+    FatalErrorIn("const volScalarField& noPyrolysis::T() const")
+        << "T field not available for " << type() << abort(FatalError);
+    return volScalarField::null();
+}
+
+
+const tmp<volScalarField> noPyrolysis::Cp() const
+{
+    FatalErrorIn("const tmp<volScalarField>& noPyrolysis::Cp() const")
+        << "Cp field not available for " << type() << abort(FatalError);
+
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "noPyrolysis::Cp",
+                time().timeName(),
+                primaryMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            primaryMesh(),
+            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0.0)
+        )
+    );
+}
+
+
+const volScalarField& noPyrolysis::kappa() const
+{
+    FatalErrorIn("const volScalarField& noPyrolysis::kappa() const")
+        << "kappa field not available for " << type() << abort(FatalError);
+    return volScalarField::null();
+}
+
+
+const volScalarField& noPyrolysis::K() const
+{
+    FatalErrorIn("const volScalarField& noPyrolysis::K() const")
+        << "K field not available for " << type() << abort(FatalError);
+    return volScalarField::null();
+}
+
+
+const surfaceScalarField& noPyrolysis::phiGas() const
+{
+    FatalErrorIn("const volScalarField& noPyrolysis::phiGas() const")
+        << "phiGas field not available for " << type() << abort(FatalError);
+    return surfaceScalarField::null();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/pyrolysisModels/noPyrolysis/noPyrolysis.H

@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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::noPyrolysis
+
+Description
+    Dummy surface pyrolysis model for 'none'
+
+SourceFiles
+    noPyrolysis.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noPyrolysis_H
+#define noPyrolysis_H
+
+#include "pyrolysisModel.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace pyrolysisModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class noPyrolysis Declaration
+\*---------------------------------------------------------------------------*/
+
+class noPyrolysis
+:
+    public pyrolysisModel
+{
+private:
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        noPyrolysis(const noPyrolysis&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const noPyrolysis&);
+
+
+protected:
+
+    // Protected member functions
+
+        //- Read control parameters from dictionary
+        virtual bool read();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+
+    // Constructors
+
+        //- Construct from type name and mesh
+        noPyrolysis(const word& modelType, const fvMesh& mesh);
+
+
+    //- Destructor
+    virtual ~noPyrolysis();
+
+
+    // Member Functions
+
+        // Fields
+
+            //- Return density [kg/m3]
+            virtual const volScalarField& rho() const;
+
+            //- Return const temperature [K]
+            virtual const volScalarField& T() const;
+
+            //- Return specific heat capacity [J/kg/K]
+            virtual const tmp<volScalarField> Cp() const;
+
+            //- Return the region absorptivity [1/m]
+            virtual const volScalarField& kappa() const;
+
+            //- Return the region thermal conductivity [W/m/k]
+            virtual const volScalarField& K() const;
+
+            //- Return the total gas mass flux to primary region [kg/m2/s]
+            virtual const surfaceScalarField& phiGas() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace pyrolysisModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModel.C

@@ -0,0 +1,182 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pyrolysisModel.H"
+#include "fvMesh.H"
+#include "directMappedFieldFvPatchField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace pyrolysisModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(pyrolysisModel, 0);
+defineRunTimeSelectionTable(pyrolysisModel, mesh);
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+void pyrolysisModel::constructMeshObjects()
+{
+    // construct filmDelta field if coupled to film model
+    if (filmCoupled_)
+    {
+        filmDeltaPtr_.reset
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "filmDelta",
+                    time_.timeName(),
+                    regionMesh(),
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                regionMesh()
+            )
+        );
+
+        const volScalarField& filmDelta = filmDeltaPtr_();
+
+        bool foundCoupledPatch = false;
+        forAll(filmDelta.boundaryField(), patchI)
+        {
+            const fvPatchField<scalar>& fvp = filmDelta.boundaryField()[patchI];
+            if (isA<directMappedFieldFvPatchField<scalar> >(fvp))
+            {
+                foundCoupledPatch = true;
+                break;
+            }
+        }
+
+        if (!foundCoupledPatch)
+        {
+            WarningIn("void pyrolysisModels::constructMeshObjects()")
+                << "filmCoupled flag set to true, but no "
+                << directMappedFieldFvPatchField<scalar>::typeName
+                << " patches found on " << filmDelta.name() << " field"
+                << endl;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+bool pyrolysisModel::read()
+{
+    if (regionModel1D::read())
+    {
+        filmCoupled_ = readBool(coeffs_.lookup("filmCoupled"));
+        reactionDeltaMin_ =
+            coeffs_.lookupOrDefault<scalar>("reactionDeltaMin", 0.0);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+pyrolysisModel::pyrolysisModel(const fvMesh& mesh)
+:
+    regionModel1D(mesh),
+    filmCoupled_(false),
+    filmDeltaPtr_(NULL),
+    reactionDeltaMin_(0.0)
+{}
+
+
+pyrolysisModel::pyrolysisModel(const word& modelType, const fvMesh& mesh)
+:
+    regionModel1D(mesh, "pyrolysis", modelType),
+    filmCoupled_(false),
+    filmDeltaPtr_(NULL),
+    reactionDeltaMin_(0.0)
+{
+    if (active_)
+    {
+        read();
+        constructMeshObjects();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+pyrolysisModel::~pyrolysisModel()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+scalar pyrolysisModel::addMassSources
+(
+    const label patchI,
+    const label faceI
+)
+{
+    return 0.0;
+}
+
+
+void pyrolysisModel::preEvolveRegion()
+{
+    if (filmCoupled_)
+    {
+        filmDeltaPtr_->correctBoundaryConditions();
+    }
+}
+
+
+scalar pyrolysisModel::solidRegionDiffNo() const
+{
+    return VSMALL;
+}
+
+
+scalar pyrolysisModel::maxDiff() const
+{
+    return GREAT;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModel.H

@@ -0,0 +1,199 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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::pyrolysisModel
+
+Description
+
+SourceFiles
+    pyrolysisModelI.H
+    pyrolysisModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pyrolysisModel_H
+#define pyrolysisModel_H
+
+#include "runTimeSelectionTables.H"
+#include "volFieldsFwd.H"
+#include "solidChemistryModel.H"
+#include "basicSolidThermo.H"
+#include "regionModel1D.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class fvMesh;
+class Time;
+
+namespace regionModels
+{
+namespace pyrolysisModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class pyrolysisModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class pyrolysisModel
+:
+    public regionModel1D
+{
+private:
+
+    // Private Member Functions
+
+        //- Construct fields
+        void constructMeshObjects();
+
+        //- Disallow default bitwise copy construct
+        pyrolysisModel(const pyrolysisModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const pyrolysisModel&);
+
+
+protected:
+
+    // Protected Data
+
+        //- Flag to indicate whether pyrolysis region coupled to a film region
+        bool filmCoupled_;
+
+        //- Pointer to film thickness field
+        autoPtr<volScalarField> filmDeltaPtr_;
+
+        //- Film height below which reactions can occur [m]
+        scalar reactionDeltaMin_;
+
+
+    // Protected Member Functions
+
+        //- Read control parameters from dictionary
+        virtual bool read();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("pyrolysisModel");
+
+
+    // Declare runtime constructor selection table
+
+         declareRunTimeSelectionTable
+         (
+             autoPtr,
+             pyrolysisModel,
+             mesh,
+             (
+                const word& modelType,
+                const fvMesh& mesh
+             ),
+             (modelType, mesh)
+         );
+
+    // Constructors
+
+        //- Construct null from mesh
+        pyrolysisModel(const fvMesh& mesh);
+
+        //- Construct from type name and mesh
+        pyrolysisModel(const word& modelType, const fvMesh& mesh);
+
+
+    // Selectors
+
+        //- Return a reference to the selected surface film model
+        static autoPtr<pyrolysisModel> New(const fvMesh& mesh);
+
+
+    //- Destructor
+    virtual ~pyrolysisModel();
+
+
+    // Member Functions
+
+        // Access
+
+            // Fields
+
+                //- Return density [kg/m3]
+                virtual const volScalarField& rho() const = 0;
+
+                //- Return const temperature [K]
+                virtual const volScalarField& T() const = 0;
+
+                //- Return specific heat capacity [J/kg/K]
+                virtual const tmp<volScalarField> Cp() const = 0;
+
+                //- Return the region absorptivity [1/m]
+                virtual const volScalarField& kappa() const = 0;
+
+                //- Return the region thermal conductivity [W/m/k]
+                virtual const volScalarField& K() const = 0;
+
+                //- Return the total gas mass flux to primary region [kg/m2/s]
+                virtual const surfaceScalarField& phiGas() const = 0;
+
+            // Sources
+
+                //- External hook to add mass to the primary region
+                virtual scalar addMassSources
+                (
+                    const label patchI,
+                    const label faceI
+                );
+
+
+        // Evolution
+
+            //- Pre-evolve region
+            virtual void preEvolveRegion();
+
+
+        // Helper function
+
+            //- Mean diffusion number of the solid region
+            virtual scalar solidRegionDiffNo() const;
+
+            //- Return max diffusivity allowed in the solid
+            virtual scalar maxDiff() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace pyrolysisModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModelI.H

@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pyrolysisModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline const Foam::Switch&
+Foam::pyrolysisModels::pyrolysisModel::active() const
+{
+    return active_;
+}
+
+
+inline const Foam::dictionary&
+Foam::pyrolysisModels::pyrolysisModel::coeffs() const
+{
+    return coeffs_;
+}
+
+
+// ************************************************************************* //

src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModelNew.C

@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pyrolysisModel.H"
+#include "fvMesh.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace pyrolysisModels
+{
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+autoPtr<pyrolysisModel> pyrolysisModel::New(const fvMesh& mesh)
+{
+    // get model name, but do not register the dictionary
+    const word modelType
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "pyrolysisProperties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).lookup("pyrolysisModel")
+    );
+
+    Info<< "Selecting pyrolysisModel " << modelType << endl;
+
+    meshConstructorTable::iterator cstrIter =
+        meshConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == meshConstructorTablePtr_->end())
+    {
+        FatalErrorIn("pyrolysisModel::New(const fvMesh&)")
+            << "Unknown pyrolysisModel type " << modelType
+            << nl << nl << "Valid pyrolisisModel types are:" << nl
+            << meshConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<pyrolysisModel>(cstrIter()(modelType, mesh));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C

@@ -0,0 +1,615 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2011 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "reactingOneDim.H"
+#include "addToRunTimeSelectionTable.H"
+#include "directMappedPatchBase.H"
+#include "mapDistribute.H"
+#include "zeroGradientFvPatchFields.H"
+#include "surfaceInterpolate.H"
+#include "fvm.H"
+#include "fvcDiv.H"
+#include "fvcVolumeIntegrate.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace pyrolysisModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(reactingOneDim, 0);
+
+addToRunTimeSelectionTable(pyrolysisModel, reactingOneDim, mesh);
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+bool reactingOneDim::read()
+{
+    if (pyrolysisModel::read())
+    {
+        const dictionary& solution = this->solution().subDict("SIMPLE");
+        solution.lookup("nNonOrthCorr") >> nNonOrthCorr_;
+        time_.controlDict().lookup("maxDi") >> maxDiff_;
+
+        coeffs().lookup("radFluxName") >> primaryRadFluxName_;
+        coeffs().lookup("minimumDelta") >> minimumDelta_;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void reactingOneDim::updateQr()
+{
+    // Retrieve field from coupled region using mapped boundary conditions
+    QrCoupled_.correctBoundaryConditions();
+
+    // Update local Qr from coupled Qr field
+    Qr_ == dimensionedScalar("zero", Qr_.dimensions(), 0.0);
+    forAll(intCoupledPatchIDs_, i)
+    {
+        const label patchI = intCoupledPatchIDs_[i];
+
+        scalarField& Qrp = Qr_.boundaryField()[patchI];
+
+        // Qr is negative going out the solid
+        // If the surface is emitting the radiative flux is set to zero
+        Qrp = max(Qrp, scalar(0.0));
+    }
+
+    // Propagate Qr through 1-D regions
+    forAll(intCoupledPatchIDs_, i)
+    {
+        const label patchI = intCoupledPatchIDs_[i];
+
+        const scalarField& Qrp = Qr_.boundaryField()[patchI];
+        const vectorField& Cf = regionMesh().Cf().boundaryField()[patchI];
+
+        forAll(Qrp, faceI)
+        {
+            const scalar Qr0 = Qrp[faceI];
+            point Cf0 = Cf[faceI];
+            const labelList& cells = boundaryFaceCells_[faceI];
+            scalar kappaInt = 0.0;
+            forAll(cells, k)
+            {
+                const label cellI = cells[k];
+                const point& Cf1 = regionMesh().cellCentres()[cellI];
+                const scalar delta = mag(Cf1 - Cf0);
+                kappaInt += kappa_[cellI]*delta;
+                Qr_[cellI] = Qr0*exp(-kappaInt);
+                Cf0 = Cf1;
+            }
+        }
+    }
+
+    Qr_.correctBoundaryConditions();
+}
+
+
+void reactingOneDim::updatePhiGas()
+{
+    phiHsGas_ ==  dimensionedScalar("zero", phiHsGas_.dimensions(), 0.0);
+    phiGas_ == dimensionedScalar("zero", phiGas_.dimensions(), 0.0);
+
+    const speciesTable& gasTable = solidChemistry_->gasTable();
+
+    forAll(gasTable, gasI)
+    {
+        tmp<volScalarField> tHsiGas = solidChemistry_->gasHs(T_, gasI);
+        tmp<volScalarField> tRRiGas = solidChemistry_->RRg(gasI);
+
+        const volScalarField& HsiGas = tHsiGas();
+        const volScalarField& RRiGas = tRRiGas();
+
+        const surfaceScalarField HsiGasf(fvc::interpolate(HsiGas));
+        const surfaceScalarField RRiGasf(fvc::interpolate(RRiGas));
+
+        forAll(intCoupledPatchIDs_, i)
+        {
+            const label patchI = intCoupledPatchIDs_[i];
+            const scalarField& phiGasp = phiHsGas_.boundaryField()[patchI];
+
+            forAll(phiGasp, faceI)
+            {
+                const labelList& cells = boundaryFaceCells_[faceI];
+                scalar massInt = 0.0;
+                forAllReverse(cells, k)
+                {
+                    const label cellI = cells[k];
+                    massInt += RRiGas[cellI]*regionMesh().V()[cellI];
+                    phiHsGas_[cellI] += massInt*HsiGas[cellI];
+                }
+
+                phiGas_.boundaryField()[patchI][faceI] += massInt;
+
+                if (debug)
+                {
+                    Info<< " Gas : " << gasTable[gasI]
+                        << " on patch : " << patchI
+                        << " mass produced at face(local) : "
+                        <<  faceI
+                        << " is : " << massInt
+                        << " [kg/s] " << endl;
+                }
+            }
+        }
+        tHsiGas().clear();
+    }
+}
+
+
+void reactingOneDim::updateFields()
+{
+    updateQr();
+
+    updatePhiGas();
+}
+
+
+void reactingOneDim::updateMesh(const scalarField& mass0)
+{
+    if (!moveMesh_)
+    {
+        return;
+    }
+
+    const scalarField newV(mass0/rho_);
+
+    Info<< "Initial/final volumes = " << gSum(regionMesh().V()) << ", "
+        << gSum(newV) << " [m3]" << endl;
+
+    // move the mesh
+    const labelList moveMap = moveMesh(regionMesh().V() - newV, minimumDelta_);
+
+    // flag any cells that have not moved as non-reacting
+    forAll(moveMap, i)
+    {
+        if (moveMap[i] == 0)
+        {
+            solidChemistry_->setCellReacting(i, false);
+        }
+    }
+}
+
+
+void reactingOneDim::solveContinuity()
+{
+    if (debug)
+    {
+        Info<< "reactingOneDim::solveContinuity()" << endl;
+    }
+
+    solve
+    (
+        fvm::ddt(rho_)
+     ==
+      - solidChemistry_->RRg()
+    );
+}
+
+
+void reactingOneDim::solveSpeciesMass()
+{
+    if (debug)
+    {
+        Info<< "reactingOneDim::solveSpeciesMass()" << endl;
+    }
+
+    volScalarField Yt(0.0*Ys_[0]);
+
+    for (label i=0; i<Ys_.size()-1; i++)
+    {
+        volScalarField& Yi = Ys_[i];
+
+        fvScalarMatrix YiEqn
+        (
+            fvm::ddt(rho_, Yi)
+         ==
+            solidChemistry_->RRs(i)
+        );
+
+        if (moveMesh_)
+        {
+            surfaceScalarField phiRhoMesh
+            (
+                fvc::interpolate(Yi*rho_)*regionMesh().phi()
+            );
+
+            YiEqn -= fvc::div(phiRhoMesh);
+        }
+
+        YiEqn.solve(regionMesh().solver("Yi"));
+        Yi.max(0.0);
+        Yt += Yi;
+    }
+
+    Ys_[Ys_.size() - 1] = 1.0 - Yt;
+}
+
+
+void reactingOneDim::solveEnergy()
+{
+    if (debug)
+    {
+        Info<< "reactingOneDim::solveEnergy()" << endl;
+    }
+
+    const volScalarField rhoCp(rho_*solidThermo_.Cp());
+
+    const surfaceScalarField phiQr(fvc::interpolate(Qr_)*nMagSf());
+
+    const surfaceScalarField phiGas(fvc::interpolate(phiHsGas_));
+
+    fvScalarMatrix TEqn
+    (
+        fvm::ddt(rhoCp, T_)
+      - fvm::laplacian(K_, T_)
+     ==
+        chemistrySh_
+      + fvc::div(phiQr)
+      + fvc::div(phiGas)
+    );
+
+    if (moveMesh_)
+    {
+        surfaceScalarField phiMesh
+        (
+            fvc::interpolate(rhoCp*T_)*regionMesh().phi()
+        );
+
+        TEqn -= fvc::div(phiMesh);
+    }
+
+    TEqn.relax();
+    TEqn.solve();
+
+    Info<< "pyrolysis min/max(T) = " << min(T_).value() << ", "
+        << max(T_).value() << endl;
+}
+
+
+void reactingOneDim::calculateMassTransfer()
+{
+    totalGasMassFlux_ = 0;
+    forAll(intCoupledPatchIDs_, i)
+    {
+        const label patchI = intCoupledPatchIDs_[i];
+        totalGasMassFlux_ += gSum(phiGas_.boundaryField()[patchI]);
+    }
+
+    if (infoOutput_)
+    {
+        totalHeatRR_ = fvc::domainIntegrate(chemistrySh_);
+
+        addedGasMass_ +=
+            fvc::domainIntegrate(solidChemistry_->RRg())*time_.deltaT();
+        lostSolidMass_ +=
+            fvc::domainIntegrate(solidChemistry_->RRs())*time_.deltaT();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
+:
+    pyrolysisModel(modelType, mesh),
+    solidChemistry_(solidChemistryModel::New(regionMesh())),
+    solidThermo_(solidChemistry_->solidThermo()),
+    kappa_(solidThermo_.kappa()),
+    K_(solidThermo_.K()),
+    rho_(solidThermo_.rho()),
+    Ys_(solidThermo_.composition().Y()),
+    T_(solidThermo_.T()),
+    primaryRadFluxName_(coeffs().lookupOrDefault<word>("radFluxName", "Qr")),
+    nNonOrthCorr_(-1),
+    maxDiff_(10),
+    minimumDelta_(1e-4),
+
+    phiGas_
+    (
+        IOobject
+        (
+            "phiGas",
+            time().timeName(),
+            regionMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar("zero", dimMass/dimTime, 0.0)
+    ),
+
+    phiHsGas_
+    (
+        IOobject
+        (
+            "phiHsGas",
+            time().timeName(),
+            regionMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar("zero", dimEnergy/dimTime, 0.0)
+    ),
+
+    chemistrySh_
+    (
+        IOobject
+        (
+            "chemistrySh",
+            time().timeName(),
+            regionMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
+    ),
+
+    QrCoupled_
+    (
+        IOobject
+        (
+            primaryRadFluxName_,
+            time().timeName(),
+            regionMesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh()
+    ),
+
+    Qr_
+    (
+        IOobject
+        (
+            "QrPyr",
+            time().timeName(),
+            regionMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar("zero", dimEnergy/dimArea/dimTime, 0.0),
+        zeroGradientFvPatchVectorField::typeName
+    ),
+
+    lostSolidMass_(dimensionedScalar("zero", dimMass, 0.0)),
+    addedGasMass_(dimensionedScalar("zero", dimMass, 0.0)),
+    totalGasMassFlux_(0.0),
+    totalHeatRR_(dimensionedScalar("zero", dimEnergy/dimTime, 0.0))
+{
+    if (active_)
+    {
+        read();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+reactingOneDim::~reactingOneDim()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+scalar reactingOneDim::addMassSources(const label patchI, const label faceI)
+{
+    label index = 0;
+    forAll(primaryPatchIDs_, i)
+    {
+        if (primaryPatchIDs_[i] == patchI)
+        {
+            index = i;
+            break;
+        }
+    }
+
+    const label localPatchId =  intCoupledPatchIDs_[index];
+
+    const scalar massAdded = phiGas_.boundaryField()[localPatchId][faceI];
+
+    if (debug)
+    {
+        Info<< "\nPyrolysis region: " << type() << "added mass : "
+            << massAdded << endl;
+    }
+
+    return massAdded;
+}
+
+
+scalar reactingOneDim::solidRegionDiffNo() const
+{
+    scalar DiNum = 0.0;
+    scalar meanDiNum = 0.0;
+    if (regionMesh().nInternalFaces() > 0)
+    {
+        surfaceScalarField KrhoCpbyDelta
+        (
+            regionMesh().surfaceInterpolation::deltaCoeffs()
+          * fvc::interpolate(K_)
+          / fvc::interpolate(Cp()*rho_)
+        );
+
+        DiNum = max(KrhoCpbyDelta.internalField())*time_.deltaTValue();
+
+        meanDiNum = average(KrhoCpbyDelta.internalField())*time().deltaTValue();
+    }
+
+    return DiNum;
+}
+
+
+scalar reactingOneDim::maxDiff() const
+{
+    return maxDiff_;
+}
+
+
+const volScalarField& reactingOneDim::rho() const
+{
+    return rho_;
+}
+
+
+const volScalarField& reactingOneDim::T() const
+{
+    return T_;
+}
+
+
+const tmp<volScalarField> reactingOneDim::Cp() const
+{
+    return solidThermo_.Cp();
+}
+
+
+const volScalarField& reactingOneDim::kappa() const
+{
+    return kappa_;
+}
+
+
+const volScalarField& reactingOneDim::K() const
+{
+    return K_;
+}
+
+
+const surfaceScalarField& reactingOneDim::phiGas() const
+{
+    return phiGas_;
+}
+
+
+void reactingOneDim::preEvolveRegion()
+{
+    pyrolysisModel::preEvolveRegion();
+
+    // Initialise all cells as able to react
+    forAll(T_, cellI)
+    {
+        solidChemistry_->setCellReacting(cellI, true);
+    }
+
+    // De-activate reactions if pyrolysis region coupled to (valid) film
+    if (filmCoupled_)
+    {
+        const volScalarField& filmDelta = filmDeltaPtr_();
+
+        forAll(intCoupledPatchIDs_, i)
+        {
+            const label patchI = intCoupledPatchIDs_[i];
+            const scalarField& filmDeltap = filmDelta.boundaryField()[patchI];
+
+            forAll(filmDeltap, faceI)
+            {
+                const scalar filmDelta0 = filmDeltap[faceI];
+                if (filmDelta0 > reactionDeltaMin_)
+                {
+                    const labelList& cells = boundaryFaceCells_[faceI];
+
+                    // TODO: only limit cell adjacent to film?
+                    //solidChemistry_->setCellNoReacting(cells[0])
+
+                    // Propagate flag through 1-D region
+                    forAll(cells, k)
+                    {
+                        solidChemistry_->setCellReacting(cells[k], false);
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+void reactingOneDim::evolveRegion()
+{
+    const scalarField mass0 = rho_*regionMesh().V();
+
+    solidChemistry_->solve
+    (
+        time().value() - time().deltaTValue(),
+        time().deltaTValue()
+    );
+
+    solveContinuity();
+
+    updateMesh(mass0);
+
+    chemistrySh_ = solidChemistry_->Sh()();
+
+    updateFields();
+
+    solveSpeciesMass();
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr_; nonOrth++)
+    {
+        solveEnergy();
+    }
+
+    calculateMassTransfer();
+
+    solidThermo_.correct();
+}
+
+
+void reactingOneDim::info() const
+{
+    Info<< "\nPyrolysis: " << type() << endl;
+
+    Info<< indent << "Total gas mass produced  [kg] = "
+        << addedGasMass_.value() << nl
+        << indent << "Total solid mass lost    [kg] = "
+        << lostSolidMass_.value() << nl
+        << indent << "Total pyrolysis gases  [kg/s] = "
+        << totalGasMassFlux_ << nl
+        << indent << "Total heat release rate [J/s] = "
+        << totalHeatRR_.value() << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace regionModels
+} // End namespace pyrolysisModels
+
+// ************************************************************************* //