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

Conflicts:
	doc/Doxygen/css/doxygen.css

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

@@ -4,7 +4,8 @@
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
       + turb.divDevRhoReff(U)
-      + fvOptions(rho, U)
+       ==
+        fvOptions(rho, U)
     );
 
     UEqn().relax();

applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C

@@ -77,8 +77,6 @@ int main(int argc, char *argv[])
 
         #include "setrDeltaT.H"
 
-        tmp<surfaceScalarField> tphiAlpha;
-
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
@@ -91,8 +89,6 @@ int main(int argc, char *argv[])
                 #define LTSSOLVE
                 #include "alphaEqnSubCycle.H"
                 #undef LTSSOLVE
-
-                interface.correct();
             }
 
             turbulence->correct();

applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C

@@ -81,8 +81,6 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        tmp<surfaceScalarField> tphiAlpha;
-
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {

applications/solvers/multiphase/interFoam/alphaEqn.H

@@ -2,11 +2,30 @@
     word alphaScheme("div(phi,alpha)");
     word alpharScheme("div(phirb,alpha)");
 
-    surfaceScalarField phic(mag(phi/mesh.magSf()));
+    // Standard face-flux compression coefficient
-    phic = min(interface.cAlpha()*phic, max(phic));
+    surfaceScalarField phic(interface.cAlpha()*mag(phi/mesh.magSf()));
-    surfaceScalarField phir(phic*interface.nHatf());
+
+    // Add the optional isotropic compression contribution
+    if (icAlpha > 0)
+    {
+        phic *= (1.0 - icAlpha);
+        phic += (interface.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
+    }
+
+    // Do not compress interface at non-coupled boundary faces
+    // (inlets, outlets etc.)
+    forAll(phic.boundaryField(), patchi)
+    {
+        fvsPatchScalarField& phicp = phic.boundaryField()[patchi];
+
+        if (!phicp.coupled())
+        {
+            phicp == 0;
+        }
+    }
+
+    tmp<surfaceScalarField> tphiAlpha;
 
-    //***HGW if (pimple.firstIter() && MULESCorr)
     if (MULESCorr)
     {
         fvScalarMatrix alpha1Eqn
@@ -32,12 +51,37 @@
             << "  Max(alpha1) = " << max(alpha1).value()
             << endl;
 
-        tphiAlpha = alpha1Eqn.flux();
+        tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux());
+        tphiAlpha = tmp<surfaceScalarField>
+        (
+            new surfaceScalarField(tphiAlphaUD())
+        );
+
+        if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
+        {
+            Info<< "Applying the previous iteration compression flux" << endl;
+            #ifdef LTSSOLVE
+            MULES::LTScorrect(alpha1, tphiAlpha(), tphiAlphaCorr0(), 1, 0);
+            #else
+            MULES::correct(alpha1, tphiAlpha(), tphiAlphaCorr0(), 1, 0);
+            #endif
+
+            tphiAlpha() += tphiAlphaCorr0();
+        }
+
+        // Cache the upwind-flux
+        tphiAlphaCorr0 = tphiAlphaUD;
+
+        alpha2 = 1.0 - alpha1;
+
+        interface.correct();
     }
 
     for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
     {
-        tmp<surfaceScalarField> tphiAlpha0
+        surfaceScalarField phir(phic*interface.nHatf());
+
+        tmp<surfaceScalarField> tphiAlphaUn
         (
             fvc::flux
             (
@@ -47,7 +91,7 @@
             )
           + fvc::flux
             (
-                -fvc::flux(-phir, alpha2, alpharScheme),
+               -fvc::flux(-phir, alpha2, alpharScheme),
                 alpha1,
                 alpharScheme
             )
@@ -55,17 +99,29 @@
 
         if (MULESCorr)
         {
-            tphiAlpha0() -= tphiAlpha();
+            tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - tphiAlpha());
+            volScalarField alpha10(alpha1);
+
             #ifdef LTSSOLVE
-            MULES::LTScorrect(alpha1, tphiAlpha0(), 1, 0);
+            MULES::LTScorrect(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
             #else
-            MULES::correct(alpha1, tphiAlpha0(), 1, 0);
+            MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
             #endif
-            tphiAlpha() += tphiAlpha0();
+
+            // Under-relax the correction for more than 3 correctors
+            if (aCorr < 3)
+            {
+                tphiAlpha() += tphiAlphaCorr();
+            }
+            else
+            {
+                alpha1 = 0.5*alpha1 + 0.5*alpha10;
+                tphiAlpha() += 0.5*tphiAlphaCorr();
+            }
         }
         else
         {
-            tphiAlpha = tphiAlpha0;
+            tphiAlpha = tphiAlphaUn;
 
             #ifdef LTSSOLVE
             MULES::explicitLTSSolve(alpha1, phi, tphiAlpha(), 1, 0);
@@ -75,10 +131,17 @@
         }
 
         alpha2 = 1.0 - alpha1;
+
+        interface.correct();
     }
 
     rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2;
 
+    if (alphaApplyPrevCorr && MULESCorr)
+    {
+        tphiAlphaCorr0 = tphiAlpha() - tphiAlphaCorr0;
+    }
+
     Info<< "Phase-1 volume fraction = "
         << alpha1.weightedAverage(mesh.Vsc()).value()
         << "  Min(alpha1) = " << min(alpha1).value()

applications/solvers/multiphase/interFoam/createFields.H

@@ -135,3 +135,6 @@
     }
 
     fv::IOoptionList fvOptions(mesh);
+
+
+    tmp<surfaceScalarField> tphiAlphaCorr0;

applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C

@@ -77,8 +77,6 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        tmp<surfaceScalarField> tphiAlpha;
-
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {

applications/solvers/multiphase/interFoam/interFoam.C

@@ -80,8 +80,6 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        tmp<surfaceScalarField> tphiAlpha;
-
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {

applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C

@@ -74,8 +74,6 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        tmp<surfaceScalarField> tphiAlpha;
-
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {

applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C

@@ -83,8 +83,6 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        tmp<surfaceScalarField> tphiAlpha;
-
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H

@@ -71,6 +71,7 @@
             (
                 geometricOneField(),
                 alpha1,
+                tphiAlpha(),
                 tphiAlphaCorr(),
                 vDotvmcAlphal,
                 (

doc/Doxygen/css/doxygen.css

@@ -11,49 +11,49 @@ body, table, div, p, dl {
 /* @group Heading Levels */
 
 h1 {
-	text-align: center;
+        text-align: center;
-	font-size: 150%;
+        font-size: 150%;
 }
 
 h2 {
-	font-size: 120%;
+        font-size: 120%;
 }
 
 h3 {
-	font-size: 100%;
+        font-size: 100%;
 }
 
 dt {
-	font-weight: bold;
+        font-weight: bold;
 }
 
 div.multicol {
-	-moz-column-gap: 1em;
+        -moz-column-gap: 1em;
-	-webkit-column-gap: 1em;
+        -webkit-column-gap: 1em;
-	-moz-column-count: 3;
+        -moz-column-count: 3;
-	-webkit-column-count: 3;
+        -webkit-column-count: 3;
 }
 
 p.startli, p.startdd, p.starttd {
-	margin-top: 2px;
+        margin-top: 2px;
 }
 
 p.endli {
-	margin-bottom: 0px;
+        margin-bottom: 0px;
 }
 
 p.enddd {
-	margin-bottom: 4px;
+        margin-bottom: 4px;
 }
 
 p.endtd {
-	margin-bottom: 2px;
+        margin-bottom: 2px;
 }
 
 /* @end */
 
 caption {
-	font-weight: bold;
+        font-weight: bold;
 }
 
 span.legend {
@@ -67,47 +67,47 @@ h3.version {
 }
 
 div.qindex, div.navtab{
-	background-color: #e8eef2;
+        background-color: #e8eef2;
-	border: 1px solid #84b0c7;
+        border: 1px solid #84b0c7;
-	text-align: center;
+        text-align: center;
-	margin: 2px;
+        margin: 2px;
-	padding: 2px;
+        padding: 2px;
 }
 
 div.qindex, div.navpath {
-	width: 100%;
+        width: 100%;
-	line-height: 140%;
+        line-height: 140%;
 }
 
 div.navtab {
-	margin-right: 15px;
+        margin-right: 15px;
 }
 
 /* @group Link Styling */
 
 a {
-	color: #153788;
+        color: #153788;
-	font-weight: normal;
+        font-weight: normal;
-	text-decoration: none;
+        text-decoration: none;
 }
 
 .contents a:visited {
-	color: #1b77c5;
+        color: #1b77c5;
 }
 
 a:hover {
-	text-decoration: underline;
+        text-decoration: underline;
 }
 
 a.qindex {
-	font-weight: bold;
+        font-weight: bold;
 }
 
 a.qindexHL {
-	font-weight: bold;
+        font-weight: bold;
-	background-color: #6666cc;
+        background-color: #6666cc;
-	color: #ffffff;
+        color: #ffffff;
-	border: 1px double #9295C2;
+        border: 1px double #9295C2;
 }
 
 .contents a.qindexHL:visited {
@@ -115,7 +115,7 @@ a.qindexHL {
 }
 
 a.el {
-	font-weight: bold;
+        font-weight: bold;
 }
 
 a.elRef {
@@ -134,7 +134,7 @@ a.codeRef {
 /* @end */
 
 dl.el {
-	margin-left: -1cm;
+        margin-left: -1cm;
 }
 
 .fragment {
@@ -143,64 +143,64 @@ dl.el {
 }
 
 pre.fragment {
-	border: 1px solid #CCCCCC;
+        border: 1px solid #CCCCCC;
-	background-color: #f5f5f5;
+        background-color: #f5f5f5;
-	padding: 4px 6px;
+        padding: 4px 6px;
-	margin: 4px 8px 4px 2px;
+        margin: 4px 8px 4px 2px;
-	overflow: auto;
+        overflow: auto;
-	word-wrap: break-word;
+        word-wrap: break-word;
-	font-size:  9pt;
+        font-size:  9pt;
-	line-height: 125%;
+        line-height: 125%;
 }
 
 div.ah {
-	background-color: black;
+        background-color: black;
-	font-weight: bold;
+        font-weight: bold;
-	color: #ffffff;
+        color: #ffffff;
-	margin-bottom: 3px;
+        margin-bottom: 3px;
-	margin-top: 3px
+        margin-top: 3px
 }
 
 div.groupHeader {
-	margin-left: 16px;
+        margin-left: 16px;
-	margin-top: 12px;
+        margin-top: 12px;
-	margin-bottom: 6px;
+        margin-bottom: 6px;
-	font-weight: bold;
+        font-weight: bold;
 }
 
 div.groupText {
-	margin-left: 16px;
+        margin-left: 16px;
-	font-style: italic;
+        font-style: italic;
 }
 
 body {
-	background: white;
+        background: white;
-	color: black;
+        color: black;
-	margin-right: 20px;
+        margin-right: 20px;
-	margin-left: 20px;
+        margin-left: 20px;
 }
 
 td.indexkey {
-	background-color: #e8eef2;
+        background-color: #e8eef2;
-	font-weight: bold;
+        font-weight: bold;
-	border: 1px solid #CCCCCC;
+        border: 1px solid #CCCCCC;
-	margin: 2px 0px 2px 0;
+        margin: 2px 0px 2px 0;
-	padding: 2px 10px;
+        padding: 2px 10px;
 }
 
 td.indexvalue {
-	background-color: #e8eef2;
+        background-color: #e8eef2;
-	border: 1px solid #CCCCCC;
+        border: 1px solid #CCCCCC;
-	padding: 2px 10px;
+        padding: 2px 10px;
-	margin: 2px 0px;
+        margin: 2px 0px;
 }
 
 tr.memlist {
-	background-color: #f0f0f0;
+        background-color: #f0f0f0;
 }
 
 p.formulaDsp {
-	text-align: center;
+        text-align: center;
 }
 
 img.formulaDsp {
@@ -208,23 +208,23 @@ img.formulaDsp {
 }
 
 img.formulaInl {
-	vertical-align: middle;
+        vertical-align: middle;
 }
 
 div.center {
-	text-align: center;
+        text-align: center;
         margin-top: 0px;
         margin-bottom: 0px;
         padding: 0px;
 }
 
 div.center img {
-	border: 0px;
+        border: 0px;
 }
 
 img.footer {
-	border: 0px;
+        border: 0px;
-	vertical-align: middle;
+        vertical-align: middle;
 }
 
 /* @group Code Colorization */
@@ -287,45 +287,45 @@ span.vhdllogic {
 /* @end */
 
 .search {
-	color: #003399;
+        color: #003399;
-	font-weight: bold;
+        font-weight: bold;
 }
 
 form.search {
-	margin-bottom: 0px;
+        margin-bottom: 0px;
-	margin-top: 0px;
+        margin-top: 0px;
 }
 
 input.search {
-	font-size: 75%;
+        font-size: 75%;
-	color: #000080;
+        color: #000080;
-	font-weight: normal;
+        font-weight: normal;
-	background-color: #e8eef2;
+        background-color: #e8eef2;
 }
 
 td.tiny {
-	font-size: 75%;
+        font-size: 75%;
 }
 
 .dirtab {
-	padding: 4px;
+        padding: 4px;
-	border-collapse: collapse;
+        border-collapse: collapse;
-	border: 1px solid #84b0c7;
+        border: 1px solid #84b0c7;
 }
 
 th.dirtab {
-	background: #e8eef2;
+        background: #e8eef2;
-	font-weight: bold;
+        font-weight: bold;
 }
 
 hr {
-	height: 0px;
+        height: 0px;
-	border: none;
+        border: none;
-	border-top: 1px solid #666;
+        border-top: 1px solid #666;
 }
 
 hr.footer {
-	height: 1px;
+        height: 1px;
 }
 
 /* @group Member Descriptions */
@@ -333,19 +333,19 @@ hr.footer {
 .mdescLeft, .mdescRight,
 .memItemLeft, .memItemRight,
 .memTemplItemLeft, .memTemplItemRight, .memTemplParams {
-	background-color: #FAFAFA;
+        background-color: #FAFAFA;
-	border: none;
+        border: none;
-	margin: 4px;
+        margin: 4px;
-	padding: 1px 0 0 8px;
+        padding: 1px 0 0 8px;
 }
 
 .mdescLeft, .mdescRight {
-	padding: 0px 8px 4px 8px;
+        padding: 0px 8px 4px 8px;
-	color: #555;
+        color: #555;
 }
 
 .memItemLeft, .memItemRight, .memTemplParams {
-	border-top: 1px solid #ccc;
+        border-top: 1px solid #ccc;
 }
 
 .memItemLeft, .memTemplItemLeft {
@@ -353,7 +353,7 @@ hr.footer {
 }
 
 .memTemplParams {
-	color: #606060;
+        color: #606060;
         white-space: nowrap;
 }
 
@@ -364,24 +364,24 @@ hr.footer {
 /* Styles for detailed member documentation */
 
 .memtemplate {
-	font-size: 80%;
+        font-size: 80%;
-	color: #606060;
+        color: #606060;
-	font-weight: normal;
+        font-weight: normal;
-	margin-left: 3px;
+        margin-left: 3px;
 }
 
 .memnav {
-	background-color: #e8eef2;
+        background-color: #e8eef2;
-	border: 1px solid #84b0c7;
+        border: 1px solid #84b0c7;
-	text-align: center;
+        text-align: center;
-	margin: 2px;
+        margin: 2px;
-	margin-right: 15px;
+        margin-right: 15px;
-	padding: 2px;
+        padding: 2px;
 }
 
 .memitem {
-	padding: 0;
+        padding: 0;
-	margin-bottom: 10px;
+        margin-bottom: 10px;
 }
 
 .memname {
@@ -428,19 +428,19 @@ hr.footer {
 }
 
 .paramkey {
-	text-align: right;
+        text-align: right;
 }
 
 .paramtype {
-	white-space: nowrap;
+        white-space: nowrap;
 }
 
 .paramname {
-	color: #602020;
+        color: #602020;
-	white-space: nowrap;
+        white-space: nowrap;
 }
 .paramname em {
-	font-style: normal;
+        font-style: normal;
 }
 
 /* @end */
@@ -450,21 +450,21 @@ hr.footer {
 /* for the tree view */
 
 .ftvtree {
-	font-family: sans-serif;
+        font-family: sans-serif;
-	margin: 0.5em;
+        margin: 0.5em;
 }
 
 /* these are for tree view when used as main index */
 
 .directory {
-	font-size: 9pt;
+        font-size: 9pt;
-	font-weight: bold;
+        font-weight: bold;
 }
 
 .directory h3 {
-	margin: 0px;
+        margin: 0px;
-	margin-top: 1em;
+        margin-top: 1em;
-	font-size: 11pt;
+        font-size: 11pt;
 }
 
 /*
@@ -476,86 +476,86 @@ proper pixel height of your image.
 
 /*
 .directory h3.swap {
-	height: 61px;
+        height: 61px;
-	background-repeat: no-repeat;
+        background-repeat: no-repeat;
-	background-image: url("yourimage.gif");
+        background-image: url("yourimage.gif");
 }
 .directory h3.swap span {
-	display: none;
+        display: none;
 }
 */
 
 .directory > h3 {
-	margin-top: 0;
+        margin-top: 0;
 }
 
 .directory p {
-	margin: 0px;
+        margin: 0px;
-	white-space: nowrap;
+        white-space: nowrap;
 }
 
 .directory div {
-	display: none;
+        display: none;
-	margin: 0px;
+        margin: 0px;
 }
 
 .directory img {
-	vertical-align: -30%;
+        vertical-align: -30%;
 }
 
 /* these are for tree view when not used as main index */
 
 .directory-alt {
-	font-size: 100%;
+        font-size: 100%;
-	font-weight: bold;
+        font-weight: bold;
 }
 
 .directory-alt h3 {
-	margin: 0px;
+        margin: 0px;
-	margin-top: 1em;
+        margin-top: 1em;
-	font-size: 11pt;
+        font-size: 11pt;
 }
 
 .directory-alt > h3 {
-	margin-top: 0;
+        margin-top: 0;
 }
 
 .directory-alt p {
-	margin: 0px;
+        margin: 0px;
-	white-space: nowrap;
+        white-space: nowrap;
 }
 
 .directory-alt div {
-	display: none;
+        display: none;
-	margin: 0px;
+        margin: 0px;
 }
 
 .directory-alt img {
-	vertical-align: -30%;
+        vertical-align: -30%;
 }
 
 /* @end */
 
 address {
-	font-style: normal;
+        font-style: normal;
-	color: #333;
+        color: #333;
 }
 
 table.doxtable {
-	border-collapse:collapse;
+        border-collapse:collapse;
 }
 
 table.doxtable td, table.doxtable th {
-	border: 1px solid #153788;
+        border: 1px solid #153788;
-	padding: 3px 7px 2px;
+        padding: 3px 7px 2px;
 }
 
 table.doxtable th {
-	background-color: #254798;
+        background-color: #254798;
-	color: #FFFFFF;
+        color: #FFFFFF;
-	font-size: 110%;
+        font-size: 110%;
-	padding-bottom: 4px;
+        padding-bottom: 4px;
-	padding-top: 5px;
+        padding-top: 5px;
-	text-align:left;
+        text-align:left;
 }
 

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

@@ -110,7 +110,8 @@ Foam::OutputFilterFunctionObject<OutputFilter>::OutputFilterFunctionObject
     storeFilter_(true),
     timeStart_(-VGREAT),
     timeEnd_(VGREAT),
-    outputControl_(t, dict)
+    outputControl_(t, dict, "output"),
+    evaluateControl_(t, dict, "evaluate")
 {
     readDict();
 }
@@ -159,7 +160,10 @@ bool Foam::OutputFilterFunctionObject<OutputFilter>::execute
             allocateFilter();
         }
 
-        ptr_->execute();
+        if (evaluateControl_.output())
+        {
+            ptr_->execute();
+        }
 
         if (forceWrite || outputControl_.output())
         {
@@ -241,14 +245,14 @@ bool Foam::OutputFilterFunctionObject<OutputFilter>::adjustTimeStep()
 
         // function objects modify deltaT inside nStepsToStartTimeChange range
         // NOTE: Potential problem if two function objects dump inside the same
-        //interval
+        // interval
         if (nSteps < nStepsToStartTimeChange_)
         {
             label nStepsToNextWrite = label(nSteps) + 1;
 
             scalar newDeltaT = timeToNextWrite/nStepsToNextWrite;
 
-            //Adjust time step
+            // Adjust time step
             if (newDeltaT < deltaT)
             {
                 deltaT = max(newDeltaT, 0.2*deltaT);

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

@@ -98,6 +98,9 @@ class OutputFilterFunctionObject
         //- Output controls
         outputFilterOutputControl outputControl_;
 
+        //- Evaluate controls
+        outputFilterOutputControl evaluateControl_;
+
         //- Pointer to the output filter
         autoPtr<OutputFilter> ptr_;
 

src/OpenFOAM/db/functionObjects/functionObjectFile/functionObjectFile.C

@@ -31,9 +31,18 @@ License
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 const Foam::word Foam::functionObjectFile::outputPrefix = "postProcessing";
+Foam::label Foam::functionObjectFile::addChars = 7;
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+void Foam::functionObjectFile::initStream(Ostream& os) const
+{
+    os.setf(ios_base::scientific, ios_base::floatfield);
+//    os.precision(IOstream::defaultPrecision());
+    os.width(charWidth());
+}
+
+
 Foam::fileName Foam::functionObjectFile::baseFileDir() const
 {
     fileName baseDir = obr_.time().path();
@@ -96,6 +105,8 @@ void Foam::functionObjectFile::createFiles()
 
                 filePtrs_.set(i, new OFstream(outputDir/(fName + ".dat")));
 
+                initStream(filePtrs_[i]);
+
                 writeFileHeader(i);
 
                 i++;
@@ -147,6 +158,12 @@ void Foam::functionObjectFile::resetName(const word& name)
 }
 
 
+Foam::Omanip<int> Foam::functionObjectFile::valueWidth(const label offset) const
+{
+    return setw(IOstream::defaultPrecision() + addChars + offset);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::functionObjectFile::functionObjectFile
@@ -289,4 +306,42 @@ Foam::OFstream& Foam::functionObjectFile::file(const label i)
 }
 
 
+Foam::label Foam::functionObjectFile::charWidth() const
+{
+    return IOstream::defaultPrecision() + addChars;
+}
+
+
+void Foam::functionObjectFile::writeCommented
+(
+    Ostream& os,
+    const string& str
+) const
+{
+    os  << setw(1) << "#" << setw(1) << ' '
+        << setw(charWidth() - 2) << str.c_str();
+}
+
+
+void Foam::functionObjectFile::writeTabbed
+(
+    Ostream& os,
+    const string& str
+) const
+{
+    os  << tab << setw(charWidth()) << str.c_str();
+}
+
+
+void Foam::functionObjectFile::writeHeader
+(
+    Ostream& os,
+    const string& str
+) const
+{
+    os  << setw(1) << "#" << setw(1) << ' '
+        << setf(ios_base::left) << setw(charWidth() - 2) << str.c_str() << nl;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/db/functionObjects/functionObjectFile/functionObjectFile.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,6 +43,7 @@ SourceFiles
 #include "OFstream.H"
 #include "PtrList.H"
 #include "HashSet.H"
+#include "IOmanip.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -77,6 +78,9 @@ protected:
 
     // Protected Member Functions
 
+        //- Initialise the output stream for writing
+        virtual void initStream(Ostream& os) const;
+
         //- Return the base directory for output
         virtual fileName baseFileDir() const;
 
@@ -98,6 +102,9 @@ protected:
         //- Reset the list of names to a single name entry
         virtual void resetName(const word& name);
 
+        //- Return the value width when writing to stream with optional offset
+        virtual Omanip<int> valueWidth(const label offset = 0) const;
+
         //- Disallow default bitwise copy construct
         functionObjectFile(const functionObjectFile&);
 
@@ -110,6 +117,9 @@ public:
     //- Folder prefix
     static const word outputPrefix;
 
+    //- Additional characters for writing
+    static label addChars;
+
     // Constructors
 
         //- Construct null
@@ -149,6 +159,39 @@ public:
 
         //- Return file 'i'
         OFstream& file(const label i);
+
+        //- Write a commented string to stream
+        void writeCommented
+        (
+            Ostream& os,
+            const string& str
+        ) const;
+
+        //- Write a tabbed string to stream
+        void writeTabbed
+        (
+            Ostream& os,
+            const string& str
+        ) const;
+
+        //- Write a commented header to stream
+        void writeHeader
+        (
+            Ostream& os,
+            const string& str
+        ) const;
+
+        //- Write a (commented) header property and value pair
+        template<class Type>
+        void writeHeaderValue
+        (
+            Ostream& os,
+            const string& property,
+            const Type& value
+        ) const;
+
+        //- Return width of character stream output
+        label charWidth() const;
 };
 
 
@@ -158,6 +201,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+    #include "functionObjectFileTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/OpenFOAM/db/functionObjects/functionObjectFile/functionObjectFileTemplates.C

@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "OStringStream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::functionObjectFile::writeHeaderValue
+(
+    Ostream& os,
+    const string& property,
+    const Type& value
+) const
+{
+    os  << setw(1) << '#' << setw(1) << ' '
+        << setw(charWidth() - 2) << setf(ios_base::left) << property.c_str()
+        << setw(1) << ':' << setw(1) << ' ' << value << nl;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/db/functionObjects/outputFilterOutputControl/outputFilterOutputControl.C

@@ -31,7 +31,7 @@ License
 namespace Foam
 {
     template<>
-    const char* NamedEnum<outputFilterOutputControl::outputControls, 6>::
+    const char* NamedEnum<outputFilterOutputControl::outputControls, 7>::
     names[] =
     {
         "timeStep",
@@ -39,11 +39,12 @@ namespace Foam
         "adjustableTime",
         "runTime",
         "clockTime",
-        "cpuTime"
+        "cpuTime",
+        "none"
     };
 }
 
-const Foam::NamedEnum<Foam::outputFilterOutputControl::outputControls, 6>
+const Foam::NamedEnum<Foam::outputFilterOutputControl::outputControls, 7>
     Foam::outputFilterOutputControl::outputControlNames_;
 
 
@@ -52,10 +53,12 @@ const Foam::NamedEnum<Foam::outputFilterOutputControl::outputControls, 6>
 Foam::outputFilterOutputControl::outputFilterOutputControl
 (
     const Time& t,
-    const dictionary& dict
+    const dictionary& dict,
+    const word& prefix
 )
 :
     time_(t),
+    prefix_(prefix),
     outputControl_(ocTimeStep),
     outputInterval_(0),
     outputTimeLastDump_(0),
@@ -75,9 +78,12 @@ Foam::outputFilterOutputControl::~outputFilterOutputControl()
 
 void Foam::outputFilterOutputControl::read(const dictionary& dict)
 {
-    if (dict.found("outputControl"))
+    const word controlName(prefix_ + "Control");
+    const word intervalName(prefix_ + "Interval");
+
+    if (dict.found(controlName))
     {
-        outputControl_ = outputControlNames_.read(dict.lookup("outputControl"));
+        outputControl_ = outputControlNames_.read(dict.lookup(controlName));
     }
     else
     {
@@ -88,13 +94,13 @@ void Foam::outputFilterOutputControl::read(const dictionary& dict)
     {
         case ocTimeStep:
         {
-            outputInterval_ = dict.lookupOrDefault<label>("outputInterval", 0);
+            outputInterval_ = dict.lookupOrDefault<label>(intervalName, 0);
             break;
         }
 
         case ocOutputTime:
         {
-            outputInterval_ = dict.lookupOrDefault<label>("outputInterval", 1);
+            outputInterval_ = dict.lookupOrDefault<label>(intervalName, 1);
             break;
         }
 
@@ -190,6 +196,11 @@ bool Foam::outputFilterOutputControl::output()
             break;
         }
 
+        case ocNone:
+        {
+            return false;
+        }
+
         default:
         {
             // this error should not actually be possible

src/OpenFOAM/db/functionObjects/outputFilterOutputControl/outputFilterOutputControl.H

@@ -61,7 +61,8 @@ public:
         ocAdjustableTime, /*!< Adjust time step for dumping */
         ocRunTime,        /*!< run time for dumping */
         ocClockTime,      /*!< clock time for dumping */
-        ocCpuTime         /*!< cpu time for dumping */
+        ocCpuTime,        /*!< cpu time for dumping */
+        ocNone            /*!< no output */
     };
 
 
@@ -72,8 +73,11 @@ private:
         //- Time object
         const Time& time_;
 
+        //- Prefix
+        const word prefix_;
+
         //- String representation of outputControls enums
-        static const NamedEnum<outputControls, 6> outputControlNames_;
+        static const NamedEnum<outputControls, 7> outputControlNames_;
 
         //- Type of output
         outputControls outputControl_;
@@ -101,7 +105,12 @@ public:
     // Constructors
 
         //- Construct from Time object and dictionary
-        outputFilterOutputControl(const Time&, const dictionary&);
+        outputFilterOutputControl
+        (
+            const Time&,
+            const dictionary&,
+            const word& prefix
+        );
 
 
     //- Destructor

src/OpenFOAM/primitives/functions/DataEntry/CSV/CSV.H

@@ -32,7 +32,7 @@ Description
         <entryName>   csvFile;
         csvFileCoeffs
         {
-            nHeaderLines        4;
+            nHeaderLine         4;
             refColumn           0;          // reference column index
             componentColumns    (1 2 3);    // component column indices
             separator           ",";        // optional (defaults to ",")

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C

@@ -3185,7 +3185,8 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
 
     forAll(refineCell, cellI)
     {
-        if (refineCell.get(cellI))
+//        if (refineCell.get(cellI))
+        if (refineCell[cellI])
         {
             nRefined++;
         }
@@ -3196,7 +3197,8 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
 
     forAll(refineCell, cellI)
     {
-        if (refineCell.get(cellI))
+//        if (refineCell.get(cellI))
+        if (refineCell[cellI])
         {
             newCellsToRefine[nRefined++] = cellI;
         }

src/finiteVolume/Make/files

@@ -189,9 +189,11 @@ $(derivedFvPatchFields)/uniformJumpAMI/uniformJumpAMIFvPatchFields.C
 $(derivedFvPatchFields)/uniformTotalPressure/uniformTotalPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/variableHeightFlowRate/variableHeightFlowRateFvPatchField.C
 $(derivedFvPatchFields)/variableHeightFlowRateInletVelocity/variableHeightFlowRateInletVelocityFvPatchVectorField.C
+$(derivedFvPatchFields)/outletPhaseMeanVelocity/outletPhaseMeanVelocityFvPatchVectorField.C
 $(derivedFvPatchFields)/waveTransmissive/waveTransmissiveFvPatchFields.C
 $(derivedFvPatchFields)/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
 $(derivedFvPatchFields)/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.C
+$(derivedFvPatchFields)/prghPressure/prghPressureFvPatchScalarField.C
 
 fvsPatchFields = fields/fvsPatchFields
 $(fvsPatchFields)/fvsPatchField/fvsPatchFields.C
@@ -376,6 +378,7 @@ $(laplacianSchemes)/gaussLaplacianScheme/gaussLaplacianSchemes.C
 
 finiteVolume/fvc/fvcMeshPhi.C
 finiteVolume/fvc/fvcSmooth/fvcSmooth.C
+finiteVolume/fvc/fvcReconstructMag.C
 
 general = cfdTools/general
 $(general)/findRefCell/findRefCell.C

src/finiteVolume/cfdTools/general/include/alphaControls.H

@@ -10,3 +10,17 @@ bool alphaOuterCorrectors
 (
     alphaControls.lookupOrDefault<Switch>("alphaOuterCorrectors", false)
 );
+
+// Apply the compression correction from the previous iteration
+// Improves efficiency for steady-simulations but can only be applied
+// once the alpha field is reasonably steady, i.e. fully developed
+bool alphaApplyPrevCorr
+(
+    alphaControls.lookupOrDefault<Switch>("alphaApplyPrevCorr", false)
+);
+
+// Isotropic compression coefficient
+scalar icAlpha
+(
+    alphaControls.lookupOrDefault<scalar>("icAlpha", 0)
+);

src/finiteVolume/fields/fvPatchFields/derived/outletPhaseMeanVelocity/outletPhaseMeanVelocityFvPatchVectorField.C

@@ -0,0 +1,188 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "outletPhaseMeanVelocityFvPatchVectorField.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::outletPhaseMeanVelocityFvPatchVectorField
+::outletPhaseMeanVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    mixedFvPatchField<vector>(p, iF),
+    Umean_(0),
+    alphaName_("none")
+{
+    refValue() = vector::zero;
+    refGrad() = vector::zero;
+    valueFraction() = 0.0;
+}
+
+
+Foam::outletPhaseMeanVelocityFvPatchVectorField
+::outletPhaseMeanVelocityFvPatchVectorField
+(
+    const outletPhaseMeanVelocityFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchField<vector>(ptf, p, iF, mapper),
+    Umean_(ptf.Umean_),
+    alphaName_(ptf.alphaName_)
+{}
+
+
+Foam::outletPhaseMeanVelocityFvPatchVectorField
+::outletPhaseMeanVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchField<vector>(p, iF),
+    Umean_(readScalar(dict.lookup("Umean"))),
+    alphaName_(dict.lookup("alpha"))
+{
+    refValue() = vector::zero;
+    refGrad() = vector::zero;
+    valueFraction() = 0.0;
+
+    if (dict.found("value"))
+    {
+        fvPatchVectorField::operator=
+        (
+            vectorField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        fvPatchVectorField::operator=(patchInternalField());
+    }
+}
+
+
+Foam::outletPhaseMeanVelocityFvPatchVectorField
+::outletPhaseMeanVelocityFvPatchVectorField
+(
+    const outletPhaseMeanVelocityFvPatchVectorField& ptf
+)
+:
+    mixedFvPatchField<vector>(ptf),
+    Umean_(ptf.Umean_),
+    alphaName_(ptf.alphaName_)
+{}
+
+
+Foam::outletPhaseMeanVelocityFvPatchVectorField
+::outletPhaseMeanVelocityFvPatchVectorField
+(
+    const outletPhaseMeanVelocityFvPatchVectorField& ptf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    mixedFvPatchField<vector>(ptf, iF),
+    Umean_(ptf.Umean_),
+    alphaName_(ptf.alphaName_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::outletPhaseMeanVelocityFvPatchVectorField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    scalarField alphap =
+        patch().lookupPatchField<volScalarField, scalar>(alphaName_);
+
+    alphap = max(alphap, scalar(0));
+    alphap = min(alphap, scalar(1));
+
+    // Get the patchInternalField (zero-gradient field)
+    vectorField Uzg(patchInternalField());
+
+    // Calculate the phase mean zero-gradient velocity
+    scalar Uzgmean =
+        gSum(alphap*(patch().Sf() & Uzg))
+       /gSum(alphap*patch().magSf());
+
+    // Set the refValue and valueFraction to adjust the boundary field
+    // such that the phase mean is Umean_
+    if (Uzgmean >= Umean_)
+    {
+        refValue() = vector::zero;
+        valueFraction() = 1.0 - Umean_/Uzgmean;
+    }
+    else
+    {
+        refValue() = (Umean_ + Uzgmean)*patch().nf();
+        valueFraction() = 1.0 - Uzgmean/Umean_;
+    }
+
+    mixedFvPatchField<vector>::updateCoeffs();
+}
+
+
+void Foam::outletPhaseMeanVelocityFvPatchVectorField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchField<vector>::write(os);
+
+    os.writeKeyword("Umean") << Umean_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("alpha") << alphaName_
+        << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+   makePatchTypeField
+   (
+       fvPatchVectorField,
+       outletPhaseMeanVelocityFvPatchVectorField
+   );
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/outletPhaseMeanVelocity/outletPhaseMeanVelocityFvPatchVectorField.H

@@ -0,0 +1,197 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::outletPhaseMeanVelocityFvPatchVectorField
+
+Group
+    grpOutletBoundaryConditions
+
+Description
+    This boundary condition adjusts the velocity for the given phase to achieve
+    the specified mean thus causing the phase-fraction to adjust according to
+    the mass flow rate.
+
+    Typical usage is as the outlet condition for a towing-tank ship simulation
+    to maintain the outlet water level at the level as the inlet.
+
+    \heading Patch usage
+    \table
+        Property     | Description             | Required    | Default value
+        Umean        | mean velocity normal to the boundary [m/s] | yes |
+        alpha        | phase-fraction field    | yes |
+    \endtable
+
+    Example of the boundary condition specification:
+    \verbatim
+    myPatch
+    {
+        type            outletPhaseMeanVelocity;
+        Umean           1.2;
+        alpha           alpha.water;
+        value           uniform (1.2 0 0);
+    }
+    \endverbatim
+
+SeeAlso
+    Foam::mixedFvPatchField
+    Foam::variableHeightFlowRateInletVelocityFvPatchVectorField
+
+SourceFiles
+    outletPhaseMeanVelocityFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef outletPhaseMeanVelocityFvPatchVectorField_H
+#define outletPhaseMeanVelocityFvPatchVectorField_H
+
+#include "mixedFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+    Class outletPhaseMeanVelocityFvPatchVectorField Declaration
+\*---------------------------------------------------------------------------*/
+
+class outletPhaseMeanVelocityFvPatchVectorField
+:
+    public mixedFvPatchVectorField
+{
+    // Private data
+
+        //- Inlet integral flow rate
+        scalar Umean_;
+
+        //- Name of the phase-fraction field
+        word alphaName_;
+
+
+public:
+
+   //- Runtime type information
+   TypeName("outletPhaseMeanVelocity");
+
+
+   // Constructors
+
+        //- Construct from patch and internal field
+        outletPhaseMeanVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        outletPhaseMeanVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  outletPhaseMeanVelocityFvPatchVectorField
+        //  onto a new patch
+        outletPhaseMeanVelocityFvPatchVectorField
+        (
+            const outletPhaseMeanVelocityFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        outletPhaseMeanVelocityFvPatchVectorField
+        (
+            const outletPhaseMeanVelocityFvPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new outletPhaseMeanVelocityFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        outletPhaseMeanVelocityFvPatchVectorField
+        (
+            const outletPhaseMeanVelocityFvPatchVectorField&,
+            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 outletPhaseMeanVelocityFvPatchVectorField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+
+    // Member functions
+
+        // Access
+
+            //- Return the flux
+            scalar Umean() const
+            {
+                return Umean_;
+            }
+
+            //- Return reference to the flux to allow adjustment
+            scalar& Umean()
+            {
+                return Umean_;
+            }
+
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/pressureInletOutletParSlipVelocity/pressureInletOutletParSlipVelocityFvPatchVectorField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,8 +41,7 @@ pressureInletOutletParSlipVelocityFvPatchVectorField
 :
     mixedFvPatchVectorField(p, iF),
     phiName_("phi"),
-    rhoName_("rho"),
+    rhoName_("rho")
-    UName_("U")
 {
     refValue() = *this;
     refGrad() = vector::zero;
@@ -61,8 +60,7 @@ pressureInletOutletParSlipVelocityFvPatchVectorField
 :
     mixedFvPatchVectorField(ptf, p, iF, mapper),
     phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_),
+    rhoName_(ptf.rhoName_)
-    UName_(ptf.UName_)
 {}
 
 
@@ -76,8 +74,7 @@ pressureInletOutletParSlipVelocityFvPatchVectorField
 :
     mixedFvPatchVectorField(p, iF),
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
-    UName_(dict.lookupOrDefault<word>("U", "U"))
 {
     fvPatchVectorField::operator=(vectorField("value", dict, p.size()));
     refValue() = *this;
@@ -94,8 +91,7 @@ pressureInletOutletParSlipVelocityFvPatchVectorField
 :
     mixedFvPatchVectorField(pivpvf),
     phiName_(pivpvf.phiName_),
-    rhoName_(pivpvf.rhoName_),
+    rhoName_(pivpvf.rhoName_)
-    UName_(pivpvf.UName_)
 {}
 
 
@@ -108,8 +104,7 @@ pressureInletOutletParSlipVelocityFvPatchVectorField
 :
     mixedFvPatchVectorField(pivpvf, iF),
     phiName_(pivpvf.phiName_),
-    rhoName_(pivpvf.rhoName_),
+    rhoName_(pivpvf.rhoName_)
-    UName_(pivpvf.UName_)
 {}
 
 
@@ -122,8 +117,6 @@ void Foam::pressureInletOutletParSlipVelocityFvPatchVectorField::updateCoeffs()
         return;
     }
 
-    const label patchI = patch().index();
-
     const surfaceScalarField& phi =
         db().lookupObject<surfaceScalarField>(phiName_);
 
@@ -131,23 +124,22 @@ void Foam::pressureInletOutletParSlipVelocityFvPatchVectorField::updateCoeffs()
         patch().patchField<surfaceScalarField, scalar>(phi);
 
     tmp<vectorField> n = patch().nf();
-    const Field<scalar>& magS = patch().magSf();
+    const Field<scalar>& magSf = patch().magSf();
 
-    const volVectorField& U = db().lookupObject<volVectorField>(UName_);
+    // Get the tangential component from the internalField (zero-gradient)
-
+    vectorField Ut(patchInternalField());
-    vectorField Uc(U.boundaryField()[patchI].patchInternalField());
+    Ut -= n()*(Ut & n());
-    Uc -= n()*(Uc & n());
 
     if (phi.dimensions() == dimVelocity*dimArea)
     {
-        refValue() = Uc + n*phip/magS;
+        refValue() = Ut + n*phip/magSf;
     }
     else if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
     {
         const fvPatchField<scalar>& rhop =
             patch().lookupPatchField<volScalarField, scalar>(rhoName_);
 
-        refValue() = Uc + n*phip/(rhop*magS);
+        refValue() = Ut + n*phip/(rhop*magSf);
     }
     else
     {
@@ -176,7 +168,6 @@ void Foam::pressureInletOutletParSlipVelocityFvPatchVectorField::write
     fvPatchVectorField::write(os);
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
-    writeEntryIfDifferent<word>(os, "U", "U", UName_);
     writeEntry("value", os);
 }
 
@@ -188,11 +179,7 @@ void Foam::pressureInletOutletParSlipVelocityFvPatchVectorField::operator=
     const fvPatchField<vector>& pvf
 )
 {
-    fvPatchField<vector>::operator=
+    fvPatchField<vector>::operator=(pvf);
-    (
-        valueFraction()*(patch().nf()*(patch().nf() & pvf))
-      + (1 - valueFraction())*pvf
-    );
 }
 
 

src/finiteVolume/fields/fvPatchFields/derived/pressureInletOutletParSlipVelocity/pressureInletOutletParSlipVelocityFvPatchVectorField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,7 +41,6 @@ Description
         Property     | Description             | Required    | Default value
         phi          | flux field name         | no          | phi
         rho          | density field name      | no          | rho
-        U            | velocity field name     | no          | U
     \endtable
 
     Example of the boundary condition specification:
@@ -95,9 +94,6 @@ class pressureInletOutletParSlipVelocityFvPatchVectorField
         //- Density field name
         word rhoName_;
 
-        //- Velocity field name
-        word UName_;
-
 
 public:
 

src/finiteVolume/fields/fvPatchFields/derived/pressureInletOutletVelocity/pressureInletOutletVelocityFvPatchVectorField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,10 +25,10 @@ License
 
 #include "pressureInletOutletVelocityFvPatchVectorField.H"
 #include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
 #include "volFields.H"
 #include "surfaceFields.H"
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::pressureInletOutletVelocityFvPatchVectorField::
@@ -38,12 +38,14 @@ pressureInletOutletVelocityFvPatchVectorField
     const DimensionedField<vector, volMesh>& iF
 )
 :
-    directionMixedFvPatchVectorField(p, iF),
+    mixedFvPatchVectorField(p, iF),
-    phiName_("phi")
+    phiName_("phi"),
+    rhoName_("rho"),
+    applyTangentialVelocity_(false)
 {
-    refValue() = vector::zero;
+    refValue() = *this;
     refGrad() = vector::zero;
-    valueFraction() = symmTensor::zero;
+    valueFraction() = 0.0;
 }
 
 
@@ -56,10 +58,12 @@ pressureInletOutletVelocityFvPatchVectorField
     const fvPatchFieldMapper& mapper
 )
 :
-    directionMixedFvPatchVectorField(ptf, p, iF, mapper),
+    mixedFvPatchVectorField(ptf, p, iF, mapper),
-    phiName_(ptf.phiName_)
+    phiName_(ptf.phiName_),
+    rhoName_(ptf.rhoName_),
+    applyTangentialVelocity_(ptf.applyTangentialVelocity_)
 {
-    if (ptf.tangentialVelocity_.size())
+    if (applyTangentialVelocity_)
     {
         tangentialVelocity_ = mapper(ptf.tangentialVelocity_);
     }
@@ -74,25 +78,26 @@ pressureInletOutletVelocityFvPatchVectorField
     const dictionary& dict
 )
 :
-    directionMixedFvPatchVectorField(p, iF),
+    mixedFvPatchVectorField(p, iF),
-    phiName_(dict.lookupOrDefault<word>("phi", "phi"))
+    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
+    applyTangentialVelocity_(false)
 {
     fvPatchVectorField::operator=(vectorField("value", dict, p.size()));
 
     if (dict.found("tangentialVelocity"))
     {
+        applyTangentialVelocity_ = true;
+
         setTangentialVelocity
         (
             vectorField("tangentialVelocity", dict, p.size())
         );
     }
-    else
-    {
-        refValue() = vector::zero;
-    }
 
+    refValue() = *this;
     refGrad() = vector::zero;
-    valueFraction() = symmTensor::zero;
+    valueFraction() = 0.0;
 }
 
 
@@ -102,9 +107,11 @@ pressureInletOutletVelocityFvPatchVectorField
     const pressureInletOutletVelocityFvPatchVectorField& pivpvf
 )
 :
-    directionMixedFvPatchVectorField(pivpvf),
+    mixedFvPatchVectorField(pivpvf),
     phiName_(pivpvf.phiName_),
-    tangentialVelocity_(pivpvf.tangentialVelocity_)
+    rhoName_(pivpvf.rhoName_),
+    tangentialVelocity_(pivpvf.tangentialVelocity_),
+    applyTangentialVelocity_(pivpvf.applyTangentialVelocity_)
 {}
 
 
@@ -115,9 +122,11 @@ pressureInletOutletVelocityFvPatchVectorField
     const DimensionedField<vector, volMesh>& iF
 )
 :
-    directionMixedFvPatchVectorField(pivpvf, iF),
+    mixedFvPatchVectorField(pivpvf, iF),
     phiName_(pivpvf.phiName_),
-    tangentialVelocity_(pivpvf.tangentialVelocity_)
+    rhoName_(pivpvf.rhoName_),
+    tangentialVelocity_(pivpvf.tangentialVelocity_),
+    applyTangentialVelocity_(pivpvf.applyTangentialVelocity_)
 {}
 
 
@@ -126,9 +135,10 @@ pressureInletOutletVelocityFvPatchVectorField
 void Foam::pressureInletOutletVelocityFvPatchVectorField::
 setTangentialVelocity(const vectorField& tangentialVelocity)
 {
+    applyTangentialVelocity_ = true;
     tangentialVelocity_ = tangentialVelocity;
-    const vectorField n(patch().nf());
+    vectorField n(patch().nf());
-    refValue() = tangentialVelocity_ - n*(n & tangentialVelocity_);
+    tangentialVelocity_ -= n*(n & tangentialVelocity_);
 }
 
 
@@ -137,8 +147,8 @@ void Foam::pressureInletOutletVelocityFvPatchVectorField::autoMap
     const fvPatchFieldMapper& m
 )
 {
-    directionMixedFvPatchVectorField::autoMap(m);
+    mixedFvPatchVectorField::autoMap(m);
-    if (tangentialVelocity_.size())
+    if (applyTangentialVelocity_)
     {
         tangentialVelocity_.autoMap(m);
     }
@@ -151,9 +161,9 @@ void Foam::pressureInletOutletVelocityFvPatchVectorField::rmap
     const labelList& addr
 )
 {
-    directionMixedFvPatchVectorField::rmap(ptf, addr);
+    mixedFvPatchVectorField::rmap(ptf, addr);
 
-    if (tangentialVelocity_.size())
+    if (applyTangentialVelocity_)
     {
         const pressureInletOutletVelocityFvPatchVectorField& tiptf =
             refCast<const pressureInletOutletVelocityFvPatchVectorField>(ptf);
@@ -170,25 +180,67 @@ void Foam::pressureInletOutletVelocityFvPatchVectorField::updateCoeffs()
         return;
     }
 
+    const surfaceScalarField& phi =
+        db().lookupObject<surfaceScalarField>(phiName_);
+
     const fvsPatchField<scalar>& phip =
-        patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);
+        patch().patchField<surfaceScalarField, scalar>(phi);
 
-    valueFraction() = neg(phip)*(I - sqr(patch().nf()));
+    vectorField n(patch().nf());
+    const Field<scalar>& magSf = patch().magSf();
 
-    directionMixedFvPatchVectorField::updateCoeffs();
+    if (phi.dimensions() == dimVelocity*dimArea)
-    directionMixedFvPatchVectorField::evaluate();
+    {
+        refValue() = n*phip/magSf;
+    }
+    else if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
+    {
+        const fvPatchField<scalar>& rhop =
+            patch().lookupPatchField<volScalarField, scalar>(rhoName_);
+
+        refValue() = n*phip/(rhop*magSf);
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "pressureInletOutletVelocityFvPatchVectorField::"
+            "updateCoeffs()"
+        )   << "dimensions of phi are not correct" << nl
+            << "    on patch " << this->patch().name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalError);
+    }
+
+    if (applyTangentialVelocity_)
+    {
+        // Adjust the tangential velocity to conserve kinetic energy
+        // of the entrained fluid
+        // scalarField magSqrUt(magSqr(tangentialVelocity_));
+        // scalarField scale
+        // (
+        //     sqrt(max(magSqrUt - magSqr(refValue()), scalar(0))/magSqrUt)
+        // );
+        // refValue() += scale*tangentialVelocity_;
+        refValue() += tangentialVelocity_;
+    }
+
+    valueFraction() = 1.0 - pos(phip);
+
+    mixedFvPatchVectorField::updateCoeffs();
 }
 
 
 void Foam::pressureInletOutletVelocityFvPatchVectorField::write
 (
     Ostream& os
-)
+) const
-const
 {
     fvPatchVectorField::write(os);
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
-    if (tangentialVelocity_.size())
+    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    if (applyTangentialVelocity_)
     {
         tangentialVelocity_.writeEntry("tangentialVelocity", os);
     }
@@ -203,9 +255,7 @@ void Foam::pressureInletOutletVelocityFvPatchVectorField::operator=
     const fvPatchField<vector>& pvf
 )
 {
-    tmp<vectorField> normalValue = transform(valueFraction(), refValue());
+    fvPatchField<vector>::operator=(pvf);
-    tmp<vectorField> transformGradValue = transform(I - valueFraction(), pvf);
-    fvPatchField<vector>::operator=(normalValue + transformGradValue);
 }
 
 

src/finiteVolume/fields/fvPatchFields/derived/pressureInletOutletVelocity/pressureInletOutletVelocityFvPatchVectorField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,21 +25,21 @@ Class
     Foam::pressureInletOutletVelocityFvPatchVectorField
 
 Group
-    grpInletletBoundaryConditions grpOutletBoundaryConditions
+    grpInletBoundaryConditions grpOutletBoundaryConditions
 
 Description
     This velocity inlet/outlet boundary condition is applied to pressure
     boundaries where the pressure is specified.  A zero-gradient condition is
     applied for outflow (as defined by the flux); for inflow, the velocity is
-    obtained from the patch-face normal component of the internal-cell value.
+    obtained from the patch-face normal component of the internal-cell value and
-
+    the tangential patch velocity can be optionally specified.
-    The tangential patch velocity can be optionally specified.
 
     \heading Patch usage
 
     \table
         Property     | Description             | Required    | Default value
         phi          | flux field name         | no          | phi
+        rho          | density field name      | no          | rho
         tangentialVelocity | tangential velocity field | no  |
     \endtable
 
@@ -49,6 +49,7 @@ Description
     {
         type            pressureInletOutletVelocity;
         phi             phi;
+        rho             rho;
         tangentialVelocity uniform (0 0 0);
         value           uniform 0;
     }
@@ -69,7 +70,7 @@ SourceFiles
 #define pressureInletOutletVelocityFvPatchVectorField_H
 
 #include "fvPatchFields.H"
-#include "directionMixedFvPatchFields.H"
+#include "mixedFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -77,21 +78,29 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-       Class pressureInletOutletVelocityFvPatchVectorField Declaration
+    Class pressureInletOutletVelocityFvPatchVectorField Declaration
 \*---------------------------------------------------------------------------*/
 
 class pressureInletOutletVelocityFvPatchVectorField
 :
-    public directionMixedFvPatchVectorField
+    public mixedFvPatchVectorField
 {
     // Private data
 
         //- Flux field name
         word phiName_;
 
+        //- Density field name
+        word rhoName_;
+
+
+protected:
+
         //- Optional tangential velocity component
         vectorField tangentialVelocity_;
 
+        bool applyTangentialVelocity_;
+
 
 public:
 
@@ -108,16 +117,9 @@ public:
             const DimensionedField<vector, volMesh>&
         );
 
-        //- Construct from patch, internal field and dictionary
-        pressureInletOutletVelocityFvPatchVectorField
-        (
-            const fvPatch&,
-            const DimensionedField<vector, volMesh>&,
-            const dictionary&
-        );
-
         //- Construct by mapping given
-        //  pressureInletOutletVelocityFvPatchVectorField onto a new patch
+        //  pressureInletOutletVelocityFvPatchVectorField
+        //  onto a new patch
         pressureInletOutletVelocityFvPatchVectorField
         (
             const pressureInletOutletVelocityFvPatchVectorField&,
@@ -126,6 +128,14 @@ public:
             const fvPatchFieldMapper&
         );
 
+        //- Construct from patch, internal field and dictionary
+        pressureInletOutletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
         //- Construct as copy
         pressureInletOutletVelocityFvPatchVectorField
         (
@@ -137,7 +147,10 @@ public:
         {
             return tmp<fvPatchVectorField>
             (
-                new pressureInletOutletVelocityFvPatchVectorField(*this)
+                new pressureInletOutletVelocityFvPatchVectorField
+                (
+                    *this
+                )
             );
         }
 
@@ -156,7 +169,11 @@ public:
         {
             return tmp<fvPatchVectorField>
             (
-                new pressureInletOutletVelocityFvPatchVectorField(*this, iF)
+                new pressureInletOutletVelocityFvPatchVectorField
+                (
+                    *this,
+                    iF
+                )
             );
         }
 
@@ -177,6 +194,18 @@ public:
                 return phiName_;
             }
 
+            //- Return the name of rho
+            const word& rhoName() const
+            {
+                return rhoName_;
+            }
+
+            //- Return reference to the name of rho to allow adjustment
+            word& rhoName()
+            {
+                return rhoName_;
+            }
+
             //- Return the tangential velocity
             const vectorField& tangentialVelocity() const
             {

src/finiteVolume/fields/fvPatchFields/derived/prghPressure/prghPressureFvPatchScalarField.C

@@ -0,0 +1,185 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "prghPressureFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "uniformDimensionedFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::prghPressureFvPatchScalarField::
+prghPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    rhoName_("rho"),
+    p_(p.size(), 0.0)
+{}
+
+
+Foam::prghPressureFvPatchScalarField::
+prghPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    rhoName_(dict.lookupOrDefault<word>("rhoName", "rho")),
+    p_("p", dict, p.size())
+{
+    if (dict.found("value"))
+    {
+        fvPatchScalarField::operator=
+        (
+            scalarField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        fvPatchField<scalar>::operator=(p_);
+    }
+}
+
+
+Foam::prghPressureFvPatchScalarField::
+prghPressureFvPatchScalarField
+(
+    const prghPressureFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    rhoName_(ptf.rhoName_),
+    p_(ptf.p_, mapper)
+{}
+
+
+Foam::prghPressureFvPatchScalarField::
+prghPressureFvPatchScalarField
+(
+    const prghPressureFvPatchScalarField& ptf
+)
+:
+    fixedValueFvPatchScalarField(ptf),
+    rhoName_(ptf.rhoName_),
+    p_(ptf.p_)
+{}
+
+
+Foam::prghPressureFvPatchScalarField::
+prghPressureFvPatchScalarField
+(
+    const prghPressureFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(ptf, iF),
+    rhoName_(ptf.rhoName_),
+    p_(ptf.p_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::prghPressureFvPatchScalarField::autoMap
+(
+    const fvPatchFieldMapper& m
+)
+{
+    fixedValueFvPatchScalarField::autoMap(m);
+    p_.autoMap(m);
+}
+
+
+void Foam::prghPressureFvPatchScalarField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList& addr
+)
+{
+    fixedValueFvPatchScalarField::rmap(ptf, addr);
+
+    const prghPressureFvPatchScalarField& tiptf =
+        refCast<const prghPressureFvPatchScalarField>(ptf);
+
+    p_.rmap(tiptf.p_, addr);
+}
+
+
+void Foam::prghPressureFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const scalarField& rhop = patch().lookupPatchField<volScalarField, scalar>
+    (
+        rhoName_
+    );
+
+    const uniformDimensionedVectorField& g =
+        db().lookupObject<uniformDimensionedVectorField>("g");
+
+    operator==(p_ - rhop*((g.value() & patch().Cf())));
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::prghPressureFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    if (rhoName_ != "rho")
+    {
+        os.writeKeyword("rhoName")
+            << rhoName_ << token::END_STATEMENT << nl;
+    }
+    p_.writeEntry("p", os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        prghPressureFvPatchScalarField
+    );
+}
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/prghPressure/prghPressureFvPatchScalarField.H

@@ -0,0 +1,236 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::prghPressureFvPatchScalarField
+
+Group
+    grpGenericBoundaryConditions
+
+Description
+    This boundary condition provides static pressure condition for p_rgh,
+    calculated as:
+
+        \f[
+            p_rgh = p - \rho g h
+        \f]
+
+    where
+    \vartable
+        p_rgh   | Pseudo hydrostatic pressure [Pa]
+        p       | Static pressure [Pa]
+        h       | Height in the opposite direction to gravity
+        \rho    | density
+        g       | acceleration due to gravity [m/s2]
+    \endtable
+
+    \heading Patch usage
+
+    \table
+        Property     | Description             | Required    | Default value
+        rhoName      | rho field name          | no          | rho
+        p            | static pressure         | yes         |
+    \endtable
+
+    Example of the boundary condition specification:
+    \verbatim
+    myPatch
+    {
+        type            prghPressure;
+        rhoName         rho;
+        p               uniform 0;
+        value           uniform 0; // optional initial value
+    }
+    \endverbatim
+
+SeeAlso
+    Foam::fixedValueFvPatchScalarField
+
+SourceFiles
+    prghPressureFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef prghPressureFvPatchScalarField_H
+#define prghPressureFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+         Class prghPressureFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class prghPressureFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+
+protected:
+
+    // Protected data
+
+        //- Name of phase-fraction field
+        word rhoName_;
+
+        //- Static pressure
+        scalarField p_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("prghPressure");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        prghPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        prghPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  prghPressureFvPatchScalarField onto a new patch
+        prghPressureFvPatchScalarField
+        (
+            const prghPressureFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        prghPressureFvPatchScalarField
+        (
+            const prghPressureFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField >
+            (
+                new prghPressureFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        prghPressureFvPatchScalarField
+        (
+            const prghPressureFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new prghPressureFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Access
+
+            //- Return the rhoName
+            const word& rhoName() const
+            {
+                return rhoName_;
+            }
+
+            //- Return reference to the rhoName to allow adjustment
+            word& rhoName()
+            {
+                return rhoName_;
+            }
+
+            //- Return the static pressure
+            const scalarField& p() const
+            {
+                return p_;
+            }
+
+            //- Return reference to the static pressure to allow adjustment
+            scalarField& p()
+            {
+                return p_;
+            }
+
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
+
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/rotatingPressureInletOutletVelocity/rotatingPressureInletOutletVelocityFvPatchVectorField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,17 +33,17 @@ License
 void Foam::rotatingPressureInletOutletVelocityFvPatchVectorField::
 calcTangentialVelocity()
 {
+    applyTangentialVelocity_ = true;
+
     const scalar t = this->db().time().timeOutputValue();
     vector om = omega_->value(t);
 
     vector axisHat = om/mag(om);
-    const vectorField tangentialVelocity
+    tangentialVelocity_ =
-    (
+        (-om) ^ (patch().Cf() - axisHat*(axisHat & patch().Cf()));
-        (-om) ^ (patch().Cf() - axisHat*(axisHat & patch().Cf()))
-    );
 
     const vectorField n(patch().nf());
-    refValue() = tangentialVelocity - n*(n & tangentialVelocity);
+    tangentialVelocity_ -= n*(n & tangentialVelocity_);
 }
 
 
@@ -72,9 +72,7 @@ rotatingPressureInletOutletVelocityFvPatchVectorField
 :
     pressureInletOutletVelocityFvPatchVectorField(ptf, p, iF, mapper),
     omega_(ptf.omega_().clone().ptr())
-{
+{}
-    calcTangentialVelocity();
-}
 
 
 Foam::rotatingPressureInletOutletVelocityFvPatchVectorField::
@@ -87,9 +85,7 @@ rotatingPressureInletOutletVelocityFvPatchVectorField
 :
     pressureInletOutletVelocityFvPatchVectorField(p, iF, dict),
     omega_(DataEntry<vector>::New("omega", dict))
-{
+{}
-    calcTangentialVelocity();
-}
 
 
 Foam::rotatingPressureInletOutletVelocityFvPatchVectorField::
@@ -100,9 +96,7 @@ rotatingPressureInletOutletVelocityFvPatchVectorField
 :
     pressureInletOutletVelocityFvPatchVectorField(rppvf),
     omega_(rppvf.omega_().clone().ptr())
-{
+{}
-    calcTangentialVelocity();
-}
 
 
 Foam::rotatingPressureInletOutletVelocityFvPatchVectorField::
@@ -114,13 +108,24 @@ rotatingPressureInletOutletVelocityFvPatchVectorField
 :
     pressureInletOutletVelocityFvPatchVectorField(rppvf, iF),
     omega_(rppvf.omega_().clone().ptr())
-{
+{}
-    calcTangentialVelocity();
-}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::rotatingPressureInletOutletVelocityFvPatchVectorField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    calcTangentialVelocity();
+
+    pressureInletOutletVelocityFvPatchVectorField::updateCoeffs();
+}
+
+
 void Foam::rotatingPressureInletOutletVelocityFvPatchVectorField::write
 (
     Ostream& os

src/finiteVolume/fields/fvPatchFields/derived/rotatingPressureInletOutletVelocity/rotatingPressureInletOutletVelocityFvPatchVectorField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -177,6 +177,9 @@ public:
 
     // Member functions
 
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
         //- Write
         virtual void write(Ostream&) const;
 };

src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.H

@@ -33,47 +33,47 @@ Description
 
     1. incompressible subsonic:
         \f[
-            p_T = p_0 + 0.5 |U|^2
+            p_p = p_0 - 0.5 |U|^2
         \f]
         where
         \vartable
-            p_T     | incompressible total pressure [m2/s2]
+            p_p     | incompressible pressure at patch [m2/s2]
-            p_0     | incompressible reference pressure [m2/s2]
+            p_0     | incompressible total pressure [m2/s2]
             U       | velocity
         \endvartable
 
     2. compressible subsonic:
         \f[
-            p_T = p_0 + 0.5 \rho |U|^2
+            p_p = p_0 - 0.5 \rho |U|^2
         \f]
         where
         \vartable
-            p_T     | total pressure [Pa]
+            p_p     | pressure at patch [Pa]
-            p_0     | reference pressure [Pa]
+            p_0     | total pressure [Pa]
             \rho    | density [kg/m3]
             U       | velocity
         \endvartable
 
     3. compressible transonic (\gamma <= 1):
         \f[
-            p_T = \frac{p_0}{1 + 0.5 \psi |U|^2}
+            p_p = \frac{p_0}{1 + 0.5 \psi |U|^2}
         \f]
         where
         \vartable
-            p_T     | total pressure [Pa]
+            p_p     | pressure at patch [Pa]
-            p_0     | reference pressure [Pa]
+            p_0     | total pressure [Pa]
             G       | coefficient given by \f$\frac{\gamma}{1-\gamma}\f$
         \endvartable
 
     4. compressible supersonic (\gamma > 1):
         \f[
-            p_T = \frac{p_0}{(1 + 0.5 \psi G |U|^2)^{\frac{1}{G}}}
+            p_p = \frac{p_0}{(1 + 0.5 \psi G |U|^2)^{\frac{1}{G}}}
         \f]
         where
         \vartable
+            p_p     | pressure at patch [Pa]
+            p_0     | total pressure [Pa]
             \gamma  | ratio of specific heats (Cp/Cv)
-            p_T     | total pressure [Pa]
-            p_0     | reference pressure [Pa]
             \psi    | compressibility [m2/s2]
             G       | coefficient given by \f$\frac{\gamma}{1-\gamma}\f$
         \endvartable
@@ -98,7 +98,7 @@ Description
         rho          | density field name      | no          | none
         psi          | compressibility field name | no       | none
         gamma        | ratio of specific heats (Cp/Cv) | yes |
-        p0           | static pressure reference | yes       |
+        p0           | total pressure          | yes       |
     \endtable
 
     Example of the boundary condition specification:

src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C

@@ -25,6 +25,9 @@ License
 
 #include "fvcReconstruct.H"
 #include "fvMesh.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "fvcSurfaceIntegrate.H"
 #include "zeroGradientFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/finiteVolume/finiteVolume/fvc/fvcReconstruct.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,6 +69,9 @@ namespace fvc
     (
         const tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >&
     );
+
+    tmp<volScalarField> reconstructMag(const surfaceScalarField&);
+    tmp<volScalarField> reconstructMag(const tmp<surfaceScalarField>&);
 }
 
 

src/finiteVolume/finiteVolume/fvc/fvcReconstructMag.C

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvcReconstruct.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "zeroGradientFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fvc
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+tmp<volScalarField> reconstructMag(const surfaceScalarField& ssf)
+{
+    const fvMesh& mesh = ssf.mesh();
+
+    const labelUList& owner = mesh.owner();
+    const labelUList& neighbour = mesh.neighbour();
+
+    const volVectorField& C = mesh.C();
+    const surfaceVectorField& Cf = mesh.Cf();
+    const surfaceVectorField& Sf = mesh.Sf();
+    const surfaceScalarField& magSf = mesh.magSf();
+
+    tmp<volScalarField> treconField
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "reconstruct("+ssf.name()+')',
+                ssf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar
+            (
+                "0",
+                ssf.dimensions()/dimArea,
+                scalar(0)
+            ),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    scalarField& rf = treconField();
+
+    forAll(owner, facei)
+    {
+        label own = owner[facei];
+        label nei = neighbour[facei];
+
+        rf[own] += (Sf[facei] & (Cf[facei] - C[own]))*ssf[facei]/magSf[facei];
+        rf[nei] -= (Sf[facei] & (Cf[facei] - C[nei]))*ssf[facei]/magSf[facei];
+    }
+
+    const surfaceScalarField::GeometricBoundaryField& bsf = ssf.boundaryField();
+
+    forAll(bsf, patchi)
+    {
+        const fvsPatchScalarField& psf = bsf[patchi];
+
+        const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
+        const vectorField& pCf = Cf.boundaryField()[patchi];
+        const vectorField& pSf = Sf.boundaryField()[patchi];
+        const scalarField& pMagSf = magSf.boundaryField()[patchi];
+
+        forAll(pOwner, pFacei)
+        {
+            label own = pOwner[pFacei];
+            rf[own] +=
+                (pSf[pFacei] & (pCf[pFacei] - C[own]))
+               *psf[pFacei]/pMagSf[pFacei];
+        }
+    }
+
+    rf /= mesh.V();
+
+    treconField().correctBoundaryConditions();
+
+    return treconField;
+}
+
+
+tmp<volScalarField> reconstructMag(const tmp<surfaceScalarField>& tssf)
+{
+    tmp<volScalarField> tvf
+    (
+        fvc::reconstructMag(tssf())
+    );
+    tssf.clear();
+    return tvf;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fvc
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -30,6 +30,7 @@ License
 void Foam::MULES::correct
 (
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiPsiCorr,
     const scalar psiMax,
     const scalar psiMin
@@ -39,6 +40,7 @@ void Foam::MULES::correct
     (
         geometricOneField(),
         psi,
+        phi,
         phiPsiCorr,
         zeroField(), zeroField(),
         psiMax, psiMin
@@ -49,6 +51,7 @@ void Foam::MULES::correct
 void Foam::MULES::LTScorrect
 (
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiPsiCorr,
     const scalar psiMax,
     const scalar psiMin
@@ -58,6 +61,7 @@ void Foam::MULES::LTScorrect
     (
         geometricOneField(),
         psi,
+        phi,
         phiPsiCorr,
         zeroField(), zeroField(),
         psiMax, psiMin

src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H

@@ -66,6 +66,7 @@ void correct
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su
@@ -76,6 +77,7 @@ void correct
 (
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -86,6 +88,7 @@ void correct
 void correct
 (
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const scalar psiMax,
     const scalar psiMin
@@ -96,6 +99,7 @@ void LTScorrect
 (
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -106,6 +110,7 @@ void LTScorrect
 void LTScorrect
 (
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const scalar psiMax,
     const scalar psiMin
@@ -119,6 +124,7 @@ void limiterCorr
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
+    const surfaceScalarField& phi,
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -133,6 +139,7 @@ void limitCorr
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,

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

@@ -37,6 +37,7 @@ void Foam::MULES::correct
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su
@@ -77,6 +78,7 @@ void Foam::MULES::correct
 (
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -94,8 +96,18 @@ void Foam::MULES::correct
         readLabel(MULEScontrols.lookup("nLimiterIter"))
     );
 
-    limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
+    limitCorr
-    correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
+    (
+        rDeltaT,
+        rho,
+        psi,
+        phi,
+        phiCorr,
+        Sp, Su,
+        psiMax, psiMin,
+        nLimiterIter
+    );
+    correct(rDeltaT, rho, psi, phi, phiCorr, Sp, Su);
 }
 
 
@@ -104,6 +116,7 @@ void Foam::MULES::LTScorrect
 (
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -123,8 +136,18 @@ void Foam::MULES::LTScorrect
         readLabel(MULEScontrols.lookup("nLimiterIter"))
     );
 
-    limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
+    limitCorr
-    correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
+    (
+        rDeltaT,
+        rho,
+        psi,
+        phi,
+        phiCorr,
+        Sp, Su,
+        psiMax, psiMin,
+        nLimiterIter
+    );
+    correct(rDeltaT, rho, psi, phi, phiCorr, Sp, Su);
 }
 
 
@@ -135,6 +158,7 @@ void Foam::MULES::limiterCorr
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
+    const surfaceScalarField& phi,
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -153,6 +177,9 @@ void Foam::MULES::limiterCorr
     tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
     const scalarField& V = tVsc();
 
+    const surfaceScalarField::GeometricBoundaryField& phiBf =
+        phi.boundaryField();
+
     const scalarField& phiCorrIf = phiCorr;
     const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
         phiCorr.boundaryField();
@@ -408,12 +435,12 @@ void Foam::MULES::limiterCorr
             {
                 const labelList& pFaceCells =
                     mesh.boundary()[patchi].faceCells();
-                const scalarField& phiCorrPf = phiCorrBf[patchi];
+                const scalarField& phiPf = phiBf[patchi];
 
                 forAll(lambdaPf, pFacei)
                 {
                     // Limit outlet faces only
-                    if (phiCorrPf[pFacei] > SMALL*SMALL)
+                    if (phiPf[pFacei] > SMALL*SMALL)
                     {
                         label pfCelli = pFaceCells[pFacei];
 
@@ -443,6 +470,7 @@ void Foam::MULES::limitCorr
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -478,6 +506,7 @@ void Foam::MULES::limitCorr
         rDeltaT,
         rho,
         psi,
+        phi,
         phiCorr,
         Sp,
         Su,

src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C

@@ -173,10 +173,55 @@ void Foam::fv::interRegionExplicitPorositySource::addSup
 
     if (eqn.dimensions() == dimForce)
     {
-        const volScalarField& rhoNbr =
+        volScalarField rhoNbr
-            nbrMesh.lookupObject<volScalarField>(rhoName_);
+        (
-        const volScalarField& muNbr =
+            IOobject
-            nbrMesh.lookupObject<volScalarField>(muName_);
+            (
+                "rho:UNbr",
+                nbrMesh.time().timeName(),
+                nbrMesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            nbrMesh,
+            dimensionedScalar("zero", dimDensity, 0.0)
+        );
+
+        volScalarField muNbr
+        (
+            IOobject
+            (
+                "mu:UNbr",
+                nbrMesh.time().timeName(),
+                nbrMesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            nbrMesh,
+            dimensionedScalar("zero", dimViscosity, 0.0)
+        );
+
+        const volScalarField& rho =
+            mesh_.lookupObject<volScalarField>(rhoName_);
+
+        const volScalarField& mu =
+            mesh_.lookupObject<volScalarField>(muName_);
+
+        // map local rho onto neighbour region
+        meshInterp().mapSrcToTgt
+        (
+            rho.internalField(),
+            plusEqOp<scalar>(),
+            rhoNbr.internalField()
+        );
+
+        // map local mu onto neighbour region
+        meshInterp().mapSrcToTgt
+        (
+            mu.internalField(),
+            plusEqOp<scalar>(),
+            muNbr.internalField()
+        );
 
         porosityPtr_->addResistance(nbrEqn, rhoNbr, muNbr);
     }

src/lagrangian/basic/particle/particleTemplates.C

@@ -424,7 +424,7 @@ Foam::scalar Foam::particle::trackToFace
                 (
                     position_,
                     endPosition,
-                    triI,
+                    tI,
                     tetAreas[tI],
                     tetPlaneBasePtIs[tI],
                     cellI_,

src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C

@@ -90,16 +90,16 @@ void Foam::cyclicACMIPolyPatch::resetAMI
             AMIPatchToPatchInterpolation::imPartialFaceAreaWeight
         );
 
-        const scalarField& srcWeightSum = AMI().srcWeightsSum();
+        srcMask_ =
+            min(1.0 - tolerance_, max(tolerance_, AMI().srcWeightsSum()));
+
+        tgtMask_ =
+            min(1.0 - tolerance_, max(tolerance_, AMI().tgtWeightsSum()));
 
-        // set patch face areas based on sum of AMI weights per face
         forAll(Sf, faceI)
         {
-            scalar w = srcWeightSum[faceI];
+            Sf[faceI] *= srcMask_[faceI];
-            w = min(1.0 - tolerance_, max(tolerance_, w));
+            noSf[faceI] *= 1.0 - srcMask_[faceI];
-
-            Sf[faceI] *= w;
-            noSf[faceI] *= 1.0 - w;
         }
 
         setNeighbourFaceAreas();
@@ -116,8 +116,6 @@ void Foam::cyclicACMIPolyPatch::setNeighbourFaceAreas() const
         refCast<const cyclicACMIPolyPatch>(this->neighbPatch());
     const polyPatch& pp = cp.nonOverlapPatch();
 
-    const scalarField& tgtWeightSum = AMI().tgtWeightsSum();
-
     const vectorField& faceAreas0 = cp.faceAreas0();
 
     vectorField::subField Sf = cp.faceAreas();
@@ -125,11 +123,8 @@ void Foam::cyclicACMIPolyPatch::setNeighbourFaceAreas() const
 
     forAll(Sf, faceI)
     {
-        scalar w = tgtWeightSum[faceI];
+        Sf[faceI] = tgtMask_[faceI]*faceAreas0[faceI];
-        w = min(1.0 - tolerance_, max(tolerance_, w));
+        noSf[faceI] = (1.0 - tgtMask_[faceI])*faceAreas0[faceI];
-
-        Sf[faceI] = w*faceAreas0[faceI];
-        noSf[faceI] = (1.0 - w)*faceAreas0[faceI];
     }
 }
 
@@ -184,6 +179,18 @@ void Foam::cyclicACMIPolyPatch::clearGeom()
 }
 
 
+const Foam::scalarField& Foam::cyclicACMIPolyPatch::srcMask() const
+{
+    return srcMask_;
+}
+
+
+const Foam::scalarField& Foam::cyclicACMIPolyPatch::tgtMask() const
+{
+    return tgtMask_;
+}
+
+
 // * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * //
 
 Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
@@ -201,6 +208,8 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     faceAreas0_(),
     nonOverlapPatchName_(word::null),
     nonOverlapPatchID_(-1),
+    srcMask_(),
+    tgtMask_(),
     updated_(false)
 {
     // Non-overlapping patch might not be valid yet so cannot determine
@@ -221,6 +230,8 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     faceAreas0_(),
     nonOverlapPatchName_(dict.lookup("nonOverlapPatch")),
     nonOverlapPatchID_(-1),
+    srcMask_(),
+    tgtMask_(),
     updated_(false)
 {
     if (nonOverlapPatchName_ == name)
@@ -256,6 +267,8 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     faceAreas0_(),
     nonOverlapPatchName_(pp.nonOverlapPatchName_),
     nonOverlapPatchID_(-1),
+    srcMask_(),
+    tgtMask_(),
     updated_(false)
 {
     // Non-overlapping patch might not be valid yet so cannot determine
@@ -278,6 +291,8 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     faceAreas0_(),
     nonOverlapPatchName_(nonOverlapPatchName),
     nonOverlapPatchID_(-1),
+    srcMask_(),
+    tgtMask_(),
     updated_(false)
 {
     if (nonOverlapPatchName_ == name())
@@ -314,6 +329,8 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
     faceAreas0_(),
     nonOverlapPatchName_(pp.nonOverlapPatchName_),
     nonOverlapPatchID_(-1),
+    srcMask_(),
+    tgtMask_(),
     updated_(false)
 {}
 
@@ -326,6 +343,13 @@ Foam::cyclicACMIPolyPatch::~cyclicACMIPolyPatch()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+const Foam::cyclicACMIPolyPatch& Foam::cyclicACMIPolyPatch::neighbPatch() const
+{
+    const polyPatch& pp = this->boundaryMesh()[neighbPatchID()];
+    return refCast<const cyclicACMIPolyPatch>(pp);
+}
+
+
 Foam::label Foam::cyclicACMIPolyPatch::nonOverlapPatchID() const
 {
     if (nonOverlapPatchID_ == -1)

src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.H

@@ -66,6 +66,12 @@ private:
         //- Index of non-overlapping patch
         mutable label nonOverlapPatchID_;
 
+        //- Mask/weighting for source patch
+        mutable scalarField srcMask_;
+
+        //- Mask/weighting for target patch
+        mutable scalarField tgtMask_;
+
         //- Flag to indicate that AMI has been updated
         mutable bool updated_;
 
@@ -111,6 +117,12 @@ protected:
         //- Clear geometry
         virtual void clearGeom();
 
+        //- Return the mask/weighting for the source patch
+        virtual const scalarField& srcMask() const;
+
+        //- Return the mask/weighting for the target patch
+        virtual const scalarField& tgtMask() const;
+
 
 public:
 
@@ -245,6 +257,9 @@ public:
             //- Return access to the original patch face areas
             inline const vectorField& faceAreas0() const;
 
+            //- Return a reference to the neighbour patch
+            virtual const cyclicACMIPolyPatch& neighbPatch() const;
+
             //- Non-overlapping patch name
             inline const word& nonOverlapPatchName() const;
 

src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchI.H

@@ -71,11 +71,11 @@ inline const Foam::scalarField& Foam::cyclicACMIPolyPatch::mask() const
 {
     if (owner())
     {
-        return AMI().srcWeightsSum();
+        return srcMask_;
     }
     else
     {
-        return neighbPatch().AMI().tgtWeightsSum();
+        return neighbPatch().tgtMask();
     }
 }
 

src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchTemplates.C

@@ -34,7 +34,7 @@ Foam::tmp<Foam::Field<Type> > Foam::cyclicACMIPolyPatch::interpolate
 {
     if (owner())
     {
-        const scalarField& w = AMI().srcWeightsSum();
+        const scalarField& w = srcMask_;
 
         return
             w*AMI().interpolateToSource(fldCouple)
@@ -42,7 +42,7 @@ Foam::tmp<Foam::Field<Type> > Foam::cyclicACMIPolyPatch::interpolate
     }
     else
     {
-        const scalarField& w = neighbPatch().AMI().tgtWeightsSum();
+        const scalarField& w = neighbPatch().tgtMask();
 
         return
             w*neighbPatch().AMI().interpolateToTarget(fldCouple)
@@ -73,14 +73,14 @@ void Foam::cyclicACMIPolyPatch::interpolate
 {
     if (owner())
     {
-        const scalarField& w = AMI().srcWeightsSum();
+        const scalarField& w = srcMask_;
 
         AMI().interpolateToSource(fldCouple, cop, result);
         result = w*result + (1.0 - w)*fldNonOverlap;
     }
     else
     {
-        const scalarField& w = neighbPatch().AMI().tgtWeightsSum();
+        const scalarField& w = neighbPatch().tgtMask();
 
         neighbPatch().AMI().interpolateToTarget(fldCouple, cop, result);
         result = w*result + (1.0 - w)*fldNonOverlap;

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C

@@ -642,6 +642,13 @@ bool Foam::cyclicAMIPolyPatch::owner() const
 }
 
 
+const Foam::cyclicAMIPolyPatch& Foam::cyclicAMIPolyPatch::neighbPatch() const
+{
+    const polyPatch& pp = this->boundaryMesh()[neighbPatchID()];
+    return refCast<const cyclicAMIPolyPatch>(pp);
+}
+
+
 const Foam::autoPtr<Foam::searchableSurface>&
 Foam::cyclicAMIPolyPatch::surfPtr() const
 {

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H

@@ -284,7 +284,7 @@ public:
             virtual bool owner() const;
 
             //- Return a reference to the neighbour patch
-            inline const cyclicAMIPolyPatch& neighbPatch() const;
+            virtual const cyclicAMIPolyPatch& neighbPatch() const;
 
             //- Return a reference to the projection surface
             const autoPtr<searchableSurface>& surfPtr() const;

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchI.H

@@ -38,14 +38,6 @@ inline const Foam::word& Foam::cyclicAMIPolyPatch::neighbPatchName() const
 }
 
 
-inline const Foam::cyclicAMIPolyPatch&
-Foam::cyclicAMIPolyPatch::neighbPatch() const
-{
-    const polyPatch& pp = this->boundaryMesh()[neighbPatchID()];
-    return refCast<const cyclicAMIPolyPatch>(pp);
-}
-
-
 inline const Foam::vector& Foam::cyclicAMIPolyPatch::rotationAxis() const
 {
     return rotationAxis_;

src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C

@@ -203,7 +203,7 @@ void Foam::mappedPatchBase::findSamples
     {
         case NEARESTCELL:
         {
-            if (samplePatch().size() && samplePatch() != "none")
+            if (samplePatch_.size() && samplePatch_ != "none")
             {
                 FatalErrorIn
                 (
@@ -495,7 +495,6 @@ void Foam::mappedPatchBase::findSamples
 void Foam::mappedPatchBase::calcMapping() const
 {
     static bool hasWarned = false;
-
     if (mapPtr_.valid())
     {
         FatalErrorIn("mappedPatchBase::calcMapping() const")
@@ -509,10 +508,8 @@ void Foam::mappedPatchBase::calcMapping() const
     // Get offsetted points
     const pointField offsettedPoints(samplePoints(patchPoints()));
 
-
+    // Do a sanity check - am I sampling my own patch?
-    // Do a sanity check
+    // This only makes sense for a non-zero offset.
-    // Am I sampling my own patch? This only makes sense for a non-zero
-    // offset.
     bool sampleMyself =
     (
         mode_ == NEARESTPATCHFACE
@@ -550,7 +547,6 @@ void Foam::mappedPatchBase::calcMapping() const
             << "offsetMode_:" << offsetModeNames_[offsetMode_] << endl;
     }
 
-
     // Get global list of all samples and the processor and face they come from.
     pointField samples;
     labelList patchFaceProcs;
@@ -565,7 +561,6 @@ void Foam::mappedPatchBase::calcMapping() const
         patchFc
     );
 
-
     // Find processor and cell/face samples are in and actual location.
     labelList sampleProcs;
     labelList sampleIndices;
@@ -641,7 +636,6 @@ void Foam::mappedPatchBase::calcMapping() const
         }
     }
 
-
     // Now we have all the data we need:
     // - where sample originates from (so destination when mapping):
     //   patchFaces, patchFaceProcs.
@@ -688,7 +682,6 @@ void Foam::mappedPatchBase::calcMapping() const
         }
     }
 
-
     // Determine schedule.
     mapPtr_.reset(new mapDistribute(sampleProcs, patchFaceProcs));
 
@@ -806,38 +799,34 @@ void Foam::mappedPatchBase::calcAMI() const
     }
 
     AMIPtr_.clear();
-/*
-    const polyPatch& nbr = samplePolyPatch();
-
-//    pointField nbrPoints(offsettedPoints());
-    pointField nbrPoints(nbr.localPoints());
 
     if (debug)
     {
+        const polyPatch& nbr = samplePolyPatch();
+
+        pointField nbrPoints(nbr.localPoints());
+
         OFstream os(patch_.name() + "_neighbourPatch-org.obj");
         meshTools::writeOBJ(os, samplePolyPatch().localFaces(), nbrPoints);
-    }
 
-    // transform neighbour patch to local system
+        // transform neighbour patch to local system
-    primitivePatch nbrPatch0
+        primitivePatch nbrPatch0
-    (
-        SubList<face>
         (
-            nbr.localFaces(),
+            SubList<face>
-            nbr.size()
+            (
-        ),
+                nbr.localFaces(),
-        nbrPoints
+                nbr.size()
-    );
+            ),
+            nbrPoints
+        );
 
-    if (debug)
-    {
         OFstream osN(patch_.name() + "_neighbourPatch-trans.obj");
         meshTools::writeOBJ(osN, nbrPatch0, nbrPoints);
 
         OFstream osO(patch_.name() + "_ownerPatch.obj");
         meshTools::writeOBJ(osO, patch_.localFaces(), patch_.localPoints());
     }
-*/
+
     // Construct/apply AMI interpolation to determine addressing and weights
     AMIPtr_.reset
     (

src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H

@@ -392,6 +392,7 @@ public:
             //- Get the patch on the region
             const polyPatch& samplePolyPatch() const;
 
+
         // Helpers
 
             //- Get the sample points

src/postProcessing/functionObjects/cloud/cloudInfo/cloudInfo.C

@@ -39,7 +39,10 @@ defineTypeNameAndDebug(cloudInfo, 0);
 
 void Foam::cloudInfo::writeFileHeader(const label i)
 {
-    file(i) << "# Time" << tab << "nParcels" << tab << "mass" << endl;
+    writeHeader(file(), "Cloud information");
+    writeCommented(file(), "Time");
+    writeTabbed(file(), "nParcels");
+    writeTabbed(file(), "mass");
 }
 
 

src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C

@@ -335,7 +335,6 @@ void Foam::fieldAverage::write()
 {
     if (active_)
     {
-        calcAverages();
         writeAverages();
         writeAveragingProperties();
 

src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransform.C

@@ -96,7 +96,20 @@ void Foam::fieldCoordinateSystemTransform::read(const dictionary& dict)
 
 void Foam::fieldCoordinateSystemTransform::execute()
 {
-    // Do nothing
+    if (active_)
+    {
+        Info<< type() << " " << name_ << " output:" << nl;
+
+        forAll(fieldSet_, fieldI)
+        {
+            // If necessary load field
+            transform<scalar>(fieldSet_[fieldI]);
+            transform<vector>(fieldSet_[fieldI]);
+            transform<sphericalTensor>(fieldSet_[fieldI]);
+            transform<symmTensor>(fieldSet_[fieldI]);
+            transform<tensor>(fieldSet_[fieldI]);
+        }
+    }
 }
 
 
@@ -114,16 +127,23 @@ void Foam::fieldCoordinateSystemTransform::timeSet()
 
 void Foam::fieldCoordinateSystemTransform::write()
 {
-    Info<< type() << " " << name_ << " output:" << nl;
+    if (active_)
-
-    forAll(fieldSet_, fieldI)
     {
-        // If necessary load field
+        Info<< type() << " " << name_ << " output:" << nl;
-        transform<scalar>(fieldSet_[fieldI]);
+
-        transform<vector>(fieldSet_[fieldI]);
+        forAll(fieldSet_, fieldI)
-        transform<sphericalTensor>(fieldSet_[fieldI]);
+        {
-        transform<symmTensor>(fieldSet_[fieldI]);
+            const word fieldName = fieldSet_[fieldI] + ":Transformed";
-        transform<tensor>(fieldSet_[fieldI]);
+
+            const regIOobject& field =
+                obr_.lookupObject<regIOobject>(fieldName);
+
+            Info<< "    writing field " << field.name() << nl;
+
+            field.write();
+        }
+
+        Info<< endl;
     }
 }
 

src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransformTemplates.C

@@ -37,7 +37,7 @@ void Foam::fieldCoordinateSystemTransform::transformField
     const Type& field
 ) const
 {
-    const word& fieldName = field.name() + "Transformed";
+    const word& fieldName = field.name() + ":Transformed";
 
     if (!obr_.foundObject<Type>(fieldName))
     {

src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.C

@@ -112,20 +112,23 @@ void Foam::fieldMinMax::read(const dictionary& dict)
 
 void Foam::fieldMinMax::writeFileHeader(const label i)
 {
-    file()
+    writeHeader(file(), "Field minima and maxima");
-        << "# Time" << token::TAB << "field" << token::TAB
+    writeCommented(file(), "Time");
-        << "min" << token::TAB << "position(min)";
+    writeTabbed(file(), "field");
+    writeTabbed(file(), "min");
+    writeTabbed(file(), "position(min)");
 
     if (Pstream::parRun())
     {
-        file() << token::TAB << "proc";
+        writeTabbed(file(), "processor");
     }
 
-    file() << token::TAB << "max" << token::TAB << "position(max)";
+    writeTabbed(file(), "max");
+    writeTabbed(file(), "position(max)");
 
     if (Pstream::parRun())
     {
-        file() << token::TAB << "proc";
+        writeTabbed(file(), "processor");
     }
 
     file() << endl;

src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -111,21 +111,25 @@ void Foam::fieldMinMax::calcMinMaxFields
                     scalar maxValue = maxVs[maxI];
                     const vector& maxC = maxCs[maxI];
 
+                    file()<< obr_.time().value();
+                    writeTabbed(file(), fieldName);
+
                     file()
-                        << obr_.time().value() << token::TAB
+                        << token::TAB << minValue
-                        << fieldName << token::TAB
+                        << token::TAB << minC;
-                        << minValue << token::TAB << minC;
 
                     if (Pstream::parRun())
                     {
-                        file() << token::TAB << minI;
+                        file()<< token::TAB << minI;
                     }
 
-                    file() << token::TAB << maxValue << token::TAB << maxC;
+                    file()
+                        << token::TAB << maxValue
+                        << token::TAB << maxC;
 
                     if (Pstream::parRun())
                     {
-                        file() << token::TAB << maxI;
+                        file()<< token::TAB << maxI;
                     }
 
                     file() << endl;
@@ -212,21 +216,25 @@ void Foam::fieldMinMax::calcMinMaxFields
                     Type maxValue = maxVs[maxI];
                     const vector& maxC = maxCs[maxI];
 
+                    file()<< obr_.time().value();
+                    writeTabbed(file(), fieldName);
+
                     file()
-                        << obr_.time().value() << token::TAB
+                        << token::TAB << minValue
-                        << fieldName << token::TAB
+                        << token::TAB << minC;
-                        << minValue << token::TAB << minC;
 
                     if (Pstream::parRun())
                     {
-                        file() << token::TAB << minI;
+                        file()<< token::TAB << minI;
                     }
 
-                    file() << token::TAB << maxValue << token::TAB << maxC;
+                    file()
+                        << token::TAB << maxValue
+                        << token::TAB << maxC;
 
                     if (Pstream::parRun())
                     {
-                        file() << token::TAB << maxI;
+                        file()<< token::TAB << maxI;
                     }
 
                     file() << endl;

src/postProcessing/functionObjects/field/fieldValues/fieldValueDelta/fieldValueDelta.C

@@ -94,11 +94,10 @@ void Foam::fieldValues::fieldValueDelta::writeFileHeader(const label i)
 
     Ostream& os = file();
 
-    os  << "# Source1   : " << source1Ptr_->name() << nl
+    writeHeaderValue(os, "Source1", source1Ptr_->name());
-        << "# Source2   : " << source2Ptr_->name() << nl
+    writeHeaderValue(os, "Source2", source2Ptr_->name());
-        << "# Operation : " << operationTypeNames_[operation_] << nl;
+    writeHeaderValue(os, "Operation", operationTypeNames_[operation_]);
-
+    writeCommented(os, "Time");
-    os  << "# Time";
 
     forAll(commonFields, i)
     {
@@ -156,7 +155,7 @@ void Foam::fieldValues::fieldValueDelta::write()
 
     if (Pstream::master())
     {
-        file()<< obr_.time().timeName();
+        file()<< obr_.time().value();
     }
 
     if (log_)

src/postProcessing/functionObjects/field/nearWallFields/nearWallFields.C

@@ -328,32 +328,8 @@ void Foam::nearWallFields::execute()
     {
         Info<< "nearWallFields:execute()" << endl;
     }
-}
 
 
-void Foam::nearWallFields::end()
-{
-    if (debug)
-    {
-        Info<< "nearWallFields:end()" << endl;
-    }
-}
-
-
-void Foam::nearWallFields::timeSet()
-{
-    // Do nothing
-}
-
-
-void Foam::nearWallFields::write()
-{
-    if (debug)
-    {
-        Info<< "nearWallFields:write()" << endl;
-    }
-
-    // Do nothing
     if (active_)
     {
         if
@@ -380,7 +356,7 @@ void Foam::nearWallFields::write()
 
         Info<< type() << " " << name_ << " output:" << nl;
 
-        Info<< "    Writing sampled fields to " << obr_.time().timeName()
+        Info<< "    Sampling fields fields to " << obr_.time().timeName()
             << endl;
 
         sampleFields(vsf_);
@@ -388,8 +364,37 @@ void Foam::nearWallFields::write()
         sampleFields(vSpheretf_);
         sampleFields(vSymmtf_);
         sampleFields(vtf_);
+    }
+}
+
+
+void Foam::nearWallFields::end()
+{
+    if (debug)
+    {
+        Info<< "nearWallFields:end()" << endl;
+    }
+}
+
+
+void Foam::nearWallFields::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::nearWallFields::write()
+{
+    if (debug)
+    {
+        Info<< "nearWallFields:write()" << endl;
+    }
+
+    if (active_)
+    {
+        Info<< "    Writing sampled fields to " << obr_.time().timeName()
+            << endl;
 
-        // Write fields
         forAll(vsf_, i)
         {
             vsf_[i].write();

src/postProcessing/functionObjects/field/processorField/processorField.C

@@ -53,6 +53,27 @@ Foam::processorField::processorField
     if (isA<fvMesh>(obr_))
     {
         read(dict);
+
+        const fvMesh& mesh = refCast<const fvMesh>(obr_);
+
+        volScalarField* procFieldPtr
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "processorID",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedScalar("0", dimless, 0.0)
+            )
+        );
+
+        mesh.objectRegistry::store(procFieldPtr);
     }
     else
     {
@@ -88,7 +109,14 @@ void Foam::processorField::read(const dictionary& dict)
 
 void Foam::processorField::execute()
 {
-    // Do nothing
+    if (active_)
+    {
+        const volScalarField& procField =
+            obr_.lookupObject<volScalarField>("processorID");
+
+        const_cast<volScalarField&>(procField) ==
+            dimensionedScalar("procI", dimless, Pstream::myProcNo());
+    }
 }
 
 
@@ -108,20 +136,8 @@ void Foam::processorField::write()
 {
     if (active_)
     {
-        const fvMesh& mesh = refCast<const fvMesh>(obr_);
+        const volScalarField& procField =
-        volScalarField procField
+            obr_.lookupObject<volScalarField>("processorID");
-        (
-            IOobject
-            (
-                "processorID",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionedScalar("procI", dimless, Pstream::myProcNo())
-        );
 
         procField.write();
     }

src/postProcessing/functionObjects/field/processorField/processorFieldTemplates.C

@@ -1,158 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "processorField.H"
-#include "volFields.H"
-#include "surfaceFields.H"
-#include "Time.H"
-#include "transformGeometricField.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type>
-void Foam::processorField::transformField
-(
-    const Type& field
-) const
-{
-    const word& fieldName = field.name() + "Transformed";
-
-    dimensionedTensor R("R", field.dimensions(), coordSys_.R());
-
-    if (obr_.foundObject<Type>(fieldName))
-    {
-        Type& transField =
-            const_cast<Type&>(obr_.lookupObject<Type>(fieldName));
-
-        transField == field;
-
-        forAll(field, i)
-        {
-            Foam::transform(transField, R, transField);
-        }
-
-        transField.write();
-    }
-    else
-    {
-        Type& transField = obr_.store
-        (
-            new Type
-            (
-                IOobject
-                (
-                    fieldName,
-                    obr_.time().timeName(),
-                    obr_,
-                    IOobject::READ_IF_PRESENT,
-                    IOobject::NO_WRITE
-                ),
-                field
-            )
-        );
-
-        forAll(field, i)
-        {
-            Foam::transform(transField, R, transField);
-        }
-
-        transField.write();
-    }
-}
-
-
-template<class Type>
-void Foam::processorField::transform
-(
-    const word& fieldName
-) const
-{
-    typedef GeometricField<Type, fvPatchField, volMesh> vfType;
-    typedef GeometricField<Type, fvsPatchField, surfaceMesh> sfType;
-
-    if (obr_.foundObject<vfType>(fieldName))
-    {
-        if (debug)
-        {
-            Info<< type() << ": Field " << fieldName << " already in database"
-                << endl;
-        }
-
-        transformField<vfType>(obr_.lookupObject<vfType>(fieldName));
-    }
-    else if (obr_.foundObject<sfType>(fieldName))
-    {
-        if (debug)
-        {
-            Info<< type() << ": Field " << fieldName << " already in database"
-                << endl;
-        }
-
-        transformField<sfType>(obr_.lookupObject<sfType>(fieldName));
-    }
-    else
-    {
-        IOobject fieldHeader
-        (
-            fieldName,
-            obr_.time().timeName(),
-            obr_,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        );
-
-        if
-        (
-            fieldHeader.headerOk()
-         && fieldHeader.headerClassName() == vfType::typeName
-        )
-        {
-            if (debug)
-            {
-                Info<< type() << ": Field " << fieldName << " read from file"
-                    << endl;
-            }
-
-            transformField<vfType>(obr_.lookupObject<vfType>(fieldName));
-        }
-        else if
-        (
-            fieldHeader.headerOk()
-         && fieldHeader.headerClassName() == sfType::typeName
-        )
-        {
-            if (debug)
-            {
-                Info<< type() << ": Field " << fieldName << " read from file"
-                    << endl;
-            }
-
-            transformField<sfType>(obr_.lookupObject<sfType>(fieldName));
-        }
-    }
-}
-
-
-// ************************************************************************* //

src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H

@@ -47,7 +47,7 @@ Description
     - alpha_liquidCore : alpha with outside liquid core set to 0
     - alpha_background : alpha with outside background set to 0.
 
-    Histogram:
+    %Histogram:
     - determine histogram of diameter (given minDiameter, maxDiameter, nBins)
     - write graph of number of droplets per bin
     - write graph of sum, average and deviation of droplet volume per bin

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

@@ -45,34 +45,65 @@ void Foam::forceCoeffs::writeFileHeader(const label i)
     {
         // force coeff data
 
+        writeHeader(file(i), "Force coefficients");
+        writeHeaderValue(file(i), "liftDir", liftDir_);
+        writeHeaderValue(file(i), "dragDir", dragDir_);
+        writeHeaderValue(file(i), "pitchAxis", pitchAxis_);
+        writeHeaderValue(file(i), "magUInf", magUInf_);
+        writeHeaderValue(file(i), "lRef", lRef_);
+        writeHeaderValue(file(i), "Aref", Aref_);
+        writeHeaderValue(file(i), "CofR", coordSys_.origin());
+        writeCommented(file(i), "Time");
+        writeTabbed(file(i), "Cm");
+        writeTabbed(file(i), "Cd");
+        writeTabbed(file(i), "Cl");
+        writeTabbed(file(i), "Cl(f)");
+        writeTabbed(file(i), "Cl(r)");
         file(i)
-            << "# liftDir   : " << liftDir_ << nl
+            << tab << "Cm" << tab << "Cd" << tab << "Cl" << tab << "Cl(f)"
-            << "# dragDir   : " << dragDir_ << nl
+            << tab << "Cl(r)";
-            << "# pitchAxis : " << pitchAxis_ << nl
-            << "# magUInf   : " << magUInf_ << nl
-            << "# lRef      : " << lRef_ << nl
-            << "# Aref      : " << Aref_ << nl
-            << "# CofR      : " << coordSys_.origin() << nl
-            << "# Time" << tab << "Cm" << tab << "Cd" << tab << "Cl" << tab
-            << "Cl(f)" << tab << "Cl(r)";
     }
     else if (i == 1)
     {
         // bin coeff data
 
-        file(i)
+        writeHeader(file(i), "Force coefficient bins");
-            << "# bins      : " << nBin_ << nl
+        writeHeaderValue(file(i), "bins", nBin_);
-            << "# start     : " << binMin_ << nl
+        writeHeaderValue(file(i), "start", binMin_);
-            << "# delta     : " << binDx_ << nl
+        writeHeaderValue(file(i), "delta", binDx_);
-            << "# direction : " << binDir_ << nl
+        writeHeaderValue(file(i), "direction", binDir_);
-            << "# Time";
+
+        vectorField binPoints(nBin_);
+        writeCommented(file(i), "x co-ords  :");
+        forAll(binPoints, pointI)
+        {
+            binPoints[pointI] = (binMin_ + (pointI + 1)*binDx_)*binDir_;
+            file(i) << tab << binPoints[pointI].x();
+        }
+        file(i) << nl;
+
+        writeCommented(file(i), "y co-ords  :");
+        forAll(binPoints, pointI)
+        {
+            file(i) << tab << binPoints[pointI].y();
+        }
+        file(i) << nl;
+
+        writeCommented(file(i), "z co-ords  :");
+        forAll(binPoints, pointI)
+        {
+            file(i) << tab << binPoints[pointI].z();
+        }
+        file(i) << nl;
+
+        writeCommented(file(i), "Time");
 
         for (label j = 0; j < nBin_; j++)
         {
             const word jn('[' + Foam::name(j) + ']');
-
+            writeTabbed(file(i), "Cm" + jn);
-            file(i)
+            writeTabbed(file(i), "Cd" + jn);
-                << tab << "Cm" << jn << tab << "Cd" << jn << tab << "Cl" << jn;
+            writeTabbed(file(i), "Cl" + jn);
         }
     }
     else
@@ -193,9 +224,8 @@ void Foam::forceCoeffs::write()
         scalar Clr = Cl/2.0 - Cm;
 
         file(0)
-            << obr_.time().value() << tab
+            << obr_.time().value() << tab << Cm << tab  << Cd
-            << Cm << tab << Cd << tab << Cl << tab << Clf << tab << Clr
+            << tab << Cl << tab << Clf << tab << Clr << endl;
-            << endl;
 
         if (log_)
         {

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

@@ -73,11 +73,13 @@ void Foam::forces::writeFileHeader(const label i)
     {
         // force data
 
+        writeHeader(file(i), "Forces");
+        writeHeaderValue(file(i), "CofR", coordSys_.origin());
+        writeCommented(file(i), "Time");
+
         file(i)
-            << "# CofR      : " << coordSys_.origin() << nl
+            << "forces[pressure,viscous,porous] "
-            << "# Time" << tab
+            << "moment[pressure,viscous,porous]";
-            << "forces(pressure,viscous,porous) "
-            << "moment(pressure,viscous,porous)";
 
         if (localSystem_)
         {
@@ -91,32 +93,54 @@ void Foam::forces::writeFileHeader(const label i)
     {
         // bin data
 
-        file(i)
+        writeHeader(file(i), "Force bins");
-            << "# bins      : " << nBin_ << nl
+        writeHeaderValue(file(i), "bins", nBin_);
-            << "# start     : " << binMin_ << nl
+        writeHeaderValue(file(i), "start", binMin_);
-            << "# delta     : " << binDx_ << nl
+        writeHeaderValue(file(i), "delta", binDx_);
-            << "# direction : " << binDir_ << nl
+        writeHeaderValue(file(i), "direction", binDir_);
-            << "# Time";
+
+        vectorField binPoints(nBin_);
+        writeCommented(file(i), "x co-ords  :");
+        forAll(binPoints, pointI)
+        {
+            binPoints[pointI] = (binMin_ + (pointI + 1)*binDx_)*binDir_;
+            file(i) << tab << binPoints[pointI].x();
+        }
+        file(i) << nl;
+
+        writeCommented(file(i), "y co-ords  :");
+        forAll(binPoints, pointI)
+        {
+            file(i) << tab << binPoints[pointI].y();
+        }
+        file(i) << nl;
+
+        writeCommented(file(i), "z co-ords  :");
+        forAll(binPoints, pointI)
+        {
+            file(i) << tab << binPoints[pointI].z();
+        }
+        file(i) << nl;
+
+        writeCommented(file(i), "Time");
 
         for (label j = 0; j < nBin_; j++)
         {
             const word jn('[' + Foam::name(j) + ']');
+            const word f("forces" + jn + "[pressure,viscous,porous]");
+            const word m("moments" + jn + "[pressure,viscous,porous]");
 
-            file(i)
+            file(i)<< tab << f << tab << m;
-                << tab
-                << "forces" << jn << "(pressure,viscous,porous) "
-                << "moment" << jn << "(pressure,viscous,porous)";
         }
         if (localSystem_)
         {
             for (label j = 0; j < nBin_; j++)
             {
                 const word jn('[' + Foam::name(j) + ']');
+                const word f("localForces" + jn + "[pressure,viscous,porous]");
+                const word m("localMoments" + jn + "[pressure,viscous,porous]");
 
-                file(i)
+                file(i)<< tab << f << tab << m;
-                    << tab
-                    << "localForces" << jn << "(pressure,viscous,porous) "
-                    << "localMoments" << jn << "(pressure,viscous,porous)";
             }
         }
     }
@@ -370,28 +394,24 @@ void Foam::forces::writeForces()
     if (log_)
     {
         Info<< type() << " " << name_ << " output:" << nl
-            << "    forces(pressure,viscous,porous) = ("
+            << "    sum of forces:" << nl
-            << sum(force_[0]) << ","
+            << "        pressure : " << sum(force_[0]) << nl
-            << sum(force_[1]) << ","
+            << "        viscous  : " << sum(force_[1]) << nl
-            << sum(force_[2]) << ")" << nl
+            << "        porous   : " << sum(force_[2]) << nl
-            << "    moment(pressure,viscous,porous) = ("
+            << "    sum of moments:" << nl
-            << sum(moment_[0]) << ","
+            << "        pressure : " << sum(moment_[0]) << nl
-            << sum(moment_[1]) << ","
+            << "        viscous  : " << sum(moment_[1]) << nl
-            << sum(moment_[2]) << ")"
+            << "        porous   : " << sum(moment_[2])
-            << nl;
+            << endl;
     }
 
-    file(0) << obr_.time().value() << tab
+    file(0) << obr_.time().value() << tab << setw(1) << '['
-        << "("
+        << sum(force_[0]) << setw(1) << ','
-        << sum(force_[0]) << ","
+        << sum(force_[1]) << setw(1) << ","
-        << sum(force_[1]) << ","
+        << sum(force_[2]) << setw(3) << "] ["
-        << sum(force_[2])
+        << sum(moment_[0]) << setw(1) << ","
-        << ") "
+        << sum(moment_[1]) << setw(1) << ","
-        << "("
+        << sum(moment_[2]) << setw(1) << "]"
-        << sum(moment_[0]) << ","
-        << sum(moment_[1]) << ","
-        << sum(moment_[2])
-        << ")"
         << endl;
 
     if (localSystem_)
@@ -403,17 +423,13 @@ void Foam::forces::writeForces()
         vectorField localMomentT(coordSys_.localVector(moment_[1]));
         vectorField localMomentP(coordSys_.localVector(moment_[2]));
 
-        file(0) << obr_.time().value() << tab
+        file(0) << obr_.time().value() << tab << setw(1) << "["
-            << "("
+            << sum(localForceN) << setw(1) << ","
-            << sum(localForceN) << ","
+            << sum(localForceT) << setw(1) << ","
-            << sum(localForceT) << ","
+            << sum(localForceP) << setw(3) << "] ["
-            << sum(localForceP)
+            << sum(localMomentN) << setw(1) << ","
-            << ") "
+            << sum(localMomentT) << setw(1) << ","
-            << "("
+            << sum(localMomentP) << setw(1) << "]"
-            << sum(localMomentN) << ","
-            << sum(localMomentT) << ","
-            << sum(localMomentP)
-            << ")"
             << endl;
     }
 }
@@ -448,9 +464,13 @@ void Foam::forces::writeBins()
     forAll(f[0], i)
     {
         file(1)
-            << tab
+            << tab << setw(1) << "["
-            << "(" << f[0][i] << "," << f[1][i] << "," << f[2][i] << ") "
+            << f[0][i] << setw(1) << ","
-            << "(" << m[0][i] << "," << m[1][i] << "," << m[2][i] << ")";
+            << f[1][i] << setw(1) << ","
+            << f[2][i] << setw(3) << "] ["
+            << m[0][i] << setw(1) << ","
+            << m[1][i] << setw(1) << ","
+            << m[2][i] << setw(1) << "]";
     }
 
     if (localSystem_)
@@ -480,9 +500,13 @@ void Foam::forces::writeBins()
         forAll(lf[0], i)
         {
             file(1)
-                << tab
+                << tab << setw(1) << "["
-                << "(" << lf[0][i] << "," << lf[1][i] << "," << lf[2][i] << ") "
+                << lf[0][i] << setw(1) << ","
-                << "(" << lm[0][i] << "," << lm[1][i] << "," << lm[2][i] << ")";
+                << lf[1][i] << setw(1) << ","
+                << lf[2][i] << setw(3) << "] ["
+                << lm[0][i] << setw(1) << ","
+                << lm[1][i] << setw(1) << ","
+                << lm[2][i] << setw(1) << "]";
         }
     }
 

src/postProcessing/functionObjects/fvTools/calcFvcDiv/calcFvcDiv.C

@@ -125,29 +125,16 @@ void Foam::calcFvcDiv::read(const dictionary& dict)
     {
         dict.lookup("fieldName") >> fieldName_;
         dict.lookup("resultName") >> resultName_;
+
+        if (resultName_ == "none")
+        {
+            resultName_ = "fvc::div(" + fieldName_ + ")";
+        }
     }
 }
 
 
 void Foam::calcFvcDiv::execute()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::calcFvcDiv::end()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::calcFvcDiv::timeSet()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::calcFvcDiv::write()
 {
     if (active_)
     {
@@ -165,4 +152,34 @@ void Foam::calcFvcDiv::write()
 }
 
 
+void Foam::calcFvcDiv::end()
+{
+    // Do nothing
+}
+
+
+void Foam::calcFvcDiv::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::calcFvcDiv::write()
+{
+    if (active_)
+    {
+        if (obr_.foundObject<regIOobject>(resultName_))
+        {
+            const regIOobject& field =
+                obr_.lookupObject<regIOobject>(resultName_);
+
+            Info<< type() << " " << name_ << " output:" << nl
+                << "    writing field " << field.name() << nl << endl;
+
+            field.write();
+        }
+    }
+}
+
+
 // ************************************************************************* //

src/postProcessing/functionObjects/fvTools/calcFvcDiv/calcFvcDivTemplates.C

@@ -38,25 +38,14 @@ void Foam::calcFvcDiv::calcDiv
 {
     const fvMesh& mesh = refCast<const fvMesh>(obr_);
 
-    word divName = resultName;
-    if (divName == "none")
-    {
-        divName = "fvc::div(" + fieldName + ")";
-    }
-
     if (mesh.foundObject<FieldType>(fieldName))
     {
         const FieldType& vf = mesh.lookupObject<FieldType>(fieldName);
 
-        volScalarField& field = divField(divName, vf.dimensions());
+        volScalarField& field = divField(resultName, vf.dimensions());
 
         field = fvc::div(vf);
 
-        Info<< type() << " " << name_ << " output:" << nl
-            << "    writing field " << field.name() << nl << endl;
-
-        field.write();
-
         processed = true;
     }
 }

src/postProcessing/functionObjects/fvTools/calcFvcGrad/calcFvcGrad.C

@@ -87,29 +87,16 @@ void Foam::calcFvcGrad::read(const dictionary& dict)
     {
         dict.lookup("fieldName") >> fieldName_;
         dict.lookup("resultName") >> resultName_;
+
+        if (resultName_ == "none")
+        {
+            resultName_ = "fvc::grad(" + fieldName_ + ")";
+        }
     }
 }
 
 
 void Foam::calcFvcGrad::execute()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::calcFvcGrad::end()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::calcFvcGrad::timeSet()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::calcFvcGrad::write()
 {
     if (active_)
     {
@@ -127,4 +114,34 @@ void Foam::calcFvcGrad::write()
 }
 
 
+void Foam::calcFvcGrad::end()
+{
+    // Do nothing
+}
+
+
+void Foam::calcFvcGrad::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::calcFvcGrad::write()
+{
+    if (active_)
+    {
+        if (obr_.foundObject<regIOobject>(resultName_))
+        {
+            const regIOobject& field =
+                obr_.lookupObject<regIOobject>(resultName_);
+
+            Info<< type() << " " << name_ << " output:" << nl
+                << "    writing field " << field.name() << nl << endl;
+
+            field.write();
+        }
+    }
+}
+
+
 // ************************************************************************* //

src/postProcessing/functionObjects/fvTools/calcFvcGrad/calcFvcGradTemplates.C

@@ -91,40 +91,25 @@ void Foam::calcFvcGrad::calcGrad
 
     const fvMesh& mesh = refCast<const fvMesh>(obr_);
 
-    word gradName = resultName;
-    if (gradName == "none")
-    {
-        gradName = "fvc::grad(" + fieldName + ")";
-    }
 
     if (mesh.foundObject<vfType>(fieldName))
     {
         const vfType& vf = mesh.lookupObject<vfType>(fieldName);
 
-        vfGradType& field = gradField<Type>(gradName, vf.dimensions());
+        vfGradType& field = gradField<Type>(resultName, vf.dimensions());
 
         field = fvc::grad(vf);
 
-        Info<< type() << " output:" << nl
-            << "    writing " << field.name() << " field" << nl << endl;
-
-        field.write();
-
         processed = true;
     }
     else if (mesh.foundObject<sfType>(fieldName))
     {
         const sfType& sf = mesh.lookupObject<sfType>(fieldName);
 
-        vfGradType& field = gradField<Type>(gradName, sf.dimensions());
+        vfGradType& field = gradField<Type>(resultName, sf.dimensions());
 
         field = fvc::grad(sf);
 
-        Info<< type() << " " << name_ << " output:" << nl
-            << "    writing field " << field.name() << nl << endl;
-
-        field.write();
-
         processed = true;
     }
 }

src/postProcessing/functionObjects/fvTools/calcMag/calcMag.C

@@ -87,29 +87,16 @@ void Foam::calcMag::read(const dictionary& dict)
     {
         dict.lookup("fieldName") >> fieldName_;
         dict.lookup("resultName") >> resultName_;
+
+        if (resultName_ == "none")
+        {
+            resultName_ = "mag(" + fieldName_ + ")";
+        }
     }
 }
 
 
 void Foam::calcMag::execute()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::calcMag::end()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::calcMag::timeSet()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::calcMag::write()
 {
     if (active_)
     {
@@ -130,4 +117,34 @@ void Foam::calcMag::write()
 }
 
 
+void Foam::calcMag::end()
+{
+    // Do nothing
+}
+
+
+void Foam::calcMag::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::calcMag::write()
+{
+    if (active_)
+    {
+        if (obr_.foundObject<regIOobject>(resultName_))
+        {
+            const regIOobject& field =
+                obr_.lookupObject<regIOobject>(resultName_);
+
+            Info<< type() << " " << name_ << " output:" << nl
+                << "    writing field " << field.name() << nl << endl;
+
+            field.write();
+        }
+    }
+}
+
+
 // ************************************************************************* //

src/postProcessing/functionObjects/fvTools/calcMag/calcMagTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,26 +80,15 @@ void Foam::calcMag::calc
 
     const fvMesh& mesh = refCast<const fvMesh>(obr_);
 
-    word magName = resultName;
-    if (magName == "none")
-    {
-        magName = "mag(" + fieldName + ")";
-    }
-
     if (mesh.foundObject<vfType>(fieldName))
     {
         const vfType& vf = mesh.lookupObject<vfType>(fieldName);
 
         volScalarField& field =
-            magField<volScalarField>(magName, vf.dimensions());
+            magField<volScalarField>(resultName_, vf.dimensions());
 
         field = mag(vf);
 
-        Info<< type() << " output:" << nl
-            << "    writing " << field.name() << " field" << nl << endl;
-
-        field.write();
-
         processed = true;
     }
     else if (mesh.foundObject<sfType>(fieldName))
@@ -107,15 +96,10 @@ void Foam::calcMag::calc
         const sfType& sf = mesh.lookupObject<sfType>(fieldName);
 
         surfaceScalarField& field =
-            magField<surfaceScalarField>(magName, sf.dimensions());
+            magField<surfaceScalarField>(resultName_, sf.dimensions());
 
         field = mag(sf);
 
-        Info<< type() << " output:" << nl
-            << "    writing " << field.name() << " field" << nl << endl;
-
-        field.write();
-
         processed = true;
     }
 }

src/postProcessing/functionObjects/utilities/CourantNo/CourantNo.C

@@ -154,24 +154,6 @@ void Foam::CourantNo::read(const dictionary& dict)
 
 
 void Foam::CourantNo::execute()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::CourantNo::end()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::CourantNo::timeSet()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::CourantNo::write()
 {
     if (active_)
     {
@@ -197,6 +179,28 @@ void Foam::CourantNo::write()
         iField = 0.5*sumPhi/mesh.V().field()*mesh.time().deltaTValue();
 
         CourantNo.correctBoundaryConditions();
+    }
+}
+
+
+void Foam::CourantNo::end()
+{
+    // Do nothing
+}
+
+
+void Foam::CourantNo::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::CourantNo::write()
+{
+    if (active_)
+    {
+        const volScalarField& CourantNo =
+            obr_.lookupObject<volScalarField>(type());
 
         Info<< type() << " " << name_ << " output:" << nl
             << "    writing field " << CourantNo.name() << nl

src/postProcessing/functionObjects/utilities/DESModelRegions/DESModelRegions.C

@@ -42,9 +42,11 @@ defineTypeNameAndDebug(DESModelRegions, 0);
 
 void Foam::DESModelRegions::writeFileHeader(const label i)
 {
-    file() << "# DES model region coverage (% volume)" << nl
+    writeHeader(file(), "DES model region coverage (% volume)");
-        << "# time " << token::TAB << "LES" << token::TAB << "RAS"
+
-        << endl;
+    writeCommented(file(), "Time");
+    writeTabbed(file(), "LES");
+    writeTabbed(file(), "RAS");
 }
 
 
@@ -127,24 +129,6 @@ void Foam::DESModelRegions::read(const dictionary& dict)
 
 
 void Foam::DESModelRegions::execute()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::DESModelRegions::end()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::DESModelRegions::timeSet()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::DESModelRegions::write()
 {
     typedef incompressible::turbulenceModel icoModel;
     typedef incompressible::DESModel icoDESModel;
@@ -206,19 +190,17 @@ void Foam::DESModelRegions::write()
 
             if (Pstream::master() && log_)
             {
-                file() << obr_.time().timeName() << token::TAB
+                file() << obr_.time().value()
-                    << prc << token::TAB << 100.0 - prc << endl;
+                    << token::TAB << prc
+                    << token::TAB << 100.0 - prc
+                    << endl;
             }
 
             if (log_)
             {
                 Info<< "    LES = " << prc << " % (volume)" << nl
-                    << "    RAS = " << 100.0 - prc << " % (volume)" << nl
+                    << "    RAS = " << 100.0 - prc << " % (volume)" << endl;
-                    << "    writing field " << DESModelRegions.name() << nl
-                    << endl;
             }
-
-            DESModelRegions.write();
         }
         else
         {
@@ -232,4 +214,33 @@ void Foam::DESModelRegions::write()
 }
 
 
+void Foam::DESModelRegions::end()
+{
+    // Do nothing
+}
+
+
+void Foam::DESModelRegions::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::DESModelRegions::write()
+{
+    if (log_)
+    {
+        const volScalarField& DESModelRegions =
+            obr_.lookupObject<volScalarField>(type());
+
+
+        Info<< type() << " " << name_ <<  " output:" << nl
+            << "    writing field " << DESModelRegions.name() << nl
+            << endl;
+
+        DESModelRegions.write();
+    }
+}
+
+
 // ************************************************************************* //

src/postProcessing/functionObjects/utilities/Lambda2/Lambda2.C

@@ -33,7 +33,7 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(Lambda2, 0);
+    defineTypeNameAndDebug(Lambda2, 0);
 }
 
 
@@ -115,24 +115,6 @@ void Foam::Lambda2::read(const dictionary& dict)
 
 
 void Foam::Lambda2::execute()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::Lambda2::end()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::Lambda2::timeSet()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::Lambda2::write()
 {
     if (active_)
     {
@@ -156,6 +138,28 @@ void Foam::Lambda2::write()
             );
 
         Lambda2 = -eigenValues(SSplusWW)().component(vector::Y);
+    }
+}
+
+
+void Foam::Lambda2::end()
+{
+    // Do nothing
+}
+
+
+void Foam::Lambda2::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::Lambda2::write()
+{
+    if (active_)
+    {
+        const volScalarField& Lambda2 =
+            obr_.lookupObject<volScalarField>(type());
 
         Info<< type() << " " << name_ << " output:" << nl
             << "    writing field " << Lambda2.name() << nl

src/postProcessing/functionObjects/utilities/Peclet/Peclet.C

@@ -35,7 +35,7 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(Peclet, 0);
+    defineTypeNameAndDebug(Peclet, 0);
 }
 
 
@@ -119,23 +119,6 @@ void Foam::Peclet::read(const dictionary& dict)
 
 
 void Foam::Peclet::execute()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::Peclet::end()
-{
-    // Do nothing - only valid on write
-}
-
-void Foam::Peclet::timeSet()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::Peclet::write()
 {
     typedef compressible::turbulenceModel cmpTurbModel;
     typedef incompressible::turbulenceModel icoTurbModel;
@@ -208,6 +191,27 @@ void Foam::Peclet::write()
                *mesh.surfaceInterpolation::deltaCoeffs()
                *fvc::interpolate(nuEff)
             );
+    }
+}
+
+
+void Foam::Peclet::end()
+{
+    // Do nothing - only valid on write
+}
+
+void Foam::Peclet::timeSet()
+{
+    // Do nothing - only valid on write
+}
+
+
+void Foam::Peclet::write()
+{
+    if (active_)
+    {
+        const surfaceScalarField& Peclet =
+            obr_.lookupObject<surfaceScalarField>(type());
 
         Info<< type() << " " << name_ << " output:" << nl
             << "    writing field " << Peclet.name() << nl

src/postProcessing/functionObjects/utilities/Q/Q.C

@@ -114,24 +114,6 @@ void Foam::Q::read(const dictionary& dict)
 
 
 void Foam::Q::execute()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::Q::end()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::Q::timeSet()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::Q::write()
 {
     if (active_)
     {
@@ -149,6 +131,28 @@ void Foam::Q::write()
             );
 
         Q = 0.5*(sqr(tr(gradU)) - tr(((gradU) & (gradU))));
+    }
+}
+
+
+void Foam::Q::end()
+{
+    // Do nothing
+}
+
+
+void Foam::Q::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::Q::write()
+{
+    if (active_)
+    {
+        const volScalarField& Q =
+            obr_.lookupObject<volScalarField>(type());
 
         Info<< type() << " " << name_ << " output:" << nl
             << "    writing field " << Q.name() << nl

src/postProcessing/functionObjects/utilities/pressureTools/pressureTools.C

@@ -31,7 +31,7 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(pressureTools, 0);
+    defineTypeNameAndDebug(pressureTools, 0);
 }
 
 
@@ -291,24 +291,6 @@ void Foam::pressureTools::read(const dictionary& dict)
 
 
 void Foam::pressureTools::execute()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::pressureTools::end()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::pressureTools::timeSet()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::pressureTools::write()
 {
     if (active_)
     {
@@ -321,6 +303,28 @@ void Foam::pressureTools::write()
             );
 
         pResult == convertToCoeff(rhoScale(p)*p + pDyn(p) + pRef());
+    }
+}
+
+
+void Foam::pressureTools::end()
+{
+    // Do nothing
+}
+
+
+void Foam::pressureTools::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::pressureTools::write()
+{
+    if (active_)
+    {
+        const volScalarField& pResult =
+            obr_.lookupObject<volScalarField>(pName());
 
         Info<< type() << " " << name_ << " output:" << nl
             << "    writing field " << pResult.name() << nl

src/postProcessing/functionObjects/utilities/wallShearStress/wallShearStress.C

@@ -43,9 +43,11 @@ defineTypeNameAndDebug(wallShearStress, 0);
 void Foam::wallShearStress::writeFileHeader(const label i)
 {
     // Add headers to output data
-    file() << "# Wall shear stress" << nl
+    writeHeader(file(), "Wall shear stress");
-        << "# time " << token::TAB << "patch" << token::TAB
+    writeCommented(file(), "Time");
-        << "min" << token::TAB << "max" << endl;
+    writeTabbed(file(), "patch");
+    writeTabbed(file(), "min");
+    writeTabbed(file(), "max");
 }
 
 
@@ -73,9 +75,11 @@ void Foam::wallShearStress::calcShearStress
 
         if (Pstream::master())
         {
-            file() << mesh.time().timeName() << token::TAB
+            file() << mesh.time().value()
-                << pp.name() << token::TAB << minSsp
+                << token::TAB << pp.name()
-                << token::TAB << maxSsp << endl;
+                << token::TAB << minSsp
+                << token::TAB << maxSsp
+                << endl;
         }
 
         if (log_)
@@ -220,24 +224,6 @@ void Foam::wallShearStress::read(const dictionary& dict)
 
 
 void Foam::wallShearStress::execute()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::wallShearStress::end()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::wallShearStress::timeSet()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::wallShearStress::write()
 {
     typedef compressible::turbulenceModel cmpModel;
     typedef incompressible::turbulenceModel icoModel;
@@ -282,8 +268,31 @@ void Foam::wallShearStress::write()
                 << "database" << exit(FatalError);
         }
 
-
         calcShearStress(mesh, Reff(), wallShearStress);
+    }
+}
+
+
+void Foam::wallShearStress::end()
+{
+    // Do nothing
+}
+
+
+void Foam::wallShearStress::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::wallShearStress::write()
+{
+    if (active_)
+    {
+        functionObjectFile::write();
+
+        const volVectorField& wallShearStress =
+            obr_.lookupObject<volVectorField>(type());
 
         if (log_)
         {

src/postProcessing/functionObjects/utilities/yPlusLES/yPlusLES.C

@@ -36,7 +36,7 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(yPlusLES, 0);
+    defineTypeNameAndDebug(yPlusLES, 0);
 }
 
 
@@ -44,10 +44,13 @@ defineTypeNameAndDebug(yPlusLES, 0);
 
 void Foam::yPlusLES::writeFileHeader(const label i)
 {
-    file() << "# y+ (LES)" << nl
+    writeHeader(file(), "y+ (LES)");
-        << "# time " << token::TAB << "patch" << token::TAB
+
-        << "min" << token::TAB << "max" << token::TAB << "average"
+    writeCommented(file(), "Time");
-        << endl;
+    writeTabbed(file(), "patch");
+    writeTabbed(file(), "min");
+    writeTabbed(file(), "max");
+    writeTabbed(file(), "average");
 }
 
 
@@ -100,10 +103,12 @@ void Foam::yPlusLES::calcIncompressibleYPlus
 
             if (Pstream::master())
             {
-                file() << obr_.time().value() << token::TAB
+                file() << obr_.time().value()
-                    << currPatch.name() << token::TAB
+                    << token::TAB << currPatch.name()
-                    << minYp << token::TAB << maxYp << token::TAB
+                    << token::TAB << minYp
-                    << avgYp << endl;
+                    << token::TAB << maxYp
+                    << token::TAB << avgYp
+                    << endl;
             }
         }
     }
@@ -166,10 +171,12 @@ void Foam::yPlusLES::calcCompressibleYPlus
 
             if (Pstream::master())
             {
-                file() << obr_.time().value() << token::TAB
+                file() << obr_.time().value()
-                    << currPatch.name() << token::TAB
+                    << token::TAB << currPatch.name()
-                    << minYp << token::TAB << maxYp << token::TAB
+                    << token::TAB << minYp
-                    << avgYp << endl;
+                    << token::TAB << maxYp
+                    << token::TAB << avgYp
+                    << endl;
             }
          }
     }
@@ -261,24 +268,6 @@ void Foam::yPlusLES::read(const dictionary& dict)
 
 
 void Foam::yPlusLES::execute()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::yPlusLES::end()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::yPlusLES::timeSet()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::yPlusLES::write()
 {
     if (active_)
     {
@@ -310,6 +299,30 @@ void Foam::yPlusLES::write()
         {
             calcIncompressibleYPlus(mesh, U, yPlusLES);
         }
+    }
+}
+
+
+void Foam::yPlusLES::end()
+{
+    // Do nothing
+}
+
+
+void Foam::yPlusLES::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::yPlusLES::write()
+{
+    if (active_)
+    {
+        functionObjectFile::write();
+
+        const volScalarField& yPlusLES =
+            obr_.lookupObject<volScalarField>(type());
 
         if (log_)
         {

src/postProcessing/functionObjects/utilities/yPlusRAS/yPlusRAS.C

@@ -46,10 +46,13 @@ defineTypeNameAndDebug(yPlusRAS, 0);
 
 void Foam::yPlusRAS::writeFileHeader(const label i)
 {
-    file() << "# y+ (RAS)" << nl
+    writeHeader(file(), "y+ (RAS)");
-        << "# time " << token::TAB << "patch" << token::TAB
+
-        << "min" << token::TAB << "max" << token::TAB << "average"
+    writeCommented(file(), "Time");
-        << endl;
+    writeTabbed(file(), "patch");
+    writeTabbed(file(), "min");
+    writeTabbed(file(), "max");
+    writeTabbed(file(), "average");
 }
 
 
@@ -95,10 +98,12 @@ void Foam::yPlusRAS::calcIncompressibleYPlus
 
             if (Pstream::master())
             {
-                file() << obr_.time().value() << token::TAB
+                file() << obr_.time().value()
-                    << nutPw.patch().name() << token::TAB
+                    << token::TAB << nutPw.patch().name()
-                    << minYp << token::TAB << maxYp << token::TAB
+                    << token::TAB << minYp
-                    << avgYp << endl;
+                    << token::TAB << maxYp
+                    << token::TAB << avgYp
+                    << endl;
             }
         }
     }
@@ -153,10 +158,12 @@ void Foam::yPlusRAS::calcCompressibleYPlus
 
             if (Pstream::master())
             {
-                file() << obr_.time().value() << token::TAB
+                file() << obr_.time().value()
-                    << mutPw.patch().name() << token::TAB
+                    << token::TAB << mutPw.patch().name()
-                    << minYp << token::TAB << maxYp << token::TAB
+                    << token::TAB << minYp
-                    << avgYp << endl;
+                    << token::TAB << maxYp
+                    << token::TAB << avgYp
+                    << endl;
             }
         }
     }
@@ -248,24 +255,6 @@ void Foam::yPlusRAS::read(const dictionary& dict)
 
 
 void Foam::yPlusRAS::execute()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::yPlusRAS::end()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::yPlusRAS::timeSet()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::yPlusRAS::write()
 {
     if (active_)
     {
@@ -295,6 +284,30 @@ void Foam::yPlusRAS::write()
         {
             calcIncompressibleYPlus(mesh, yPlusRAS);
         }
+    }
+}
+
+
+void Foam::yPlusRAS::end()
+{
+    // Do nothing
+}
+
+
+void Foam::yPlusRAS::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::yPlusRAS::write()
+{
+    if (active_)
+    {
+        functionObjectFile::write();
+
+        const volScalarField& yPlusRAS =
+            obr_.lookupObject<volScalarField>(type());
 
         if (log_)
         {

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C

@@ -440,6 +440,7 @@ kinematicSingleLayer::kinematicSingleLayer
     cumulativeContErr_(0.0),
 
     deltaSmall_("deltaSmall", dimLength, SMALL),
+    deltaCoLimit_(solution().lookupOrDefault("deltaCoLimit", 1e-4)),
 
     rho_
     (
@@ -898,15 +899,17 @@ scalar kinematicSingleLayer::CourantNumber() const
 
     if (regionMesh().nInternalFaces() > 0)
     {
-        const scalarField sumPhi(fvc::surfaceSum(mag(phi_)));
+        const scalarField sumPhi
+        (
+            fvc::surfaceSum(mag(phi_))().internalField()
+          / (deltaRho_.internalField() + ROOTVSMALL)
+        );
 
-        const scalarField& V = regionMesh().V();
+        forAll(delta_, i)
-
-        forAll(deltaRho_, i)
         {
-            if (deltaRho_[i] > SMALL)
+            if (delta_[i] > deltaCoLimit_)
             {
-                CoNum = max(CoNum, sumPhi[i]/deltaRho_[i]/V[i]);
+                CoNum = max(CoNum, sumPhi[i]/(delta_[i]*magSf()[i]));
             }
         }
 

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H

@@ -100,6 +100,9 @@ protected:
             //- Small delta
             const dimensionedScalar deltaSmall_;
 
+            //- Film thickness above which Courant number calculation in valid
+            scalar deltaCoLimit_;
+
 
         // Thermo properties
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,11 +45,35 @@ addToRunTimeSelectionTable(injectionModel, removeInjection, dictionary);
 removeInjection::removeInjection
 (
     const surfaceFilmModel& owner,
-    const dictionary&
+    const dictionary& dict
 )
 :
-    injectionModel(owner)
+    injectionModel(type(), owner, dict),
-{}
+    deltaStable_(coeffs_.lookupOrDefault<scalar>("deltaStable", 0.0)),
+    mask_(owner.regionMesh().nCells(), -1)
+{
+    wordReList patches;
+    if (coeffs_.readIfPresent("patches", patches))
+    {
+        Info<< "        applying to patches:" << nl;
+        const polyBoundaryMesh& pbm = owner.regionMesh().boundaryMesh();
+        const labelHashSet patchSet = pbm.patchSet(patches);
+
+        forAllConstIter(labelHashSet, patchSet, iter)
+        {
+            label patchI = iter.key();
+            const polyPatch& pp = pbm[patchI];
+            UIndirectList<scalar>(mask_, pp.faceCells()) = 1.0;
+
+            Info<< "            " << pp.name() << endl;
+        }
+    }
+    else
+    {
+        Info<< "            applying to all patches" << endl;
+        mask_ = 1.0;
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -64,11 +88,23 @@ void removeInjection::correct
 (
     scalarField& availableMass,
     scalarField& massToInject,
-    scalarField&
+    scalarField& diameterToInject
 )
 {
-    massToInject = availableMass;
+    const scalarField& delta = owner().delta();
-    availableMass = 0.0;
+    const scalarField& rho = owner().rho();
+    const scalarField& magSf = owner().magSf();
+
+    forAll(delta, cellI)
+    {
+        if (mask_[cellI] > 0)
+        {
+            scalar ddelta = max(0.0, delta[cellI] - deltaStable_);
+            scalar dMass = ddelta*rho[cellI]*magSf[cellI];
+            massToInject[cellI] += dMass;
+            availableMass[cellI] -= dMass;
+        }
+    }
 }
 
 

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/removeInjection/removeInjection.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,6 +65,16 @@ private:
         void operator=(const removeInjection&);
 
 
+protected:
+
+        //- Stable film thickness - mass only removed if thickness execeeds
+        //  this threhold value
+        scalar deltaStable_;
+
+        //- Mask per cell to indicate whether mass can be removed
+        scalarField mask_;
+
+
 public:
 
     //- Runtime type information

src/sampling/sampledSurface/sampledPatch/sampledPatch.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -248,6 +248,7 @@ Foam::tmp<Foam::vectorField> Foam::sampledPatch::sample
     return sampleField(vField);
 }
 
+
 Foam::tmp<Foam::sphericalTensorField> Foam::sampledPatch::sample
 (
     const volSphericalTensorField& vField
@@ -292,6 +293,7 @@ Foam::tmp<Foam::vectorField> Foam::sampledPatch::sample
     return sampleField(sField);
 }
 
+
 Foam::tmp<Foam::sphericalTensorField> Foam::sampledPatch::sample
 (
     const surfaceSphericalTensorField& sField
@@ -318,6 +320,7 @@ Foam::tmp<Foam::tensorField> Foam::sampledPatch::sample
     return sampleField(sField);
 }
 
+
 Foam::tmp<Foam::scalarField> Foam::sampledPatch::interpolate
 (
     const interpolation<scalar>& interpolator
@@ -335,6 +338,7 @@ Foam::tmp<Foam::vectorField> Foam::sampledPatch::interpolate
     return interpolateField(interpolator);
 }
 
+
 Foam::tmp<Foam::sphericalTensorField> Foam::sampledPatch::interpolate
 (
     const interpolation<sphericalTensor>& interpolator

src/sampling/sampledSurface/sampledPatch/sampledPatch.H

@@ -53,7 +53,8 @@ class sampledPatch
     public sampledSurface
 {
     //- Private typedefs for convenience
-        typedef MeshedSurface<face> MeshStorage;
+    typedef MeshedSurface<face> MeshStorage;
+
 
     // Private data
 
@@ -78,16 +79,17 @@ class sampledPatch
         //- Start indices (in patchFaceLabels_) of patches
         labelList patchStart_;
 
+
     // Private Member Functions
 
-        //- sample field on faces
+        //- Sample field on faces
         template<class Type>
         tmp<Field<Type> > sampleField
         (
             const GeometricField<Type, fvPatchField, volMesh>& vField
         ) const;
 
-        //- sample surface field on faces
+        //- Sample surface field on faces
         template<class Type>
         tmp<Field<Type> > sampleField
         (
@@ -95,10 +97,9 @@ class sampledPatch
         ) const;
 
         template<class Type>
-        tmp<Field<Type> >
+        tmp<Field<Type> > interpolateField(const interpolation<Type>&) const;
-        interpolateField(const interpolation<Type>&) const;
 
-        //- remap action on triangulation or cleanup
+        //- Re-map action on triangulation or cleanup
         virtual void remapFaces(const labelUList& faceMap);
 
 
@@ -180,96 +181,99 @@ public:
         }
 
 
-        //- sample field on surface
+        // Sample
-        virtual tmp<scalarField> sample
+            //- Sample field on surface
-        (
+            virtual tmp<scalarField> sample
-            const volScalarField&
+            (
-        ) const;
+                const volScalarField&
+            ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
-        virtual tmp<vectorField> sample
+            virtual tmp<vectorField> sample
-        (
+            (
-            const volVectorField&
+                const volVectorField&
-        ) const;
+            ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
-        virtual tmp<sphericalTensorField> sample
+            virtual tmp<sphericalTensorField> sample
-        (
+            (
-            const volSphericalTensorField&
+                const volSphericalTensorField&
-        ) const;
+            ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
-        virtual tmp<symmTensorField> sample
+            virtual tmp<symmTensorField> sample
-        (
+            (
-            const volSymmTensorField&
+                const volSymmTensorField&
-        ) const;
+            ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
-        virtual tmp<tensorField> sample
+            virtual tmp<tensorField> sample
-        (
+            (
-            const volTensorField&
+                const volTensorField&
-        ) const;
+            ) const;
 
-        //- Surface sample field on surface
+            //- Surface sample field on surface
-        virtual tmp<scalarField> sample
+            virtual tmp<scalarField> sample
-        (
+            (
-            const surfaceScalarField&
+                const surfaceScalarField&
-        ) const;
+            ) const;
 
-        //- Surface Sample field on surface
+            //- Surface Sample field on surface
-        virtual tmp<vectorField> sample
+            virtual tmp<vectorField> sample
-        (
+            (
-            const surfaceVectorField&
+                const surfaceVectorField&
-        ) const;
+            ) const;
 
-        //- Surface sample field on surface
+            //- Surface sample field on surface
-        virtual tmp<sphericalTensorField> sample
+            virtual tmp<sphericalTensorField> sample
-        (
+            (
-            const surfaceSphericalTensorField&
+                const surfaceSphericalTensorField&
-        ) const;
+            ) const;
 
-        //- Surface sample field on surface
+            //- Surface sample field on surface
-        virtual tmp<symmTensorField> sample
+            virtual tmp<symmTensorField> sample
-        (
+            (
-            const surfaceSymmTensorField&
+                const surfaceSymmTensorField&
-        ) const;
+            ) const;
 
-        //- Surface sample field on surface
+            //- Surface sample field on surface
-        virtual tmp<tensorField> sample
+            virtual tmp<tensorField> sample
-        (
+            (
-            const surfaceTensorField&
+                const surfaceTensorField&
-        ) const;
+            ) const;
-
-        //- interpolate field on surface
-        virtual tmp<scalarField> interpolate
-        (
-            const interpolation<scalar>&
-        ) const;
 
 
-        //- interpolate field on surface
+        // Interpolate
-        virtual tmp<vectorField> interpolate
-        (
-            const interpolation<vector>&
-        ) const;
 
-        //- interpolate field on surface
+            //- Interpolate field on surface
-        virtual tmp<sphericalTensorField> interpolate
+            virtual tmp<scalarField> interpolate
-        (
+            (
-            const interpolation<sphericalTensor>&
+                const interpolation<scalar>&
-        ) const;
+            ) const;
 
-        //- interpolate field on surface
+            //- Interpolate field on surface
-        virtual tmp<symmTensorField> interpolate
+            virtual tmp<vectorField> interpolate
-        (
+            (
-            const interpolation<symmTensor>&
+                const interpolation<vector>&
-        ) const;
+            ) const;
 
-        //- interpolate field on surface
+            //- Interpolate field on surface
-        virtual tmp<tensorField> interpolate
+            virtual tmp<sphericalTensorField> interpolate
-        (
+            (
-            const interpolation<tensor>&
+                const interpolation<sphericalTensor>&
-        ) const;
+            ) const;
+
+            //- Interpolate field on surface
+            virtual tmp<symmTensorField> interpolate
+            (
+                const interpolation<symmTensor>&
+            ) const;
+
+            //- Interpolate field on surface
+            virtual tmp<tensorField> interpolate
+            (
+                const interpolation<tensor>&
+            ) const;
 
         //- Write
         virtual void print(Ostream&) const;

src/sampling/sampledSurface/sampledPatchInternalField/sampledPatchInternalField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -204,6 +204,7 @@ Foam::tmp<Foam::vectorField> Foam::sampledPatchInternalField::interpolate
     return interpolateField(interpolator);
 }
 
+
 Foam::tmp<Foam::sphericalTensorField>
 Foam::sampledPatchInternalField::interpolate
 (

src/sampling/sampledSurface/sampledPatchInternalField/sampledPatchInternalField.H

@@ -49,7 +49,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                       Class sampledPatchInternalField Declaration
+                  Class sampledPatchInternalField Declaration
 \*---------------------------------------------------------------------------*/
 
 class sampledPatchInternalField
@@ -74,6 +74,7 @@ class sampledPatchInternalField
         template<class Type>
         tmp<Field<Type> > interpolateField(const interpolation<Type>&) const;
 
+
 public:
 
     //- Runtime type information
@@ -97,67 +98,71 @@ public:
 
     // Member Functions
 
-        //- sample field on surface
+        // Sample
-        virtual tmp<scalarField> sample
-        (
-            const volScalarField&
-        ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
-        virtual tmp<vectorField> sample
+            virtual tmp<scalarField> sample
-        (
+            (
-            const volVectorField&
+                const volScalarField&
-        ) const;
+            ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
-        virtual tmp<sphericalTensorField> sample
+            virtual tmp<vectorField> sample
-        (
+            (
-            const volSphericalTensorField&
+                const volVectorField&
-        ) const;
+            ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
-        virtual tmp<symmTensorField> sample
+            virtual tmp<sphericalTensorField> sample
-        (
+            (
-            const volSymmTensorField&
+                const volSphericalTensorField&
-        ) const;
+            ) const;
 
-        //- sample field on surface
+            //- Sample field on surface
-        virtual tmp<tensorField> sample
+            virtual tmp<symmTensorField> sample
-        (
+            (
-            const volTensorField&
+                const volSymmTensorField&
-        ) const;
+            ) const;
+
+            //- Sample field on surface
+            virtual tmp<tensorField> sample
+            (
+                const volTensorField&
+            ) const;
 
 
-        //- interpolate field on surface
+        // Interpolate
-        virtual tmp<scalarField> interpolate
-        (
-            const interpolation<scalar>&
-        ) const;
 
+            //- Interpolate field on surface
+            virtual tmp<scalarField> interpolate
+            (
+                const interpolation<scalar>&
+            ) const;
 
-        //- interpolate field on surface
+            //- Interpolate field on surface
-        virtual tmp<vectorField> interpolate
+            virtual tmp<vectorField> interpolate
-        (
+            (
-            const interpolation<vector>&
+                const interpolation<vector>&
-        ) const;
+            ) const;
 
-        //- interpolate field on surface
+            //- Interpolate field on surface
-        virtual tmp<sphericalTensorField> interpolate
+            virtual tmp<sphericalTensorField> interpolate
-        (
+            (
-            const interpolation<sphericalTensor>&
+                const interpolation<sphericalTensor>&
-        ) const;
+            ) const;
 
-        //- interpolate field on surface
+            //- Interpolate field on surface
-        virtual tmp<symmTensorField> interpolate
+            virtual tmp<symmTensorField> interpolate
-        (
+            (
-            const interpolation<symmTensor>&
+                const interpolation<symmTensor>&
-        ) const;
+            ) const;
+
+            //- Interpolate field on surface
+            virtual tmp<tensorField> interpolate
+            (
+                const interpolation<tensor>&
+            ) const;
 
-        //- interpolate field on surface
-        virtual tmp<tensorField> interpolate
-        (
-            const interpolation<tensor>&
-        ) const;
 
         //- Write
         virtual void print(Ostream&) const;

src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C

@@ -113,9 +113,9 @@ void Foam::sampledSurfaces::sampleAndWrite
 )
 {
     // interpolator for this field
-    autoPtr< interpolation<Type> > interpolator;
+    autoPtr<interpolation<Type> > interpolatorPtr;
 
-    const word& fieldName   = vField.name();
+    const word& fieldName = vField.name();
     const fileName outputDir = outputPath_/vField.time().timeName();
 
     forAll(*this, surfI)
@@ -126,16 +126,16 @@ void Foam::sampledSurfaces::sampleAndWrite
 
         if (s.interpolate())
         {
-            if (interpolator.empty())
+            if (interpolatorPtr.empty())
             {
-                interpolator = interpolation<Type>::New
+                interpolatorPtr = interpolation<Type>::New
                 (
                     interpolationScheme_,
                     vField
                 );
             }
 
-            values = s.interpolate(interpolator());
+            values = s.interpolate(interpolatorPtr());
         }
         else
         {

src/thermophysicalModels/radiationModels/Make/files

@@ -27,6 +27,11 @@ submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionE
 submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
 submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C
 
+/* Soot model */
+submodels/sootModel/sootModel/sootModel.C
+submodels/sootModel/sootModel/sootModelNew.C
+submodels/sootModel/mixtureFractionSoot/mixtureFractionSoots.C
+submodels/sootModel/noSoot/noSoot.C
 
 /* Boundary conditions */
 derivedFvPatchFields/MarshakRadiation/MarshakRadiationFvPatchScalarField.C

src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C

@@ -337,6 +337,17 @@ Foam::radiation::greyMeanAbsorptionEmission::ECont(const label bandI) const
             }
         }
     }
+    else
+    {
+        WarningIn
+        (
+            "tmp<volScalarField>"
+            "radiation::greyMeanAbsorptionEmission::ECont"
+            "("
+                "const label"
+            ") const"
+        ) << "dQ field not found in mesh" << endl;
+    }
 
     return E;
 }

src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -255,14 +255,39 @@ Foam::radiation::wideBandAbsorptionEmission::ECont(const label bandI) const
         )
     );
 
-    if (mesh().foundObject<volScalarField>("hrr"))
+    if (mesh().foundObject<volScalarField>("dQ"))
     {
-        const volScalarField& hrr = mesh().lookupObject<volScalarField>("hrr");
+        const volScalarField& dQ = mesh().lookupObject<volScalarField>("dQ");
-        E().internalField() =
+
-            iEhrrCoeffs_[bandI]
+        if (dQ.dimensions() == dimEnergy/dimTime)
-           *hrr.internalField()
+        {
-           *(iBands_[bandI][1] - iBands_[bandI][0])
+            E().internalField() =
-           /totalWaveLength_;
+                iEhrrCoeffs_[bandI]
+               *dQ.internalField()
+               *(iBands_[bandI][1] - iBands_[bandI][0])
+               /totalWaveLength_
+               /mesh_.V();
+        }
+        else if (dQ.dimensions() == dimEnergy/dimTime/dimVolume)
+        {
+            E().internalField() =
+                iEhrrCoeffs_[bandI]
+               *dQ.internalField()
+               *(iBands_[bandI][1] - iBands_[bandI][0])
+               /totalWaveLength_;
+        }
+        else
+        {
+            WarningIn
+            (
+                "tmp<volScalarField>"
+                "radiation::wideBandAbsorptionEmission::ECont"
+                "("
+                    "const label"
+                ") const"
+            )
+                << "Incompatible dimensions for dQ field" << endl;
+        }
     }
 
     return E;
@@ -289,9 +314,8 @@ void Foam::radiation::wideBandAbsorptionEmission::correct
 
     for (label j=0; j<nBands_; j++)
     {
-        Info<< "Calculating absorption in band: " << j << endl;
         aLambda[j].internalField() = this->a(j);
-        Info<< "Calculated absorption in band: " << j << endl;
+
         a.internalField() +=
             aLambda[j].internalField()
            *(iBands_[j][1] - iBands_[j][0])

src/thermophysicalModels/radiationModels/submodels/sootModel/mixtureFractionSoot/mixtureFractionSoot.C

@@ -0,0 +1,166 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "mixtureFractionSoot.H"
+#include "singleStepReactingMixture.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+
+template<class ThermoType>
+const Foam::singleStepReactingMixture<ThermoType>&
+Foam::radiation::mixtureFractionSoot<ThermoType>::checkThermo
+(
+    const fluidThermo& thermo
+)
+{
+    if (isA<singleStepReactingMixture<ThermoType> >(thermo))
+    {
+        return dynamic_cast<const singleStepReactingMixture<ThermoType>& >
+        (
+            thermo
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "template<class ThermoType> "
+            "Foam::radiation::mixtureFractionSoot "
+            "("
+                "const dictionary&, "
+                "const fvMesh&"
+            ")"
+        )
+            << "Inconsistent thermo package for " << thermo.type()
+            << "Please select a thermo package based on "
+            << "singleStepReactingMixture" << exit(FatalError);
+
+        return dynamic_cast<const singleStepReactingMixture<ThermoType>& >
+        (
+            thermo
+        );
+    }
+
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ThermoType>
+Foam::radiation::mixtureFractionSoot<ThermoType>::mixtureFractionSoot
+(
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const word& modelType
+)
+:
+    sootModel(dict, mesh, modelType),
+    soot_
+    (
+        IOobject
+        (
+            "soot",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    coeffsDict_(dict.subOrEmptyDict(modelType + "Coeffs")),
+    nuSoot_(readScalar(coeffsDict_.lookup("nuSoot"))),
+    Wsoot_(readScalar(coeffsDict_.lookup("Wsoot"))),
+    sootMax_(-1),
+    mappingFieldName_
+    (
+        coeffsDict_.lookupOrDefault<word>("mappingFieldName", "none")
+    ),
+    mapFieldMax_(1),
+    thermo_(mesh.lookupObject<fluidThermo>("thermophysicalProperties")),
+    mixture_(checkThermo(thermo_))
+{
+    const Reaction<ThermoType>& reaction = mixture_.operator[](0);
+
+    const scalarList& specieStoichCoeffs(mixture_.specieStoichCoeffs());
+
+    scalar totalMol = 0.0;
+    forAll(reaction.rhs(), i)
+    {
+        label specieI = reaction.rhs()[i].index;
+        totalMol += mag(specieStoichCoeffs[specieI]);
+    }
+
+    totalMol += nuSoot_;
+
+    scalarList Xi(reaction.rhs().size());
+
+    scalar Wm = 0.0;
+    forAll(reaction.rhs(), i)
+    {
+        const label specieI = reaction.rhs()[i].index;
+        Xi[i] = mag(specieStoichCoeffs[specieI])/totalMol;
+        Wm += Xi[i]*mixture_.speciesData()[specieI].W();
+    }
+
+    const scalar XSoot = nuSoot_/totalMol;
+    Wm += XSoot*Wsoot_;
+
+    sootMax_ = XSoot*Wsoot_/Wm;
+
+    Info << "Maximum soot mass concentrations: " << sootMax_ << nl;
+
+    if (mappingFieldName_ == "none")
+    {
+        const label index = reaction.rhs()[0].index;
+        mappingFieldName_ = mixture_.Y(index).name();
+    }
+
+    const label mapFieldIndex = mixture_.species()[mappingFieldName_];
+
+    mapFieldMax_ = mixture_.Yprod0()[mapFieldIndex];
+
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class ThermoType>
+Foam::radiation::mixtureFractionSoot<ThermoType>::~mixtureFractionSoot()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ThermoType>
+void Foam::radiation::mixtureFractionSoot<ThermoType>::correct()
+{
+    const volScalarField& mapField =
+        mesh_.lookupObject<volScalarField>(mappingFieldName_);
+
+    soot_ = sootMax_*(mapField/mapFieldMax_);
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/radiationModels/submodels/sootModel/mixtureFractionSoot/mixtureFractionSoot.H

@@ -0,0 +1,179 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::radiation::mixtureFractionSoot
+
+Description
+    This soot model is purely an state model. The ammount of soot produced is
+    determined by a single step chemistry as :
+
+        nuf Fuel + nuOx Ox = nuP P + nuSoot soot
+
+    nuSoot is prescribed by the user.
+
+    The single step chemistry used is read from the combustion.
+    The soot is not considered into the thermodynamics of the system and it
+    is not considered as an extra specie in the solver.
+
+    The spacial distribution is given by the normalization of the first product
+    on the rhs of the reaction by default or it can be added as input.
+
+    The input dictionary reads like in the radiationProperties dictionary:
+
+    sootModel mixtureFractionSoot<gasHThermoPhysics>;
+
+    mixtureFractionSootCoeffs
+    {
+        nuSoot              0.015;
+        Wsoot               12;
+        mappingFieldName    P;
+    }
+
+SourceFiles
+    mixtureFractionSoot.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mixtureFractionSoot_H
+#define mixtureFractionSoot_H
+
+#include "interpolationLookUpTable.H"
+#include "sootModel.H"
+#include "HashTable.H"
+
+#include "fluidThermo.H"
+#include "singleStepReactingMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class mixtureFractionSoot Declaration
+\*---------------------------------------------------------------------------*/
+template<class ThermoType>
+class mixtureFractionSoot
+:
+    public sootModel
+{
+
+    // Static functions
+
+        //- Check mixture in thermo
+        static const singleStepReactingMixture<ThermoType>& checkThermo
+        (
+            const fluidThermo&
+        );
+
+
+    // Private data
+
+         //- Soot mass fraction
+        volScalarField soot_;
+
+        //- Soot model dictionary
+        dictionary coeffsDict_;
+
+        //- Soot yield
+        scalar nuSoot_;
+
+        //- Soot molecular weight
+        scalar Wsoot_;
+
+        //- Maximum soot mass concentration at stoichiometric
+        scalar sootMax_;
+
+        //- Name of the field mapping the soot
+        word mappingFieldName_;
+
+        //- Maximum value of the map field
+        scalar mapFieldMax_;
+
+        //- Thermo package
+        const fluidThermo& thermo_;
+
+        //- Auto Ptr to singleStepReactingMixture
+        const singleStepReactingMixture<ThermoType>& mixture_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("mixtureFractionSoot");
+
+
+    // Constructors
+
+        //- Construct from components
+        mixtureFractionSoot
+        (
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const word& modelType
+        );
+
+
+    //- Destructor
+    virtual ~mixtureFractionSoot();
+
+
+    // Member Functions
+
+        // Edit
+
+            //- Main update/correction routine
+            virtual void correct();
+
+
+        // Access
+
+            //- Return Ysoot
+            const volScalarField& soot() const
+            {
+                return soot_;
+            }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace radiation
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "mixtureFractionSoot.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiationModels/submodels/sootModel/mixtureFractionSoot/mixtureFractionSoots.C

@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+
+#include "mixtureFractionSoot.H"
+#include "makeSootTypes.H"
+#include "thermoPhysicsTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+makeSootTypesThermo(mixtureFractionSoot, gasHThermoPhysics);
+makeSootTypesThermo(mixtureFractionSoot, gasEThermoPhysics);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace radiation
+
+// ************************************************************************* //

src/thermophysicalModels/radiationModels/submodels/sootModel/noSoot/noSoot.C

@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "noSoot.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace radiation
+    {
+        defineTypeNameAndDebug(noSoot, 0);
+
+        addToRunTimeSelectionTable
+        (
+            sootModel,
+            noSoot,
+            dictionary
+        );
+    }
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::noSoot::noSoot
+(
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const word& modelType
+)
+:
+    sootModel(dict, mesh, modelType)
+{}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+Foam::radiation::noSoot::~noSoot()
+{}
+
+
+void Foam::radiation::noSoot::correct()
+{
+    //Do nothing
+}
+
+const Foam::volScalarField& Foam::radiation::noSoot::soot() const
+{
+    notImplemented
+    (
+        "Foam::volScalarField& Foam::radiation::noSoot::soot() const"
+    );
+    return tmp<volScalarField>(NULL);
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/radiationModels/submodels/sootModel/noSoot/noSoot.H

@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::radiation::noSoot
+
+Description
+    noSoot
+
+
+SourceFiles
+    noSoot.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noSoot_H
+#define noSoot_H
+
+
+#include "sootModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class noSoot Declaration
+\*---------------------------------------------------------------------------*/
+class noSoot
+:
+    public sootModel
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+
+    // Constructors
+
+        //- Construct from components
+        noSoot(const dictionary& dict, const fvMesh& mesh, const word&);
+
+
+    //- Destructor
+    virtual ~noSoot();
+
+
+    // Member Functions
+
+        // Edit
+
+            //- Main update/correction routine
+            void correct();
+
+        // Access
+
+            //- Return Ysoot
+            const volScalarField& soot() const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace radiation
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiationModels/submodels/sootModel/sootModel/makeSootTypes.H

@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeSootTypes_H
+#define makeSootTypes_H
+
+#include "addToRunTimeSelectionTable.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeSootTypesThermo(sootModelType, Thermo)                            \
+                                                                              \
+    typedef sootModelType<Thermo> sootModelType##Thermo;                      \
+                                                                              \
+    defineTemplateTypeNameAndDebugWithName                                    \
+    (                                                                         \
+        sootModelType##Thermo,                                                \
+        #sootModelType"<"#Thermo">",                                          \
+        0                                                                     \
+    );                                                                        \
+                                                                              \
+    addToRunTimeSelectionTable                                                \
+    (                                                                         \
+        sootModel,                                                            \
+        sootModelType##Thermo,                                                \
+        dictionary                                                            \
+    );
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiationModels/submodels/sootModel/sootModel/sootModel.C

@@ -0,0 +1,62 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sootModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+    defineTypeNameAndDebug(sootModel, 0);
+    defineRunTimeSelectionTable(sootModel, dictionary);
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::sootModel::sootModel
+(
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const word& modelType
+)
+:
+    dict_(dict),
+    mesh_(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor    * * * * * * * * * * * * * * //
+
+Foam::radiation::sootModel::~sootModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiationModels/submodels/sootModel/sootModel/sootModel.H

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::radiation::sootModel
+
+Description
+    Base class for soor models
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef radiationsootModel_H
+#define radiationsootModel_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "volFields.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace radiation
+{
+
+
+/*---------------------------------------------------------------------------*\
+                            Class sootModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class sootModel
+{
+
+protected:
+
+    // Protected data
+
+        //- Radiation model dictionary
+        const dictionary dict_;
+
+        //- Reference to the fvMesh
+        const fvMesh& mesh_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("sootModel");
+
+    //- Declare runtime constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            sootModel,
+            dictionary,
+            (
+                const dictionary& dict,
+                const fvMesh& mesh,
+                const word& modelType
+            ),
+            (dict, mesh, modelType)
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        sootModel
+        (
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const word& modelType
+        );
+
+
+    //- Selector
+    static autoPtr<sootModel> New
+    (
+        const dictionary& dict,
+        const fvMesh& mesh
+    );
+
+
+    //- Destructor
+    virtual ~sootModel();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Reference to the mesh
+            inline const fvMesh& mesh() const
+            {
+                return mesh_;
+            }
+
+            //- Reference to the dictionary
+            inline const dictionary& dict() const
+            {
+                return dict_;
+            }
+
+
+            // Member Functions
+
+                // Edit
+
+                    //- Main update/correction routine
+                    virtual void correct() = 0;
+
+
+                // Access
+
+                    //- Return const reference to soot
+                    virtual const volScalarField& soot() const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace radiation
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiationModels/submodels/sootModel/sootModel/sootModelNew.C

@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "sootModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::radiation::sootModel>
+Foam::radiation::sootModel::New
+(
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+{
+    word modelType("none");
+
+    if (dict.found("sootModel"))
+    {
+        dict.lookup("sootModel") >> modelType;
+
+        Info<< "Selecting sootModel " << modelType << endl;
+    }
+
+    dictionaryConstructorTable::iterator cstrIter =
+            dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "sootModel::New(const dictionary&, const fvMesh&)"
+        )   << "Unknown sootModel type "
+            << modelType << nl << nl
+            << "Valid sootModel types are :" << nl
+            << dictionaryConstructorTablePtr_->sortedToc() << exit(FatalError);
+    }
+
+    const label tempOpen = modelType.find('<');
+
+    const word className = modelType(0, tempOpen);
+
+    return autoPtr<sootModel>(cstrIter()(dict, mesh, className));
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C

@@ -281,9 +281,8 @@ void Foam::chemkinReader::addPressureDependentReaction
                     thirdBodyEfficiencies(speciesTable_, efficiencies)
                 )
             );
+            break;
         }
-        break;
-
         case Troe:
         {
             scalarList TroeCoeffs
@@ -337,9 +336,8 @@ void Foam::chemkinReader::addPressureDependentReaction
                     thirdBodyEfficiencies(speciesTable_, efficiencies)
                 )
             );
+            break;
         }
-        break;
-
         case SRI:
         {
             scalarList SRICoeffs
@@ -395,27 +393,16 @@ void Foam::chemkinReader::addPressureDependentReaction
                     thirdBodyEfficiencies(speciesTable_, efficiencies)
                 )
             );
+            break;
         }
-        break;
-
         default:
         {
-            if (fofType < 4)
+            FatalErrorIn("chemkinReader::addPressureDependentReaction")
-            {
+                << "Fall-off function type "
-                FatalErrorIn("chemkinReader::addPressureDependentReaction")
+                << fallOffFunctionNames[fofType]
-                    << "Fall-off function type "
+                << " on line " << lineNo_-1
-                    << fallOffFunctionNames[fofType]
+                << " not implemented"
-                    << " on line " << lineNo_-1
+                << exit(FatalError);
-                    << " not implemented"
-                    << exit(FatalError);
-            }
-            else
-            {
-                FatalErrorIn("chemkinReader::addPressureDependentReaction")
-                    << "Unknown fall-off function type " << fofType
-                    << " on line " << lineNo_-1
-                    << exit(FatalError);
-            }
         }
     }
 }
@@ -537,9 +524,8 @@ void Foam::chemkinReader::addReaction
                     )
                 );
             }
+            break;
         }
-        break;
-
         case thirdBodyArrhenius:
         {
             if (rType == nonEquilibriumReversible)
@@ -597,9 +583,8 @@ void Foam::chemkinReader::addReaction
                     )
                 );
             }
+            break;
         }
-        break;
-
         case unimolecularFallOff:
         {
             addPressureDependentReaction<FallOffReactionRate>
@@ -616,9 +601,8 @@ void Foam::chemkinReader::addReaction
                 Afactor,
                 RR
             );
+            break;
         }
-        break;
-
         case chemicallyActivatedBimolecular:
         {
             addPressureDependentReaction<ChemicallyActivatedReactionRate>
@@ -635,9 +619,8 @@ void Foam::chemkinReader::addReaction
                 Afactor/concFactor,
                 RR
             );
+            break;
         }
-        break;
-
         case LandauTeller:
         {
             const scalarList& LandauTellerCoeffs =
@@ -705,9 +688,8 @@ void Foam::chemkinReader::addReaction
                     )
                 );
             }
+            break;
         }
-        break;
-
         case Janev:
         {
             const scalarList& JanevCoeffs =
@@ -727,9 +709,8 @@ void Foam::chemkinReader::addReaction
                     FixedList<scalar, 9>(JanevCoeffs)
                 )
             );
+            break;
         }
-        break;
-
         case powerSeries:
         {
             const scalarList& powerSeriesCoeffs =
@@ -749,35 +730,23 @@ void Foam::chemkinReader::addReaction
                     FixedList<scalar, 4>(powerSeriesCoeffs)
                 )
             );
+            break;
         }
-        break;
-
         case unknownReactionRateType:
         {
             FatalErrorIn("chemkinReader::addReaction")
                 << "Internal error on line " << lineNo_-1
                 << ": reaction rate type has not been set"
                 << exit(FatalError);
+            break;
         }
-        break;
-
         default:
         {
-            if (rrType < 9)
+            FatalErrorIn("chemkinReader::addReaction")
-            {
+                << "Reaction rate type " << reactionRateTypeNames[rrType]
-                FatalErrorIn("chemkinReader::addReaction")
+                << " on line " << lineNo_-1
-                    << "Reaction rate type " << reactionRateTypeNames[rrType]
+                << " not implemented"
-                    << " on line " << lineNo_-1
+                << exit(FatalError);
-                    << " not implemented"
-                    << exit(FatalError);
-            }
-            else
-            {
-                FatalErrorIn("chemkinReader::addReaction")
-                    << "Unknown reaction rate type " << rrType
-                    << " on line " << lineNo_-1
-                    << exit(FatalError);
-            }
         }
     }