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

src/OpenFOAM/db/functionObjects/outputFilterOutputControl/outputFilterOutputControl.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,11 +30,8 @@ License
 namespace Foam
 {
     template<>
-    const char* Foam::NamedEnum
+    const char* NamedEnum<outputFilterOutputControl::outputControls, 2>::
-    <
+    names[] =
-        Foam::outputFilterOutputControl::outputControls,
-        2
-    >::names[] =
     {
         "timeStep",
         "outputTime"

src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.C

@@ -41,8 +41,18 @@ namespace Foam
         fixedTemperatureSource,
         dictionary
     );
+
+    template<>
+    const char* NamedEnum<fixedTemperatureSource::temperatureMode, 2>::names[] =
+    {
+        "constant",
+        "lookup"
+    };
 }
 
+const Foam::NamedEnum<Foam::fixedTemperatureSource::temperatureMode, 2>
+    Foam::fixedTemperatureSource::temperatureModeNames_;
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -55,8 +65,29 @@ Foam::fixedTemperatureSource::fixedTemperatureSource
 )
 :
     basicSource(name, modelType, dict, mesh),
-    T_(readScalar(coeffs_.lookup("temperature")))
+    mode_(temperatureModeNames_.read(coeffs_.lookup("mode"))),
+    Tconstant_(0.0),
+    TName_("T")
 {
+    switch (mode_)
+    {
+        case tmConstant:
+        {
+            coeffs_.lookup("temperature") >> Tconstant_;
+            break;
+        }
+        case tmLookup:
+        {
+            TName_ = coeffs_.lookupOrDefault<word>("TName", "T");
+            break;
+        }
+        default:
+        {
+            // error handling done by NamedEnum
+        }
+    }
+
+
     fieldNames_.setSize(1, "energy");
     applied_.setSize(1, false);
 }
@@ -81,8 +112,31 @@ void Foam::fixedTemperatureSource::setValue
 
     if (eqn.psi().name() == thermo.he().name())
     {
-        const scalarField Tfield(cells_.size(), T_);
+        switch (mode_)
-        eqn.setValues(cells_, thermo.he(thermo.p(), Tfield, cells_));
+        {
+            case tmConstant:
+            {
+                scalarField Tconst(cells_.size(), Tconstant_);
+                eqn.setValues(cells_, thermo.he(thermo.p(), Tconst, cells_));
+
+                break;
+            }
+            case tmLookup:
+            {
+                const volScalarField& T =
+                    mesh().lookupObject<volScalarField>(TName_);
+
+                scalarField Tlookup(T, cells_);
+                eqn.setValues(cells_, thermo.he(thermo.p(), Tlookup, cells_));
+
+                break;
+            }
+            default:
+            {
+                // error handling done by NamedEnum
+            }
+        }
+
     }
 }
 
@@ -98,7 +152,8 @@ bool Foam::fixedTemperatureSource::read(const dictionary& dict)
 {
     if (basicSource::read(dict))
     {
-        coeffs_.readIfPresent("T", T_);
+        coeffs_.readIfPresent("temperature", Tconstant_);
+        coeffs_.readIfPresent("TName", TName_);
 
         return true;
     }

src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.H

@@ -33,7 +33,13 @@ Description
         fixedTemperatureSourceCoeffs
         {
             fieldNames      (h e hs);   // names of fields to apply source
+            mode            constant;   // constant or lookup
+
+            // constant option
             temperature     500;        // fixed temperature [K]
+
+            // loolup option
+            // TName        T;          // optional temperature field name
         }
 
 
@@ -46,6 +52,7 @@ SourceFiles
 #define fixedTemperatureSource_H
 
 #include "basicSource.H"
+#include "NamedEnum.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,12 +68,32 @@ class fixedTemperatureSource
     public basicSource
 {
 
+public:
+
+    //- Temperature mode
+    enum temperatureMode
+    {
+        tmConstant,
+        tmLookup
+    };
+
+
+    //- String representation of temperatureMode enums
+    static const NamedEnum<temperatureMode, 2> temperatureModeNames_;
+
+
 protected:
 
     // Protected data
 
-        //- Fixed temperature [K]
+        //- Operation mode
-        scalar T_;
+        temperatureMode mode_;
+
+        //- Constant temperature [K]
+        scalar Tconstant_;
+
+        //- Temperature field name
+        word TName_;
 
 
 private:

src/finiteVolume/fields/fvPatchFields/derived/fixedJump/fixedJumpFvPatchField.C

@@ -63,7 +63,19 @@ Foam::fixedJumpFvPatchField<Type>::fixedJumpFvPatchField
 :
     jumpCyclicFvPatchField<Type>(p, iF),
     jump_("jump", dict, p.size())
-{}
+{
+    if (dict.found("value"))
+    {
+        fvPatchField<Type>::operator=
+        (
+            Field<Type>("value", dict, p.size())
+        );
+    }
+    else
+    {
+        this->evaluate(Pstream::blocking);
+    }
+}
 
 
 template<class Type>

src/finiteVolume/fields/fvPatchFields/derived/fixedJumpAMI/fixedJumpAMIFvPatchField.C

@@ -63,7 +63,19 @@ Foam::fixedJumpAMIFvPatchField<Type>::fixedJumpAMIFvPatchField
 :
     jumpCyclicAMIFvPatchField<Type>(p, iF),
     jump_("jump", dict, p.size())
-{}
+{
+    if (dict.found("value"))
+    {
+        fvPatchField<Type>::operator=
+        (
+            Field<Type>("value", dict, p.size())
+        );
+    }
+    else
+    {
+        this->evaluate(Pstream::blocking);
+    }
+}
 
 
 template<class Type>

src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C

@@ -53,7 +53,7 @@ timeVaryingMappedFixedValueFvPatchField
     endSampleTime_(-1),
     endSampledValues_(0),
     endAverage_(pTraits<Type>::zero),
-    offSet_()
+    offset_()
 {}
 
 
@@ -79,7 +79,7 @@ timeVaryingMappedFixedValueFvPatchField
     endSampleTime_(-1),
     endSampledValues_(0),
     endAverage_(pTraits<Type>::zero),
-    offSet_()
+    offset_()
 {}
 
 
@@ -104,7 +104,7 @@ timeVaryingMappedFixedValueFvPatchField
     endSampleTime_(-1),
     endSampledValues_(0),
     endAverage_(pTraits<Type>::zero),
-    offSet_(DataEntry<Type>::New("offSet", dict))
+    offset_(DataEntry<Type>::New("offset", dict))
 {
     dict.readIfPresent("fieldTableName", fieldTableName_);
 
@@ -138,7 +138,7 @@ timeVaryingMappedFixedValueFvPatchField
     endSampleTime_(ptf.endSampleTime_),
     endSampledValues_(ptf.endSampledValues_),
     endAverage_(ptf.endAverage_),
-    offSet_(ptf.offSet_().clone().ptr())
+    offset_(ptf.offset_().clone().ptr())
 {}
 
 
@@ -163,7 +163,7 @@ timeVaryingMappedFixedValueFvPatchField
     endSampleTime_(ptf.endSampleTime_),
     endSampledValues_(ptf.endSampledValues_),
     endAverage_(ptf.endAverage_),
-    offSet_(ptf.offSet_().clone().ptr())
+    offset_(ptf.offset_().clone().ptr())
 {}
 
 
@@ -487,7 +487,7 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::updateCoeffs()
 
     // apply offset to mapped values
     const scalar t = this->db().time().timeOutputValue();
-    this->operator==(*this + offSet_->value(t));
+    this->operator==(*this + offset_->value(t));
 
     if (debug)
     {
@@ -512,7 +512,7 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::write(Ostream& os) const
             << token::END_STATEMENT << nl;
     }
 
-    offSet_->writeData(os);
+    offset_->writeData(os);
 
     this->writeEntry("value", os);
 }

src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H

@@ -132,7 +132,7 @@ class timeVaryingMappedFixedValueFvPatchField
         Type endAverage_;
 
         //- Time varying offset values to interpolated data
-        autoPtr<DataEntry<Type> > offSet_;
+        autoPtr<DataEntry<Type> > offset_;
 
 
 public:

src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.C

@@ -73,6 +73,10 @@ Foam::uniformJumpAMIFvPatchField<Type>::uniformJumpAMIFvPatchField
     {
         fvPatchField<Type>::operator=(Field<Type>("value", dict, p.size()));
     }
+    else
+    {
+        this->evaluate(Pstream::blocking);
+    }
 }
 
 

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C

@@ -326,7 +326,7 @@ bool Foam::KinematicParcel<ParcelType>::move
 
         p.age() += dt;
 
-        td.cloud().functions().postMove(p, cellI, dt);
+        td.cloud().functions().postMove(p, cellI, dt, td.keepParticle);
     }
 
     return td.keepParticle;
@@ -340,7 +340,7 @@ void Foam::KinematicParcel<ParcelType>::hitFace(TrackData& td)
     typename TrackData::cloudType::parcelType& p =
         static_cast<typename TrackData::cloudType::parcelType&>(*this);
 
-    td.cloud().functions().postFace(p, p.face());
+    td.cloud().functions().postFace(p, p.face(), td.keepParticle);
 }
 
 
@@ -364,7 +364,14 @@ bool Foam::KinematicParcel<ParcelType>::hitPatch
         static_cast<typename TrackData::cloudType::parcelType&>(*this);
 
     // Invoke post-processing model
-    td.cloud().functions().postPatch(p, pp, trackFraction, tetIs);
+    td.cloud().functions().postPatch
+    (
+        p,
+        pp,
+        trackFraction,
+        tetIs,
+        td.keepParticle
+    );
 
     // Invoke surface film model
     if (td.cloud().surfaceFilm().transferParcel(p, pp, td.keepParticle))

src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,6 +29,7 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "FacePostProcessing.H"
+#include "ParticleCollector.H"
 #include "ParticleErosion.H"
 #include "ParticleTracks.H"
 #include "ParticleTrap.H"
@@ -42,6 +43,7 @@ License
     makeCloudFunctionObject(CloudType);                                       \
                                                                               \
     makeCloudFunctionObjectType(FacePostProcessing, CloudType);               \
+    makeCloudFunctionObjectType(ParticleCollector, CloudType);                \
     makeCloudFunctionObjectType(ParticleErosion, CloudType);                  \
     makeCloudFunctionObjectType(ParticleTracks, CloudType);                   \
     makeCloudFunctionObjectType(ParticleTrap, CloudType);                     \

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

@@ -96,7 +96,8 @@ void Foam::CloudFunctionObject<CloudType>::postMove
 (
     const typename CloudType::parcelType&,
     const label,
-    const scalar
+    const scalar,
+    bool&
 )
 {
     // do nothing
@@ -109,7 +110,8 @@ void Foam::CloudFunctionObject<CloudType>::postPatch
     const typename CloudType::parcelType&,
     const polyPatch&,
     const scalar,
-    const tetIndices&
+    const tetIndices&,
+    bool&
 )
 {
     // do nothing
@@ -120,7 +122,8 @@ template<class CloudType>
 void Foam::CloudFunctionObject<CloudType>::postFace
 (
     const typename CloudType::parcelType&,
-    const label
+    const label,
+    bool&
 )
 {
     // do nothing

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

@@ -46,6 +46,9 @@ SourceFiles
 namespace Foam
 {
 
+class polyPatch;
+class tetIndices;
+
 /*---------------------------------------------------------------------------*\
                     Class CloudFunctionObject Declaration
 \*---------------------------------------------------------------------------*/
@@ -134,7 +137,8 @@ public:
             (
                 const typename CloudType::parcelType& p,
                 const label cellI,
-                const scalar dt
+                const scalar dt,
+                bool& keepParticle
             );
 
             //- Post-patch hook
@@ -143,14 +147,16 @@ public:
                 const typename CloudType::parcelType& p,
                 const polyPatch& pp,
                 const scalar trackFraction,
-                const tetIndices& testIs
+                const tetIndices& testIs,
+                bool& keepParticle
             );
 
             //- Post-face hook
             virtual void postFace
             (
                 const typename CloudType::parcelType& p,
-                const label faceI
+                const label faceI,
+                bool& keepParticle
             );
 };
 

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

@@ -132,12 +132,18 @@ void Foam::CloudFunctionObjectList<CloudType>::postMove
 (
     const typename CloudType::parcelType& p,
     const label cellI,
-    const scalar dt
+    const scalar dt,
+    bool& keepParticle
 )
 {
     forAll(*this, i)
     {
-        this->operator[](i).postMove(p, cellI, dt);
+        this->operator[](i).postMove(p, cellI, dt, keepParticle);
+
+        if (!keepParticle)
+        {
+            return;
+        }
     }
 }
 
@@ -148,12 +154,25 @@ void Foam::CloudFunctionObjectList<CloudType>::postPatch
     const typename CloudType::parcelType& p,
     const polyPatch& pp,
     const scalar trackFraction,
-    const tetIndices& tetIs
+    const tetIndices& tetIs,
+    bool& keepParticle
 )
 {
     forAll(*this, i)
     {
-        this->operator[](i).postPatch(p, pp, trackFraction, tetIs);
+        this->operator[](i).postPatch
+        (
+            p,
+            pp,
+            trackFraction,
+            tetIs,
+            keepParticle
+        );
+
+        if (!keepParticle)
+        {
+            return;
+        }
     }
 }
 
@@ -162,12 +181,18 @@ template<class CloudType>
 void Foam::CloudFunctionObjectList<CloudType>::postFace
 (
     const typename CloudType::parcelType& p,
-    const label faceI
+    const label faceI,
+    bool& keepParticle
 )
 {
     forAll(*this, i)
     {
-        this->operator[](i).postFace(p, faceI);
+        this->operator[](i).postFace(p, faceI, keepParticle);
+
+        if (!keepParticle)
+        {
+            return;
+        }
     }
 }
 

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

@@ -114,7 +114,8 @@ public:
             (
                 const typename CloudType::parcelType& p,
                 const label cellI,
-                const scalar dt
+                const scalar dt,
+                bool& keepParticle
             );
 
             //- Post-patch hook
@@ -123,14 +124,16 @@ public:
                 const typename CloudType::parcelType& p,
                 const polyPatch& pp,
                 const scalar trackFraction,
-                const tetIndices& tetIs
+                const tetIndices& tetIs,
+                bool& keepParticle
             );
 
             //- Post-face hook
             virtual void postFace
             (
                 const typename CloudType::parcelType& p,
-                const label faceI
+                const label faceI,
+                bool& keepParticle
             );
 };
 

src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -133,7 +133,7 @@ void Foam::FacePostProcessing<CloudType>::write()
         {
             OFstream& os = outputFilePtr_[zoneI];
             os  << time.timeName() << token::TAB << sumMassTotal << token::TAB
-                <<  sumMassFlowRate<< endl;
+                << sumMassFlowRate<< endl;
         }
     }
 
@@ -376,7 +376,8 @@ template<class CloudType>
 void Foam::FacePostProcessing<CloudType>::postFace
 (
     const parcelType& p,
-    const label faceI
+    const label faceI,
+    bool&
 )
 {
     if

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

@@ -160,7 +160,8 @@ public:
             virtual void postFace
             (
                 const parcelType& p,
-                const label faceI
+                const label faceI,
+                bool& keepParticle
             );
 };
 

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

@@ -0,0 +1,677 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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 "ParticleCollector.H"
+#include "Pstream.H"
+#include "surfaceWriter.H"
+#include "unitConversion.H"
+#include "Random.H"
+#include "triangle.H"
+#include "cloud.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::ParticleCollector<CloudType>::makeLogFile
+(
+    const faceList& faces,
+    const Field<point>& points,
+    const Field<scalar>& area
+)
+{
+    // Create the output file if not already created
+    if (log_)
+    {
+        if (debug)
+        {
+            Info<< "Creating output file" << endl;
+        }
+
+        if (Pstream::master())
+        {
+            const fileName logDir = outputDir_/this->owner().time().timeName();
+
+            // Create directory if does not exist
+            mkDir(logDir);
+
+            // Open new file at start up
+            outputFilePtr_.reset
+            (
+                new OFstream(logDir/(type() + ".dat"))
+            );
+
+            outputFilePtr_()
+                << "# Source     : " << type() << nl
+                << "# Total area : " << sum(area) << nl
+                << "# Time";
+
+            forAll(faces, i)
+            {
+                word id = Foam::name(i);
+
+                outputFilePtr_()
+                    << tab << "area[" << id << "]"
+                    << tab << "mass[" << id << "]"
+                    << tab << "massFlowRate[" << id << "]"
+                    << endl;
+            }
+        }
+    }
+}
+
+
+template<class CloudType>
+void Foam::ParticleCollector<CloudType>::initPolygons()
+{
+    mode_ = mtPolygon;
+
+    List<Field<point> > polygons(this->coeffDict().lookup("polygons"));
+    label nPoints = 0;
+    forAll(polygons, polyI)
+    {
+        label np = polygons[polyI].size();
+        if (np < 3)
+        {
+            FatalIOErrorIn
+            (
+                "Foam::ParticleCollector<CloudType>::initPolygons()",
+                this->coeffDict()
+            )
+                << "polygons must consist of at least 3 points"
+                << exit(FatalIOError);
+        }
+
+        nPoints += np;
+    }
+
+    label pointOffset = 0;
+    points_.setSize(nPoints);
+    faces_.setSize(polygons.size());
+    faceTris_.setSize(polygons.size());
+    area_.setSize(polygons.size());
+    forAll(faces_, faceI)
+    {
+        const Field<point>& polyPoints = polygons[faceI];
+        face f(identity(polyPoints.size()) + pointOffset);
+        UIndirectList<point>(points_, f) = polyPoints;
+        area_[faceI] = f.mag(points_);
+
+        DynamicList<face> tris;
+        f.triangles(points_, tris);
+        faceTris_[faceI].transfer(tris);
+
+        faces_[faceI].transfer(f);
+
+        pointOffset += polyPoints.size();
+    }
+}
+
+
+template<class CloudType>
+void Foam::ParticleCollector<CloudType>::initConcentricCircles()
+{
+    mode_ = mtConcentricCircle;
+
+    vector origin(this->coeffDict().lookup("origin"));
+
+    radius_ = this->coeffDict().lookup("radius");
+    nSector_ = readLabel(this->coeffDict().lookup("nSector"));
+ 
+    label nS = nSector_;
+
+    vector refDir;
+    if (nSector_ > 1)
+    {
+        refDir = this->coeffDict().lookup("refDir");
+        refDir -= normal_*(normal_ & refDir);
+        refDir /= mag(refDir);
+    }
+    else
+    {
+        // set 4 quadrants for single sector cases
+        nS = 4;
+
+        vector tangent = vector::zero;
+        scalar magTangent = 0.0;
+
+        Random rnd(1234);
+        while (magTangent < SMALL)
+        {
+            vector v = rnd.vector01();
+
+            tangent = v - (v & normal_)*normal_;
+            magTangent = mag(tangent);
+        }
+
+        refDir = tangent/magTangent;
+    }
+
+    scalar dTheta = 5.0;
+    scalar dThetaSector = 360.0/scalar(nS);
+    label intervalPerSector = max(1, ceil(dThetaSector/dTheta));
+    dTheta = dThetaSector/scalar(intervalPerSector);
+
+    label nPointPerSector = intervalPerSector + 1;
+
+    label nPointPerRadius = nS*(nPointPerSector - 1);
+    label nPoint = radius_.size()*nPointPerRadius;
+    label nFace = radius_.size()*nS;
+
+    // add origin
+    nPoint++;
+
+    points_.setSize(nPoint);
+    faces_.setSize(nFace);
+    area_.setSize(nFace);
+
+    coordSys_ = cylindricalCS("coordSys", origin, normal_, refDir, false);
+
+    List<label> ptIDs(identity(nPointPerRadius));
+
+    points_[0] = origin;
+
+    // points
+    forAll(radius_, radI)
+    {
+        label pointOffset = radI*nPointPerRadius + 1;
+
+        for (label i = 0; i < nPointPerRadius; i++)
+        {
+            label pI = i + pointOffset;
+            point pCyl(radius_[radI], degToRad(i*dTheta), 0.0);
+            points_[pI] = coordSys_.globalPosition(pCyl);
+        }
+    }
+
+    // faces
+    DynamicList<label> facePts(2*nPointPerSector);
+    forAll(radius_, radI)
+    {
+        if (radI == 0)
+        {
+            for (label secI = 0; secI < nS; secI++)
+            {
+                facePts.clear();
+
+                // append origin point
+                facePts.append(0);
+
+                for (label ptI = 0; ptI < nPointPerSector; ptI++)
+                {
+                    label i = ptI + secI*(nPointPerSector - 1);
+                    label id = ptIDs.fcIndex(i - 1) + 1;
+                    facePts.append(id);
+                }
+
+                label faceI = secI + radI*nS;
+
+                faces_[faceI] = face(facePts);
+                area_[faceI] = faces_[faceI].mag(points_);
+            }
+        }
+        else
+        {
+            for (label secI = 0; secI < nS; secI++)
+            {
+                facePts.clear();
+
+                label offset = (radI - 1)*nPointPerRadius + 1;
+
+                for (label ptI = 0; ptI < nPointPerSector; ptI++)
+                {
+                    label i = ptI + secI*(nPointPerSector - 1);
+                    label id = offset + ptIDs.fcIndex(i - 1);
+                    facePts.append(id);
+                }
+                for (label ptI = nPointPerSector-1; ptI >= 0; ptI--)
+                {
+                    label i = ptI + secI*(nPointPerSector - 1);
+                    label id = offset + nPointPerRadius + ptIDs.fcIndex(i - 1);
+                    facePts.append(id);
+                }
+
+                label faceI = secI + radI*nS;
+
+                faces_[faceI] = face(facePts);
+                area_[faceI] = faces_[faceI].mag(points_);
+            }
+        }
+    }
+}
+
+
+template<class CloudType>
+Foam::label Foam::ParticleCollector<CloudType>::collectParcelPolygon
+(
+    const point& position,
+    const vector& U
+) const
+{
+    scalar dt = this->owner().db().time().deltaTValue();
+
+    point end(position + dt*U);
+
+    label dummyNearType = -1;
+    label dummyNearLabel = -1;
+
+    forAll(faces_, faceI)
+    {
+        const label facePoint0 = faces_[faceI][0];
+
+        const point p0 = points_[facePoint0];
+
+        const scalar d1 = normal_ & (position - p0);
+        const scalar d2 = normal_ & (end - p0);
+
+        if (sign(d1) == sign(d2))
+        {
+            // did not cross polygon plane
+            continue;
+        }
+
+        // intersection point
+        point pIntersect = position + (d1/(d1 - d2))*dt*U;
+
+        const List<face>& tris = faceTris_[faceI];
+
+        // identify if point is within poly bounds
+        forAll(tris, triI)
+        {
+            const face& tri = tris[triI];
+            triPointRef t
+            (
+                points_[tri[0]],
+                points_[tri[1]],
+                points_[tri[2]]
+            );
+
+            if (t.classify(pIntersect, dummyNearType, dummyNearLabel))
+            {
+                return faceI;
+            }
+        }
+    }
+
+    return -1;
+}
+
+
+template<class CloudType>
+Foam::label Foam::ParticleCollector<CloudType>::collectParcelConcentricCircles
+(
+    const point& position,
+    const vector& U
+) const
+{
+    label secI = -1;
+
+    scalar dt = this->owner().db().time().deltaTValue();
+
+    point end(position + dt*U);
+
+    const scalar d1 = normal_ & (position - coordSys_.origin());
+    const scalar d2 = normal_ & (end - coordSys_.origin());
+
+    if (sign(d1) == sign(d2))
+    {
+        // did not cross plane
+        return secI;
+    }
+
+    // intersection point in cylindrical co-ordinate system
+    point pCyl = coordSys_.localPosition(position + (d1/(d1 - d2))*dt*U);
+
+    scalar r = pCyl[0];
+
+    if (r < radius_.last())
+    {
+        label radI = 0;
+        while (r > radius_[radI])
+        {
+            radI++;
+        }
+
+        if (nSector_ == 1)
+        {
+            secI = 4*radI;
+        }
+        else
+        {
+            scalar theta = pCyl[1] + constant::mathematical::pi;
+
+            secI =
+                nSector_*radI
+              + floor
+                (
+                    scalar(nSector_)*theta/constant::mathematical::twoPi
+                );
+        }
+    }
+
+    return secI;
+}
+
+
+template<class CloudType>
+void Foam::ParticleCollector<CloudType>::write()
+{
+    const fvMesh& mesh = this->owner().mesh();
+    const Time& time = mesh.time();
+    scalar timeNew = time.value();
+    scalar timeElapsed = timeNew - timeOld_;
+
+    totalTime_ += timeElapsed;
+
+    const scalar alpha = (totalTime_ - timeElapsed)/totalTime_;
+    const scalar beta = timeElapsed/totalTime_;
+
+    forAll(faces_, faceI)
+    {
+        massFlowRate_[faceI] =
+            alpha*massFlowRate_[faceI] + beta*mass_[faceI]/timeElapsed;
+        massTotal_[faceI] += mass_[faceI];
+    }
+
+    const label procI = Pstream::myProcNo();
+
+    Info<< type() << " output:" << nl;
+
+    if (outputFilePtr_.valid())
+    {
+        outputFilePtr_() << time.timeName();
+    }
+
+
+    Field<scalar> faceMassTotal(mass_.size());
+    Field<scalar> faceMassFlowRate(massFlowRate_.size());
+    
+    forAll(faces_, faceI)
+    {
+        scalarList allProcMass(Pstream::nProcs());
+        allProcMass[procI] = massTotal_[faceI];
+        Pstream::gatherList(allProcMass);
+        faceMassTotal[faceI] = sum(allProcMass);
+
+        scalarList allProcMassFlowRate(Pstream::nProcs());
+        allProcMassFlowRate[procI] = massFlowRate_[faceI];
+        Pstream::gatherList(allProcMassFlowRate);
+        faceMassFlowRate[faceI] = sum(allProcMassFlowRate);
+
+        Info<< "    face " << faceI
+            << ": total mass = " << faceMassTotal[faceI]
+            << "; average mass flow rate = " << faceMassFlowRate[faceI]
+            << nl;
+
+        if (outputFilePtr_.valid())
+        {
+            outputFilePtr_()
+                << tab << area_[faceI]
+                << tab << faceMassTotal[faceI]
+                << tab << faceMassFlowRate[faceI]
+                << endl;
+        }
+    }
+
+    Info<< endl;
+
+
+    if (surfaceFormat_ != "none")
+    {
+        if (Pstream::master())
+        {
+            autoPtr<surfaceWriter> writer(surfaceWriter::New(surfaceFormat_));
+
+            writer->write
+            (
+                outputDir_/time.timeName(),
+                "collector",
+                points_,
+                faces_,
+                "massTotal",
+                faceMassTotal,
+                false
+            );
+
+            writer->write
+            (
+                outputDir_/time.timeName(),
+                "collector",
+                points_,
+                faces_,
+                "massFlowRate",
+                faceMassFlowRate,
+                false
+            );
+        }
+    }
+
+
+    if (resetOnWrite_)
+    {
+        forAll(faces_, faceI)
+        {
+            massFlowRate_[faceI] = 0.0;
+        }
+        timeOld_ = timeNew;
+        totalTime_ = 0.0;
+    }
+
+    forAll(faces_, faceI)
+    {
+        mass_[faceI] = 0.0;
+    }
+
+    // writeProperties();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::ParticleCollector<CloudType>::ParticleCollector
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    CloudFunctionObject<CloudType>(dict, owner, typeName),
+    mode_(mtUnknown),
+    parcelType_(this->coeffDict().lookupOrDefault("parcelType", -1)),
+    removeCollected_(this->coeffDict().lookup("removeCollected")),
+    points_(),
+    faces_(),
+    faceTris_(),
+    nSector_(0),
+    radius_(),
+    coordSys_(false),
+    normal_(this->coeffDict().lookup("normal")),
+    negateParcelsOppositeNormal_
+    (
+        readBool(this->coeffDict().lookup("negateParcelsOppositeNormal"))
+    ),
+    surfaceFormat_(this->coeffDict().lookup("surfaceFormat")),
+    resetOnWrite_(this->coeffDict().lookup("resetOnWrite")),
+    totalTime_(0.0),
+    mass_(),
+    massTotal_(),
+    massFlowRate_(),
+    log_(this->coeffDict().lookup("log")),
+    outputFilePtr_(),
+    outputDir_(owner.mesh().time().path()),
+    timeOld_(owner.mesh().time().value())
+{
+    if (Pstream::parRun())
+    {
+        // Put in undecomposed case (Note: gives problems for
+        // distributed data running)
+        outputDir_ =
+            outputDir_/".."/"postProcessing"/cloud::prefix/owner.name();
+    }
+    else
+    {
+        outputDir_ =
+            outputDir_/"postProcessing"/cloud::prefix/owner.name();
+    }
+
+    normal_ /= mag(normal_);
+
+    word mode(this->coeffDict().lookup("mode"));
+    if (mode == "polygon")
+    {
+        initPolygons();
+    }
+    else if (mode == "concentricCircle")
+    {
+        initConcentricCircles();
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "Foam::ParticleCollector<CloudType>::ParticleCollector"
+            "("
+                "const dictionary& dict,"
+                "CloudType& owner"
+            ")"
+        )
+            << "Unknown mode " << mode << ".  Available options are "
+            << "polygon and concentricCircle" << exit(FatalError);
+    }
+
+    mass_.setSize(faces_.size(), 0.0);
+    massTotal_.setSize(faces_.size(), 0.0);
+    massFlowRate_.setSize(faces_.size(), 0.0);
+
+    makeLogFile(faces_, points_, area_);
+
+    // readProperties(); AND initialise mass... fields
+}
+
+
+template<class CloudType>
+Foam::ParticleCollector<CloudType>::ParticleCollector
+(
+    const ParticleCollector<CloudType>& pc
+)
+:
+    CloudFunctionObject<CloudType>(pc),
+    mode_(pc.mode_),
+    parcelType_(pc.parcelType_),
+    points_(pc.points_),
+    faces_(pc.faces_),
+    faceTris_(pc.faceTris_),
+    nSector_(pc.nSector_),
+    radius_(pc.radius_),
+    coordSys_(pc.coordSys_),
+    normal_(pc.normal_),
+    negateParcelsOppositeNormal_(pc.negateParcelsOppositeNormal_),
+    surfaceFormat_(pc.surfaceFormat_),
+    resetOnWrite_(pc.resetOnWrite_),
+    totalTime_(pc.totalTime_),
+    mass_(pc.mass_),
+    massTotal_(pc.massTotal_),
+    massFlowRate_(pc.massFlowRate_),
+    log_(pc.log_),
+    outputFilePtr_(),
+    outputDir_(pc.outputDir_),
+    timeOld_(0.0)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::ParticleCollector<CloudType>::~ParticleCollector()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::ParticleCollector<CloudType>::postMove
+(
+    const parcelType& p,
+    const label cellI,
+    const scalar dt,
+    bool& keepParticle
+)
+{
+    if ((parcelType_ != -1) && (parcelType_ != p.typeId()))
+    {
+        return;
+    }
+
+    label faceI = -1;
+
+    switch (mode_)
+    {
+        case mtPolygon:
+        {
+            faceI = collectParcelPolygon(p.position(), p.U());
+            break;
+        }
+        case mtConcentricCircle:
+        {
+            faceI = collectParcelConcentricCircles(p.position(), p.U());
+            break;
+        }
+        default:
+        {
+        }
+    }
+
+    if (faceI != -1)
+    {
+        scalar m = p.nParticle()*p.mass();
+
+        if (negateParcelsOppositeNormal_)
+        {
+            vector Uhat = p.U();
+            Uhat /= mag(Uhat) + ROOTVSMALL;
+            if ((Uhat & normal_) < 0)
+            {
+                m *= -1.0;
+            }
+        }
+
+        // add mass contribution
+        mass_[faceI] += m;
+
+        if (nSector_ == 1)
+        {
+            mass_[faceI + 1] += m;
+            mass_[faceI + 2] += m;
+            mass_[faceI + 3] += m;
+        }
+
+        if (removeCollected_)
+        {
+            keepParticle = false;
+        }
+    }
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,259 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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::ParticleCollector
+
+Description
+    Function object to collect the parcel mass- and mass flow rate over a
+    set of polygons.  The polygons are defined as lists of points.  If a
+    parcel is 'collected', it is subsequently flagged to be removed from the
+    domain.
+
+SourceFiles
+    ParticleCollector.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ParticleCollector_H
+#define ParticleCollector_H
+
+#include "CloudFunctionObject.H"
+#include "cylindricalCS.H"
+#include "face.H"
+#include "Switch.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class ParticleCollector Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class ParticleCollector
+:
+    public CloudFunctionObject<CloudType>
+{
+public:
+
+    enum modeType
+    {
+        mtPolygon,
+        mtConcentricCircle,
+        mtUnknown
+    };
+
+
+private:
+
+    // Private Data
+
+        // Typedefs
+
+            //- Convenience typedef for parcel type
+            typedef typename CloudType::parcelType parcelType;
+
+        //- Collector mode type
+        modeType mode_;
+
+        //- Index of parcel types to collect (-1 by default = all particles)
+        const label parcelType_;
+
+        //- Flag to remove collected particles
+        Switch removeCollected_;
+
+        //- List of points
+        Field<point> points_;
+
+        //- List of faces
+        List<face> faces_;
+
+
+        // Polygon collector
+
+            //- Triangulation of faces
+            List<List<face> > faceTris_;
+
+        // Concentric circles collector
+
+            //- Number of sectors per circle
+            label nSector_;
+
+            //- List of radii
+            List<scalar> radius_;
+
+            //- Cylindrical co-ordinate system
+            cylindricalCS coordSys_;
+
+
+        //- Face areas
+        Field<scalar> area_;
+
+        //- Polygon normal vector
+        vector normal_;
+
+        //- Remove mass of parcel travelling in opposite direction to normal_
+        bool negateParcelsOppositeNormal_;
+
+        //- Surface output format
+        const word surfaceFormat_;
+
+        //- Flag to indicate whether data should be reset/cleared on writing
+        Switch resetOnWrite_;
+
+        //- Total time
+        scalar totalTime_;
+
+        //- Mass storage
+        List<scalar> mass_;
+
+        //- Mass total storage
+        List<scalar> massTotal_;
+
+        //- Mass flow rate storage
+        List<scalar> massFlowRate_;
+
+        //- Flag to indicate whether data should be written to file
+        Switch log_;
+
+        //- Output file pointer
+        autoPtr<OFstream> outputFilePtr_;
+
+        //- Output directory
+        fileName outputDir_;
+
+        //- Last calculation time
+        scalar timeOld_;
+
+
+    // Private Member Functions
+
+        //- Helper function to create log files
+        void makeLogFile
+        (
+            const faceList& faces,
+            const Field<point>& points,
+            const Field<scalar>& area
+        );
+
+        //- Initialise polygon collectors
+        void initPolygons();
+
+        //- Initialise concentric circle collectors
+        void initConcentricCircles();
+
+        //- Collect parcels in polygon collectors
+        label collectParcelPolygon
+        (
+            const point& position,  
+            const vector& U
+        ) const;
+
+        //- Collect parcels in concentric circle collectors
+        label collectParcelConcentricCircles
+        (
+            const point& position,  
+            const vector& U
+        ) const;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Write post-processing info
+        void write();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("particleCollector");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        ParticleCollector(const dictionary& dict, CloudType& owner);
+
+        //- Construct copy
+        ParticleCollector(const ParticleCollector<CloudType>& pc);
+
+        //- Construct and return a clone
+        virtual autoPtr<CloudFunctionObject<CloudType> > clone() const
+        {
+            return autoPtr<CloudFunctionObject<CloudType> >
+            (
+                new ParticleCollector<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~ParticleCollector();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return const access to the reset on write flag
+            inline const Switch& resetOnWrite() const;
+
+
+        // Evaluation
+
+            //- Post-move hook
+            virtual void postMove
+            (
+                const parcelType& p,
+                const label cellI,
+                const scalar dt,
+                bool& keepParticle
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "ParticleCollectorI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "ParticleCollector.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollectorI.H

@@ -0,0 +1,34 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+inline const Foam::Switch&
+Foam::ParticleCollector<CloudType>::resetOnWrite() const
+{
+    return resetOnWrite_;
+}
+
+
+// ************************************************************************* //

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

@@ -168,7 +168,8 @@ void Foam::ParticleErosion<CloudType>::postPatch
     const parcelType& p,
     const polyPatch& pp,
     const scalar trackFraction,
-    const tetIndices& tetIs
+    const tetIndices& tetIs,
+    bool&
 )
 {
     const label patchI = pp.index();

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

@@ -128,7 +128,8 @@ public:
                 const parcelType& p,
                 const polyPatch& pp,
                 const scalar trackFraction,
-                const tetIndices& tetIs
+                const tetIndices& tetIs,
+                bool& keepParticle
             );
 };
 

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

@@ -112,7 +112,8 @@ template<class CloudType>
 void Foam::ParticleTracks<CloudType>::postFace
 (
     const parcelType& p,
-    const label
+    const label,
+    bool&
 )
 {
     if

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

@@ -143,7 +143,8 @@ public:
             virtual void postFace
             (
                 const parcelType& p,
-                const label faceI
+                const label faceI,
+                bool& keepParticle
             );
 };
 

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

@@ -104,7 +104,8 @@ void Foam::ParticleTrap<CloudType>::postMove
 (
     parcelType& p,
     const label cellI,
-    const scalar
+    const scalar,
+    bool&
 )
 {
     if (alphaPtr_->internalField()[cellI] < threshold_)

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

@@ -125,7 +125,8 @@ public:
             (
                 typename CloudType::parcelType& p,
                 const label cellI,
-                const scalar dt
+                const scalar dt,
+                bool& keepParticle
             );
 };
 

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

@@ -220,7 +220,8 @@ void Foam::PatchPostProcessing<CloudType>::postPatch
     const parcelType& p,
     const polyPatch& pp,
     const scalar,
-    const tetIndices& tetIs
+    const tetIndices& tetIs,
+    bool&
 )
 {
     const label patchI = pp.index();

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

@@ -129,7 +129,8 @@ public:
                 const parcelType& p,
                 const polyPatch& pp,
                 const scalar trackFraction,
-                const tetIndices& tetIs
+                const tetIndices& tetIs,
+                bool& keepParticle
             );
 };
 

src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -125,7 +125,8 @@ void Foam::VoidFraction<CloudType>::postMove
 (
     const parcelType& p,
     const label cellI,
-    const scalar dt
+    const scalar dt,
+    bool&
 )
 {
     volScalarField& theta = thetaPtr_();

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

@@ -115,7 +115,8 @@ public:
             (
                 const parcelType& p,
                 const label cellI,
-                const scalar dt
+                const scalar dt,
+                bool& keepParticle
             );
 };
 

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C

@@ -302,13 +302,15 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
                     if (particleHit)
                     {
-                        this->owner().functions().postFace(p, realFaceI);
+                        bool keep = true;
+                        this->owner().functions().postFace(p, realFaceI, keep);
                         this->owner().functions().postPatch
                         (
                             p,
                             mesh.boundaryMesh()[patchI],
                             1.0,
-                            p.currentTetIndices()
+                            p.currentTetIndices(),
+                            keep
                         );
                      }
                 }

src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C

@@ -398,13 +398,15 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
     // surface energy of secondary parcels [J]
     scalar ESigmaSec = 0;
 
-    // sample splash distribution to detrmine secondary parcel diameters
+    // sample splash distribution to determine secondary parcel diameters
     scalarList dNew(parcelsPerSplash_);
+    scalarList npNew(parcelsPerSplash_);
     forAll(dNew, i)
     {
         const scalar y = rndGen_.sample01<scalar>();
         dNew[i] = -dBarSplash*log(exp(-dMin/dBarSplash) - y*K);
-        ESigmaSec += sigma*p.areaS(dNew[i]);
+        npNew[i] = mRatio*np*pow3(d)/pow3(dNew[i])/parcelsPerSplash_;
+        ESigmaSec += npNew[i]*sigma*p.areaS(dNew[i]);
     }
 
     // incident kinetic energy [J]
@@ -459,7 +461,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
         // perturb new parcels towards the owner cell centre
         pPtr->position() += 0.5*rndGen_.sample01<scalar>()*(posC - posCf);
 
-        pPtr->nParticle() = mRatio*np*pow3(d)/pow3(dNew[i])/parcelsPerSplash_;
+        pPtr->nParticle() = npNew[i];
 
         pPtr->d() = dNew[i];
 

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C

@@ -90,58 +90,6 @@ bool reactingOneDim::read(const dictionary& dict)
 }
 
 
-void reactingOneDim::updateQr()
-{
-    // Retrieve field from coupled region using mapped boundary conditions
-    QrCoupled_.correctBoundaryConditions();
-
-    // Update local Qr from coupled Qr field
-    Qr_ == dimensionedScalar("zero", Qr_.dimensions(), 0.0);
-    forAll(intCoupledPatchIDs_, i)
-    {
-        const label patchI = intCoupledPatchIDs_[i];
-
-        scalarField& Qrp = Qr_.boundaryField()[patchI];
-
-        // Qr is negative going out the solid
-        // If the surface is emitting the radiative flux is set to zero
-        Qrp = max(Qrp, scalar(0.0));
-    }
-
-    const volScalarField kappaRad_(kappaRad());
-
-    // Propagate Qr through 1-D regions
-    label totalFaceId = 0;
-    forAll(intCoupledPatchIDs_, i)
-    {
-        const label patchI = intCoupledPatchIDs_[i];
-
-        const scalarField& Qrp = Qr_.boundaryField()[patchI];
-        const vectorField& Cf = regionMesh().Cf().boundaryField()[patchI];
-
-        forAll(Qrp, faceI)
-        {
-            const scalar Qr0 = Qrp[faceI];
-            point Cf0 = Cf[faceI];
-            const labelList& cells = boundaryFaceCells_[totalFaceId];
-            scalar kappaInt = 0.0;
-            forAll(cells, k)
-            {
-                const label cellI = cells[k];
-                const point& Cf1 = regionMesh().cellCentres()[cellI];
-                const scalar delta = mag(Cf1 - Cf0);
-                kappaInt += kappaRad_[cellI]*delta;
-                Qr_[cellI] = Qr0*exp(-kappaInt);
-                Cf0 = Cf1;
-            }
-            totalFaceId ++;
-        }
-    }
-
-    Qr_.correctBoundaryConditions();
-}
-
-
 void reactingOneDim::updatePhiGas()
 {
     phiHsGas_ ==  dimensionedScalar("zero", phiHsGas_.dimensions(), 0.0);
@@ -198,8 +146,6 @@ void reactingOneDim::updatePhiGas()
 
 void reactingOneDim::updateFields()
 {
-    updateQr();
-
     updatePhiGas();
 }
 
@@ -305,8 +251,6 @@ void reactingOneDim::solveEnergy()
 
     tmp<volScalarField> alpha(solidThermo_.alpha());
 
-    const surfaceScalarField phiQr(fvc::interpolate(Qr_)*nMagSf());
-
     const surfaceScalarField phiGas(fvc::interpolate(phiHsGas_));
 
     fvScalarMatrix hEqn
@@ -315,7 +259,6 @@ void reactingOneDim::solveEnergy()
       - fvm::laplacian(alpha, h_)
      ==
         chemistrySh_
-      + fvc::div(phiQr)
       + fvc::div(phiGas)
     );
 
@@ -380,7 +323,6 @@ reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
     ),
     Ys_(solidThermo_.composition().Y()),
     h_(solidThermo_.he()),
-    primaryRadFluxName_(coeffs().lookupOrDefault<word>("radFluxName", "Qr")),
     nNonOrthCorr_(-1),
     maxDiff_(10),
     minimumDelta_(1e-4),
@@ -427,34 +369,6 @@ reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
         dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
     ),
 
-    QrCoupled_
-    (
-        IOobject
-        (
-            primaryRadFluxName_,
-            time().timeName(),
-            regionMesh(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        regionMesh()
-    ),
-
-    Qr_
-    (
-        IOobject
-        (
-            "QrPyr",
-            time().timeName(),
-            regionMesh(),
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        regionMesh(),
-        dimensionedScalar("zero", dimEnergy/dimArea/dimTime, 0.0),
-        zeroGradientFvPatchVectorField::typeName
-    ),
-
     lostSolidMass_(dimensionedScalar("zero", dimMass, 0.0)),
     addedGasMass_(dimensionedScalar("zero", dimMass, 0.0)),
     totalGasMassFlux_(0.0),
@@ -492,7 +406,6 @@ reactingOneDim::reactingOneDim
     ),
     Ys_(solidThermo_.composition().Y()),
     h_(solidThermo_.he()),
-    primaryRadFluxName_(dict.lookupOrDefault<word>("radFluxName", "Qr")),
     nNonOrthCorr_(-1),
     maxDiff_(10),
     minimumDelta_(1e-4),
@@ -539,34 +452,6 @@ reactingOneDim::reactingOneDim
         dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
     ),
 
-    QrCoupled_
-    (
-        IOobject
-        (
-            primaryRadFluxName_,
-            time().timeName(),
-            regionMesh(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        regionMesh()
-    ),
-
-    Qr_
-    (
-        IOobject
-        (
-            "QrPyr",
-            time().timeName(),
-            regionMesh(),
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        regionMesh(),
-        dimensionedScalar("zero", dimEnergy/dimArea/dimTime, 0.0),
-        zeroGradientFvPatchVectorField::typeName
-    ),
-
     lostSolidMass_(dimensionedScalar("zero", dimMass, 0.0)),
     addedGasMass_(dimensionedScalar("zero", dimMass, 0.0)),
     totalGasMassFlux_(0.0),
@@ -684,36 +569,6 @@ void reactingOneDim::preEvolveRegion()
     {
         solidChemistry_->setCellReacting(cellI, true);
     }
-
-    // De-activate reactions if pyrolysis region coupled to (valid) film
-    if (filmCoupled_)
-    {
-        const volScalarField& filmDelta = filmDeltaPtr_();
-
-        forAll(intCoupledPatchIDs_, i)
-        {
-            const label patchI = intCoupledPatchIDs_[i];
-            const scalarField& filmDeltap = filmDelta.boundaryField()[patchI];
-
-            forAll(filmDeltap, faceI)
-            {
-                const scalar filmDelta0 = filmDeltap[faceI];
-                if (filmDelta0 > reactionDeltaMin_)
-                {
-                    const labelList& cells = boundaryFaceCells_[faceI];
-
-                    // TODO: only limit cell adjacent to film?
-                    //solidChemistry_->setCellNoReacting(cells[0])
-
-                    // Propagate flag through 1-D region
-                    forAll(cells, k)
-                    {
-                        solidChemistry_->setCellReacting(cells[k], false);
-                    }
-                }
-            }
-        }
-    }
 }
 
 

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -125,16 +125,6 @@ protected:
             volScalarField chemistrySh_;
 
 
-        // Source term fields
-
-            //- Coupled region radiative heat flux [W/m2]
-            //  Requires user to input mapping info for coupled patches
-            volScalarField QrCoupled_;
-
-            //- In depth radiative heat flux [W/m2]
-            volScalarField Qr_;
-
-
         // Checks
 
             //- Cumulative lost mass of the condensed phase [kg]
@@ -164,9 +154,6 @@ protected:
         //- Update/move mesh based on change in mass
         void updateMesh(const scalarField& mass0);
 
-        //- Update radiative flux in pyrolysis region
-        void updateQr();
-
         //- Update enthalpy flux for pyrolysis gases
         void updatePhiGas();
 
@@ -259,12 +246,6 @@ public:
             virtual scalar solidRegionDiffNo() const;
 
 
-        // Source fields (read/write access)
-
-            //- In depth radiative heat flux
-            inline const volScalarField& Qr() const;
-
-
        // Evolution
 
             //- Pre-evolve region

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDimI.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,11 +34,4 @@ Foam::regionModels::pyrolysisModels::reactingOneDim::nNonOrthCorr() const
 }
 
 
-inline const Foam::volScalarField&
-Foam::regionModels::pyrolysisModels::reactingOneDim::Qr() const
-{
-    return Qr_;
-}
-
-
 // ************************************************************************* //

src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisTemperatureCoupled/filmPyrolysisTemperatureCoupledFvPatchScalarField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,8 +40,7 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 :
     fixedValueFvPatchScalarField(p, iF),
     phiName_("phi"),
-    rhoName_("rho"),
+    rhoName_("rho")
-    deltaWet_(1e-6)
 {}
 
 
@@ -56,8 +55,7 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 :
     fixedValueFvPatchScalarField(ptf, p, iF, mapper),
     phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_),
+    rhoName_(ptf.rhoName_)
-    deltaWet_(ptf.deltaWet_)
 {}
 
 
@@ -71,8 +69,7 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 :
     fixedValueFvPatchScalarField(p, iF),
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
-    deltaWet_(dict.lookupOrDefault<scalar>("deltaWet", 1e-6))
 {
     fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
 }
@@ -86,8 +83,7 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 :
     fixedValueFvPatchScalarField(fptpsf),
     phiName_(fptpsf.phiName_),
-    rhoName_(fptpsf.rhoName_),
+    rhoName_(fptpsf.rhoName_)
-    deltaWet_(fptpsf.deltaWet_)
 {}
 
 
@@ -100,8 +96,7 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 :
     fixedValueFvPatchScalarField(fptpsf, iF),
     phiName_(fptpsf.phiName_),
-    rhoName_(fptpsf.rhoName_),
+    rhoName_(fptpsf.rhoName_)
-    deltaWet_(fptpsf.deltaWet_)
 {}
 
 
@@ -151,13 +146,12 @@ void Foam::filmPyrolysisTemperatureCoupledFvPatchScalarField::updateCoeffs()
 
     const label filmPatchI = filmModel.regionPatchID(patchI);
 
-    scalarField deltaFilm = filmModel.delta().boundaryField()[filmPatchI];
+    scalarField alphaFilm = filmModel.alpha().boundaryField()[filmPatchI];
-    filmModel.toPrimary(filmPatchI, deltaFilm);
+    filmModel.toPrimary(filmPatchI, alphaFilm);
 
     scalarField TFilm = filmModel.Ts().boundaryField()[filmPatchI];
     filmModel.toPrimary(filmPatchI, TFilm);
 
-
     // Retrieve pyrolysis model
     const pyrModelType& pyrModel =
         db().lookupObject<pyrModelType>("pyrolysisProperties");
@@ -168,19 +162,8 @@ void Foam::filmPyrolysisTemperatureCoupledFvPatchScalarField::updateCoeffs()
     pyrModel.toPrimary(pyrPatchI, TPyr);
 
 
-    forAll(deltaFilm, i)
+    // Evaluate temperature
-    {
+    Tp = alphaFilm*TFilm + (1.0 - alphaFilm)*TPyr;
-        if (deltaFilm[i] > deltaWet_)
-        {
-            // temperature set by film
-            Tp[i] = TFilm[i];
-        }
-        else
-        {
-            // temperature set by pyrolysis model
-            Tp[i] = TPyr[i];
-        }
-    }
 
     // Restore tag
     UPstream::msgType() = oldTag;
@@ -197,7 +180,6 @@ void Foam::filmPyrolysisTemperatureCoupledFvPatchScalarField::write
     fvPatchScalarField::write(os);
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
-    os.writeKeyword("deltaWet") << deltaWet_ << token::END_STATEMENT << nl;
     writeEntry("value", os);
 }
 

src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisTemperatureCoupled/filmPyrolysisTemperatureCoupledFvPatchScalarField.H

@@ -28,15 +28,7 @@ Description
     This boundary condition is designed to be used in conjunction with surface
     film and pyrolysis modelling.  It provides a temperature boundary condition
     for patches on the primary region based on whether the patch is seen to
-    be 'wet', specified by:
+    be 'wet', retrieved from the film alpha field.
-
-        \f[
-            delta > delta_wet
-        \f]
-
-    where
-        \var delta     = film height [m]
-        \var delta_wet = film height above which the surface is considered wet
 
     \li if the patch is wet, the temperature is set using the film temperature
     \li otherwise, it is set using pyrolysis temperature
@@ -84,9 +76,6 @@ class filmPyrolysisTemperatureCoupledFvPatchScalarField
         //- Name of density field
         word rhoName_;
 
-        //- Film height threshold beyond which it is considered 'wet' [m]
-        scalar deltaWet_;
-
 
 public:
 

src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,8 +40,7 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 :
     fixedValueFvPatchVectorField(p, iF),
     phiName_("phi"),
-    rhoName_("rho"),
+    rhoName_("rho")
-    deltaWet_(1e-6)
 {}
 
 
@@ -56,8 +55,7 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 :
     fixedValueFvPatchVectorField(ptf, p, iF, mapper),
     phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_),
+    rhoName_(ptf.rhoName_)
-    deltaWet_(ptf.deltaWet_)
 {}
 
 
@@ -71,8 +69,7 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 :
     fixedValueFvPatchVectorField(p, iF),
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
-    deltaWet_(dict.lookupOrDefault<scalar>("deltaWet", 1e-6))
 {
     fvPatchVectorField::operator=(vectorField("value", dict, p.size()));
 }
@@ -86,8 +83,7 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 :
     fixedValueFvPatchVectorField(fpvpvf),
     phiName_(fpvpvf.phiName_),
-    rhoName_(fpvpvf.rhoName_),
+    rhoName_(fpvpvf.rhoName_)
-    deltaWet_(fpvpvf.deltaWet_)
 {}
 
 
@@ -100,8 +96,7 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 :
     fixedValueFvPatchVectorField(fpvpvf, iF),
     phiName_(fpvpvf.phiName_),
-    rhoName_(fpvpvf.rhoName_),
+    rhoName_(fpvpvf.rhoName_)
-    deltaWet_(fpvpvf.deltaWet_)
 {}
 
 
@@ -154,13 +149,12 @@ void Foam::filmPyrolysisVelocityCoupledFvPatchVectorField::updateCoeffs()
 
     const label filmPatchI = filmModel.regionPatchID(patchI);
 
-    scalarField deltaFilm = filmModel.delta().boundaryField()[filmPatchI];
+    scalarField alphaFilm = filmModel.alpha().boundaryField()[filmPatchI];
-    filmModel.toPrimary(filmPatchI, deltaFilm);
+    filmModel.toPrimary(filmPatchI, alphaFilm);
 
     vectorField UFilm = filmModel.Us().boundaryField()[filmPatchI];
     filmModel.toPrimary(filmPatchI, UFilm);
 
-
     // Retrieve pyrolysis model
     const pyrModelType& pyrModel =
         db().objectRegistry::lookupObject<pyrModelType>
@@ -203,19 +197,9 @@ void Foam::filmPyrolysisVelocityCoupledFvPatchVectorField::updateCoeffs()
     const scalarField UAvePyr(-phiPyr/patch().magSf());
     const vectorField& nf = patch().nf();
 
-    forAll(deltaFilm, i)
+
-    {
+    // Evaluate velocity
-        if (deltaFilm[i] > deltaWet_)
+    Up = alphaFilm*UFilm + (1.0 - alphaFilm)*UAvePyr*nf;
-        {
-            // velocity set by film
-            Up[i] = UFilm[i];
-        }
-        else
-        {
-            // velocity set by pyrolysis model
-            Up[i] = UAvePyr[i]*nf[i];
-        }
-    }
 
     // Restore tag
     UPstream::msgType() = oldTag;
@@ -232,7 +216,6 @@ void Foam::filmPyrolysisVelocityCoupledFvPatchVectorField::write
     fvPatchVectorField::write(os);
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
-    os.writeKeyword("deltaWet") << deltaWet_ << token::END_STATEMENT << nl;
     writeEntry("value", os);
 }
 

src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.H

@@ -28,15 +28,7 @@ Description
     This boundary condition is designed to be used in conjunction with surface
     film and pyrolysis modelling.  It provides a velocity boundary condition
     for patches on the primary region based on whether the patch is seen to
-    be 'wet', specified by:
+    be 'wet', retrieved from the film alpha field.
-
-        \f[
-            delta > delta_wet
-        \f]
-
-    where
-        \var delta     = film height [m]
-        \var delta_wet = film height above which the surface is considered wet
 
     \li if the patch is wet, the velocity is set using the film velocity
     \li otherwise, it is set using pyrolysis out-gassing velocity
@@ -84,9 +76,6 @@ class filmPyrolysisVelocityCoupledFvPatchVectorField
         //- Name of density field
         word rhoName_;
 
-        //- Film height threshold beyond which it is considered 'wet'
-        scalar deltaWet_;
-
 
 public:
 

src/regionModels/regionModel/Make/files

@@ -8,4 +8,8 @@ derivedFvPatches/mappedVariableThicknessWall/mappedVariableThicknessWallFvPatch.
 
 regionProperties/regionProperties.C
 
+regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObject.C
+regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectNew.C
+regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectList.C
+
 LIB = $(FOAM_LIBBIN)/libregionModels

src/regionModels/regionModel/regionModel/regionModel.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -190,6 +190,66 @@ bool Foam::regionModels::regionModel::read(const dictionary& dict)
 }
 
 
+Foam::label Foam::regionModels::regionModel::nbrCoupledPatchID
+(
+    const regionModel& nbrRegion,
+    const label regionPatchI
+) const
+{
+    label nbrPatchI = -1;
+
+    // region
+    const fvMesh& nbrRegionMesh = nbrRegion.regionMesh();
+
+    // boundary mesh
+    const polyBoundaryMesh& nbrPbm = nbrRegionMesh.boundaryMesh();
+
+    const mappedPatchBase& mpb =
+        refCast<const mappedPatchBase>
+        (
+            regionMesh().boundaryMesh()[regionPatchI]
+        );
+
+    // sample patch name on the primary region
+    const word& primaryPatchName = mpb.samplePatch();
+
+    // find patch on nbr region that has the same sample patch name
+    forAll(nbrRegion.intCoupledPatchIDs(), j)
+    {
+        const label nbrRegionPatchI = nbrRegion.intCoupledPatchIDs()[j];
+
+        const mappedPatchBase& mpb =
+            refCast<const mappedPatchBase>(nbrPbm[nbrRegionPatchI]);
+
+        if (mpb.samplePatch() == primaryPatchName)
+        {
+            nbrPatchI = nbrRegionPatchI;
+            break;
+        }
+    }
+
+    if (nbrPatchI == -1)
+    {
+        const polyPatch& p = regionMesh().boundaryMesh()[regionPatchI];
+
+        FatalErrorIn
+        (
+            "Foam::tmp<Foam::Field<Type> > "
+            "Foam::regionModels::regionModel::nbrCoupledPatchID"
+            "("
+                "const regionModel& , "
+                "const label"
+            ") const"
+        )
+            << "Unable to find patch pair for local patch "
+            << p.name() << " and region " << nbrRegion.name()
+            << abort(FatalError);
+    }
+
+    return nbrPatchI;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::regionModels::regionModel::regionModel(const fvMesh& mesh)
@@ -214,7 +274,8 @@ Foam::regionModels::regionModel::regionModel(const fvMesh& mesh)
     coeffs_(dictionary::null),
     primaryPatchIDs_(),
     intCoupledPatchIDs_(),
-    regionName_("none")
+    regionName_("none"),
+    functions_(*this)
 {}
 
 
@@ -246,7 +307,8 @@ Foam::regionModels::regionModel::regionModel
     coeffs_(subOrEmptyDict(modelName + "Coeffs")),
     primaryPatchIDs_(),
     intCoupledPatchIDs_(),
-    regionName_(lookup("regionName"))
+    regionName_(lookup("regionName")),
+    functions_(*this, subOrEmptyDict("functions"))
 {
     if (active_)
     {
@@ -274,7 +336,7 @@ Foam::regionModels::regionModel::regionModel
     (
         IOobject
         (
-            regionType,
+            regionType + "Properties",
             mesh.time().constant(),
             mesh,
             IOobject::NO_READ,
@@ -292,7 +354,8 @@ Foam::regionModels::regionModel::regionModel
     coeffs_(dict.subOrEmptyDict(modelName + "Coeffs")),
     primaryPatchIDs_(),
     intCoupledPatchIDs_(),
-    regionName_(dict.lookup("regionName"))
+    regionName_(dict.lookup("regionName")),
+    functions_(*this, subOrEmptyDict("functions"))
 {
     if (active_)
     {
@@ -315,18 +378,6 @@ Foam::regionModels::regionModel::~regionModel()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::regionModels::regionModel::preEvolveRegion()
-{
-    // do nothing
-}
-
-
-void Foam::regionModels::regionModel::evolveRegion()
-{
-    // do nothing
-}
-
-
 void Foam::regionModels::regionModel::evolve()
 {
     if (active_)
@@ -334,15 +385,14 @@ void Foam::regionModels::regionModel::evolve()
         Info<< "\nEvolving " << modelName_ << " for region "
             << regionMesh().name() << endl;
 
-        // Update any input information
         //read();
 
-        // Pre-evolve
         preEvolveRegion();
 
-        // Increment the region equations up to the new time level
         evolveRegion();
 
+        postEvolveRegion();
+
         // Provide some feedback
         if (infoOutput_)
         {
@@ -354,6 +404,24 @@ void Foam::regionModels::regionModel::evolve()
 }
 
 
+void Foam::regionModels::regionModel::preEvolveRegion()
+{
+    functions_.preEvolveRegion();
+}
+
+
+void Foam::regionModels::regionModel::evolveRegion()
+{
+    // do nothing
+}
+
+
+void Foam::regionModels::regionModel::postEvolveRegion()
+{
+    functions_.postEvolveRegion();
+}
+
+
 void Foam::regionModels::regionModel::info() const
 {
     // do nothing

src/regionModels/regionModel/regionModel/regionModel.H

@@ -41,16 +41,13 @@ SourceFiles
 #include "labelList.H"
 #include "volFields.H"
 #include "mappedPatchBase.H"
+#include "regionModelFunctionObjectList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// Forward declaration of classes
-//class fvMesh;
-//class Time;
-
 namespace regionModels
 {
 
@@ -121,6 +118,9 @@ protected:
         //- Region name
         word regionName_;
 
+        //- Region model function objects
+        regionModelFunctionObjectList functions_;
+
 
     // Protected member functions
 
@@ -219,6 +219,46 @@ public:
 
         // Helper
 
+            //- Return the coupled patch ID paired with coupled patch
+            //  regionPatchI
+            label nbrCoupledPatchID
+            (
+                const regionModel& nbrRegion,
+                const label regionPatchI
+            ) const;
+
+            //- Map patch field from another region model to local patch
+            template<class Type>
+            tmp<Foam::Field<Type> > mapRegionPatchField
+            (
+                const regionModel& nbrRegion,
+                const label regionPatchI,
+                const label nbrPatchI,
+                const Field<Type>& nbrField,
+                const bool flip = false
+            ) const;
+
+            //- Map patch field from another region model to local patch
+            template<class Type>
+            tmp<Field<Type> > mapRegionPatchField
+            (
+                const word& regionModelName,
+                const word& fieldName,
+                const label regionPatchI,
+                const bool flip = false
+            ) const;
+
+            //- Map patch internal field from another region model to local
+            //  patch
+            template<class Type>
+            tmp<Field<Type> > mapRegionPatchInternalField
+            (
+                const word& regionModelName,
+                const word& fieldName,
+                const label regionPatchI,
+                const bool flip = false
+            ) const;
+
             //- Convert a local region field to the primary region
             template<class Type>
             void toPrimary
@@ -256,14 +296,17 @@ public:
 
         // Evolution
 
+            //- Main driver routing to evolve the region - calls other evolves
+            virtual void evolve();
+
             //- Pre-evolve region
             virtual void preEvolveRegion();
 
             //- Evolve the region
             virtual void evolveRegion();
 
-            //- Evolve the region
+            //- Post-evolve region
-            virtual void evolve();
+            virtual void postEvolveRegion();
 
 
         // I-O

src/regionModels/regionModel/regionModel/regionModelTemplates.C

@@ -23,6 +23,164 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::regionModels::regionModel::mapRegionPatchField
+(
+    const regionModel& nbrRegion,
+    const label regionPatchI,
+    const label nbrPatchI,
+    const Field<Type>& nbrField,
+    const bool flip
+) const
+{
+    const fvMesh& nbrRegionMesh = nbrRegion.regionMesh();
+
+    const polyPatch& p = regionMesh().boundaryMesh()[regionPatchI];
+
+    const polyPatch& nbrP = nbrRegionMesh.boundaryMesh()[nbrPatchI];
+
+    int oldTag = UPstream::msgType();
+    UPstream::msgType() = oldTag + 1;
+
+    AMIPatchToPatchInterpolation ami(p, nbrP, faceAreaIntersect::tmMesh, flip);
+
+    tmp<Field<Type> > tresult(ami.interpolateToSource(nbrField));
+
+    UPstream::msgType() = oldTag;
+
+    return tresult;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::regionModels::regionModel::mapRegionPatchField
+(
+    const word& regionModelName,
+    const word& fieldName,
+    const label regionPatchI,
+    const bool flip
+) const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    const regionModel& nbrRegion =
+        this->primaryMesh_.lookupObject<regionModel>(regionModelName);
+
+    const fvMesh& nbrRegionMesh = nbrRegion.regionMesh();
+
+    const polyPatch& p = regionMesh().boundaryMesh()[regionPatchI];
+
+    if (nbrRegionMesh.foundObject<fieldType>(fieldName))
+    {
+        const fieldType& nbrField =
+            nbrRegionMesh.lookupObject<fieldType>(fieldName);
+
+        const label nbrPatchI = nbrCoupledPatchID(nbrRegion, regionPatchI);
+
+        const polyPatch& nbrP = nbrRegionMesh.boundaryMesh()[nbrPatchI];
+
+        int oldTag = UPstream::msgType();
+        UPstream::msgType() = oldTag + 1;
+
+        AMIPatchToPatchInterpolation ami
+        (
+            p,
+            nbrP,
+            faceAreaIntersect::tmMesh,
+            flip
+        );
+
+        const Field<Type>& nbrFieldp = nbrField.boundaryField()[nbrPatchI];
+
+        tmp<Field<Type> > tresult(ami.interpolateToSource(nbrFieldp));
+
+        UPstream::msgType() = oldTag;
+
+        return tresult;
+    }
+    else
+    {
+        return
+            tmp<Field<Type> >
+            (
+                new Field<Type>
+                (
+                    p.size(),
+                    pTraits<Type>::zero
+                )
+            );
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::regionModels::regionModel::mapRegionPatchInternalField
+(
+    const word& regionModelName,
+    const word& fieldName,
+    const label regionPatchI,
+    const bool flip
+) const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    const regionModel& nbrRegion =
+        this->primaryMesh_.lookupObject<regionModel>(regionModelName);
+
+    const fvMesh& nbrRegionMesh = nbrRegion.regionMesh();
+
+    const polyPatch& p = regionMesh().boundaryMesh()[regionPatchI];
+
+    if (nbrRegionMesh.foundObject<fieldType>(fieldName))
+    {
+        const fieldType& nbrField =
+            nbrRegionMesh.lookupObject<fieldType>(fieldName);
+
+        const label nbrPatchI = nbrCoupledPatchID(nbrRegion, regionPatchI);
+
+        const polyPatch& nbrP = nbrRegionMesh.boundaryMesh()[nbrPatchI];
+
+        int oldTag = UPstream::msgType();
+        UPstream::msgType() = oldTag + 1;
+
+        AMIPatchToPatchInterpolation ami
+        (
+            p,
+            nbrP,
+            faceAreaIntersect::tmMesh,
+            flip
+        );
+
+        const fvPatchField<Type>& nbrFieldp =
+            nbrField.boundaryField()[nbrPatchI];
+
+        tmp<Field<Type> > tresult
+        (
+            ami.interpolateToSource(nbrFieldp.patchInternalField())
+        );
+
+        UPstream::msgType() = oldTag;
+
+        return tresult;
+    }
+    else
+    {
+        return
+            tmp<Field<Type> >
+            (
+                new Field<Type>
+                (
+                    p.size(),
+                    pTraits<Type>::zero
+                )
+            );
+    }
+}
+
+
 template<class Type>
 void Foam::regionModels::regionModel::toPrimary
 (

src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObject.C

@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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 "regionModelFunctionObject.H"
+#include "regionModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+    defineTypeNameAndDebug(regionModelFunctionObject, 0);
+    defineRunTimeSelectionTable(regionModelFunctionObject, dictionary);
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::regionModels::regionModelFunctionObject::regionModelFunctionObject
+(
+    regionModel& owner
+)
+:
+    dict_(dictionary::null),
+    owner_(owner),
+    modelType_("modelType")
+{}
+
+
+Foam::regionModels::regionModelFunctionObject::regionModelFunctionObject
+(
+    const dictionary& dict,
+    regionModel& owner,
+    const word& type
+)
+:
+    dict_(dict),
+    owner_(owner),
+    modelType_(type)
+{}
+
+
+Foam::regionModels::regionModelFunctionObject::regionModelFunctionObject
+(
+    const regionModelFunctionObject& rmfo
+)
+:
+    dict_(rmfo.dict_),
+    owner_(rmfo.owner_),
+    modelType_(rmfo.modelType_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::regionModels::regionModelFunctionObject::~regionModelFunctionObject()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::regionModels::regionModelFunctionObject::preEvolveRegion()
+{
+    // do nothing
+}
+
+
+void Foam::regionModels::regionModelFunctionObject::postEvolveRegion()
+{
+    if (owner_.regionMesh().time().outputTime())
+    {
+        write();
+    }
+}
+
+
+void Foam::regionModels::regionModelFunctionObject::write() const
+{
+    // do nothing
+}
+
+
+// ************************************************************************* //

src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObject.H

@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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::regionModelFunctionObject
+
+Description
+    Region model function object base class
+
+SourceFiles
+    regionModelFunctionObject.C
+    regionModelFunctionObjectNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef regionModelFunctionObject_H
+#define regionModelFunctionObject_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace regionModels
+{
+
+class regionModel;
+
+/*---------------------------------------------------------------------------*\
+                  Class regionModelFunctionObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class regionModelFunctionObject
+{
+protected:
+
+    // Protected data
+
+        //- Dictionary
+        dictionary dict_;
+
+        //- Reference to the region model
+        regionModel& owner_;
+
+        //- Model type name
+        word modelType_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("regionModelFunctionObject");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        regionModelFunctionObject,
+        dictionary,
+        (
+            const dictionary& dict,
+            regionModel& owner
+        ),
+        (dict, owner)
+    );
+
+
+    // Constructors
+
+        //- Construct null from owner
+        regionModelFunctionObject(regionModel& owner);
+
+        //- Construct from dictionary
+        regionModelFunctionObject
+        (
+            const dictionary& dict,
+            regionModel& owner,
+            const word& modelType
+        );
+
+        //- Construct copy
+        regionModelFunctionObject(const regionModelFunctionObject& ppm);
+
+        //- Construct and return a clone
+        virtual autoPtr<regionModelFunctionObject> clone() const
+        {
+            return autoPtr<regionModelFunctionObject>
+            (
+                new regionModelFunctionObject(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~regionModelFunctionObject();
+
+
+    //- Selector
+    static autoPtr<regionModelFunctionObject> New
+    (
+        const dictionary& dict,
+        regionModel& owner,
+        const word& modelType
+    );
+
+
+    // Member Functions
+
+        // Evaluation
+
+            //- Pre-evolve region hook
+            virtual void preEvolveRegion();
+
+            //- Post-evolve region hook
+            virtual void postEvolveRegion();
+
+    // I-O
+
+        //- write
+        virtual void write() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectList.C

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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 "regionModelFunctionObjectList.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::regionModels::regionModelFunctionObjectList::regionModelFunctionObjectList
+(
+    regionModel& owner
+)
+:
+    PtrList<regionModelFunctionObject>(),
+    owner_(owner),
+    dict_(dictionary::null)
+{}
+
+
+Foam::regionModels::regionModelFunctionObjectList::regionModelFunctionObjectList
+(
+    regionModel& owner,
+    const dictionary& dict,
+    const bool readFields
+)
+:
+    PtrList<regionModelFunctionObject>(),
+    owner_(owner),
+    dict_(dict)
+{
+    if (readFields)
+    {
+        wordList modelNames(dict.toc());
+
+        Info<< "    Selecting region model functions" << endl;
+
+        if (modelNames.size() > 0)
+        {
+            this->setSize(modelNames.size());
+
+            forAll(modelNames, i)
+            {
+                const word& modelName = modelNames[i];
+
+                this->set
+                (
+                    i,
+                    regionModelFunctionObject::New
+                    (
+                        dict,
+                        owner,
+                        modelName
+                    )
+                );
+            }
+        }
+        else
+        {
+            Info<< "    none" << endl;
+        }
+    }
+}
+
+
+Foam::regionModels::regionModelFunctionObjectList::regionModelFunctionObjectList
+(
+    const regionModelFunctionObjectList& cfol
+)
+:
+    PtrList<regionModelFunctionObject>(cfol),
+    owner_(cfol.owner_),
+    dict_(cfol.dict_)
+{}
+
+
+// * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
+
+Foam::regionModels::regionModelFunctionObjectList::
+~regionModelFunctionObjectList()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::regionModels::regionModelFunctionObjectList::preEvolveRegion()
+{
+    forAll(*this, i)
+    {
+        this->operator[](i).preEvolveRegion();
+    }
+}
+
+
+void Foam::regionModels::regionModelFunctionObjectList::postEvolveRegion()
+{
+    forAll(*this, i)
+    {
+        this->operator[](i).postEvolveRegion();
+    }
+}
+
+
+// ************************************************************************* //

src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectList.H

@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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::regionModelFunctionObjectList
+
+Description
+    List of cloud function objects
+
+SourceFiles
+    regionModelFunctionObjectListI.H
+    regionModelFunctionObjectList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef regionModelFunctionObjectList_H
+#define regionModelFunctionObjectList_H
+
+#include "PtrList.H"
+#include "regionModelFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+
+class regionModel;
+
+/*---------------------------------------------------------------------------*\
+                Class regionModelFunctionObjectList Declaration
+\*---------------------------------------------------------------------------*/
+
+class regionModelFunctionObjectList
+:
+    public PtrList<regionModelFunctionObject>
+{
+protected:
+
+    // Protected Data
+
+        //- Reference to the owner region model
+        regionModel& owner_;
+
+        //- Dictionary
+        const dictionary dict_;
+
+
+public:
+
+    // Constructors
+
+        //- Null constructor
+        regionModelFunctionObjectList(regionModel& owner);
+
+        //- Construct from mesh
+        regionModelFunctionObjectList
+        (
+            regionModel& owner,
+            const dictionary& dict,
+            const bool readFields = true
+        );
+
+        //- Construct copy
+        regionModelFunctionObjectList
+        (
+            const regionModelFunctionObjectList& rmfol
+        );
+
+
+    //- Destructor
+    virtual ~regionModelFunctionObjectList();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return const access to the cloud owner
+            inline const regionModel& owner() const;
+
+            //- Return refernce to the cloud owner
+            inline regionModel& owner();
+
+            //- Return the forces dictionary
+            inline const dictionary& dict() const;
+
+
+        // Evaluation
+
+            //- Pre-evolve hook
+            virtual void preEvolveRegion();
+
+            //- Post-evolve hook
+            virtual void postEvolveRegion();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "regionModelFunctionObjectListI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectListI.H

@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+inline const Foam::regionModels::regionModel&
+Foam::regionModels::regionModelFunctionObjectList::owner() const
+{
+    return owner_;
+}
+
+
+inline Foam::regionModels::regionModel&
+Foam::regionModels::regionModelFunctionObjectList::owner()
+{
+    return owner_;
+}
+
+
+inline const Foam::dictionary&
+Foam::regionModels::regionModelFunctionObjectList::dict() const
+{
+    return dict_;
+}
+
+
+// ************************************************************************* //

src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectNew.C

@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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 "regionModelFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::regionModels::regionModelFunctionObject>
+Foam::regionModels::regionModelFunctionObject::New
+(
+    const dictionary& dict,
+    regionModel& owner,
+    const word& modelName
+)
+{
+    const word modelType = dict.subDict(modelName).lookup("type");
+
+    Info<< "        " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "regionModelFunctionObject::New"
+            "("
+                "const dictionary&, "
+                "regionModel&, "
+                "const word&"
+            ")"
+        )   << "Unknown region model function type "
+            << modelType << nl << nl
+            << "Valid region model function types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return
+        autoPtr<regionModelFunctionObject>
+        (
+            cstrIter()
+            (
+                dict.subDict(modelName),
+                owner
+            )
+        );
+}
+
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/Make/files

@@ -38,12 +38,16 @@ $(THERMOMODELS)/heatTransferModel/mappedConvectiveHeatTransfer/mappedConvectiveH
 $(THERMOMODELS)/filmRadiationModel/filmRadiationModel/filmRadiationModel.C
 $(THERMOMODELS)/filmRadiationModel/filmRadiationModel/filmRadiationModelNew.C
 $(THERMOMODELS)/filmRadiationModel/noRadiation/noRadiation.C
+$(THERMOMODELS)/filmRadiationModel/constantRadiation/constantRadiation.C
+$(THERMOMODELS)/filmRadiationModel/primaryRadiation/primaryRadiation.C
 $(THERMOMODELS)/filmRadiationModel/standardRadiation/standardRadiation.C
 
 
 /* Boundary conditions */
 PATCHFIELDS=derivedFvPatchFields
 $(PATCHFIELDS)/filmHeightInletVelocity/filmHeightInletVelocityFvPatchVectorField.C
+$(PATCHFIELDS)/inclinedFilmNusseltHeight/inclinedFilmNusseltHeightFvPatchScalarField.C
+$(PATCHFIELDS)/inclinedFilmNusseltInletVelocity/inclinedFilmNusseltInletVelocityFvPatchVectorField.C
 
 /* Wall functions for primary region */
 WALLFUNCS=$(PATCHFIELDS)/wallFunctions

src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltHeight/inclinedFilmNusseltHeightFvPatchScalarField.C

@@ -0,0 +1,213 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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 "inclinedFilmNusseltHeightFvPatchScalarField.H"
+#include "volFields.H"
+#include "kinematicSingleLayer.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::inclinedFilmNusseltHeightFvPatchScalarField::
+inclinedFilmNusseltHeightFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    GammaMean_(),
+    a_(),
+    omega_()
+{}
+
+
+Foam::inclinedFilmNusseltHeightFvPatchScalarField::
+inclinedFilmNusseltHeightFvPatchScalarField
+(
+    const inclinedFilmNusseltHeightFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    GammaMean_(ptf.GammaMean_().clone().ptr()),
+    a_(ptf.a_().clone().ptr()),
+    omega_(ptf.omega_().clone().ptr())
+{}
+
+
+Foam::inclinedFilmNusseltHeightFvPatchScalarField::
+inclinedFilmNusseltHeightFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    GammaMean_(DataEntry<scalar>::New("GammaMean", dict)),
+    a_(DataEntry<scalar>::New("a", dict)),
+    omega_(DataEntry<scalar>::New("omega", dict))
+{
+    fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
+}
+
+
+Foam::inclinedFilmNusseltHeightFvPatchScalarField::
+inclinedFilmNusseltHeightFvPatchScalarField
+(
+    const inclinedFilmNusseltHeightFvPatchScalarField& wmfrhpsf
+)
+:
+    fixedValueFvPatchScalarField(wmfrhpsf),
+    GammaMean_(wmfrhpsf.GammaMean_().clone().ptr()),
+    a_(wmfrhpsf.a_().clone().ptr()),
+    omega_(wmfrhpsf.omega_().clone().ptr())
+{}
+
+
+Foam::inclinedFilmNusseltHeightFvPatchScalarField::
+inclinedFilmNusseltHeightFvPatchScalarField
+(
+    const inclinedFilmNusseltHeightFvPatchScalarField& wmfrhpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(wmfrhpsf, iF),
+    GammaMean_(wmfrhpsf.GammaMean_().clone().ptr()),
+    a_(wmfrhpsf.a_().clone().ptr()),
+    omega_(wmfrhpsf.omega_().clone().ptr())
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::inclinedFilmNusseltHeightFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const label patchI = patch().index();
+
+    const scalar t = db().time().timeOutputValue();
+
+    // retrieve the film region from the database
+
+    const regionModels::regionModel& region =
+        db().lookupObject<regionModels::regionModel>("surfaceFilmProperties");
+
+    const regionModels::surfaceFilmModels::kinematicSingleLayer& film =
+        dynamic_cast
+        <
+            const regionModels::surfaceFilmModels::kinematicSingleLayer&
+        >(region);
+
+    // calculate the vector tangential to the patch
+
+    const vectorField n(patch().nf());
+
+    const volVectorField& nHat = film.nHat();
+
+    const vectorField nHatp(nHat.boundaryField()[patchI].patchInternalField());
+
+    vectorField nTan(nHatp ^ n);
+    nTan /= mag(nTan) + ROOTVSMALL;
+
+    // calculate distance in patch tangential direction
+
+    const vectorField& Cf = patch().Cf();
+    scalarField d(nTan & Cf);
+
+    // calculate the wavy film height
+
+    const scalar GMean = GammaMean_->value(t);
+    const scalar a = a_->value(t);
+    const scalar omega = omega_->value(t);
+
+    const scalarField G(GMean + a*sin(omega*constant::mathematical::twoPi*d));
+
+    const volScalarField& mu = film.mu();
+    const scalarField mup(mu.boundaryField()[patchI].patchInternalField());
+
+    const volScalarField& rho = film.rho();
+    const scalarField rhop(rho.boundaryField()[patchI].patchInternalField());
+
+    const scalarField Re(max(G, 0.0)/mup);
+
+    // TODO: currently re-evaluating the entire gTan field to return this patch
+    const scalarField gTan(film.gTan()().boundaryField()[patchI] & n);
+
+    if (max(mag(gTan)) < SMALL)
+    {
+        WarningIn
+        (
+            "void Foam::inclinedFilmNusseltHeightFvPatchScalarField::"
+            "updateCoeffs()"
+        )
+            << "Tangential gravity component is zero.  This boundary condition "
+            << "is designed to operate on patches inclined with respect to "
+            << "gravity"
+            << endl;
+    }
+
+    operator==
+    (
+        pow(3.0*sqr(mup/rhop)/(gTan + ROOTVSMALL), 0.333)*pow(Re, 0.333)
+    );
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::inclinedFilmNusseltHeightFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fixedValueFvPatchScalarField::write(os);
+    GammaMean_->writeData(os);
+    a_->writeData(os);
+    omega_->writeData(os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        inclinedFilmNusseltHeightFvPatchScalarField
+    );
+}
+
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltHeight/inclinedFilmNusseltHeightFvPatchScalarField.H

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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::inclinedFilmNusseltHeightFvPatchScalarField
+
+Description
+    Film height boundary condition for inclined films that imposes a
+    sinusoidal perturbation on top of a mean flow rate, where the height is
+    calculated using the Nusselt solution.
+
+SourceFiles
+    inclinedFilmNusseltHeightFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef inclinedFilmNusseltHeightFvPatchScalarField_H
+#define inclinedFilmNusseltHeightFvPatchScalarField_H
+
+#include "fvPatchFields.H"
+#include "fixedValueFvPatchFields.H"
+#include "DataEntry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+         Class inclinedFilmNusseltHeightFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class inclinedFilmNusseltHeightFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+    // Private data
+
+        //- Mean mass flow rate per unit length [kg/s/m]
+        autoPtr<DataEntry<scalar> > GammaMean_;
+
+        //- Perturbation amplitude [m]
+        autoPtr<DataEntry<scalar> > a_;
+
+        //- Perturbation frequency [rad/s/m]
+        autoPtr<DataEntry<scalar> > omega_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("inclinedFilmNusseltHeight");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        inclinedFilmNusseltHeightFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        inclinedFilmNusseltHeightFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given 
+        // inclinedFilmNusseltHeightFvPatchScalarField onto a new patch
+        inclinedFilmNusseltHeightFvPatchScalarField
+        (
+            const inclinedFilmNusseltHeightFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        inclinedFilmNusseltHeightFvPatchScalarField
+        (
+            const inclinedFilmNusseltHeightFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new inclinedFilmNusseltHeightFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        inclinedFilmNusseltHeightFvPatchScalarField
+        (
+            const inclinedFilmNusseltHeightFvPatchScalarField&,
+            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 inclinedFilmNusseltHeightFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltInletVelocity/inclinedFilmNusseltInletVelocityFvPatchVectorField.C

@@ -0,0 +1,210 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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 "inclinedFilmNusseltInletVelocityFvPatchVectorField.H"
+#include "volFields.H"
+#include "kinematicSingleLayer.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::
+inclinedFilmNusseltInletVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchVectorField(p, iF),
+    GammaMean_(),
+    a_(),
+    omega_()
+{}
+
+
+Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::
+inclinedFilmNusseltInletVelocityFvPatchVectorField
+(
+    const inclinedFilmNusseltInletVelocityFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchVectorField(ptf, p, iF, mapper),
+    GammaMean_(ptf.GammaMean_().clone().ptr()),
+    a_(ptf.a_().clone().ptr()),
+    omega_(ptf.omega_().clone().ptr())
+{}
+
+
+Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::
+inclinedFilmNusseltInletVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchVectorField(p, iF),
+    GammaMean_(DataEntry<scalar>::New("GammaMean", dict)),
+    a_(DataEntry<scalar>::New("a", dict)),
+    omega_(DataEntry<scalar>::New("omega", dict))
+{
+    fvPatchVectorField::operator=(vectorField("value", dict, p.size()));
+}
+
+
+Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::
+inclinedFilmNusseltInletVelocityFvPatchVectorField
+(
+    const inclinedFilmNusseltInletVelocityFvPatchVectorField& fmfrpvf
+)
+:
+    fixedValueFvPatchVectorField(fmfrpvf),
+    GammaMean_(fmfrpvf.GammaMean_().clone().ptr()),
+    a_(fmfrpvf.a_().clone().ptr()),
+    omega_(fmfrpvf.omega_().clone().ptr())
+{}
+
+
+Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::
+inclinedFilmNusseltInletVelocityFvPatchVectorField
+(
+    const inclinedFilmNusseltInletVelocityFvPatchVectorField& fmfrpvf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchVectorField(fmfrpvf, iF),
+    GammaMean_(fmfrpvf.GammaMean_().clone().ptr()),
+    a_(fmfrpvf.a_().clone().ptr()),
+    omega_(fmfrpvf.omega_().clone().ptr())
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const label patchI = patch().index();
+
+    const scalar t = db().time().timeOutputValue();
+
+    // retrieve the film region from the database
+
+    const regionModels::regionModel& region =
+        db().lookupObject<regionModels::regionModel>("surfaceFilmProperties");
+
+    const regionModels::surfaceFilmModels::kinematicSingleLayer& film =
+        dynamic_cast
+        <
+            const regionModels::surfaceFilmModels::kinematicSingleLayer&
+        >(region);
+
+    // calculate the vector tangential to the patch
+
+    const vectorField n(patch().nf());
+
+    const volVectorField& nHat = film.nHat();
+
+    const vectorField nHatp(nHat.boundaryField()[patchI].patchInternalField());
+
+    vectorField nTan(nHatp ^ n);
+    nTan /= mag(nTan) + ROOTVSMALL;
+
+    // calculate distance in patch tangential direction
+
+    const vectorField& Cf = patch().Cf();
+    scalarField d(nTan & Cf);
+
+    // calculate the wavy film height
+
+    const scalar GMean = GammaMean_->value(t);
+    const scalar a = a_->value(t);
+    const scalar omega = omega_->value(t);
+
+    const scalarField G(GMean + a*sin(omega*constant::mathematical::twoPi*d));
+
+    const volScalarField& mu = film.mu();
+    const scalarField mup(mu.boundaryField()[patchI].patchInternalField());
+
+    const volScalarField& rho = film.rho();
+    const scalarField rhop(rho.boundaryField()[patchI].patchInternalField());
+
+    const scalarField Re(max(G, 0.0)/mup);
+
+    // TODO: currently re-evaluating the entire gTan field to return this patch
+    const scalarField gTan(film.gTan()().boundaryField()[patchI] & n);
+
+    if (max(mag(gTan)) < SMALL)
+    {
+        WarningIn
+        (
+            "void Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::"
+            "updateCoeffs()"
+        )
+            << "Tangential gravity component is zero.  This boundary condition "
+            << "is designed to operate on patches inclined with respect to "
+            << "gravity"
+            << endl;
+    }
+
+    operator==(n*pow(gTan*mup/(3.0*rhop), 0.333)*pow(Re, 0.666));
+
+    fixedValueFvPatchVectorField::updateCoeffs();
+}
+
+
+void Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchVectorField::write(os);
+    GammaMean_->writeData(os);
+    a_->writeData(os);
+    omega_->writeData(os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchVectorField,
+        inclinedFilmNusseltInletVelocityFvPatchVectorField
+    );
+}
+
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltInletVelocity/inclinedFilmNusseltInletVelocityFvPatchVectorField.H

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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::inclinedFilmNusseltInletVelocityFvPatchVectorField
+
+Description
+    Film velocity boundary condition for inclined films that imposes a
+    sinusoidal perturbation on top of a mean flow rate, where the velocity is
+    calculated using the Nusselt solution.
+
+SourceFiles
+    inclinedFilmNusseltInletVelocityFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef inclinedFilmNusseltInletVelocityFvPatchVectorField_H
+#define inclinedFilmNusseltInletVelocityFvPatchVectorField_H
+
+#include "fvPatchFields.H"
+#include "fixedValueFvPatchFields.H"
+#include "DataEntry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+     Class inclinedFilmNusseltInletVelocityFvPatchVectorField Declaration
+\*---------------------------------------------------------------------------*/
+
+class inclinedFilmNusseltInletVelocityFvPatchVectorField
+:
+    public fixedValueFvPatchVectorField
+{
+    // Private data
+
+        //- Mean mass flow rate per unit length [kg/s/m]
+        autoPtr<DataEntry<scalar> > GammaMean_;
+
+        //- Perturbation amplitude [m]
+        autoPtr<DataEntry<scalar> > a_;
+
+        //- Perturbation frequency [rad/s/m]
+        autoPtr<DataEntry<scalar> > omega_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("inclinedFilmNusseltInletVelocity");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        inclinedFilmNusseltInletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        inclinedFilmNusseltInletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        // inclinedFilmNusseltInletVelocityFvPatchVectorField onto a new patch
+        inclinedFilmNusseltInletVelocityFvPatchVectorField
+        (
+            const inclinedFilmNusseltInletVelocityFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        inclinedFilmNusseltInletVelocityFvPatchVectorField
+        (
+            const inclinedFilmNusseltInletVelocityFvPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new inclinedFilmNusseltInletVelocityFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        inclinedFilmNusseltInletVelocityFvPatchVectorField
+        (
+            const inclinedFilmNusseltInletVelocityFvPatchVectorField&,
+            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 inclinedFilmNusseltInletVelocityFvPatchVectorField
+                (
+                    *this, iF
+                )
+            );
+        }
+
+
+    // Member functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C

@@ -61,7 +61,7 @@ bool kinematicSingleLayer::read()
     {
         const dictionary& solution = this->solution().subDict("PISO");
         solution.lookup("momentumPredictor") >> momentumPredictor_;
-        solution.lookup("nOuterCorr") >> nOuterCorr_;
+        solution.readIfPresent("nOuterCorr", nOuterCorr_);
         solution.lookup("nCorr") >> nCorr_;
         solution.lookup("nNonOrthCorr") >> nNonOrthCorr_;
 
@@ -197,6 +197,12 @@ tmp<volScalarField> kinematicSingleLayer::pp()
 }
 
 
+void kinematicSingleLayer::correctAlpha()
+{
+    alpha_ == pos(delta_ - dimensionedScalar("SMALL", dimLength, SMALL));
+}
+
+
 void kinematicSingleLayer::updateSubmodels()
 {
     if (debug)
@@ -276,8 +282,8 @@ void kinematicSingleLayer::updateSurfaceVelocities()
     Uw_ -= nHat()*(Uw_ & nHat());
     Uw_.correctBoundaryConditions();
 
-    // TODO: apply quadratic profile to determine surface velocity
+    // apply quadratic profile to surface velocity
-    Us_ = U_;
+    Us_ = 2.0*U_;
     Us_.correctBoundaryConditions();
 }
 
@@ -439,7 +445,7 @@ kinematicSingleLayer::kinematicSingleLayer
     surfaceFilmModel(modelType, mesh, g),
 
     momentumPredictor_(solution().subDict("PISO").lookup("momentumPredictor")),
-    nOuterCorr_(readLabel(solution().subDict("PISO").lookup("nOuterCorr"))),
+    nOuterCorr_(solution().subDict("PISO").lookupOrDefault("nOuterCorr", 1)),
     nCorr_(readLabel(solution().subDict("PISO").lookup("nCorr"))),
     nNonOrthCorr_
     (
@@ -503,6 +509,19 @@ kinematicSingleLayer::kinematicSingleLayer
         ),
         regionMesh()
     ),
+    alpha_
+    (
+        IOobject
+        (
+            "alpha",
+            time().timeName(),
+            regionMesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar("zero", dimless, 0.0)
+    ),
     U_
     (
         IOobject
@@ -817,6 +836,8 @@ void kinematicSingleLayer::preEvolveRegion()
         Info<< "kinematicSingleLayer::preEvolveRegion()" << endl;
     }
 
+    surfaceFilmModel::preEvolveRegion();
+
     transferPrimaryRegionThermoFields();
 
     correctThermoFields();
@@ -838,6 +859,8 @@ void kinematicSingleLayer::evolveRegion()
         Info<< "kinematicSingleLayer::evolveRegion()" << endl;
     }
 
+    correctAlpha();
+
     updateSubmodels();
 
     // Solve continuity for deltaRho_
@@ -846,7 +869,7 @@ void kinematicSingleLayer::evolveRegion()
     // Implicit pressure source coefficient - constant
     tmp<volScalarField> tpp(this->pp());
 
-    for (int oCorr=0; oCorr<nOuterCorr_; oCorr++)
+    for (int oCorr=1; oCorr<=nOuterCorr_; oCorr++)
     {
         // Explicit pressure source contribution - varies with delta_
         tmp<volScalarField> tpu(this->pu());
@@ -1032,7 +1055,9 @@ void kinematicSingleLayer::info() const
         << indent << "min/max(mag(U))    = " << min(mag(U_)).value() << ", "
         << max(mag(U_)).value() << nl
         << indent << "min/max(delta)     = " << min(delta_).value() << ", "
-        << max(delta_).value() << nl;
+        << max(delta_).value() << nl
+        << indent << "coverage           = "
+        << gSum(alpha_.internalField()*magSf())/gSum(magSf()) <<  nl;
 
     injection_.info(Info);
 }

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -113,6 +113,9 @@ protected:
             //- Film thickness / [m]
             volScalarField delta_;
 
+            //- Film coverage indicator, 1 = covered, 0 = uncovered / []
+            volScalarField alpha_;
+
             //- Velocity - mean / [m/s]
             volVectorField U_;
 
@@ -221,12 +224,15 @@ protected:
         //- Transfer source fields from the primary region to the film region
         virtual void transferPrimaryRegionSourceFields();
 
-        // Explicit pressure source contribution
+        //- Explicit pressure source contribution
         virtual tmp<volScalarField> pu();
 
-        // Implicit pressure source coefficient
+        //- Implicit pressure source coefficient
         virtual tmp<volScalarField> pp();
 
+        //- Correct film coverage field
+        virtual void correctAlpha();
+
         //- Update the film sub-models
         virtual void updateSubmodels();
 
@@ -323,6 +329,9 @@ public:
             //- Return const access to the film thickness / [m]
             inline const volScalarField& delta() const;
 
+            //- Return the film coverage, 1 = covered, 0 = uncovered / []
+            inline const volScalarField& alpha() const;
+
             //- Return the film velocity [m/s]
             virtual const volVectorField& U() const;
 

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,6 +79,12 @@ inline const volScalarField& kinematicSingleLayer::delta() const
 }
 
 
+inline const volScalarField& kinematicSingleLayer::alpha() const
+{
+    return alpha_;
+}
+
+
 inline volVectorField& kinematicSingleLayer::USpPrimary()
 {
     return USpPrimary_;
@@ -173,7 +179,8 @@ inline tmp<volScalarField> kinematicSingleLayer::netMass() const
 {
     dimensionedScalar d0("SMALL", dimLength, ROOTVSMALL);
     return
-        fvc::surfaceSum(phi_/(fvc::interpolate(delta_) + d0))*time().deltaT()
+        fvc::surfaceSum(pos(phi_)*phi_/(fvc::interpolate(delta_) + d0))
+       *time().deltaT()
       + rho_*delta_*magSf();
 }
 

src/regionModels/surfaceFilmModels/noFilm/noFilm.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,10 +102,10 @@ const volScalarField& noFilm::delta() const
 }
 
 
-const volScalarField& noFilm::sigma() const
+const volScalarField& noFilm::alpha() const
 {
-    FatalErrorIn("const volScalarField& noFilm::sigma() const")
+    FatalErrorIn("const volScalarField& noFilm::alpha() const")
-        << "sigma field not available for " << type() << abort(FatalError);
+        << "alpha field not available for " << type() << abort(FatalError);
 
     return volScalarField::null();
 }
@@ -192,6 +192,15 @@ const volScalarField& noFilm::kappa() const
 }
 
 
+const volScalarField& noFilm::sigma() const
+{
+    FatalErrorIn("const volScalarField& noFilm::sigma() const")
+        << "sigma field not available for " << type() << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
 tmp<volScalarField> noFilm::primaryMassTrans() const
 {
     return tmp<volScalarField>

src/regionModels/surfaceFilmModels/noFilm/noFilm.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,8 +115,8 @@ public:
             //- Return the film thickness [m]
             virtual const volScalarField& delta() const;
 
-            //- Return const access to the surface tension / [m/s2]
+            //- Return the film coverage, 1 = covered, 0 = uncovered / []
-            inline const volScalarField& sigma() const;
+            virtual const volScalarField& alpha() const;
 
             //- Return the film velocity [m/s]
             virtual const volVectorField& U() const;
@@ -145,6 +145,9 @@ public:
             //- Return the film thermal conductivity [W/m/K]
             virtual const volScalarField& kappa() const;
 
+            //- Return const access to the surface tension / [m/s2]
+            inline const volScalarField& sigma() const;
+
 
             // Transfer fields - to the primary region
 

src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,7 +52,6 @@ contactAngleForce::contactAngleForce
 )
 :
     force(typeName, owner, dict),
-    deltaWet_(readScalar(coeffs_.lookup("deltaWet"))),
     Ccf_(readScalar(coeffs_.lookup("Ccf"))),
     rndGen_(label(0), -1),
     distribution_
@@ -82,7 +81,7 @@ tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
         (
             IOobject
             (
-                "contactForce",
+                typeName + ".contactForce",
                 owner_.time().timeName(),
                 owner_.regionMesh(),
                 IOobject::NO_READ,
@@ -100,17 +99,11 @@ tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
 
     const scalarField& magSf = owner_.magSf();
 
-    const volScalarField& delta = owner_.delta();
+    const volScalarField& alpha = owner_.alpha();
     const volScalarField& sigma = owner_.sigma();
 
-    volScalarField alpha
-    (
-        "alpha",
-        pos(delta - dimensionedScalar("deltaWet", dimLength, deltaWet_))
-    );
     volVectorField gradAlpha(fvc::grad(alpha));
 
-
     scalarField nHits(owner_.regionMesh().nCells(), 0.0);
 
     forAll(nbr, faceI)
@@ -119,19 +112,17 @@ tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
         const label cellN = nbr[faceI];
 
         label cellI = -1;
-        if ((delta[cellO] > deltaWet_) && (delta[cellN] < deltaWet_))
+        if ((alpha[cellO] > 0.5) && (alpha[cellN] < 0.5))
         {
             cellI = cellO;
         }
-        else if ((delta[cellO] < deltaWet_) && (delta[cellN] > deltaWet_))
+        else if ((alpha[cellO] < 0.5) && (alpha[cellN] > 0.5))
         {
             cellI = cellN;
         }
 
         if (cellI != -1)
         {
-//            const scalar dx = Foam::sqrt(magSf[cellI]);
-            // bit of a cheat, but ok for regular meshes
             const scalar dx = owner_.regionMesh().deltaCoeffs()[faceI];
             const vector n =
                 gradAlpha[cellI]/(mag(gradAlpha[cellI]) + ROOTVSMALL);
@@ -141,17 +132,17 @@ tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
         }
     }
 
-    forAll(delta.boundaryField(), patchI)
+    forAll(alpha.boundaryField(), patchI)
     {
-        const fvPatchField<scalar>& df = delta.boundaryField()[patchI];
+        const fvPatchField<scalar>& alphaf = alpha.boundaryField()[patchI];
-        const scalarField& dx = df.patch().deltaCoeffs();
+        const scalarField& dx = alphaf.patch().deltaCoeffs();
-        const labelUList& faceCells = df.patch().faceCells();
+        const labelUList& faceCells = alphaf.patch().faceCells();
 
-        forAll(df, faceI)
+        forAll(alphaf, faceI)
         {
             label cellO = faceCells[faceI];
 
-            if ((delta[cellO] > deltaWet_) && (df[faceI] < deltaWet_))
+            if ((alpha[cellO] > 0.5) && (alphaf[faceI] < 0.5))
             {
                 const vector n =
                     gradAlpha[cellO]/(mag(gradAlpha[cellO]) + ROOTVSMALL);

src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -60,9 +60,6 @@ private:
 
     // Private Data
 
-        //- Threshold film thickness beyon which the film is 'wet'
-        scalar deltaWet_;
-
         //- Coefficient applied to the contact angle force
         scalar Ccf_;
 

src/regionModels/surfaceFilmModels/submodels/kinematic/force/surfaceShearForce/surfaceShearForce.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,11 +80,11 @@ tmp<fvVectorMatrix> surfaceShearForce::correct(volVectorField& U)
     tmp<volScalarField> tCs;
 
     typedef compressible::turbulenceModel turbModel;
-    if (film.primaryMesh().foundObject<turbModel>("turbulenceProperties"))
+    if (film.primaryMesh().foundObject<turbModel>("turbulenceModel"))
     {
         // local reference to turbulence model
         const turbModel& turb =
-            film.primaryMesh().lookupObject<turbModel>("turbulenceProperties");
+            film.primaryMesh().lookupObject<turbModel>("turbulenceModel");
 
         // calculate and store the stress on the primary region
         const volSymmTensorField primaryReff(turb.devRhoReff());
@@ -117,13 +117,23 @@ tmp<fvVectorMatrix> surfaceShearForce::correct(volVectorField& U)
         Reff.correctBoundaryConditions();
 
         dimensionedScalar U0("SMALL", U.dimensions(), SMALL);
-        tCs = Cf_*mag(-film.nHat() & Reff)/(mag(Up - U) + U0);
+        volVectorField UHat("UHat", (Up - U)/(mag(Up - U) + U0));
+
+        // shear stress tangential to the film
+        volVectorField tauTan
+        (
+            "tauTan",
+            UHat & (Reff + film.nHat()*(-film.nHat() & Reff))
+        );
+
+        // note: Cf_ 'should' be 1 in this case
+        tCs = Cf_*mag(tauTan)/(mag(Up - U) + U0);
     }
     else
     {
         // laminar case - employ simple coeff-based model
-        const volScalarField& rho = film.rho();
+        const volScalarField& rhop = film.rhoPrimary();
-        tCs = Cf_*rho*mag(Up - U);
+        tCs = Cf_*rhop*mag(Up - U);
     }
 
     dimensionedScalar d0("SMALL", delta.dimensions(), SMALL);
@@ -131,7 +141,7 @@ tmp<fvVectorMatrix> surfaceShearForce::correct(volVectorField& U)
     // linear coeffs to apply to velocity
     const volScalarField& Cs = tCs();
     volScalarField Cw("Cw", mu/(0.3333*(delta + d0)));
-    Cw.min(1.0e+06);
+    Cw.min(5000.0);
 
     return
     (

src/regionModels/surfaceFilmModels/submodels/subModelBase.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,7 +51,7 @@ subModelBase::subModelBase
 )
 :
     owner_(owner),
-    coeffs_(dict.subDict(type + "Coeffs"))
+    coeffs_(dict.subOrEmptyDict(type + "Coeffs"))
 {}
 
 

src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.C

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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 "constantRadiation.H"
+#include "volFields.H"
+#include "zeroGradientFvPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(constantRadiation, 0);
+
+addToRunTimeSelectionTable
+(
+    filmRadiationModel,
+    constantRadiation,
+    dictionary
+);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+constantRadiation::constantRadiation
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    filmRadiationModel(typeName, owner, dict),
+    QrConst_
+    (
+        IOobject
+        (
+            typeName + "::QrConst",
+            owner.time().timeName(),
+            owner.regionMesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        owner.regionMesh()
+    ),
+    mask_
+    (
+        IOobject
+        (
+            typeName + "::mask",
+            owner.time().timeName(),
+            owner.regionMesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        owner.regionMesh(),
+        dimensionedScalar("one", dimless, 1.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    timeStart_(readScalar(coeffs_.lookup("timeStart"))),
+    duration_(readScalar(coeffs_.lookup("duration")))
+{
+    mask_ = pos(mask_);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+constantRadiation::~constantRadiation()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void constantRadiation::correct()
+{}
+
+
+tmp<volScalarField> constantRadiation::Shs()
+{
+    tmp<volScalarField> tShs
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                typeName + "::Shs",
+                owner().time().timeName(),
+                owner().regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner().regionMesh(),
+            dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    const scalar time = owner().time().value();
+
+    if ((time >= timeStart_) && (time <= timeStart_ + duration_))
+    {
+        scalarField& Shs = tShs();
+        const scalarField& Qr = QrConst_.internalField();
+        const scalarField& alpha = owner_.alpha().internalField();
+
+        Shs = mask_*Qr*alpha;
+    }
+
+    return tShs;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.H

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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::constantRadiation
+
+Description
+    Film constant radiation model.  The constant radiative flux is specified
+    by the user, and operated over a time interval defined by a start time and
+    duration.  In addition, a mask can be applied to shield the film from the
+    radiation.
+
+SourceFiles
+    constantRadiation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constantRadiation_H
+#define constantRadiation_H
+
+#include "filmRadiationModel.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class constantRadiation Declaration
+\*---------------------------------------------------------------------------*/
+
+class constantRadiation
+:
+    public filmRadiationModel
+{
+private:
+
+    // Private data
+
+        //- Constant radiative flux [kg/s3]
+        volScalarField QrConst_;
+
+        //- Radiation mask
+        volScalarField mask_;
+
+        //- Time start [s]
+        const scalar timeStart_;
+
+        //- Duration [s]
+        const scalar duration_;
+
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        constantRadiation(const constantRadiation&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const constantRadiation&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constantRadiation");
+
+
+    // Constructors
+
+        //- Construct from surface film model and dictionary
+        constantRadiation
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~constantRadiation();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct();
+
+            //- Return the radiation sensible enthalpy source
+            //  Also updates QrNet
+            virtual tmp<volScalarField> Shs();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C

@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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 "primaryRadiation.H"
+#include "volFields.H"
+#include "zeroGradientFvPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(primaryRadiation, 0);
+
+addToRunTimeSelectionTable
+(
+    filmRadiationModel,
+    primaryRadiation,
+    dictionary
+);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+primaryRadiation::primaryRadiation
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    filmRadiationModel(typeName, owner, dict),
+    QinPrimary_
+    (
+        IOobject
+        (
+            "Qin", // same name as Qin on primary region to enable mapping
+            owner.time().timeName(),
+            owner.regionMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        owner.regionMesh(),
+        dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
+        owner.mappedPushedFieldPatchTypes<scalar>()
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+primaryRadiation::~primaryRadiation()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void primaryRadiation::correct()
+{
+    // Transfer Qin from primary region
+    QinPrimary_.correctBoundaryConditions();
+}
+
+
+tmp<volScalarField> primaryRadiation::Shs()
+{
+    tmp<volScalarField> tShs
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                typeName + "::Shs",
+                owner().time().timeName(),
+                owner().regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner().regionMesh(),
+            dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    scalarField& Shs = tShs();
+    const scalarField& QinP = QinPrimary_.internalField();
+    const scalarField& alpha = owner_.alpha().internalField();
+
+    Shs = QinP*alpha;
+
+    return tShs;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.H

@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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::primaryRadiation
+
+Description
+    Radiation model whereby the radiative heat flux is mapped from the primary
+    region
+
+SourceFiles
+    primaryRadiation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef primaryRadiation_H
+#define primaryRadiation_H
+
+#include "filmRadiationModel.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class primaryRadiation Declaration
+\*---------------------------------------------------------------------------*/
+
+class primaryRadiation
+:
+    public filmRadiationModel
+{
+private:
+
+    // Private data
+
+        //- Incident radiative flux mapped from  the primary region / [kg/s3]
+        volScalarField QinPrimary_;
+
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        primaryRadiation(const primaryRadiation&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const primaryRadiation&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("primaryRadiation");
+
+
+    // Constructors
+
+        //- Construct from surface film model and dictionary
+        primaryRadiation
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~primaryRadiation();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct();
+
+            //- Return the radiation sensible enthalpy source
+            //  Also updates QrNet
+            virtual tmp<volScalarField> Shs();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C

@@ -85,8 +85,6 @@ standardRadiation::standardRadiation
         dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
         zeroGradientFvPatchScalarField::typeName
     ),
-    delta_(owner.delta()),
-    deltaMin_(readScalar(coeffs_.lookup("deltaMin"))),
     beta_(readScalar(coeffs_.lookup("beta"))),
     kappaBar_(readScalar(coeffs_.lookup("kappaBar")))
 {}
@@ -129,9 +127,10 @@ tmp<volScalarField> standardRadiation::Shs()
 
     scalarField& Shs = tShs();
     const scalarField& QinP = QinPrimary_.internalField();
-    const scalarField& delta = delta_.internalField();
+    const scalarField& delta = owner_.delta().internalField();
+    const scalarField& alpha = owner_.alpha().internalField();
 
-    Shs = beta_*(QinP*pos(delta - deltaMin_))*(1.0 - exp(-kappaBar_*delta));
+    Shs = beta_*QinP*alpha*(1.0 - exp(-kappaBar_*delta));
 
     // Update net Qr on local region
     QrNet_.internalField() = QinP - Shs;

src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.H

@@ -65,15 +65,9 @@ private:
         //- Remaining radiative flux after removing local contribution
         volScalarField QrNet_;
 
-        //- Reference to the film thickness field / [m]
-        const volScalarField& delta_;
-
 
         // Model coefficients
 
-            //- Minimum thickness to apply radiation model
-            scalar deltaMin_;
-
             //- Beta coefficient
             scalar beta_;
 

src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,7 +77,6 @@ standardPhaseChange::standardPhaseChange
 )
 :
     phaseChangeModel(typeName, owner, dict),
-    Tb_(readScalar(coeffs_.lookup("Tb"))),
     deltaMin_(readScalar(coeffs_.lookup("deltaMin"))),
     L_(readScalar(coeffs_.lookup("L"))),
     TbFactor_(coeffs_.lookupOrDefault<scalar>("TbFactor", 1.1))
@@ -113,14 +112,10 @@ void standardPhaseChange::correctModel
     const scalarField& YInf = film.YPrimary()[vapId];
     const scalarField& pInf = film.pPrimary();
     const scalarField& T = film.T();
-    const scalarField& Tw = film.Tw();
     const scalarField& rho = film.rho();
-    const scalarField& TInf = film.TPrimary();
     const scalarField& rhoInf = film.rhoPrimary();
     const scalarField& muInf = film.muPrimary();
     const scalarField& magSf = film.magSf();
-    const scalarField hInf(film.htcs().h());
-    const scalarField hFilm(film.htcw().h());
     const vectorField dU(film.UPrimary() - film.Us());
     const scalarField limMass
     (
@@ -134,8 +129,11 @@ void standardPhaseChange::correctModel
             // cell pressure [Pa]
             const scalar pc = pInf[cellI];
 
+            // calculate the boiling temperature
+            const scalar Tb = liq.pvInvert(pc);
+
             // local temperature - impose lower limit of 200 K for stability
-            const scalar Tloc = min(TbFactor_*Tb_, max(200.0, T[cellI]));
+            const scalar Tloc = min(TbFactor_*Tb, max(200.0, T[cellI]));
 
             // saturation pressure [Pa]
             const scalar pSat = liq.pv(pc, Tloc);
@@ -147,15 +145,10 @@ void standardPhaseChange::correctModel
             if (pSat >= 0.95*pc)
             {
                 // boiling
-                const scalar qDotInf = hInf[cellI]*(TInf[cellI] - T[cellI]);
-                const scalar qDotFilm = hFilm[cellI]*(T[cellI] - Tw[cellI]);
-
                 const scalar Cp = liq.Cp(pc, Tloc);
-                const scalar Tcorr = max(0.0, T[cellI] - Tb_);
+                const scalar Tcorr = max(0.0, T[cellI] - Tb);
                 const scalar qCorr = limMass[cellI]*Cp*(Tcorr);
-                dMass[cellI] =
+                dMass[cellI] = qCorr/hVap;
-                    dt*magSf[cellI]/hVap*(qDotInf + qDotFilm)
-                  + qCorr/hVap;
             }
             else
             {

src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,9 +69,6 @@ protected:
 
     // Protected data
 
-        //- Boiling temperature / [K]
-        const scalar Tb_;
-
         //- Minimum film height for model to be active
         const scalar deltaMin_;
 

src/regionModels/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -186,6 +186,9 @@ public:
             //- Return the film thickness [m]
             virtual const volScalarField& delta() const = 0;
 
+            //- Return the film coverage, 1 = covered, 0 = uncovered / []
+            virtual const volScalarField& alpha() const = 0;
+
             //- Return the film velocity [m/s]
             virtual const volVectorField& U() const = 0;
 

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C

@@ -227,6 +227,33 @@ void thermoSingleLayer::transferPrimaryRegionSourceFields()
 }
 
 
+void thermoSingleLayer::correctAlpha()
+{
+    if (hydrophilic_)
+    {
+        const scalar hydrophilicDry = hydrophilicDryScale_*deltaWet_;
+        const scalar hydrophilicWet = hydrophilicWetScale_*deltaWet_;
+
+        forAll(alpha_, i)
+        {
+            if ((alpha_[i] < 0.5) && (delta_[i] > hydrophilicDry))
+            {
+                alpha_[i] = 1.0;
+            }
+            else if ((alpha_[i] > 0.5) && (delta_[i] < hydrophilicWet))
+            {
+                alpha_[i] = 0.0;
+            }
+        }
+    }
+    else
+    {
+        alpha_ ==
+            pos(delta_ - dimensionedScalar("deltaWet", dimLength, deltaWet_));
+    }
+}
+
+
 void thermoSingleLayer::updateSubmodels()
 {
     if (debug)
@@ -290,7 +317,8 @@ void thermoSingleLayer::solveEnergy()
       - hsSp_
       + q(hs_)
       + radiation_->Shs()
-      - fvm::SuSp(rhoSp_, hs_)
+//      - fvm::SuSp(rhoSp_, hs_)
+      - rhoSp_*hs_
     );
 
     correctThermoFields();
@@ -426,6 +454,11 @@ thermoSingleLayer::thermoSingleLayer
         zeroGradientFvPatchScalarField::typeName
     ),
 
+    deltaWet_(readScalar(coeffs_.lookup("deltaWet"))),
+    hydrophilic_(readBool(coeffs_.lookup("hydrophilic"))),
+    hydrophilicDryScale_(0.0),
+    hydrophilicWetScale_(0.0),
+
     hsSp_
     (
         IOobject
@@ -481,8 +514,20 @@ thermoSingleLayer::thermoSingleLayer
         heatTransferModel::New(*this, coeffs().subDict("lowerSurfaceModels"))
     ),
     phaseChange_(phaseChangeModel::New(*this, coeffs())),
-    radiation_(filmRadiationModel::New(*this, coeffs()))
+    radiation_(filmRadiationModel::New(*this, coeffs())),
+    Tmin_(-VGREAT),
+    Tmax_(VGREAT)
 {
+    if (coeffs().readIfPresent("Tmin", Tmin_))
+    {
+        Info<< "    limiting minimum temperature to " << Tmin_ << endl;
+    }
+
+    if (coeffs().readIfPresent("Tmax", Tmax_))
+    {
+        Info<< "    limiting maximum temperature to " << Tmax_ << endl;
+    }
+
     if (thermo_.hasMultiComponentCarrier())
     {
         YPrimary_.setSize(thermo_.carrier().species().size());
@@ -510,6 +555,12 @@ thermoSingleLayer::thermoSingleLayer
         }
     }
 
+    if (hydrophilic_)
+    {
+        coeffs_.lookup("hydrophilicDryScale") >> hydrophilicDryScale_;
+        coeffs_.lookup("hydrophilicWetScale") >> hydrophilicWetScale_;
+    }
+
     if (readFields)
     {
         transferPrimaryRegionThermoFields();
@@ -584,12 +635,14 @@ void thermoSingleLayer::evolveRegion()
         Info<< "thermoSingleLayer::evolveRegion()" << endl;
     }
 
+    correctAlpha();
+
     updateSubmodels();
 
     // Solve continuity for deltaRho_
     solveContinuity();
 
-    for (int oCorr=0; oCorr<nOuterCorr_; oCorr++)
+    for (int oCorr=1; oCorr<=nOuterCorr_; oCorr++)
     {
         // Explicit pressure source contribution
         tmp<volScalarField> tpu(this->pu());

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -123,6 +123,22 @@ protected:
                 volScalarField primaryEnergyPCTrans_;
 
 
+        //- Threshold film thickness beyond which the film is considered 'wet'
+        scalar deltaWet_;
+
+
+        // Hyprophilic/phobic properties
+
+            //- Activation flag
+            bool hydrophilic_;
+
+            //- Length scale applied to deltaWet_ for dry faces, typically 0.5
+            scalar hydrophilicDryScale_;
+
+            //- Length scale applied to deltaWet_ for wet faces, typically 0.001
+            scalar hydrophilicWetScale_;
+
+
         // Source term fields
 
             // Film region - registered to the film region mesh
@@ -166,6 +182,14 @@ protected:
             autoPtr<filmRadiationModel> radiation_;
 
 
+        // Limits
+
+            //- Minimum temperature limit (optional)
+            scalar Tmin_;
+
+            //- Maximum temperature limit (optional)
+            scalar Tmax_;
+
 
     // Protected member functions
 
@@ -190,6 +214,9 @@ protected:
         //- Transfer source fields from the primary region to the film region
         virtual void transferPrimaryRegionSourceFields();
 
+        //- Correct film coverage field
+        virtual void correctAlpha();
+
         //- Update the film sub-models
         virtual void updateSubmodels();
 
@@ -342,6 +369,15 @@ public:
             inline const filmRadiationModel& radiation() const;
 
 
+        // Derived fields (calculated on-the-fly)
+
+            //- Return the convective heat energy from film to wall
+            inline tmp<scalarField> Qconvw(const label patchI) const;
+
+            //- Return the convective heat energy from primary region to film
+            inline tmp<scalarField> Qconvp(const label patchI) const;
+
+
         // Evolution
 
             //- Pre-evolve film hook

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "thermoSingleLayer.H"
+#include "heatTransferModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -88,7 +89,7 @@ inline tmp<volScalarField> thermoSingleLayer::T
     const volScalarField& hs
 ) const
 {
-    return tmp<volScalarField>
+    tmp<volScalarField> tT
     (
         new volScalarField
         (
@@ -104,6 +105,11 @@ inline tmp<volScalarField> thermoSingleLayer::T
             zeroGradientFvPatchScalarField::typeName
         )
     );
+
+    tT().min(Tmax_);
+    tT().max(Tmin_);
+
+    return tT;
 }
 
 
@@ -155,6 +161,26 @@ inline const filmRadiationModel& thermoSingleLayer::radiation() const
 }
 
 
+inline tmp<scalarField> thermoSingleLayer::Qconvw(const label patchI) const
+{
+    const scalarField htc(htcw_->h()().boundaryField()[patchI]);
+    const scalarField& Tp = T_.boundaryField()[patchI];
+    const scalarField& Twp = Tw_.boundaryField()[patchI];
+
+    return htc*(Tp - Twp);
+}
+
+
+inline tmp<scalarField> thermoSingleLayer::Qconvp(const label patchI) const
+{
+    const scalarField htc(htcs_->h()().boundaryField()[patchI]);
+    const scalarField& Tp = T_.boundaryField()[patchI];
+    const scalarField& Tpp = TPrimary_.boundaryField()[patchI];
+
+    return htc*(Tp - Tpp);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace surfaceFilmModels

src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJump/energyJumpFvPatchScalarField.C

@@ -68,6 +68,10 @@ Foam::energyJumpFvPatchScalarField::energyJumpFvPatchScalarField
             scalarField("value", dict, p.size())
         );
     }
+    else
+    {
+        evaluate(Pstream::blocking);
+    }
 }
 
 

src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJumpAMI/energyJumpAMIFvPatchScalarField.C

@@ -68,6 +68,10 @@ Foam::energyJumpAMIFvPatchScalarField::energyJumpAMIFvPatchScalarField
             scalarField("value", dict, p.size())
         );
     }
+    else
+    {
+        evaluate(Pstream::blocking);
+    }
 }
 
 

src/thermophysicalModels/basic/heThermo/heThermo.C

@@ -35,6 +35,39 @@ License
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+template<class BasicThermo, class MixtureType>
+Foam::wordList Foam::heThermo<BasicThermo, MixtureType>::heBoundaryBaseTypes()
+{
+    const volScalarField::GeometricBoundaryField& tbf =
+        this->T_.boundaryField();
+
+    wordList hbt(tbf.size(), word::null);
+
+    forAll(tbf, patchi)
+    {
+        if (isA<fixedJumpFvPatchScalarField>(tbf[patchi]))
+        {
+            const fixedJumpFvPatchScalarField& pf =
+                dynamic_cast<const fixedJumpFvPatchScalarField&>(tbf[patchi]);
+
+            hbt[patchi] = pf.interfaceFieldType();
+        }
+        else if (isA<fixedJumpAMIFvPatchScalarField>(tbf[patchi]))
+        {
+            const fixedJumpAMIFvPatchScalarField& pf =
+                dynamic_cast<const fixedJumpAMIFvPatchScalarField&>
+                (
+                    tbf[patchi]
+                );
+
+            hbt[patchi] = pf.interfaceFieldType();
+        }
+    }
+
+    return hbt;
+}
+
+
 template<class BasicThermo, class MixtureType>
 Foam::wordList Foam::heThermo<BasicThermo, MixtureType>::heBoundaryTypes()
 {
@@ -149,7 +182,8 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo
         ),
         mesh,
         dimEnergy/dimMass,
-        this->heBoundaryTypes()
+        this->heBoundaryTypes(),
+        this->heBoundaryBaseTypes()
     )
 {
     init();
@@ -179,7 +213,8 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo
         ),
         mesh,
         dimEnergy/dimMass,
-        this->heBoundaryTypes()
+        this->heBoundaryTypes(),
+        this->heBoundaryBaseTypes()
     )
 {
     init();

src/thermophysicalModels/basic/heThermo/heThermo.H

@@ -68,6 +68,10 @@ protected:
             //  by interrogating the temperature field boundary types
             wordList heBoundaryTypes();
 
+            //- Return the enthalpy/internal energy field boundary base types
+            //  by interrogating the temperature field boundary types
+            wordList heBoundaryBaseTypes();
+
             //- Correct the enthalpy/internal energy field boundaries
             void heBoundaryCorrection(volScalarField& he);
 

tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/MRFProperties

@@ -19,6 +19,8 @@ zone1
 {
     cellZone    rotor;
 
+    active      on;
+
     // Fixed patches (by default they 'move' with the MRF zone)
     nonRotatingPatches ();
 

tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicitFixedCoeff/constant/sourcesProperties

@@ -26,6 +26,7 @@ source1
 
     fixedTemperatureSourceCoeffs
     {
+        mode            constant;
         temperature     350;
     }
 }

tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/sourcesProperties

@@ -26,6 +26,7 @@ source1
 
     fixedTemperatureSourceCoeffs
     {
+        mode            constant;
         temperature     350;
     }
 }

tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/sourcesProperties

@@ -26,6 +26,7 @@ source1
 
     fixedTemperatureSourceCoeffs
     {
+        mode            constant;
         temperature     2000;
     }
 }

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/surfaceFilmProperties

@@ -26,6 +26,9 @@ thermoSingleLayerCoeffs
     thermoModel singleComponent;
     liquid      H2O;
 
+    deltaWet    1e-4;
+    hydrophilic no;
+
     forces
     (
         surfaceShear

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/system/wallFilmRegion/fvSolution

@@ -37,8 +37,7 @@ solvers
 PISO
 {
     momentumPredictor true;
-    nOuterCorr      1;
+    nCorr           1;
-    nCorr           2;
     nNonOrthCorr    0;
 }
 

tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/surfaceFilmProperties

@@ -28,6 +28,9 @@ thermoSingleLayerCoeffs
 
     liquid      H2O;
 
+    deltaWet    1e-4;
+    hydrophilic no;
+
     radiationModel  none;
 
     upperSurfaceModels
@@ -65,7 +68,6 @@ thermoSingleLayerCoeffs
 
     contactAngleCoeffs
     {
-        deltaWet        1e-4;
         Ccf             0.085;
         contactAngleDistribution
         {

tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/system/wallFilmRegion.org/fvSolution

@@ -39,8 +39,7 @@ solvers
 PISO
 {
     momentumPredictor true;
-    nOuterCorr      1;
+    nCorr           1;
-    nCorr           2;
     nNonOrthCorr    0;
 }
 

tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/constant/surfaceFilmProperties

@@ -43,7 +43,6 @@ kinematicSingleLayerCoeffs
 
     contactAngleCoeffs
     {
-        deltaWet        1e-4;
         Ccf             1;
         contactAngleDistribution
         {

tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/system/wallFilmRegion/fvSolution

@@ -37,8 +37,7 @@ solvers
 PISO
 {
     momentumPredictor  true;
-    nOuterCorr      1;
+    nCorr           1;
-    nCorr           2;
     nNonOrthCorr    0;
 }
 

tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/surfaceFilmProperties

@@ -27,6 +27,9 @@ thermoSingleLayerCoeffs
 
     liquid      H2O;
 
+    deltaWet    1e-4;
+    hydrophilic no;
+
     forces
     (
         surfaceShear

tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/system/wallFilmRegion/fvSolution

@@ -37,8 +37,7 @@ solvers
 PISO
 {
     momentumPredictor true;
-    nOuterCorr      1;
+    nCorr           1;
-    nCorr           2;
     nNonOrthCorr    0;
 }