Merge commit 'OpenCFD/master' into olesenm

applications/utilities/mesh/generation/blockMesh/curvedEdges/arcEdge.C

@@ -30,8 +30,7 @@ Description
 
 #include "arcEdge.H"
 #include "mathematicalConstants.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -40,8 +39,7 @@ namespace Foam
     defineTypeNameAndDebug(arcEdge, 0);
 
     // Add the curvedEdge constructor functions to the hash tables
-    curvedEdge::addIstreamConstructorToTable<arcEdge>
-        addArcEdgeIstreamConstructorToTable_;
+    addToRunTimeSelectionTable(curvedEdge, arcEdge, Istream);
 }
 
 

applications/utilities/mesh/generation/blockMesh/curvedEdges/curvedEdge.C

@@ -41,26 +41,7 @@ namespace Foam
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 defineTypeNameAndDebug(curvedEdge, 0);
-
-// Define the constructor function hash tables
-HashTable<curvedEdge::IstreamConstructorPtr_>*
-    curvedEdge::IstreamConstructorTablePtr_;
-
-
-// Hash table Constructor called from the table add functions.
-
-void curvedEdge::constructTables()
-{
-    static bool constructed = false;
-
-    if (!constructed)
-    {
-        curvedEdge::IstreamConstructorTablePtr_
-            = new HashTable<curvedEdge::IstreamConstructorPtr_>;
-
-        constructed = true;
-    }
-}
+defineRunTimeSelectionTable(curvedEdge, Istream);
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -117,10 +98,11 @@ autoPtr<curvedEdge> curvedEdge::New(const pointField& points, Istream& is)
 
     word curvedEdgeType(is);
 
-    HashTable<IstreamConstructorPtr_>::iterator curvedEdgeConstructorIter =
-        IstreamConstructorTablePtr_->find(curvedEdgeType);
+    IstreamConstructorTable::iterator cstrIter =
+        IstreamConstructorTablePtr_
+            ->find(curvedEdgeType);
 
-    if (curvedEdgeConstructorIter == IstreamConstructorTablePtr_->end())
+    if (cstrIter == IstreamConstructorTablePtr_->end())
     {
         FatalErrorIn("curvedEdge::New(const pointField&, Istream&)")
             << "Unknown curvedEdge type " << curvedEdgeType << endl << endl
@@ -129,7 +111,7 @@ autoPtr<curvedEdge> curvedEdge::New(const pointField& points, Istream& is)
             << abort(FatalError);
     }
 
-    return autoPtr<curvedEdge>(curvedEdgeConstructorIter()(points, is));
+    return autoPtr<curvedEdge>(cstrIter()(points, is));
 }
 
 
@@ -177,7 +159,6 @@ Ostream& operator<<(Ostream& os, const curvedEdge& p)
 }
 
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

applications/utilities/mesh/generation/blockMesh/curvedEdges/curvedEdge.H

@@ -63,51 +63,23 @@ protected:
 
 public:
 
-    // Constructor Hash tables
-
-        //- Construct from Istream function pointer type
-        typedef autoPtr<curvedEdge> (*IstreamConstructorPtr_)
-            (const pointField&, Istream&);
-
-        //- Construct from Istream function pointer table pointer
-        static HashTable<IstreamConstructorPtr_>*
-            IstreamConstructorTablePtr_;
+    //- Runtime type information
+    TypeName("curvedEdge");
 
 
-    // Hash table constructor classes and functions
+    // Declare run-time constructor selection tables
 
-        //- Hash table Constructor.
-        //  Must be called from the table add functions below.
-        static void constructTables();
-
-
-        //- Class to add constructor from Istream to Hash table
-        template<class curvedEdgeType>
-        class addIstreamConstructorToTable
-        {
-        public:
-
-            static autoPtr<curvedEdge> New
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            curvedEdge,
+            Istream,
             (
                 const pointField& points,
                 Istream& is
-            )
-            {
-                return autoPtr<curvedEdge>(new curvedEdgeType(points, is));
-            }
-
-            addIstreamConstructorToTable()
-            {
-                curvedEdge::constructTables();
-
-                curvedEdge::IstreamConstructorTablePtr_
-                    ->insert(curvedEdgeType::typeName, New);
-            }
-        };
-
-
-    //- Runtime type information
-    TypeName("curvedEdge");
+            ),
+            (points, is)
+        );
 
 
     // Constructors

applications/utilities/mesh/generation/blockMesh/curvedEdges/polySplineEdge.C

@@ -26,6 +26,7 @@ License
 
 #include "polySplineEdge.H"
 #include "BSpline.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -34,8 +35,7 @@ namespace Foam
     defineTypeNameAndDebug(polySplineEdge, 0);
 
     // Add the curvedEdge constructor functions to the hash tables
-    curvedEdge::addIstreamConstructorToTable<polySplineEdge>
-        addPolySplineEdgeIstreamConstructorToTable_;
+    addToRunTimeSelectionTable(curvedEdge, polySplineEdge, Istream);
 }
 
 

applications/utilities/mesh/generation/blockMesh/curvedEdges/polySplineEdge.H

@@ -92,7 +92,7 @@ public:
         );
 
         //- Construct from Istream setting pointsList
-        polySplineEdge(const pointField& points,Istream& is);
+        polySplineEdge(const pointField& points, Istream& is);
 
 
     // Destructor

applications/utilities/mesh/generation/blockMesh/curvedEdges/simpleSplineEdge.C

@@ -27,9 +27,8 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "error.H"
-
 #include "simpleSplineEdge.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -39,10 +38,8 @@ namespace Foam
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 defineTypeNameAndDebug(simpleSplineEdge, 0);
+addToRunTimeSelectionTable(curvedEdge, simpleSplineEdge, Istream);
 
-// Add the curvedEdge constructor functions to the hash tables
-curvedEdge::addIstreamConstructorToTable<simpleSplineEdge>
-    addSimpleSplineEdgeIstreamConstructorToTable_;
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 

applications/utilities/postProcessing/sampling/probeLocations/probeLocations.C

@@ -39,10 +39,12 @@ using namespace Foam;
 int main(int argc, char *argv[])
 {
     timeSelector::addOptions();
+#   include "addRegionOption.H"
+
 #   include "setRootCase.H"
 #   include "createTime.H"
     instantList timeDirs = timeSelector::select0(runTime, args);
-#   include "createMesh.H"
+#   include "createNamedMesh.H"
 
     IOprobes sniff
     (

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C

@@ -143,7 +143,7 @@ void Foam::processorPolyPatch::initGeometry()
         (
             Pstream::blocking,
             neighbProcNo(),
-            3*(sizeof(label) + size()*sizeof(vector) + sizeof(float))
+            3*(sizeof(label) + size()*sizeof(vector) + sizeof(scalar))
         );
 
         toNeighbProc
@@ -163,7 +163,7 @@ void Foam::processorPolyPatch::calcGeometry()
             (
                 Pstream::blocking,
                 neighbProcNo(),
-                3*(sizeof(label) + size()*sizeof(vector) + sizeof(float))
+                3*(sizeof(label) + size()*sizeof(vector) + sizeof(scalar))
             );
             fromNeighbProc
                 >> neighbFaceCentres_

src/fvMotionSolver/Make/files

@@ -1,5 +1,6 @@
 fvMotionSolvers/fvMotionSolver/fvMotionSolver.C
 fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
+fvMotionSolvers/displacement/displacementFvMotionSolver/displacementFvMotionSolver.C
 fvMotionSolvers/displacement/interpolation/displacementInterpolationFvMotionSolver.C
 fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
 fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C

src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C

@@ -55,25 +55,10 @@ namespace Foam
 Foam::displacementSBRStressFvMotionSolver::displacementSBRStressFvMotionSolver
 (
     const polyMesh& mesh,
-    Istream&
+    Istream& is
 )
 :
-    fvMotionSolver(mesh),
-    points0_
-    (
-        pointIOField
-        (
-            IOobject
-            (
-                "points",
-                time().constant(),
-                polyMesh::meshSubDir,
-                mesh,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE
-            )
-        )
-    ),
+    displacementFvMotionSolver(mesh, is),
     pointDisplacement_
     (
         IOobject
@@ -132,7 +117,7 @@ Foam::displacementSBRStressFvMotionSolver::curPoints() const
 
     tmp<pointField> tcurPoints
     (
-        points0_ + pointDisplacement_.internalField()
+        points0() + pointDisplacement_.internalField()
     );
 
     twoDCorrectPoints(tcurPoints());
@@ -208,63 +193,7 @@ void Foam::displacementSBRStressFvMotionSolver::updateMesh
     const mapPolyMesh& mpm
 )
 {
-    fvMotionSolver::updateMesh(mpm);
-
-    // Map points0_
-    // Map points0_. Bit special since we somehow have to come up with
-    // a sensible points0 position for introduced points.
-    // Find out scaling between points0 and current points
-
-    // Get the new points either from the map or the mesh
-    const pointField& points =
-    (
-        mpm.hasMotionPoints()
-      ? mpm.preMotionPoints()
-      : fvMesh_.points()
-    );
-
-    // Note: boundBox does reduce
-    const vector span0 = boundBox(points0_).span();
-    const vector span  = boundBox(points).span();
-
-    vector scaleFactors(cmptDivide(span0, span));
-
-    pointField newPoints0(mpm.pointMap().size());
-
-    forAll(newPoints0, pointI)
-    {
-        label oldPointI = mpm.pointMap()[pointI];
-
-        if (oldPointI >= 0)
-        {
-            label masterPointI = mpm.reversePointMap()[oldPointI];
-
-            if (masterPointI == pointI)
-            {
-                newPoints0[pointI] = points0_[oldPointI];
-            }
-            else
-            {
-                // New point. Assume motion is scaling.
-                newPoints0[pointI] = points0_[oldPointI] + cmptMultiply
-                (
-                    scaleFactors,
-                    points[pointI]-points[masterPointI]
-                );
-            }
-        }
-        else
-        {
-            FatalErrorIn
-            (
-                "displacementSBRStressFvMotionSolver::updateMesh"
-                "(const mapPolyMesh& mpm)"
-            )   << "Cannot work out coordinates of introduced vertices."
-                << " New vertex " << pointI << " at coordinate "
-                << points[pointI] << exit(FatalError);
-        }
-    }
-    points0_.transfer(newPoints0);
+    displacementFvMotionSolver::updateMesh(mpm);
 
     // Update diffusivity. Note two stage to make sure old one is de-registered
     // before creating/registering new one.

src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.H

@@ -37,7 +37,7 @@ SourceFiles
 #ifndef displacementSBRStressFvMotionSolver_H
 #define displacementSBRStressFvMotionSolver_H
 
-#include "fvMotionSolver.H"
+#include "displacementFvMotionSolver.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,13 +53,10 @@ class motionDiffusivity;
 
 class displacementSBRStressFvMotionSolver
 :
-    public fvMotionSolver
+    public displacementFvMotionSolver
 {
     // Private data
 
-        //- Reference point field
-        pointField points0_;
-
         //- Point motion field
         mutable pointVectorField pointDisplacement_;
 
@@ -105,12 +102,6 @@ public:
 
     // Member Functions
 
-        //- Return reference to the reference field
-        const pointField& points0() const
-        {
-            return points0_;
-        }
-
         //- Return reference to the point motion displacement field
         pointVectorField& pointDisplacement()
         {

src/fvMotionSolver/fvMotionSolvers/displacement/displacementFvMotionSolver/displacementFvMotionSolver.C

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "displacementFvMotionSolver.H"
+#include "addToRunTimeSelectionTable.H"
+#include "mapPolyMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+//    defineTypeNameAndDebug(displacementFvMotionSolver, 0);
+//
+//    addToRunTimeSelectionTable
+//    (
+//        fvMotionSolver,
+//        displacementFvMotionSolver,
+//        dictionary
+//    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::displacementFvMotionSolver::
+displacementFvMotionSolver
+(
+    const polyMesh& mesh,
+    Istream&
+)
+:
+    fvMotionSolver(mesh),
+    points0_
+    (
+        pointIOField
+        (
+            IOobject
+            (
+                "points",
+                mesh.time().constant(),
+                polyMesh::meshSubDir,
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        )
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::displacementFvMotionSolver::~displacementFvMotionSolver()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::displacementFvMotionSolver::updateMesh(const mapPolyMesh& mpm)
+{
+    fvMotionSolver::updateMesh(mpm);
+
+    // Map points0_. Bit special since we somehow have to come up with
+    // a sensible points0 position for introduced points.
+    // Find out scaling between points0 and current points
+
+    // Get the new points either from the map or the mesh
+    const pointField& points =
+    (
+        mpm.hasMotionPoints()
+      ? mpm.preMotionPoints()
+      : fvMesh_.points()
+    );
+
+    // Note: boundBox does reduce
+    const vector span0 = boundBox(points0_).span();
+    const vector span  = boundBox(points).span();
+
+    vector scaleFactors(cmptDivide(span0, span));
+
+    pointField newPoints0(mpm.pointMap().size());
+
+    forAll(newPoints0, pointI)
+    {
+        label oldPointI = mpm.pointMap()[pointI];
+
+        if (oldPointI >= 0)
+        {
+            label masterPointI = mpm.reversePointMap()[oldPointI];
+
+            if (masterPointI == pointI)
+            {
+                newPoints0[pointI] = points0_[oldPointI];
+            }
+            else
+            {
+                // New point. Assume motion is scaling.
+                newPoints0[pointI] = points0_[oldPointI] + cmptMultiply
+                (
+                    scaleFactors,
+                    points[pointI]-points[masterPointI]
+                );
+            }
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "displacementLaplacianFvMotionSolver::updateMesh"
+                "(const mapPolyMesh& mpm)"
+            )   << "Cannot work out coordinates of introduced vertices."
+                << " New vertex " << pointI << " at coordinate "
+                << points[pointI] << exit(FatalError);
+        }
+    }
+    points0_.transfer(newPoints0);
+}
+
+
+// ************************************************************************* //

src/fvMotionSolver/fvMotionSolvers/displacement/displacementFvMotionSolver/displacementFvMotionSolver.H

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::displacementFvMotionSolver.H
+
+Description
+    Base class for fvMotionSolvers which calculate displacement.
+
+SourceFiles
+    displacementFvMotionSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef displacementFvMotionSolver_H
+#define displacementFvMotionSolver_H
+
+#include "fvMotionSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+        Class displacementFvMotionSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class displacementFvMotionSolver
+:
+    public fvMotionSolver
+{
+    // Private data
+
+        //- Reference point field
+        pointField points0_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        displacementFvMotionSolver
+        (
+            const displacementFvMotionSolver&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const displacementFvMotionSolver&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("displacementInterpolation");
+
+
+    // Constructors
+
+        //- Construct from polyMesh and data stream
+        displacementFvMotionSolver
+        (
+            const polyMesh&,
+            Istream& msDataUnused
+        );
+
+
+    // Destructor
+
+        ~displacementFvMotionSolver();
+
+
+    // Member Functions
+
+        //- Return reference to the reference field
+        const pointField& points0() const
+        {
+            return points0_;
+        }
+
+        //- Update topology
+        virtual void updateMesh(const mapPolyMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationFvMotionSolver.C

@@ -58,26 +58,10 @@ Foam::displacementInterpolationFvMotionSolver::
 displacementInterpolationFvMotionSolver
 (
     const polyMesh& mesh,
-    Istream&
+    Istream& is
 )
 :
-    fvMotionSolver(mesh),
-    points0_
-    (
-        pointIOField
-        (
-            IOobject
-            (
-                "points",
-                mesh.time().constant(),
-                polyMesh::meshSubDir,
-                mesh,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE,
-                false
-            )
-        )
-    ),
+    displacementFvMotionSolver(mesh, is),
     dynamicMeshCoeffs_
     (
         IOdictionary
@@ -174,7 +158,7 @@ displacementInterpolationFvMotionSolver
             forAll(fz().meshPoints(), localI)
             {
                 label pointI = fz().meshPoints()[localI];
-                const scalar coord = points0_[pointI][dir];
+                const scalar coord = points0()[pointI][dir];
                 minCoord = min(minCoord, coord);
                 maxCoord = max(maxCoord, coord);
             }
@@ -198,7 +182,7 @@ displacementInterpolationFvMotionSolver
         zoneCoordinates[zoneCoordinates.size()-1] += SMALL;
 
         // Check if we have static min and max mesh bounds
-        const scalarField meshCoords = points0_.component(dir);
+        const scalarField meshCoords = points0().component(dir);
 
         scalar minCoord = gMin(meshCoords);
         scalar maxCoord = gMax(meshCoords);
@@ -288,7 +272,7 @@ displacementInterpolationFvMotionSolver
                     "displacementInterpolationFvMotionSolver::"
                     "displacementInterpolationFvMotionSolver"
                     "(const polyMesh&, Istream&)"
-                )   << "Did not find point " << points0_[pointI]
+                )   << "Did not find point " << points0()[pointI]
                     << " coordinate " << meshCoords[pointI]
                     << " in ranges " << rangeToCoord
                     << abort(FatalError);
@@ -344,18 +328,18 @@ Foam::displacementInterpolationFvMotionSolver::
 Foam::tmp<Foam::pointField>
 Foam::displacementInterpolationFvMotionSolver::curPoints() const
 {
-    if (mesh().nPoints() != points0_.size())
+    if (mesh().nPoints() != points0().size())
     {
         FatalErrorIn
         (
             "displacementInterpolationFvMotionSolver::curPoints() const"
         )   << "The number of points in the mesh seems to have changed." << endl
-            << "In constant/polyMesh there are " << points0_.size()
+            << "In constant/polyMesh there are " << points0().size()
             << " points; in the current mesh there are " << mesh().nPoints()
             << " points." << exit(FatalError);
     }
 
-    tmp<pointField> tcurPoints(new pointField(points0_));
+    tmp<pointField> tcurPoints(new pointField(points0()));
     pointField& curPoints = tcurPoints();
 
     // Interpolate the diplacement of the face zones.
@@ -413,68 +397,4 @@ Foam::displacementInterpolationFvMotionSolver::curPoints() const
 }
 
 
-void Foam::displacementInterpolationFvMotionSolver::updateMesh
-(
-    const mapPolyMesh& mpm
-)
-{
-    fvMotionSolver::updateMesh(mpm);
-
-    // Map points0_. Bit special since we somehow have to come up with
-    // a sensible points0 position for introduced points.
-    // Find out scaling between points0 and current points
-
-    // Get the new points either from the map or the mesh
-    const pointField& points =
-    (
-        mpm.hasMotionPoints()
-      ? mpm.preMotionPoints()
-      : fvMesh_.points()
-    );
-
-    // Note: boundBox does reduce
-    const vector span0 = boundBox(points0_).span();
-    const vector span  = boundBox(points).span();
-
-    vector scaleFactors(cmptDivide(span0, span));
-
-    pointField newPoints0(mpm.pointMap().size());
-
-    forAll(newPoints0, pointI)
-    {
-        label oldPointI = mpm.pointMap()[pointI];
-
-        if (oldPointI >= 0)
-        {
-            label masterPointI = mpm.reversePointMap()[oldPointI];
-
-            if (masterPointI == pointI)
-            {
-                newPoints0[pointI] = points0_[oldPointI];
-            }
-            else
-            {
-                // New point. Assume motion is scaling.
-                newPoints0[pointI] = points0_[oldPointI] + cmptMultiply
-                (
-                    scaleFactors,
-                    points[pointI]-points[masterPointI]
-                );
-            }
-        }
-        else
-        {
-            FatalErrorIn
-            (
-                "displacementLaplacianFvMotionSolver::updateMesh"
-                "(const mapPolyMesh& mpm)"
-            )   << "Cannot work out coordinates of introduced vertices."
-                << " New vertex " << pointI << " at coordinate "
-                << points[pointI] << exit(FatalError);
-        }
-    }
-    points0_.transfer(newPoints0);
-}
-
-
 // ************************************************************************* //

src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationFvMotionSolver.H

@@ -48,7 +48,7 @@ SourceFiles
 #ifndef displacementInterpolationFvMotionSolver_H
 #define displacementInterpolationFvMotionSolver_H
 
-#include "fvMotionSolver.H"
+#include "displacementFvMotionSolver.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,13 +61,10 @@ namespace Foam
 
 class displacementInterpolationFvMotionSolver
 :
-    public fvMotionSolver
+    public displacementFvMotionSolver
 {
     // Private data
 
-        //- Reference point field
-        pointField points0_;
-
         //- Additional settings for motion solver
         dictionary dynamicMeshCoeffs_;
 
@@ -130,21 +127,12 @@ public:
 
     // Member Functions
 
-        //- Return reference to the reference field
-        const pointField& points0() const
-        {
-            return points0_;
-        }
-
         //- Return point location obtained from the current motion field
         virtual tmp<pointField> curPoints() const;
 
         //- Solve for motion
         virtual void solve()
         {}
-
-        //- Update topology
-        virtual void updateMesh(const mapPolyMesh&);
 };
 
 

src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C

@@ -53,26 +53,10 @@ namespace Foam
 Foam::displacementLaplacianFvMotionSolver::displacementLaplacianFvMotionSolver
 (
     const polyMesh& mesh,
-    Istream&
+    Istream& is
 )
 :
-    fvMotionSolver(mesh),
-    points0_
-    (
-        pointIOField
-        (
-            IOobject
-            (
-                "points",
-                time().constant(),
-                polyMesh::meshSubDir,
-                mesh,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE,
-                false
-            )
-        )
-    ),
+    displacementFvMotionSolver(mesh, is),
     pointDisplacement_
     (
         IOobject
@@ -186,7 +170,7 @@ Foam::displacementLaplacianFvMotionSolver::curPoints() const
         }
 
         pointLocation_().internalField() =
-            points0_
+            points0()
           + pointDisplacement_.internalField();
 
         pointLocation_().correctBoundaryConditions();
@@ -198,7 +182,7 @@ Foam::displacementLaplacianFvMotionSolver::curPoints() const
 
             forAll(pz, i)
             {
-                pointLocation_()[pz[i]] = points0_[pz[i]];
+                pointLocation_()[pz[i]] = points0()[pz[i]];
             }
         }
 
@@ -210,7 +194,7 @@ Foam::displacementLaplacianFvMotionSolver::curPoints() const
     {
         tmp<pointField> tcurPoints
         (
-            points0_ + pointDisplacement_.internalField()
+            points0() + pointDisplacement_.internalField()
         );
 
         // Implement frozen points
@@ -220,7 +204,7 @@ Foam::displacementLaplacianFvMotionSolver::curPoints() const
 
             forAll(pz, i)
             {
-                tcurPoints()[pz[i]] = points0_[pz[i]];
+                tcurPoints()[pz[i]] = points0()[pz[i]];
             }
         }
 
@@ -257,74 +241,7 @@ void Foam::displacementLaplacianFvMotionSolver::updateMesh
     const mapPolyMesh& mpm
 )
 {
-    fvMotionSolver::updateMesh(mpm);
-
-    // Map points0_. Bit special since we somehow have to come up with
-    // a sensible points0 position for introduced points.
-    // Find out scaling between points0 and current points
-
-    // Get the new points either from the map or the mesh
-    const pointField& points =
-    (
-        mpm.hasMotionPoints()
-      ? mpm.preMotionPoints()
-      : fvMesh_.points()
-    );
-
-    // Note: boundBox does reduce
-    const vector span0 = boundBox(points0_).span();
-    const vector span  = boundBox(points).span();
-
-    vector scaleFactors(cmptDivide(span0, span));
-
-    pointField newPoints0(mpm.pointMap().size());
-
-    forAll(newPoints0, pointI)
-    {
-        label oldPointI = mpm.pointMap()[pointI];
-
-        if (oldPointI >= 0)
-        {
-            label masterPointI = mpm.reversePointMap()[oldPointI];
-
-            if (masterPointI == pointI)
-            {
-                newPoints0[pointI] = points0_[oldPointI];
-            }
-            else
-            {
-                // New point. Assume motion is scaling.
-                newPoints0[pointI] = points0_[oldPointI] + cmptMultiply
-                (
-                    scaleFactors,
-                    points[pointI]-points[masterPointI]
-                );
-            }
-        }
-        else
-        {
-            FatalErrorIn
-            (
-                "displacementLaplacianFvMotionSolver::updateMesh"
-                "(const mapPolyMesh& mpm)"
-            )   << "Cannot work out coordinates of introduced vertices."
-                << " New vertex " << pointI << " at coordinate "
-                << points[pointI] << exit(FatalError);
-        }
-    }
-    points0_.transfer(newPoints0);
-
-    if (debug & 2)
-    {
-        OFstream str(time().timePath()/"points0.obj");
-        Pout<< "displacementLaplacianFvMotionSolver :"
-            << " Writing points0_ to " << str.name() << endl;
-
-        forAll(points0_, pointI)
-        {
-            meshTools::writeOBJ(str, points0_[pointI]);
-        }
-    }
+    displacementFvMotionSolver::updateMesh(mpm);
 
     // Update diffusivity. Note two stage to make sure old one is de-registered
     // before creating/registering new one.

src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.H

@@ -37,7 +37,7 @@ SourceFiles
 #ifndef displacementLaplacianFvMotionSolver_H
 #define displacementLaplacianFvMotionSolver_H
 
-#include "fvMotionSolver.H"
+#include "displacementFvMotionSolver.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,13 +53,10 @@ class motionDiffusivity;
 
 class displacementLaplacianFvMotionSolver
 :
-    public fvMotionSolver
+    public displacementFvMotionSolver
 {
     // Private data
 
-        //- Reference point field
-        pointField points0_;
-
         //- Point motion field
         mutable pointVectorField pointDisplacement_;
 
@@ -113,13 +110,6 @@ public:
 
     // Member Functions
 
-
-        //- Return reference to the reference field
-        const pointField& points0() const
-        {
-            return points0_;
-        }
-
         //- Return reference to the point motion displacement field
         pointVectorField& pointDisplacement()
         {

src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C

@@ -29,7 +29,7 @@ License
 #include "Time.H"
 #include "transformField.H"
 #include "fvMesh.H"
-#include "displacementLaplacianFvMotionSolver.H"
+#include "displacementFvMotionSolver.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,6 +52,254 @@ const NamedEnum<surfaceSlipDisplacementPointPatchVectorField::followMode, 3>
     surfaceSlipDisplacementPointPatchVectorField::followModeNames_;
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void surfaceSlipDisplacementPointPatchVectorField::calcProjection
+(
+    vectorField& displacement
+) const
+{
+    const polyMesh& mesh = patch().boundaryMesh().mesh()();
+    const pointField& localPoints = patch().localPoints();
+    const labelList& meshPoints = patch().meshPoints();
+
+    //const scalar deltaT = mesh.time().deltaT().value();
+
+    // Construct large enough vector in direction of projectDir so
+    // we're guaranteed to hit something.
+
+    //- Per point projection vector:
+    const scalar projectLen = mag(mesh.bounds().max()-mesh.bounds().min());
+
+    // For case of fixed projection vector:
+    vector projectVec;
+    if (projectMode_ == FIXEDNORMAL)
+    {
+        vector n = projectDir_/mag(projectDir_);
+        projectVec = projectLen*n;
+    }
+
+
+    // Get fixed points (bit of a hack)
+    const pointZone* zonePtr = NULL;
+
+    if (frozenPointsZone_.size() > 0)
+    {
+        const pointZoneMesh& pZones = mesh.pointZones();
+
+        zonePtr = &pZones[pZones.findZoneID(frozenPointsZone_)];
+
+        Pout<< "surfaceSlipDisplacementPointPatchVectorField : Fixing all "
+            << zonePtr->size() << " points in pointZone " << zonePtr->name()
+            << endl;
+    }
+
+    // Get the starting locations from the motionSolver
+    const displacementFvMotionSolver& motionSolver =
+        mesh.lookupObject<displacementFvMotionSolver>
+        (
+            "dynamicMeshDict"
+        );
+    const pointField& points0 = motionSolver.points0();
+
+
+    pointField start(meshPoints.size());
+    forAll(start, i)
+    {
+        start[i] = points0[meshPoints[i]] + displacement[i];
+    }
+
+    label nNotProjected = 0;
+
+    if (projectMode_ == NEAREST)
+    {
+        List<pointIndexHit> nearest;
+        labelList hitSurfaces;
+        surfaces().findNearest
+        (
+            start,
+            scalarField(start.size(), sqr(projectLen)),
+            hitSurfaces,
+            nearest
+        );
+
+        forAll(nearest, i)
+        {
+            if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
+            {
+                // Fixed point. Reset to point0 location.
+                displacement[i] = points0[meshPoints[i]] - localPoints[i];
+            }
+            else if (nearest[i].hit())
+            {
+                displacement[i] =
+                    nearest[i].hitPoint()
+                  - points0[meshPoints[i]];
+            }
+            else
+            {
+                nNotProjected++;
+
+                if (debug)
+                {
+                    Pout<< "    point:" << meshPoints[i]
+                        << " coord:" << localPoints[i]
+                        << "  did not find any surface within " << projectLen
+                        << endl;
+                }
+            }
+        }
+    }
+    else
+    {
+        // Do tests on all points. Combine later on.
+
+        // 1. Check if already on surface
+        List<pointIndexHit> nearest;
+        {
+            labelList nearestSurface;
+            surfaces().findNearest
+            (
+                start,
+                scalarField(start.size(), sqr(SMALL)),
+                nearestSurface,
+                nearest
+            );
+        }
+
+        // 2. intersection. (combined later on with information from nearest
+        // above)
+        vectorField projectVecs(start.size(), projectVec);
+
+        if (projectMode_ == POINTNORMAL)
+        {
+            projectVecs = projectLen*patch().pointNormals();
+        }
+
+        // Knock out any wedge component
+        scalarField offset(start.size(), 0.0);
+        if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
+        {
+            forAll(offset, i)
+            {
+                offset[i] = start[i][wedgePlane_];
+                start[i][wedgePlane_] = 0;
+                projectVecs[i][wedgePlane_] = 0;
+            }
+        }
+
+        List<pointIndexHit> rightHit;
+        {
+            labelList rightSurf;
+            surfaces().findAnyIntersection
+            (
+                start,
+                start+projectVecs,
+                rightSurf,
+                rightHit
+            );
+        }
+        
+        List<pointIndexHit> leftHit;
+        {
+            labelList leftSurf;
+            surfaces().findAnyIntersection
+            (
+                start,
+                start-projectVecs,
+                leftSurf,
+                leftHit
+            );
+        }
+
+        // 3. Choose either -fixed, nearest, right, left.
+        forAll(displacement, i)
+        {
+            if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
+            {
+                // Fixed point. Reset to point0 location.
+                displacement[i] = points0[meshPoints[i]] - localPoints[i];
+            }
+            else if (nearest[i].hit())
+            {
+                // Found nearest.
+                displacement[i] =
+                    nearest[i].hitPoint()
+                  - points0[meshPoints[i]];
+            }
+            else
+            {
+                pointIndexHit interPt;
+
+                if (rightHit[i].hit())
+                {
+                    if (leftHit[i].hit())
+                    {
+                        if
+                        (
+                            magSqr(rightHit[i].hitPoint()-start[i])
+                          < magSqr(leftHit[i].hitPoint()-start[i])
+                        )
+                        {
+                            interPt = rightHit[i];
+                        }
+                        else
+                        {
+                            interPt = leftHit[i];
+                        }
+                    }
+                    else
+                    {
+                        interPt = rightHit[i];
+                    }
+                }
+                else
+                {
+                    if (leftHit[i].hit())
+                    {
+                        interPt = leftHit[i];
+                    }
+                }
+
+
+                if (interPt.hit())
+                {
+                    if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
+                    {
+                        interPt.rawPoint()[wedgePlane_] += offset[i];
+                    }
+                    displacement[i] = interPt.rawPoint()-points0[meshPoints[i]];
+                }
+                else
+                {
+                    nNotProjected++;
+
+                    if (debug)
+                    {
+                        Pout<< "    point:" << meshPoints[i]
+                            << " coord:" << localPoints[i]
+                            << "  did not find any intersection between"
+                            << " ray from " << start[i]-projectVecs[i]
+                            << " to " << start[i]+projectVecs[i] << endl;
+                    }
+                }
+            }
+        }
+    }
+
+    reduce(nNotProjected, sumOp<label>());
+
+    if (nNotProjected > 0)
+    {
+        Info<< "surfaceSlipDisplacement :"
+            << " on patch " << patch().name()
+            << " did not project " << nNotProjected
+            << " out of " << returnReduce(localPoints.size(), sumOp<label>())
+            << " points." << endl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 surfaceSlipDisplacementPointPatchVectorField::
@@ -81,7 +329,7 @@ surfaceSlipDisplacementPointPatchVectorField
     projectMode_(followModeNames_.read(dict.lookup("followMode"))),
     projectDir_(dict.lookup("projectDirection")),
     wedgePlane_(readLabel(dict.lookup("wedgePlane"))),
-    frozenPointsZone_(dict.lookup("frozenPointsZone"))
+    frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
 {}
 
 
@@ -167,229 +415,10 @@ void surfaceSlipDisplacementPointPatchVectorField::evaluate
     const Pstream::commsTypes commsType
 )
 {
-    const polyMesh& mesh = patch().boundaryMesh().mesh()();
-
-    // const scalar deltaT = mesh.time().deltaT().value();
-
-    // Construct large enough vector in direction of projectDir so
-    // we're guaranteed to hit something.
-
-    const scalar projectLen = mesh.bounds().mag();
-
-    // For case of fixed projection vector:
-    vector projectVec;
-    if (projectMode_ == FIXEDNORMAL)
-    {
-        vector n = projectDir_/mag(projectDir_);
-        projectVec = projectLen*n;
-    }
-
-    //- Per point projection vector:
-
-    const pointField& localPoints = patch().localPoints();
-    const labelList& meshPoints = patch().meshPoints();
-
     vectorField displacement(this->patchInternalField());
 
-
-    // Get fixed points (bit of a hack)
-    const pointZone* zonePtr = NULL;
-
-    if (frozenPointsZone_.size())
-    {
-        const pointZoneMesh& pZones = mesh.pointZones();
-
-        zonePtr = &pZones[pZones.findZoneID(frozenPointsZone_)];
-
-        Pout<< "surfaceSlipDisplacementPointPatchVectorField : Fixing all "
-            << zonePtr->size() << " points in pointZone " << zonePtr->name()
-            << endl;
-    }
-
-    // Get the starting locations from the motionSolver
-    const displacementLaplacianFvMotionSolver& motionSolver =
-        mesh.lookupObject<displacementLaplacianFvMotionSolver>
-        (
-            "dynamicMeshDict"
-        );
-    const pointField& points0 = motionSolver.points0();
-
-
-//XXXXXX
-
-
-    pointField start(meshPoints.size());
-    forAll(start, i)
-    {
-        start[i] = points0[meshPoints[i]] + displacement[i];
-    }
-
-    if (projectMode_ == NEAREST)
-    {
-        List<pointIndexHit> nearest;
-        labelList hitSurfaces;
-        surfaces().findNearest
-        (
-            start,
-            scalarField(start.size(), sqr(projectLen)),
-            hitSurfaces,
-            nearest
-        );
-
-        forAll(nearest, i)
-        {
-            if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
-            {
-                // Fixed point. Reset to point0 location.
-                displacement[i] = points0[meshPoints[i]] - localPoints[i];
-            }
-            else if (nearest[i].hit())
-            {
-                displacement[i] =
-                    nearest[i].hitPoint()
-                  - points0[meshPoints[i]];
-            }
-            else
-            {
-                Pout<< "    point:" << meshPoints[i]
-                    << " coord:" << localPoints[i]
-                    << "  did not find any surface within " << projectLen
-                    << endl;
-            }
-        }
-    }
-    else
-    {
-        // Do tests on all points. Combine later on.
-
-        // 1. Check if already on surface
-        List<pointIndexHit> nearest;
-        {
-            labelList nearestSurface;
-            surfaces().findNearest
-            (
-                start,
-                scalarField(start.size(), sqr(SMALL)),
-                nearestSurface,
-                nearest
-            );
-        }
-
-        // 2. intersection. (combined later on with information from nearest
-        // above)
-        vectorField projectVecs(start.size(), projectVec);
-
-        if (projectMode_ == POINTNORMAL)
-        {
-            projectVecs = projectLen*patch().pointNormals();
-        }
-
-        // Knock out any wedge component
-        scalarField offset(start.size(), 0.0);
-        if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
-        {
-            forAll(offset, i)
-            {
-                offset[i] = start[i][wedgePlane_];
-                start[i][wedgePlane_] = 0;
-                projectVecs[i][wedgePlane_] = 0;
-            }
-        }
-
-        List<pointIndexHit> rightHit;
-        {
-            labelList rightSurf;
-            surfaces().findAnyIntersection
-            (
-                start,
-                start+projectVecs,
-                rightSurf,
-                rightHit
-            );
-        }
-
-        List<pointIndexHit> leftHit;
-        {
-            labelList leftSurf;
-            surfaces().findAnyIntersection
-            (
-                start,
-                start-projectVecs,
-                leftSurf,
-                leftHit
-            );
-        }
-
-        // 3. Choose either -fixed, nearest, right, left.
-        forAll(displacement, i)
-        {
-            if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
-            {
-                // Fixed point. Reset to point0 location.
-                displacement[i] = points0[meshPoints[i]] - localPoints[i];
-            }
-            else if (nearest[i].hit())
-            {
-                // Found nearest.
-                displacement[i] =
-                    nearest[i].hitPoint()
-                  - points0[meshPoints[i]];
-            }
-            else
-            {
-                pointIndexHit interPt;
-
-                if (rightHit[i].hit())
-                {
-                    if (leftHit[i].hit())
-                    {
-                        if
-                        (
-                            magSqr(rightHit[i].hitPoint()-start[i])
-                          < magSqr(leftHit[i].hitPoint()-start[i])
-                        )
-                        {
-                            interPt = rightHit[i];
-                        }
-                        else
-                        {
-                            interPt = leftHit[i];
-                        }
-                    }
-                    else
-                    {
-                        interPt = rightHit[i];
-                    }
-                }
-                else
-                {
-                    if (leftHit[i].hit())
-                    {
-                        interPt = leftHit[i];
-                    }
-                }
-
-
-                if (interPt.hit())
-                {
-                    if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
-                    {
-                        interPt.rawPoint()[wedgePlane_] += offset[i];
-                    }
-                    displacement[i] = interPt.rawPoint()-points0[meshPoints[i]];
-                }
-                else
-                {
-                    Pout<< "    point:" << meshPoints[i]
-                        << " coord:" << localPoints[i]
-                        << "  did not find any intersection between ray from "
-                        << start[i]-projectVecs[i]
-                        << " to " << start[i]+projectVecs[i]
-                        << endl;
-                }
-            }
-        }
-    }
+    // Calculate displacement to project points onto surface
+    calcProjection(displacement);
 
     // Get internal field to insert values into
     Field<vector>& iF = const_cast<Field<vector>&>(this->internalField());
@@ -412,8 +441,11 @@ void surfaceSlipDisplacementPointPatchVectorField::write(Ostream& os) const
         << token::END_STATEMENT << nl;
     os.writeKeyword("wedgePlane") << wedgePlane_
         << token::END_STATEMENT << nl;
-    os.writeKeyword("frozenPointsZone") << frozenPointsZone_
-        << token::END_STATEMENT << nl;
+    if (frozenPointsZone_ != word::null)
+    {
+        os.writeKeyword("frozenPointsZone") << frozenPointsZone_
+            << token::END_STATEMENT << nl;
+    }
 }
 
 

src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.H

@@ -107,6 +107,9 @@ private:
 
     // Private Member Functions
 
+        //- Calculate displacement to project onto surface
+        void calcProjection(vectorField& displacement) const;
+
         //- Disallow default bitwise assignment
         void operator=(const surfaceSlipDisplacementPointPatchVectorField&);
 

src/meshTools/indexedOctree/indexedOctree.C

@@ -824,8 +824,10 @@ Foam::direction Foam::indexedOctree<Type>::getFace
 }
 
 
-// Traverse a node. If intersects a triangle return first intersection point.
-// Else return the bounxing box face hit:
+// Traverse a node. If intersects a triangle return first intersection point:
+//  hitInfo.index = index of shape
+//  hitInfo.point = point on shape
+// Else return a miss and the bounding box face hit:
 //  hitInfo.point = coordinate of intersection of ray with bounding box
 //  faceID  = index of bounding box face
 template <class Type>
@@ -918,10 +920,10 @@ void Foam::indexedOctree<Type>::traverseNode
         {
             faceID = 0;
 
-            WarningIn("indexedOctree<Type>::traverseNode")
+            FatalErrorIn("indexedOctree<Type>::traverseNode(..)")
                 << "Did not hit side of box " << subBb
                 << " with ray from " << start << " to " << end
-                << endl;
+                << abort(FatalError);
         }
         else
         {
@@ -936,10 +938,10 @@ void Foam::indexedOctree<Type>::traverseNode
         {
             faceID = 0;
 
-            WarningIn("indexedOctree<Type>::traverseNode")
+            FatalErrorIn("indexedOctree<Type>::traverseNode(..)")
                 << "Did not hit side of box " << subBb
                 << " with ray from " << start << " to " << end
-                << endl;
+                << abort(FatalError);
         }
         else
         {

src/sampling/probes/probes.C

@@ -182,15 +182,24 @@ bool Foam::probes::checkFieldTypes()
     if (Pstream::master())
     {
         fileName probeDir;
+
+        fileName probeSubDir = name_;
+
+        if (obr_.name() != polyMesh::defaultRegion)
+        {
+            probeSubDir = probeSubDir/obr_.name();
+        }
+        probeSubDir = probeSubDir/obr_.time().timeName();
+
         if (Pstream::parRun())
         {
             // Put in undecomposed case
             // (Note: gives problems for distributed data running)
-            probeDir = obr_.time().path()/".."/name_/obr_.time().timeName();
+            probeDir = obr_.time().path()/".."/probeSubDir;
         }
         else
         {
-            probeDir = obr_.time().path()/name_/obr_.time().timeName();
+            probeDir = obr_.time().path()/probeSubDir;
         }
 
         // Close the file if any fields have been removed.

src/sampling/sampledSurface/isoSurface/isoSurface.H

@@ -439,8 +439,8 @@ public:
             return triPointMergeMap_;
         }
 
-        //- Interpolates cCoords,pCoords. Takes the original fields
-        //  used to create the iso surface.
+        //- Interpolates cCoords,pCoords. Uses the references to the original
+        //  fields used to create the iso surface.
         template <class Type>
         tmp<Field<Type> > interpolate
         (

src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C

@@ -153,6 +153,16 @@ void Foam::sampledIsoSurface::getIsoFields() const
             }
         }
 
+
+        // If averaging redo the volField. Can only be done now since needs the
+        // point field.
+        if (average_)
+        {
+            storedVolFieldPtr_.reset(average(fvm, *pointFieldPtr_).ptr());
+            volFieldPtr_ = storedVolFieldPtr_.operator->();
+        }
+
+
         if (debug)
         {
             Info<< "sampledIsoSurface::getIsoField() : volField "
@@ -241,6 +251,20 @@ void Foam::sampledIsoSurface::getIsoFields() const
             pointSubFieldPtr_ = storedPointSubFieldPtr_.operator->();
         }
 
+
+
+        // If averaging redo the volField. Can only be done now since needs the
+        // point field.
+        if (average_)
+        {
+            storedVolSubFieldPtr_.reset
+            (
+                average(subFvm, *pointSubFieldPtr_).ptr()
+            );
+            volSubFieldPtr_ = storedVolSubFieldPtr_.operator->();
+        }
+
+
         if (debug)
         {
             Info<< "sampledIsoSurface::getIsoField() : volSubField "
@@ -394,68 +418,33 @@ bool Foam::sampledIsoSurface::updateGeometry() const
     surfPtr_.clear();
     facesPtr_.clear();
 
-    if (average_)
+    if (subMeshPtr_.valid())
     {
-        if (subMeshPtr_.valid())
-        {
-            surfPtr_.reset
+        surfPtr_.reset
+        (
+            new isoSurface
             (
-                new isoSurface
-                (
-                    average(subMeshPtr_().subMesh(), *pointSubFieldPtr_),
-                    *pointSubFieldPtr_,
-                    isoVal_,
-                    regularise_,
-                    mergeTol_
-                )
-            );
-        }
-        else
-        {
-            surfPtr_.reset
-            (
-                new isoSurface
-                (
-                    average(fvm, *pointFieldPtr_),
-                    *pointFieldPtr_,
-                    isoVal_,
-                    regularise_,
-                    mergeTol_
-                )
-            );
-        }
+                *volSubFieldPtr_,
+                *pointSubFieldPtr_,
+                isoVal_,
+                regularise_,
+                mergeTol_
+            )
+        );
     }
     else
     {
-        if (subMeshPtr_.valid())
-        {
-            surfPtr_.reset
+        surfPtr_.reset
+        (
+            new isoSurface
             (
-                new isoSurface
-                (
-                    *volSubFieldPtr_,
-                    *pointSubFieldPtr_,
-                    isoVal_,
-                    regularise_,
-                    mergeTol_
-                )
-            );
-        }
-        else
-        {
-            surfPtr_.reset
-            (
-                new isoSurface
-                (
-                    *volFieldPtr_,
-                    *pointFieldPtr_,
-                    //average(pointMesh::New(mesh()), *volFieldPtr_),
-                    isoVal_,
-                    regularise_,
-                    mergeTol_
-                )
-            );
-        }
+                *volFieldPtr_,
+                *pointFieldPtr_,
+                isoVal_,
+                regularise_,
+                mergeTol_
+            )
+        );
     }
 
 

src/thermophysicalModels/radiation/Make/files

@@ -1,13 +1,15 @@
 /* Radiation constants */
 radiationConstants/radiationConstants.C
 
-
 /* Radiation model */
 radiationModel/radiationModel/radiationModel.C
 radiationModel/radiationModel/newRadiationModel.C
 radiationModel/noRadiation/noRadiation.C
 radiationModel/P1/P1.C
-
+radiationModel/fvDOM/fvDOM/fvDOM.C
+radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
+radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C
+radiationModel/fvDOM/absorptionCoeffs/absorptionCoeffs.C
 
 /* Scatter model */
 submodels/scatterModel/scatterModel/scatterModel.C
@@ -21,11 +23,14 @@ submodels/absorptionEmissionModel/absorptionEmissionModel/newAbsorptionEmissionM
 submodels/absorptionEmissionModel/noAbsorptionEmission/noAbsorptionEmission.C
 submodels/absorptionEmissionModel/constantAbsorptionEmission/constantAbsorptionEmission.C
 submodels/absorptionEmissionModel/binaryAbsorptionEmission/binaryAbsorptionEmission.C
+submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
+submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
 
 
 /* Boundary conditions */
 derivedFvPatchFields/MarshakRadiation/MarshakRadiationMixedFvPatchScalarField.C
 derivedFvPatchFields/MarshakRadiationFixedT/MarshakRadiationFixedTMixedFvPatchScalarField.C
-
+derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
+derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
 
 LIB = $(FOAM_LIBBIN)/libradiation

src/thermophysicalModels/radiation/Make/options

@@ -1,6 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/OpenFOAM/lnInclude \
+    -I radiationModel/fvDOM/interpolationLookUpTable
 
 LIB_LIBS = \
     -lfiniteVolume

src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C

@@ -0,0 +1,273 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "greyDiffusiveRadiationMixedFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+
+#include "fvDOM.H"
+#include "radiationConstants.H"
+#include "mathematicalConstants.H"
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::
+greyDiffusiveRadiationMixedFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    TName_("undefinedT"),
+    emissivity_(0.0),
+    rayId_(0),
+    lambdaId_(0)
+{
+    refValue() = 0.0;
+    refGrad() = 0.0;
+    valueFraction() = 1.0;
+}
+
+
+Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::
+greyDiffusiveRadiationMixedFvPatchScalarField
+(
+    const greyDiffusiveRadiationMixedFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchScalarField(ptf, p, iF, mapper),
+    TName_(ptf.TName_),
+    emissivity_(ptf.emissivity_),
+    rayId_(ptf.rayId_),
+    lambdaId_(ptf.lambdaId_)
+{}
+
+
+Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::
+greyDiffusiveRadiationMixedFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    TName_(dict.lookup("T")),
+    emissivity_(readScalar(dict.lookup("emissivity"))),
+    rayId_(-1),
+    lambdaId_(-1)
+{
+    const scalarField& Tp =
+        patch().lookupPatchField<volScalarField, scalar>(TName_);
+
+    refValue() =
+        emissivity_*4.0*radiation::sigmaSB.value()*pow4(Tp)
+       /Foam::mathematicalConstant::pi;
+
+    refGrad() = 0.0;
+
+    if (dict.found("value"))
+    {
+        fvPatchScalarField::operator=
+        (
+            scalarField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        fvPatchScalarField::operator=(refValue());
+    }
+}
+
+
+Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::
+greyDiffusiveRadiationMixedFvPatchScalarField
+(
+    const greyDiffusiveRadiationMixedFvPatchScalarField& ptf
+)
+:
+    mixedFvPatchScalarField(ptf),
+    TName_(ptf.TName_),
+    emissivity_(ptf.emissivity_),
+    rayId_(ptf.rayId_),
+    lambdaId_(ptf.lambdaId_)
+{}
+
+
+Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::
+greyDiffusiveRadiationMixedFvPatchScalarField
+(
+    const greyDiffusiveRadiationMixedFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(ptf, iF),
+    TName_(ptf.TName_),
+    emissivity_(ptf.emissivity_),
+    rayId_(ptf.rayId_),
+    lambdaId_(ptf.lambdaId_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::
+updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    const scalarField& Tp =
+        patch().lookupPatchField<volScalarField, scalar>(TName_);
+
+    const fvDOM& dom = db().lookupObject<fvDOM>("radiationProperties");
+
+    const label patchI = patch().index();
+
+    if (dom.nLambda() == 1)
+    {
+        if (rayId_ == -1)
+        {
+            for (label rayI=0; rayI < dom.nRay(); rayI++)
+            {
+                for (label lambdaI=0; lambdaI < dom.nLambda(); lambdaI++)
+                {
+                    const volScalarField& radiationField =
+                        dom.IRayLambda(rayI, lambdaI);
+                    if
+                    (
+                        &(radiationField.internalField())
+                     == &dimensionedInternalField()
+                    )
+                    {
+                        rayId_ = rayI;
+                        lambdaId_ = lambdaI;
+                        break;
+                    }
+                }
+            }
+        }
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "Foam::radiation::"
+            "greyDiffusiveRadiationMixedFvPatchScalarField::"
+            "updateCoeffs"
+        )   << " a grey boundary condition is used with a non-grey"
+            << "absorption model"
+            << exit(FatalError);
+    }
+
+    scalarField& Iw = *this;
+    vectorField n = patch().Sf()/patch().magSf();
+
+    radiativeIntensityRay& ray =
+        const_cast<radiativeIntensityRay&>(dom.IRay(rayId_));
+
+    ray.Qr().boundaryField()[patchI] += Iw*(-n & ray.dAve());
+
+    forAll(Iw, faceI)
+    {
+        scalar Ir = 0.0;
+
+        for (label rayI=0; rayI < dom.nRay(); rayI++)
+        {
+            const vector& d = dom.IRay(rayI).d();
+
+            const scalarField& Iface =
+                dom.IRay(rayI).ILambda(lambdaId_).boundaryField()[patchI];
+
+            if ((-n[faceI] & d) < 0.0) // qin into the wall
+            {
+                const vector& dAve = dom.IRay(rayI).dAve();
+                Ir += Iface[faceI]*mag(n[faceI] & dAve);
+            }
+        }
+
+        const vector& d = dom.IRay(rayId_).d();
+
+        if ((-n[faceI] & d) > 0.) //direction out of the wall
+        {
+            refGrad()[faceI] = 0.0;
+            valueFraction()[faceI] = 1.0;
+            refValue()[faceI] =
+                (
+                    Ir*(1.0 - emissivity_)
+                  + emissivity_*radiation::sigmaSB.value()*pow4(Tp[faceI])
+                )
+               /mathematicalConstant::pi;
+
+        }
+        else //direction into the wall
+        {
+            valueFraction()[faceI] = 0.0;
+            refGrad()[faceI] = 0.0;
+            refValue()[faceI] = 0.0; //not used
+        }
+    }
+
+    mixedFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("T") << TName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("emissivity") << emissivity_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        greyDiffusiveRadiationMixedFvPatchScalarField
+    );
+}
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H

@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::greyDiffusiveRadiationMixedFvPatchScalarField
+
+Description
+    Radiation temperature specified
+
+SourceFiles
+    greyDiffusiveRadiationMixedFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef greyDiffusiveRadiationMixedFvPatchScalarField_H
+#define greyDiffusiveRadiationMixedFvPatchScalarField_H
+
+#include "mixedFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+/*---------------------------------------------------------------------------*\
+        Class greyDiffusiveRadiationMixedFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class greyDiffusiveRadiationMixedFvPatchScalarField
+:
+    public mixedFvPatchScalarField
+{
+    // Private data
+
+        //- Name of temperature field
+        word TName_;
+
+        //- Emissivity
+        scalar emissivity_;
+
+        //- Ray index
+        label rayId_;
+
+        //- Wavelength index
+        label lambdaId_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("greyDiffusiveRadiation");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        greyDiffusiveRadiationMixedFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        greyDiffusiveRadiationMixedFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given a
+        //  greyDiffusiveRadiationMixedFvPatchScalarField onto a new patch
+        greyDiffusiveRadiationMixedFvPatchScalarField
+        (
+            const greyDiffusiveRadiationMixedFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        greyDiffusiveRadiationMixedFvPatchScalarField
+        (
+            const greyDiffusiveRadiationMixedFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new greyDiffusiveRadiationMixedFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        greyDiffusiveRadiationMixedFvPatchScalarField
+        (
+            const greyDiffusiveRadiationMixedFvPatchScalarField&,
+            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 greyDiffusiveRadiationMixedFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Access
+
+            //- Return the temperature field name
+            const word& TName() const
+            {
+                return TName_;
+            }
+
+            //- Return reference to the temperature field name to allow
+            //  adjustment
+            word& TName()
+            {
+                return TName_;
+            }
+
+            //- Return the emissivity
+            const scalar& emissivity() const
+            {
+                return emissivity_;
+            }
+
+            //- Return reference to the emissivity to allow adjustment
+            scalar& emissivity()
+            {
+                return emissivity_;
+            }
+
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // I-O
+
+            //- Write
+            virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C

@@ -0,0 +1,269 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "wideBandDiffusiveRadiationMixedFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+
+#include "fvDOM.H"
+#include "wideBandAbsorptionEmission.H"
+#include "radiationConstants.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::
+wideBandDiffusiveRadiationMixedFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    TName_("undefinedT"),
+    emissivity_(0.0),
+    rayId_(0),
+    lambdaId_(0)
+{
+    refValue() = 0.0;
+    refGrad() = 0.0;
+    valueFraction() = 1.0;
+}
+
+
+Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::
+wideBandDiffusiveRadiationMixedFvPatchScalarField
+(
+    const wideBandDiffusiveRadiationMixedFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchScalarField(ptf, p, iF, mapper),
+    TName_(ptf.TName_),
+    emissivity_(ptf.emissivity_),
+    rayId_(ptf.rayId_),
+    lambdaId_(ptf.lambdaId_)
+{}
+
+
+Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::
+wideBandDiffusiveRadiationMixedFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    TName_(dict.lookup("T")),
+    emissivity_(readScalar(dict.lookup("emissivity"))),
+    rayId_(0),
+    lambdaId_(0)
+{
+    const scalarField& Tp =
+        patch().lookupPatchField<volScalarField, scalar>(TName_);
+
+    refValue() =
+        emissivity_*4.0*radiation::sigmaSB.value()*pow4(Tp)
+       /Foam::mathematicalConstant::pi;
+    refGrad() = 0.0;
+
+    if (dict.found("value"))
+    {
+        fvPatchScalarField::operator=
+        (
+            scalarField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        fvPatchScalarField::operator=(refValue());
+    }
+}
+
+
+Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::
+wideBandDiffusiveRadiationMixedFvPatchScalarField
+(
+    const wideBandDiffusiveRadiationMixedFvPatchScalarField& ptf
+)
+:
+    mixedFvPatchScalarField(ptf),
+    TName_(ptf.TName_),
+    emissivity_(ptf.emissivity_),
+    rayId_(ptf.rayId_),
+    lambdaId_(ptf.lambdaId_)
+{}
+
+
+Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::
+wideBandDiffusiveRadiationMixedFvPatchScalarField
+(
+    const wideBandDiffusiveRadiationMixedFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(ptf, iF),
+    TName_(ptf.TName_),
+    emissivity_(ptf.emissivity_),
+    rayId_(ptf.rayId_),
+    lambdaId_(ptf.lambdaId_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::
+updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    const fvDOM& dom = db().lookupObject<fvDOM>("radiationProperties");
+
+    const label patchI = patch().index();
+
+    if (dom.nLambda() > 1)
+    {
+        if (rayId_ == -1)
+        {
+            for (label rayI=0; rayI < dom.nRay(); rayI++)
+            {
+                for (label lambdaI=0; lambdaI < dom.nLambda(); lambdaI++)
+                {
+                    const volScalarField& radiationField =
+                        dom.IRayLambda(rayI, lambdaI);
+                    if
+                    (
+                        &(radiationField.internalField())
+                      ==&dimensionedInternalField()
+                    )
+                    {
+                        rayId_ = rayI;
+                        lambdaId_ = lambdaI;
+                        break;
+                    }
+                }
+            }
+        }
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "Foam::radiation::"
+            "wideBandDiffusiveRadiationMixedFvPatchScalarField::"
+            "updateCoeffs"
+        )   << " a Non-grey boundary condition is used with a grey"
+            << "absorption model" << nl << exit(FatalError);
+    }
+
+    scalarField& Iw = *this;
+    vectorField n = patch().Sf()/patch().magSf();
+
+    radiativeIntensityRay& ray =
+        const_cast<radiativeIntensityRay&>(dom.IRay(rayId_));
+
+    ray.Qr().boundaryField()[patchI] += Iw*(-n & ray.dAve());
+
+    const scalarField Eb =
+        dom.blackBody().bLambda(lambdaId_).boundaryField()[patchI];
+
+    forAll(Iw, faceI)
+    {
+        scalar Ir = 0.0;
+        for (label rayI=0; rayI < dom.nRay(); rayI++)
+        {
+            const vector& d = dom.IRay(rayI).d();
+
+            const scalarField& Iface =
+                dom.IRay(rayI).ILambda(lambdaId_).boundaryField()[patchI];
+
+            if ((-n[faceI] & d) < 0.0) // qin into the wall
+            {
+                const vector& dAve = dom.IRay(rayI).dAve();
+                Ir = Ir + Iface[faceI]*mag(n[faceI] & dAve);
+            }
+        }
+
+        const vector& d = dom.IRay(rayId_).d();
+
+        if ((-n[faceI] & d) > 0.0) //direction out of the wall
+        {
+            refGrad()[faceI] = 0.0;
+            valueFraction()[faceI] = 1.0;
+            refValue()[faceI] =
+                (
+                    Ir*(1.0 - emissivity_)
+                  + emissivity_*Eb[faceI]
+                )
+               /mathematicalConstant::pi;
+        }
+        else  //direction into the wall
+        {
+            valueFraction()[faceI] = 0.0;
+            refGrad()[faceI] = 0.0;
+            refValue()[faceI] = 0.0; //not used
+        }
+    }
+
+    mixedFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("T") << TName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("emissivity") << emissivity_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        wideBandDiffusiveRadiationMixedFvPatchScalarField
+    );
+}
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.H

@@ -0,0 +1,193 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::wideBandDiffusiveRadiationMixedFvPatchScalarField
+
+Description
+    Radiation temperature specified
+
+SourceFiles
+    wideBandDiffusiveRadiationMixedFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wideBandDiffusiveRadiationMixedFvPatchScalarField_H
+#define wideBandDiffusiveRadiationMixedFvPatchScalarField_H
+
+#include "mixedFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+/*---------------------------------------------------------------------------*\
+      Class wideBandDiffusiveRadiationMixedFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class wideBandDiffusiveRadiationMixedFvPatchScalarField
+:
+    public mixedFvPatchScalarField
+{
+    // Private data
+
+        //- Name of temperature field
+        word TName_;
+
+        //- Emissivity
+        scalar emissivity_;
+
+        //- Ray index
+        label rayId_;
+
+        //- Wavelength index
+        label lambdaId_;
+
+        //- Radiative heat flux on walls.
+        scalarField qr_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("wideBandDiffusiveRadiation");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        wideBandDiffusiveRadiationMixedFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        wideBandDiffusiveRadiationMixedFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given GreyDiffusiveRadiationMixedFvPatchField
+        //  onto a new patch
+        wideBandDiffusiveRadiationMixedFvPatchScalarField
+        (
+            const wideBandDiffusiveRadiationMixedFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        wideBandDiffusiveRadiationMixedFvPatchScalarField
+        (
+            const wideBandDiffusiveRadiationMixedFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new wideBandDiffusiveRadiationMixedFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        wideBandDiffusiveRadiationMixedFvPatchScalarField
+        (
+            const wideBandDiffusiveRadiationMixedFvPatchScalarField&,
+            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 wideBandDiffusiveRadiationMixedFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Access
+
+            //- Return the temperature field name
+            const word& TName() const
+            {
+                return TName_;
+            }
+
+            //- Return reference to the temperature field name to allow
+            //  adjustment
+            word& TName()
+            {
+                return TName_;
+            }
+
+            //- Return the emissivity
+            const scalar& emissivity() const
+            {
+                return emissivity_;
+            }
+
+            //- Return reference to the emissivity to allow adjustment
+            scalar& emissivity()
+            {
+                return emissivity_;
+            }
+
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // I-O
+
+            //- Write
+            virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/P1/P1.C

@@ -50,9 +50,9 @@ namespace Foam
     }
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::radiation::P1::P1(const volScalarField& T)
 :
     radiationModel(typeName, T),
@@ -133,12 +133,8 @@ bool Foam::radiation::P1::read()
 }
 
 
-void Foam::radiation::P1::correct()
+void Foam::radiation::P1::calculate()
 {
-    if (!radiation_)
-    {
-        return;
-    }
     a_ = absorptionEmission_->a();
     e_ = absorptionEmission_->e();
     E_ = absorptionEmission_->E();

src/thermophysicalModels/radiation/radiationModel/P1/P1.H

@@ -59,7 +59,6 @@ class P1
 :
     public radiationModel
 {
-
     // Private data
 
         //- Incident radiation / [W/m2]
@@ -97,18 +96,17 @@ public:
 
 
     // Destructor
-
-        ~P1();
+    virtual ~P1();
 
 
     // Member functions
 
         // Edit
 
-            //- Update radiationSource varible
-            void correct();
+            //- Solve radiation equation(s)
+            void calculate();
 
-            //- Read radiationProperties dictionary
+            //- Read radiation properties dictionary
             bool read();
 
 

src/thermophysicalModels/radiation/radiationModel/fvDOM/absorptionCoeffs/absorptionCoeffs.C

@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "absorptionCoeffs.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::absorptionCoeffs::absorptionCoeffs(Istream& is)
+:
+   Tcommon_(readScalar(is)),
+   Tlow_(readScalar(is)),
+   Thigh_(readScalar(is)),
+   invTemp_(readBool(is))
+{
+    for (label coefLabel=0; absorptionCoeffs::nCoeffs_; coefLabel++)
+    {
+        is >> highACoeffs_[coefLabel];
+    }
+
+    for (label coefLabel=0; absorptionCoeffs::nCoeffs_; coefLabel++)
+    {
+        is >> lowACoeffs_[coefLabel];
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
+
+Foam::radiation::absorptionCoeffs::~absorptionCoeffs()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::radiation::absorptionCoeffs::checkT(const scalar T) const
+{
+    if (T < Tlow_ || T > Thigh_)
+    {
+        FatalErrorIn
+        (
+            "absorptionCoeffs::checkT(const scalar T) const"
+        )   << "attempt to use absCoeff out of temperature range:" << nl
+            << "    " << Tlow_ << " -> " << Thigh_ << ";  T = " << T
+            << nl << abort(FatalError);
+    }
+}
+
+
+const Foam::radiation::absorptionCoeffs::coeffArray&
+Foam::radiation::absorptionCoeffs::coeffs
+(
+    const scalar T
+) const
+{
+    checkT(T);
+
+    if (T < Tcommon_)
+    {
+        return lowACoeffs_;
+    }
+    else
+    {
+        return highACoeffs_;
+    }
+}
+
+
+void Foam::radiation::absorptionCoeffs::initialise(Istream&)
+{
+    absorptionCoeffs(Istream);
+}
+
+
+void Foam::radiation::absorptionCoeffs::initialise(const dictionary& dict)
+{
+    dict.lookup("Tcommon") >> Tcommon_;
+    dict.lookup("Tlow") >> Tlow_;
+    dict.lookup("Tlow") >> Thigh_;
+    dict.lookup("invTemp") >> invTemp_;
+
+    dict.lookup("loTcoeffs") >> lowACoeffs_;
+    dict.lookup("hiTcoeffs") >> highACoeffs_;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/absorptionCoeffs/absorptionCoeffs.H

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::janafThermo
+
+Description
+    Absorption coefficients class used in greyMeanAbsorptionEmission and
+    wideBandAbsorptionEmission
+
+SourceFiles
+    absorptionCoeffs.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef absorptionCoeffs_H
+#define absorptionCoeffs_H
+
+#include "List.H"
+#include "IOstreams.H"
+#include "IOdictionary.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class absorptionCoeffs Declaration
+\*---------------------------------------------------------------------------*/
+
+class absorptionCoeffs
+{
+public:
+
+    // Public data members
+
+        static const int nCoeffs_ = 6;
+        typedef FixedList<scalar, nCoeffs_> coeffArray;
+
+
+private:
+
+    // Private data
+
+        // Temperature limits of applicability for functions
+
+            scalar Tcommon_;
+
+            scalar Tlow_;
+
+            scalar Thigh_;
+
+
+        // Polynomial using inverse temperatures
+        bool invTemp_;
+
+        coeffArray highACoeffs_;
+        coeffArray lowACoeffs_;
+
+
+    // Private member functions
+
+        //- Check given temperature is within the range of the fitted coeffs
+        void checkT(const scalar T) const;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from Istream
+        absorptionCoeffs(Istream&);
+
+        // Null constructor
+        absorptionCoeffs()
+        {}
+
+
+    // Destructor
+    ~absorptionCoeffs();
+
+
+    // Member functions
+
+        //- Return the coefficients corresponding to the given temperature
+        const coeffArray& coeffs(const scalar T) const;
+
+        // Initialise from a dictionary
+        void initialise(const dictionary&);
+
+        // Initialise from an Istream
+        void initialise(Istream&);
+
+
+    // Access Functions
+
+        inline bool invTemp() const;
+
+        inline scalar Tcommon() const;
+
+        inline scalar Tlow() const;
+
+        inline scalar Thigh() const;
+
+        inline const coeffArray& highACoeffs() const;
+
+        inline const coeffArray& lowACoeffs() const;
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace radiation
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "absorptionCoeffsI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/absorptionCoeffs/absorptionCoeffsI.H

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+inline bool Foam::radiation::absorptionCoeffs::invTemp() const
+{
+    return  invTemp_;
+}
+
+
+inline Foam::scalar Foam::radiation::absorptionCoeffs::Tcommon() const
+{
+    return  Tcommon_;
+}
+
+
+inline Foam::scalar Foam::radiation::absorptionCoeffs::Tlow() const
+{
+    return  Tlow_;
+}
+
+
+inline Foam::scalar Foam::radiation::absorptionCoeffs::Thigh() const
+{
+    return  Thigh_;
+}
+
+
+inline const Foam::radiation::absorptionCoeffs::coeffArray&
+Foam::radiation::absorptionCoeffs::highACoeffs() const
+{
+    return  highACoeffs_;
+}
+
+
+inline const Foam::radiation::absorptionCoeffs::coeffArray&
+Foam::radiation::absorptionCoeffs::lowACoeffs() const
+{
+    return  lowACoeffs_;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C

@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "blackBodyEmission.H"
+#include "dimensionedConstants.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::blackBodyEmission::blackBodyEmission
+(
+    const fileName& name,
+    const word& instance,
+    const label nLambda,
+    const volScalarField& T
+)
+:
+    blackBodyEmissiveTable_(name, instance, T.mesh()),
+    C1_("C1",dimensionSet(1, 4, 3, 0, 0, 0, 0), 3.7419e-16),
+    C2_("C2",dimensionSet(0, 1, 0, 1, 0, 0, 0), 14.388e-6),
+    bLambda_(nLambda),
+    T_(T)
+{
+    forAll(bLambda_, lambdaI)
+    {
+        bLambda_.set
+        (
+            lambdaI,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "bLambda_" + Foam::name(lambdaI) ,
+                    T.mesh().time().timeName(),
+                    T.mesh(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                radiation::sigmaSB*pow4(T)
+            )
+        );
+
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::radiation::blackBodyEmission::~blackBodyEmission()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::radiation::blackBodyEmission::fLambdaT
+(
+    const scalar lambdaT
+) const
+{
+    return  blackBodyEmissiveTable_.lookUp(lambdaT*1.0e6)[1];
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::blackBodyEmission::EbDeltaLambdaT
+(
+    const volScalarField& T,
+    const Vector2D<scalar>& band
+) const
+{
+    tmp<volScalarField> Eb
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Eb",
+                T.mesh().time().timeName(),
+                T.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            radiation::sigmaSB*pow4(T)
+        )
+    );
+
+
+    if (band == Vector2D<scalar>::one)
+    {
+        return Eb;
+    }
+    else
+    {
+        forAll(T, i)
+        {
+            scalar T1 = fLambdaT(band[1]*T[i]);
+            scalar T2 = fLambdaT(band[0]*T[i]);
+            dimensionedScalar fLambdaDelta
+            (
+                "fLambdaDelta",
+                dimless,
+                T1 - T2
+            );
+            Eb()[i] = Eb()[i]*fLambdaDelta.value();
+        }
+        return Eb;
+    }
+}
+
+
+void Foam::radiation::blackBodyEmission::correct
+(
+    const label lambdaI,
+    const Vector2D<scalar>& band
+)
+{
+    bLambda_[lambdaI] = EbDeltaLambdaT(T_, band);
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.H

@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::radiation::blackBodyEmission
+
+Description
+    Class black body emission
+
+SourceFiles
+    blackBodyEmission.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef blackModyEmission_H
+#define blackModyEmission_H
+
+#include "volFields.H"
+#include "dimensionedScalar.H"
+#include "mathematicalConstants.H"
+#include "radiationConstants.H"
+#include "interpolationLookUpTable.H"
+#include "Vector2D.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class blackBodyEmission Declaration
+\*---------------------------------------------------------------------------*/
+
+class blackBodyEmission
+{
+    // Private data
+
+        mutable interpolationLookUpTable<scalar> blackBodyEmissiveTable_;
+
+        //- Constant C1
+        const dimensionedScalar C1_;
+
+        //- Constant C2
+        const dimensionedScalar C2_;
+
+        // Ptr List for black body emission energy field for each wavelength
+        PtrList<volScalarField> bLambda_;
+
+        // Reference to the temperature field
+        const volScalarField& T_;
+
+
+    // Private member functions
+
+        scalar fLambdaT(const scalar lambdaT) const;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        blackBodyEmission
+        (
+            const fileName& name,
+            const word& instance,
+            const label nLambda,
+            const volScalarField& T
+        );
+
+
+    // Destructor
+    ~blackBodyEmission();
+
+
+    // Member functions
+
+        // Access
+
+            //- Black body spectrum
+            inline const volScalarField& bLambda(const label lambdaI) const
+            {
+                return bLambda_[lambdaI];
+            }
+
+            //- Spectral emission for the black body at T and lambda
+            inline dimensionedScalar EblambdaT
+            (
+                const dimensionedScalar& T,
+                const scalar lambda
+            ) const
+            {
+                return (C1_/(pow5(lambda)*(exp(C2_/(lambda*T)) - 1.0)));
+            }
+
+            //- Integral energy at T from lambda1 to lambda2
+            tmp<Foam::volScalarField> EbDeltaLambdaT
+            (
+                const volScalarField& T,
+                const Vector2D<scalar>& band
+            ) const;
+
+
+    // Edit
+
+        // Update black body emission
+        void correct(const label lambdaI, const Vector2D<scalar>& band);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace radiation
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmissivePower

@@ -0,0 +1,174 @@
+171
+(
+1000    0.00032
+1100    0.00091
+1200    0.00213
+1300    0.00432
+1400    0.00779
+1500    0.01285
+1600    0.01972
+1700    0.02853
+1800    0.03934
+1900    0.05210
+2000    0.06672
+2100    0.08305
+2200    0.10088
+2300    0.12002
+2400    0.14025
+2500    0.16135
+2600    0.18311
+2700    0.20535
+2800    0.22788
+2900    0.25055
+3000    0.27322
+3100    0.29576
+3200    0.31809
+3300    0.34009
+3400    0.36172
+3500    0.38290
+3600    0.40359
+3700    0.42375
+3800    0.44336
+3900    0.46240
+4000    0.48085
+4100    0.49872
+4200    0.51599
+4300    0.53267
+4400    0.54877
+4500    0.56429
+4600    0.57925
+4700    0.59366
+4800    0.60753
+4900    0.62088
+5000    0.63372
+5100    0.64606
+5200    0.65794
+5300    0.66935
+5400    0.68033
+5500    0.69087
+5600    0.70101
+5700    0.71076
+5800    0.72012
+5900    0.72913
+6000    0.73778
+6100    0.74610
+6200    0.75410
+6300    0.76180
+6400    0.76920
+6500    0.77631
+6600    0.78316
+6700    0.78975
+6800    0.79609
+6900    0.80219
+7000    0.80807
+7100    0.81373
+7200    0.81918
+7300    0.82443
+7400    0.82949
+7500    0.83436
+7600    0.83906
+7700    0.84359
+7800    0.84796
+7900    0.85218
+8000    0.85625
+8100    0.86017
+8200    0.86396
+8300    0.86762
+8400    0.87115
+8500    0.87456
+8600    0.87786
+8700    0.88105
+8800    0.88413
+8900    0.88711
+9000    0.88999
+9100   0.89277
+9200   0.89547
+9300   0.89807
+9400   0.90060
+9500   0.90304
+9600   0.90541
+9700   0.90770
+9800   0.90992
+9900   0.91207
+10000   0.91415
+10200   0.91813
+10400   0.92188
+10600   0.92540
+10800   0.92872
+11000   0.93184
+11200   0.93479
+11400   0.93758
+11600   0.94021
+11800   0.94270
+12000   0.94505
+12200   0.94728
+12400   0.94939
+12600   0.95139
+12800   0.95329
+13000   0.95509
+13200   0.95680
+13400   0.95843
+13600   0.95998
+13800   0.96145
+14000   0.96285
+14200   0.96418
+14400   0.96546
+14600   0.96667
+14800   0.96783
+15000   0.96893
+15200   0.96999
+15400   0.97100
+15600   0.97196
+15800   0.97288
+16000   0.97377
+16200   0.97461
+16400   0.97542
+16600   0.97620
+16800   0.97694
+17000   0.97765
+17200   0.97834
+17400   0.97899
+17600   0.97962
+17800   0.98023
+18000   0.98080
+18200   0.98137
+18400   0.98191
+18600   0.98243
+18900   0.98293
+19000   0.98340
+19200   0.98387
+19400   0.98431
+19600   0.98474
+19800   0.98515
+20000   0.98555
+21000   0.98735
+22000   0.98886
+23000   0.99014
+24000   0.99123
+25000   0.99217
+26000   0.99297
+27000   0.99367
+28000   0.99429
+29000   0.99482
+30000   0.99529
+31000   0.99571
+32000   0.99607
+33000   0.99640
+34000   0.99669
+35000   0.99695
+35000   0.99719
+36000   0.99740
+37000   0.99759
+38000   0.99776
+39000   0.99792
+40000   0.99806
+41000   0.99819
+42000   0.99831
+43000   0.99842
+44000   0.99851
+45000   0.99861
+46000   0.99869
+47000   0.99877
+48000   0.99884
+49000   0.99890
+)

src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C

@@ -0,0 +1,372 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvDOM.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvm.H"
+
+#include "absorptionEmissionModel.H"
+#include "scatterModel.H"
+#include "mathematicalConstants.H"
+#include "radiationConstants.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace radiation
+    {
+        defineTypeNameAndDebug(fvDOM, 0);
+
+        addToRunTimeSelectionTable
+        (
+            radiationModel,
+            fvDOM,
+            dictionary
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
+:
+    radiationModel(typeName, T),
+    G_
+    (
+        IOobject
+        (
+            "G",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("G", dimMass/pow3(dimTime), 0.0)
+    ),
+    Qr_
+    (
+        IOobject
+        (
+            "Qr",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
+    ),
+    a_
+    (
+        IOobject
+        (
+            "a",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("a", dimless/dimLength, 0.0)
+    ),
+    e_
+    (
+        IOobject
+        (
+            "e",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("a", dimless/dimLength, 0.0)
+    ),
+    E_
+    (
+        IOobject
+        (
+            "E",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
+    ),
+    nTheta_(readLabel(coeffs_.lookup("nTheta"))),
+    nPhi_(readLabel(coeffs_.lookup("nPhi"))),
+    nRay_(0),
+    nLambda_(absorptionEmission_->nBands()),
+    aLambda_(nLambda_),
+    blackBody_
+    (
+        fileName("blackBodyEmissivePower"),
+        mesh_.time().constant(),
+        nLambda_,
+        T
+    ),
+    IRay_(0),
+    convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0))
+{
+    if (mesh_.nSolutionD() == 3)    //3D
+    {
+        IRay_.setSize(4*nPhi_*nTheta_);
+        nRay_ = 4*nPhi_*nTheta_;
+        scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
+        scalar deltaTheta = mathematicalConstant::pi/nTheta_;
+        label i = 0;
+        for (label n = 1; n <= nTheta_; n++)
+        {
+            for (label m = 1; m <= 4*nPhi_; m++)
+            {
+                scalar thetai = (2.0*n - 1.0)*deltaTheta/2.0;
+                scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
+                IRay_.set
+                (
+                    i,
+                    new radiativeIntensityRay
+                    (
+                        *this,
+                        mesh_,
+                        phii,
+                        thetai,
+                        deltaPhi,
+                        deltaTheta,
+                        nLambda_,
+                        absorptionEmission_,
+                        blackBody_
+                    )
+                );
+                i++;
+            }
+        }
+    }
+    else
+    {
+        if (mesh_.nSolutionD() == 2)    //2D (X & Y)
+        {
+            scalar thetai = mathematicalConstant::pi/2.0;
+            scalar deltaTheta = mathematicalConstant::pi;
+            IRay_.setSize(4*nPhi_);
+            nRay_ = 4*nPhi_;
+            scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
+            label i = 0;
+            for (label m = 1; m <= 4*nPhi_; m++)
+            {
+                scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
+                IRay_.set
+                (
+                    i,
+                    new radiativeIntensityRay
+                    (
+                        *this,
+                        mesh_,
+                        phii,
+                        thetai,
+                        deltaPhi,
+                        deltaTheta,
+                        nLambda_,
+                        absorptionEmission_,
+                        blackBody_
+                    )
+                );
+                i++;
+            }
+        }
+        else    //1D (X)
+        {
+            scalar thetai = mathematicalConstant::pi/2.0;
+            scalar deltaTheta = mathematicalConstant::pi;
+            IRay_.setSize(2);
+            nRay_ = 2;
+            scalar deltaPhi = mathematicalConstant::pi;
+            label i = 0;
+            for (label m = 1; m <= 2; m++)
+            {
+                scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
+                IRay_.set
+                (
+                    i,
+                    new radiativeIntensityRay
+                    (
+                        *this,
+                        mesh_,
+                        phii,
+                        thetai,
+                        deltaPhi,
+                        deltaTheta,
+                        nLambda_,
+                        absorptionEmission_,
+                        blackBody_
+                    )
+                );
+                i++;
+            }
+
+        }
+    }
+
+
+    // Construct absorption field for each wavelength
+    forAll(aLambda_, lambdaI)
+    {
+        aLambda_.set
+        (
+            lambdaI,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "aLambda_" + Foam::name(lambdaI) ,
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                a_
+            )
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::radiation::fvDOM::~fvDOM()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::radiation::fvDOM::read()
+{
+    if (radiationModel::read())
+    {
+        // nothing to read
+
+//        coeffs_.lookup("nTheta") >> nTheta_;
+//        coeffs_.lookup("nPhi") >> nPhi_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void Foam::radiation::fvDOM::calculate()
+{
+    absorptionEmission_->correct(a_, aLambda_);
+
+    updateBlackBodyEmission();
+
+    scalar maxResidual = 0.0;
+    label radIter = 0;
+    do
+    {
+        radIter++;
+        forAll(IRay_, rayI)
+        {
+            maxResidual = 0.0;
+            scalar maxBandResidual = IRay_[rayI].correct();
+            maxResidual = max(maxBandResidual, maxResidual);
+        }
+
+        Info << "Radiation solver iter: " <<  radIter << endl;
+
+    } while(maxResidual > convergence_);
+
+    updateG();
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Rp",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            4.0*a_*radiation::sigmaSB //absorptionEmission_->a()
+        )
+    );
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::radiation::fvDOM::Ru() const
+{
+
+    const DimensionedField<scalar, volMesh>& G =
+        G_.dimensionedInternalField();
+    const DimensionedField<scalar, volMesh> E =
+        absorptionEmission_->ECont()().dimensionedInternalField();
+    const DimensionedField<scalar, volMesh> a =
+        a_.dimensionedInternalField(); //absorptionEmission_->aCont()()
+
+    return  a*G - 4.0*E;
+}
+
+
+void Foam::radiation::fvDOM::updateBlackBodyEmission()
+{
+    for (label j=0; j < nLambda_; j++)
+    {
+        blackBody_.correct(j, absorptionEmission_->bands(j));
+    }
+}
+
+
+void Foam::radiation::fvDOM::updateG()
+{
+    G_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
+    Qr_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
+
+    forAll(IRay_, rayI)
+    {
+        IRay_[rayI].addIntensity();
+        G_ += IRay_[rayI].I()*IRay_[rayI].omega();
+        Qr_ += IRay_[rayI].Qr();
+    }
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.H

@@ -0,0 +1,227 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::radiation::fvDOM
+
+Description
+
+    Finite Volume Discrete Ordinary Method. Solves the RTE equation for n
+    directions in a participating media, not including scatter.
+
+    Available absorption models:
+        greyMeanAbsoprtionEmission
+        wideBandAbsorptionEmission
+
+    i.e. dictionary
+    fvDOMCoeffs
+    {
+        Nphi    1;          // azimuthal angles in PI/2 on X-Y.(from Y to X)
+        Ntheta  2;          // polar angles in P1 (from Z to X-Y plane)
+        convergence 1e-4;   // convergence criteria for radiation iteration
+    }
+
+    nFlowIterPerRadIter   1; // Number of flow iterations per radiation
+                                iteration
+
+    The total number of solid angles is  4*Nphi*Ntheta.
+
+    In 1D the direction of the rays is X (Nphi and Ntheta are ignored)
+    In 2D the direction of the rays is on X-Y plane (only Nphi is considered)
+    In 3D (Nphi and Ntheta are considered)
+
+SourceFiles
+    fvDOM.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef radiationModelfvDOM_H
+#define radiationModelfvDOM_H
+
+#include "radiativeIntensityRay.H"
+#include "radiationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class fvDOM Declaration
+\*---------------------------------------------------------------------------*/
+
+class fvDOM
+:
+    public radiationModel
+{
+    // Private data
+
+        //- Incident radiation  [W/m2]
+        volScalarField G_;
+
+        //- Total radiative heat flux [W/m2]
+        volScalarField Qr_;
+
+        //- Total absorption coefficient [1/m]
+        volScalarField a_;
+
+        //- Total emission coefficient [1/m]
+        volScalarField e_;
+
+        //- Emission contribution [Kg/m/s^3]
+        volScalarField E_;
+
+        //- Number of solid angles in theta
+        label nTheta_;
+
+        //- Number of solid angles in phi
+        label nPhi_ ;
+
+        //- Total number of rays (1 per direction)
+        label nRay_;
+
+        //- Number of wavelength bands
+        label nLambda_;
+
+        //- Wavelength total absorption coefficient [1/m]
+        PtrList<volScalarField> aLambda_;
+
+        //- Black body
+        blackBodyEmission blackBody_;
+
+        //- List of pointers to radiative intensity rays
+        PtrList<radiativeIntensityRay> IRay_;
+
+        //- Convergence criterion
+        scalar convergence_;
+
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        fvDOM(const fvDOM&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const fvDOM&);
+
+        //- Update Absorption Coefficients
+//        void updateAbsorptionCoeffs(void);
+
+        //- Update nlack body emission
+        void updateBlackBodyEmission();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fvDOM");
+
+
+    // Constructors
+
+        //- Construct from components
+        fvDOM(const volScalarField& T);
+
+
+    //- Destructor
+    virtual ~fvDOM();
+
+
+    // Member functions
+
+        // Edit
+
+            //- Solve radiation equation(s)
+            void calculate();
+
+            //- Read radiation properties dictionary
+            bool read();
+
+            //- Update G and calculate total heat flux on boundary
+            void updateG();
+
+            //- Source term component (for power of T^4)
+            virtual tmp<volScalarField> Rp() const;
+
+            //- Source term component (constant)
+            virtual tmp<DimensionedField<scalar, volMesh> > Ru() const;
+
+
+        // Access
+
+            //- Ray intensity for rayI
+            inline const radiativeIntensityRay& IRay(const label rayI) const;
+
+            //- Ray intensity for rayI and lambda bandwidth
+            inline const volScalarField& IRayLambda
+            (
+                const label rayI,
+                const label lambdaI
+            ) const;
+
+            //- Number of angles in theta
+            inline label nTheta() const;
+
+            //- Number of angles in phi
+            inline label nPhi() const;
+
+            //- Number of rays
+            inline label nRay() const;
+
+            //- Number of wavelengths
+            inline label nLambda() const;
+
+            // Const access to total absorption coefficient
+            inline const volScalarField& a() const;
+
+            // Const access to wavelength total absorption coefficient
+            inline const volScalarField& aLambda(const label lambdaI) const;
+
+            // Const access to incident radiation field
+            inline const volScalarField& G() const;
+
+            //  Const access to total radiative heat flux field
+            inline const volScalarField& Qr() const;
+
+            //  Const access to black body
+            inline const blackBodyEmission& blackBody() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "fvDOMI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace radiation
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOMI.H

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+inline const Foam::radiation::radiativeIntensityRay&
+Foam::radiation::fvDOM::IRay(const label rayI) const
+{
+    return  IRay_[rayI];
+}
+
+
+inline const Foam::volScalarField&
+Foam::radiation::fvDOM::IRayLambda
+(
+    const label rayI,
+    const label lambdaI
+) const
+{
+    return IRay_[rayI].ILambda(lambdaI);
+}
+
+
+inline Foam::label Foam::radiation::fvDOM::nTheta() const
+{
+    return nTheta_;
+}
+
+
+inline Foam::label Foam::radiation::fvDOM::nPhi() const
+{
+    return nPhi_;
+}
+
+
+inline Foam::label Foam::radiation::fvDOM::nRay() const
+{
+    return nRay_;
+}
+
+
+inline Foam::label Foam::radiation::fvDOM::nLambda() const
+{
+    return nLambda_;
+}
+
+
+inline const Foam::volScalarField& Foam::radiation::fvDOM::a() const
+{
+    return a_;
+}
+
+
+inline const Foam::volScalarField& Foam::radiation::fvDOM::aLambda
+(
+    const label lambdaI
+) const
+{
+    return aLambda_[lambdaI];
+}
+
+
+inline const Foam::volScalarField& Foam::radiation::fvDOM::G() const
+{
+    return G_;
+}
+
+
+inline const Foam::volScalarField& Foam::radiation::fvDOM::Qr() const
+{
+    return Qr_;
+}
+
+
+inline const Foam::radiation::blackBodyEmission&
+Foam::radiation::fvDOM::blackBody() const
+{
+    return blackBody_;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.C

@@ -0,0 +1,521 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "IFstream.H"
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+template <class Type>
+Foam::label Foam::interpolationLookUpTable <Type>::index
+(
+    const List<scalar>& indices,
+    const bool lastDim
+) const
+{
+    label totalindex = 0;
+
+    for (int i = 0; i < dim_.size() - 1; i++)
+    {
+        label dim = 1;
+        for (int j = i + 1; j < dim_.size(); j++)
+        {
+            dim *=(dim_[j]+1);
+        }
+
+        totalindex +=
+            dim
+           *min
+            (
+                max(label((indices[i] - min_[i])/delta_[i]), 0),
+                dim_[i]
+            );
+    }
+
+    if (lastDim)
+    {
+        label iLastdim = dim_.size() - 1;
+        totalindex += Foam::min
+        (
+            max
+            (
+                label((indices[iLastdim] - min_[iLastdim])/delta_[iLastdim]),
+                0
+            ),
+            dim_[iLastdim]
+        );
+    }
+
+    return totalindex;
+}
+
+
+template <class Type>
+Foam::label Foam::interpolationLookUpTable <Type>::index
+(
+    const scalar indice
+) const
+{
+    label i = 0;
+    label totalIndex =
+        Foam::min
+        (
+            Foam::max
+            (
+                label((indice - min_[i])/delta_[i]),
+                0
+            ),
+            dim_[i]
+        );
+
+    return totalIndex;
+}
+
+
+template<class Type>
+bool Foam::interpolationLookUpTable<Type>::checkRange
+(
+    const scalar lookUpValue,
+    const label interfield
+) const
+{
+    if (lookUpValue >= min_[interfield] &&  lookUpValue <= max_[interfield])
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class Type>
+Foam::scalar Foam::interpolationLookUpTable<Type>::interpolate
+(
+    const label lo,
+    const label hi,
+    const scalar lookUpValue,
+    const label ofield,
+    const label interfield
+) const
+{
+        if
+        (
+            List<scalarField>::operator[](interfield).operator[](hi)
+         != List<scalarField>::operator[](interfield).operator[](lo)
+        )
+        {
+            scalar output
+            (
+                List<scalarField>::operator[](ofield).operator[](lo)
+              + (
+                    List<scalarField>::operator[](ofield).operator[](hi)
+                  - List<scalarField>::operator[](ofield).operator[](lo)
+                )
+               *(
+                    lookUpValue
+                  - List<scalarField>::operator[](interfield).operator[](lo)
+                )
+               /(
+                    List<scalarField>::operator[](interfield).operator[](hi)
+                  - List<scalarField>::operator[](interfield).operator[](lo)
+                )
+            );
+            return output;
+        }
+        else
+        {
+            return List<scalarField>::operator[](ofield).operator[](lo);
+        }
+}
+
+
+template<class Type>
+void Foam::interpolationLookUpTable<Type>::dimensionTable()
+{
+    min_.setSize(entries_.size());
+    dim_.setSize(entries_.size());
+    delta_.setSize(entries_.size());
+    max_.setSize(entries_.size());
+    entryIndices_.setSize(entries_.size());
+    outputIndices_.setSize(output_.size());
+    label index = 0;
+    label tableDim = 1;
+
+    forAll(entries_,i)
+    {
+        dim_[i] = readLabel(entries_[i].lookup("N"));
+        max_[i] = readScalar(entries_[i].lookup("max"));
+        min_[i] = readScalar(entries_[i].lookup("min"));
+        delta_[i] = (max_[i] - min_[i])/dim_[i];
+        tableDim *= (dim_[i] + 1);
+        fieldIndices_.insert(entries_[i].lookup("name"),index);
+        entryIndices_[i] = index;
+        index++;
+    }
+
+    forAll(output_,i)
+    {
+        fieldIndices_.insert(output_[i].lookup("name"),index);
+        outputIndices_[i] = index;
+        index++;
+    }
+
+    List<scalarField>& internal = *this;
+
+    internal.setSize(entries_.size() + output_.size());
+
+    interpOutput_.setSize(entries_.size() + output_.size());
+
+    forAll(internal, i)
+    {
+        internal[i].setSize(tableDim);
+    }
+}
+
+
+template<class Type>
+void Foam::interpolationLookUpTable<Type>::readTable
+(
+    const word& instance,
+    const fvMesh& mesh
+)
+{
+    IOdictionary control
+    (
+        IOobject
+        (
+            fileName_,
+            instance,
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    control.lookup("fields") >> entries_;
+    control.lookup("output") >> output_;
+    control.lookup("values") >> *this;
+
+    dimensionTable();
+
+    check();
+
+    if (this->size() == 0)
+    {
+        FatalErrorIn
+        (
+            "Foam::interpolationLookUpTable<Type>::readTable()"
+        )   << "table is empty" << nl
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::interpolationLookUpTable<Type>::interpolationLookUpTable()
+:
+    List<scalarField >(),
+    fileName_("fileNameIsUndefined")
+{}
+
+
+template<class Type>
+Foam::interpolationLookUpTable<Type>::interpolationLookUpTable
+(
+    const fileName& fn, const word& instance, const fvMesh& mesh
+)
+:
+    List<scalarField >(),
+    fileName_(fn),
+    dim_(0),
+    min_(0),
+    delta_(0.0),
+    max_(0.0),
+    entries_(0),
+    output_(0),
+    entryIndices_(0),
+    outputIndices_(0),
+    interpOutput_(0)
+{
+    readTable(instance, mesh);
+}
+
+
+template<class Type>
+Foam::interpolationLookUpTable<Type>::interpolationLookUpTable
+(
+     const interpolationLookUpTable& interpTable
+)
+:
+    List<scalarField >(interpTable),
+    fileName_(interpTable.fileName_),
+    entryIndices_(interpTable.entryIndices_),
+    outputIndices_(interpTable.outputIndices_),
+    dim_(interpTable.dim_),
+    min_(interpTable.min_),
+    delta_(interpTable.delta_),
+    max_(interpTable.max_),
+    entries_(0),
+    output_(0),
+    interpOutput_(interpTable.interpOutput_)
+{}
+
+
+template<class Type>
+Foam::interpolationLookUpTable<Type>::interpolationLookUpTable
+(
+    const dictionary& dict
+)
+:
+    List<scalarField >(),
+    fileName_(fileName(dict.lookup("fileName")).expand()),
+    dim_(0),
+    min_(0.0),
+    delta_(0.0),
+    max_(0.0),
+    entries_(dict.lookup("fields")),
+    output_(dict.lookup("output")),
+    entryIndices_(0),
+    outputIndices_(0),
+    fieldIndices_(0),
+    interpOutput_(0)
+{
+    dimensionTable();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::interpolationLookUpTable<Type>::check() const
+{
+//  check order in the first dimension.
+
+    scalar prevValue = List<scalarField >::operator[](0).operator[](0);
+    label dim = 1 ;
+    for (int j = 1; j < dim_.size(); j++)
+    {
+        dim *=(dim_[j]+1);
+    }
+
+    for (label i = 1; i < dim_[0]; i++)
+    {
+        label index = i*dim;
+        const scalar currValue =
+            List<scalarField >::operator[](0).operator[](index);
+
+        // avoid duplicate values (divide-by-zero error)
+        if (currValue <= prevValue)
+        {
+            FatalErrorIn
+            (
+                "Foam::interpolationLookUpTable<Type>::checkOrder() const"
+            )   << "out-of-order value: "
+                << currValue << " at index " << index << nl
+                << exit(FatalError);
+        }
+        prevValue = currValue;
+    }
+}
+
+
+template<class Type>
+void Foam::interpolationLookUpTable<Type>::write
+(
+    Ostream& os,
+    const fileName& fn,
+    const word& instance,
+    const fvMesh& mesh
+) const
+{
+    IOdictionary control
+    (
+        IOobject
+        (
+            fn,
+            instance,
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    control.writeHeader(os);
+
+    os.writeKeyword("fields");
+    os << entries_ << token::END_STATEMENT << nl;
+
+    os.writeKeyword("output");
+    os << output_ << token::END_STATEMENT << nl;
+
+    if (this->size() == 0)
+    {
+        FatalErrorIn
+        (
+            "Foam::interpolationTable<Type>::write()"
+        )   << "table is empty" << nl
+            << exit(FatalError);
+    }
+    os.writeKeyword("values");
+    os << *this << token::END_STATEMENT << nl;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::scalarField&
+Foam::interpolationLookUpTable<Type>::operator[](const label i)
+{
+    label ii = i;
+    label n  = this->size();
+
+    if (n <= 1)
+    {
+        FatalErrorIn
+        (
+            "Foam::interpolationLookUpTable<Type>::operator[]"
+            "(const label) const"
+        )   << "table has (" << n << ") columns" << nl
+            << exit(FatalError);
+    }
+    else if (ii < 0)
+    {
+        FatalErrorIn
+        (
+            "Foam::interpolationLookUpTable<Type>::operator[]"
+            "(const label) const"
+        )   << "index (" << ii << ") underflow" << nl
+            << exit(FatalError);
+    }
+    else if (ii > n)
+    {
+        FatalErrorIn
+        (
+            "Foam::interpolationLookUpTable<Type>::operator[]"
+            "(const label) const"
+        )   << "index (" << ii << ") overflow" << nl
+            << exit(FatalError);
+    }
+
+    return List<scalarField >::operator[](ii);
+}
+
+
+template<class Type>
+const Foam::scalarField&
+Foam::interpolationLookUpTable<Type>::operator[](const label i) const
+{
+    label ii = i;
+    label n  = this->size();
+
+    if (n <= 1)
+    {
+        FatalErrorIn
+        (
+            "Foam::interpolationLookUpTable<Type>::operator[]"
+            "(const label) const"
+        )   << "table has (" << n << ") columns" << nl
+            << exit(FatalError);
+    }
+    else if (ii < 0)
+    {
+        FatalErrorIn
+        (
+            "Foam::interpolationLookUpTable<Type>::operator[]"
+            "(const label) const"
+        )   << "index (" << ii << ") underflow" << nl
+            << exit(FatalError);
+    }
+
+    else if (ii > n)
+    {
+        FatalErrorIn
+        (
+            "Foam::interpolationLookUpTable<Type>::operator[]"
+            "(const label) const"
+        )   << "index (" << ii << ") overflow" << nl
+            << exit(FatalError);
+    }
+
+    return List<scalarField >::operator[](ii);
+}
+
+
+template<class Type>
+bool Foam::interpolationLookUpTable<Type>::found(const word& fieldName) const
+{
+    return fieldIndices_.found(fieldName);
+}
+
+
+template<class Type>
+const Foam::scalarList&
+Foam::interpolationLookUpTable<Type>::lookUp(const scalar retvals)
+{
+    const label lo = index(retvals);
+    findHi(lo, retvals);
+    return interpOutput_;
+}
+
+
+template<class Type>
+void Foam::interpolationLookUpTable<Type>::findHi
+(
+    const label lo,
+    const scalar retvals
+)
+{
+    forAll(outputIndices_,j)
+    {
+        scalar tmp = 0;
+        label ofield = outputIndices_[j];
+        scalar baseValue = List<scalarField>::operator[](ofield).operator[](lo);
+
+        forAll(entryIndices_,i)
+        {
+            if (checkRange(retvals,entryIndices_[i]))
+            {
+                label dim = 1;
+
+                label hi = Foam::min(lo + dim, (*this)[0].size() - 1);
+
+                tmp += interpolate(lo, hi, retvals, ofield, entryIndices_[i])
+                     - baseValue;
+            }
+            interpOutput_[entryIndices_[i]] = retvals;
+        }
+
+        tmp += baseValue;
+        interpOutput_[outputIndices_[j]] = tmp;
+    }
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.H

@@ -0,0 +1,234 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::interpolationLookUpTable
+
+Description
+    A list of lists. Interpolates based on the first dimension.
+    The values must be positive and monotonically increasing in each dimension
+
+Note
+    - Accessing an empty list results in an error.
+    - Accessing a list with a single element always returns the same value.
+
+SourceFiles
+    interpolationLookUpTable.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interpolationLookUpTable_H
+#define interpolationLookUpTable_H
+
+#include "List.H"
+#include "ListOps.H"
+#include "scalarField.H"
+#include "HashTable.H"
+#include "IOdictionary.H"
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class interpolationLookUpTable Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class interpolationLookUpTable
+:
+    public List<scalarField>
+{
+private:
+
+    // Private data
+
+        //- File name
+        fileName fileName_;
+
+        //- Table dimensions
+        List<label>  dim_;
+
+        //- Min on each dimension
+        List<scalar>  min_;
+
+        //- Deltas on each dimension
+        List<scalar>  delta_;
+
+        //- Maximum on each dimension
+        List<scalar>  max_;
+
+        //- Dictionary entries
+        List<dictionary> entries_;
+
+        //- Output dictionaries
+        List<dictionary> output_;
+
+        //- Input indices from the look up table
+        List<label> entryIndices_;
+
+        //- Output Indeces from the Look Up Table
+        List<label> outputIndices_;
+
+        //- Field names and indices
+        HashTable<label> fieldIndices_;
+
+        //- Output list containing input and interpolation values of outputs
+        List<scalar> interpOutput_;
+
+
+    // Private Member Functions
+
+        //- Read the table of data from file
+        void readTable(const word& instance, const fvMesh& mesh);
+
+        //- Dimension table from dictionaries input and output
+        void dimensionTable();
+
+        //- Find table index by scalarList
+        label index(const List<scalar>&, const bool lastDim=true) const;
+
+        //- Find table index by scalar
+        label index(const scalar) const;
+
+        //- Check range of lookup value
+        bool checkRange(const scalar, const label) const;
+
+        //- Interpolate function return an scalar
+        scalar interpolate
+        (
+            const label lo,
+            const label hi,
+            const scalar lookUpValue,
+            const label ofield,
+            const label interfield
+        ) const;
+
+        // Check list is monotonically increasing
+        void check() const;
+
+        // find hi index, interpolate and populate interpOutput_
+        void findHi(const label lo, const scalar retvals);
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        interpolationLookUpTable();
+
+        //- Construct given the name of the file containing the table of data
+        interpolationLookUpTable
+        (
+            const fileName& fn,
+            const word& instance,
+            const fvMesh& mesh
+        );
+
+         //- Construct from dictionary
+        interpolationLookUpTable(const dictionary& dict);
+
+        //- Construct copy
+        interpolationLookUpTable(const interpolationLookUpTable& interpTable);
+
+
+    // Member Functions
+
+        //- Return true if the filed exists in the table
+        bool found(const word& fieldName) const;
+
+        //- Return the output list given a single input scalar
+        const List<scalar>& lookUp(const scalar);
+
+        //- Write Look Up Table to filename.
+        void write
+        (
+            Ostream& os,
+            const fileName& fn,
+            const word& instance,
+            const fvMesh& mesh
+        ) const;
+
+
+    // Access
+
+        //- Return the index of a field by name
+        inline label findFieldIndex(const word& fieldName) const;
+
+        //- Return const access to the output dictionaries
+        inline const List<dictionary>& output() const;
+
+        //- Return const access tp the dictionary entries
+        inline const List<dictionary>& entries() const;
+
+        //- Return const access to the list of min dimensions
+        inline const List<scalar>& min() const;
+
+        //- Return const access to the list of dimensions
+        inline const List<label>& dim() const;
+
+        //- Return const access to the deltas in each dimension
+        inline const List<scalar>& delta() const;
+
+        //- Return const access to the list of max dimensions
+        inline const List<scalar>& max() const;
+
+        //- Return const access to the table name
+        inline word tableName() const;
+
+
+     // Member Operators
+
+        //- Return an element of constant List<scalar, Type>
+        const scalarField& operator[](const label) const;
+
+        //- Return an element of List<scalar, Type>
+        scalarField& operator[](const label);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "interpolationLookUpTableI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#ifdef NoRepository
+#   include "interpolationLookUpTable.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTableI.H

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+inline Foam::label
+Foam::interpolationLookUpTable<Type>::findFieldIndex
+(
+    const word& fieldName
+) const
+{
+    return fieldIndices_[fieldName];
+}
+
+
+template<class Type>
+inline const Foam::List<Foam::dictionary>&
+Foam::interpolationLookUpTable<Type>::output() const
+{
+    return output_;
+}
+
+
+template<class Type>
+inline const Foam::List<Foam::dictionary>&
+Foam::interpolationLookUpTable<Type>::entries() const
+{
+    return entries_;
+}
+
+
+template<class Type>
+inline const Foam::List<Foam::scalar>&
+Foam::interpolationLookUpTable<Type>::min() const
+{
+    return min_;
+}
+
+
+template<class Type>
+inline const Foam::List<Foam::label>&
+Foam::interpolationLookUpTable<Type>::dim() const
+{
+    return dim_;
+}
+
+
+template<class Type>
+inline const Foam::List<Foam::scalar>&
+Foam::interpolationLookUpTable<Type>::delta() const
+{
+    return delta_;
+}
+
+
+template<class Type>
+inline const Foam::List<Foam::scalar>&
+Foam::interpolationLookUpTable<Type>::max() const
+{
+    return max_;
+}
+
+
+template<class Type>
+inline Foam::word Foam::interpolationLookUpTable<Type>::tableName() const
+{
+    return fileName_.name();
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C

@@ -0,0 +1,222 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "radiativeIntensityRay.H"
+#include "fvm.H"
+#include "fvDOM.H"
+
+#include "absorptionEmissionModel.H"
+#include "scatterModel.H"
+#include "mathematicalConstants.H"
+#include "radiationConstants.H"
+#include "radiationModel.H"
+#include "Vector2D.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+Foam::label Foam::radiation::radiativeIntensityRay::rayId = 0;
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
+(
+    const fvDOM& dom,
+    const fvMesh& mesh,
+    const scalar phi,
+    const scalar theta,
+    const scalar deltaPhi,
+    const scalar deltaTheta,
+    const label nLambda,
+    const absorptionEmissionModel& absorptionEmission,
+    const blackBodyEmission& blackBody
+)
+:
+    dom_(dom),
+    mesh_(mesh),
+    absorptionEmission_(absorptionEmission),
+    blackBody_(blackBody),
+    I_
+    (
+        IOobject
+        (
+            "I" + name(rayId),
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("I", dimMass/pow3(dimTime), 0.0)
+    ),
+    Qr_
+    (
+        IOobject
+        (
+            "Qr" + name(rayId),
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
+    ),
+    d_(vector::zero),
+    dAve_(vector::zero),
+    theta_(theta),
+    phi_(phi),
+    omega_(0.0),
+    nLambda_(nLambda),
+    ILambda_(nLambda)
+{
+    scalar sinTheta = Foam::sin(theta);
+    scalar cosTheta = Foam::cos(theta);
+    scalar sinPhi = Foam::sin(phi);
+    scalar cosPhi = Foam::cos(phi);
+
+    omega_ = 2.0*sinTheta*Foam::sin(deltaTheta/2.0)*deltaPhi;
+    d_ = vector(sinTheta*sinPhi, sinTheta*cosPhi, cosTheta);
+    dAve_ = vector
+    (
+        sinPhi
+       *Foam::sin(0.5*deltaPhi)
+       *(deltaTheta - Foam::cos(2.0*theta)
+       *Foam::sin(deltaTheta)),
+        cosPhi
+       *Foam::sin(0.5*deltaPhi)
+       *(deltaTheta - Foam::cos(2.0*theta)
+       *Foam::sin(deltaTheta)),
+        0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta)
+    );
+
+    forAll(ILambda_, i)
+    {
+        IOobject IHeader
+        (
+            "ILambda_" + name(rayId) + "_" + name(i),
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        );
+
+        // check if field exists and can be read
+        if (IHeader.headerOk())
+        {
+            ILambda_.set
+            (
+                i,
+                new volScalarField(IHeader, mesh_)
+            );
+        }
+        else
+        {
+            volScalarField IDefault
+            (
+                IOobject
+                (
+                    "IDefault",
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh_
+            );
+
+            ILambda_.set
+            (
+                i,
+                new volScalarField(IHeader, IDefault)
+            );
+        }
+    }
+    rayId++;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::radiation::radiativeIntensityRay::~radiativeIntensityRay()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
+{
+    // reset boundary heat flux to zero
+    Qr_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
+
+    scalar maxResidual = -GREAT;
+
+    forAll(ILambda_, lambdaI)
+    {
+        const volScalarField& k = dom_.aLambda(lambdaI);
+
+        surfaceScalarField Ji = dAve_ & mesh_.Sf();
+
+        fvScalarMatrix IiEq
+        (
+            fvm::div(Ji, ILambda_[lambdaI], " div(Ji,Ii_h)")
+          + fvm::Sp(k*omega_, ILambda_[lambdaI])
+         ==
+            1.0/Foam::mathematicalConstant::pi
+           *(
+                k*omega_*blackBody_.bLambda(lambdaI)
+              + absorptionEmission_.ECont(lambdaI)
+            )
+        );
+
+        IiEq.relax();
+
+        scalar eqnResidual = solve
+        (
+            IiEq,
+            mesh_.solver("Ii")
+        ).initialResidual();
+
+        maxResidual = max(eqnResidual, maxResidual);
+
+    }
+
+    return maxResidual;
+}
+
+
+void Foam::radiation::radiativeIntensityRay::addIntensity()
+{
+    I_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
+
+    forAll(ILambda_, lambdaI)
+    {
+        I_ += absorptionEmission_.addIntensity(lambdaI, ILambda_[lambdaI]);
+    }
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H

@@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::radiation::radiativeIntensityRay
+
+Description
+    Radiation intensity for a ray in a given direction
+
+SourceFiles
+    radiativeIntensityRay.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef radiativeIntensityRay_H
+#define radiativeIntensityRay_H
+
+#include "absorptionEmissionModel.H"
+#include "blackBodyEmission.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+// Forward declaration of classes
+class fvDOM;
+
+/*---------------------------------------------------------------------------*\
+                    Class radiativeIntensityRay Declaration
+\*---------------------------------------------------------------------------*/
+
+class radiativeIntensityRay
+{
+    // Private data
+
+        //- Refence to the owner fvDOM object
+        const fvDOM& dom_;
+
+        //- Reference to the mesh
+        const fvMesh& mesh_;
+
+        //- Absorption/emission model
+        const absorptionEmissionModel& absorptionEmission_;
+
+        //- Black body
+        const blackBodyEmission& blackBody_;
+
+        //- Total radiative intensity / [W/m2]
+        volScalarField I_;
+
+        //- Total radiative heat flux on boundary
+        volScalarField Qr_;
+
+        //- Direction
+        vector d_;
+
+        //- Average direction vector inside the solid angle
+        vector dAve_;
+
+        //- Theta angle
+        scalar theta_;
+
+        //- Phi angle
+        scalar phi_;
+
+        //- Solid angle
+        scalar omega_;
+
+        //- Number of wavelengths/bands
+        label nLambda_;
+
+        //- List of pointers to radiative intensity fields for given wavelengths
+        PtrList<volScalarField> ILambda_;
+
+        //- Global ray id - incremented in constructor
+        static label rayId;
+
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        radiativeIntensityRay(const radiativeIntensityRay&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const radiativeIntensityRay&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct form components
+        radiativeIntensityRay
+        (
+            const fvDOM& dom,
+            const fvMesh& mesh,
+            const scalar phi,
+            const scalar theta,
+            const scalar deltaPhi,
+            const scalar deltaTheta,
+            const label lambda,
+            const absorptionEmissionModel& absEmmModel_,
+            const blackBodyEmission& blackBody
+        );
+
+
+    // Destructor
+    ~radiativeIntensityRay();
+
+
+    // Member functions
+
+        // Edit
+
+            //- Update radiative intensity on i direction
+            scalar correct();
+
+            //- Initialise the ray in i direction
+            void init
+            (
+                const scalar phi,
+                const scalar theta,
+                const scalar deltaPhi,
+                const scalar deltaTheta,
+                const scalar lambda
+            );
+
+            //- Add radiative intensities from all the bands
+            void addIntensity();
+
+
+        // Access
+
+            //- Return intensity
+            inline const volScalarField& I() const;
+
+            //- Return const access to the boundary heat flux
+            inline const volScalarField& Qr() const;
+
+            //- Return non-const access to the boundary heat flux
+            inline volScalarField& Qr();
+
+            //- Return direction
+            inline const vector& d() const;
+
+            //- Return the average vector inside the solid angle
+            inline const vector& dAve() const;
+
+            //- Return the number of bands
+            inline scalar nLambda() const;
+
+            //- Return the phi angle
+            inline scalar phi() const;
+
+            //- Return the theta angle
+            inline scalar theta() const;
+
+            //- Return the solid angle
+            inline scalar omega() const;
+
+            //- Return the radiative intensity for a given wavelength
+            inline const volScalarField& ILambda(const label lambdaI) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace radiation
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "radiativeIntensityRayI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRayI.H

@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+inline const Foam::volScalarField& Foam::radiation::radiativeIntensityRay::I() const
+{
+    return I_;
+}
+
+
+inline const Foam::volScalarField& Foam::radiation::radiativeIntensityRay::Qr() const
+{
+    return Qr_;
+}
+
+
+inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::Qr()
+{
+    return Qr_;
+}
+
+
+inline const Foam::vector& Foam::radiation::radiativeIntensityRay::d() const
+{
+    return d_;
+}
+
+
+inline const Foam::vector& Foam::radiation::radiativeIntensityRay::dAve() const
+{
+    return dAve_;
+}
+
+
+inline Foam::scalar Foam::radiation::radiativeIntensityRay::nLambda() const
+{
+    return nLambda_;
+}
+
+
+inline Foam::scalar Foam::radiation::radiativeIntensityRay::phi() const
+{
+    return phi_;
+}
+
+
+inline Foam::scalar Foam::radiation::radiativeIntensityRay::theta() const
+{
+    return theta_;
+}
+
+
+inline Foam::scalar Foam::radiation::radiativeIntensityRay::omega() const
+{
+    return omega_;
+}
+
+
+inline const Foam::volScalarField& Foam::radiation::radiativeIntensityRay::ILambda
+(
+    const label lambdaI
+) const
+{
+    return ILambda_[lambdaI];
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/noRadiation/noRadiation.C

@@ -50,10 +50,9 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::radiation::noRadiation::noRadiation(const volScalarField& T)
 :
-    radiationModel(typeName, T)
+    radiationModel(T)
 {}
 
 
@@ -71,7 +70,7 @@ bool Foam::radiation::noRadiation::read()
 }
 
 
-void Foam::radiation::noRadiation::correct()
+void Foam::radiation::noRadiation::calculate()
 {
     // Do nothing
 }

src/thermophysicalModels/radiation/radiationModel/noRadiation/noRadiation.H

@@ -54,7 +54,6 @@ class noRadiation
 :
     public radiationModel
 {
-
     // Private member functions
 
         //- Disallow default bitwise copy construct
@@ -76,17 +75,16 @@ public:
         noRadiation(const volScalarField& T);
 
 
-    // Destructor
-
-        ~noRadiation();
+    //- Destructor
+    virtual ~noRadiation();
 
 
     // Member functions
 
         // Edit
 
-            //- Update radiation source
-            void correct();
+            //- Solve radiation equation(s)
+            void calculate();
 
             //- Read radiationProperties dictionary
             bool read();

src/thermophysicalModels/radiation/radiationModel/radiationModel/newRadiationModel.C

@@ -52,8 +52,8 @@ autoPtr<radiationModel> radiationModel::New
             IOobject
             (
                 "radiationProperties",
-                T.time().constant(),
-                T.db(),
+                T.mesh().time().constant(),
+                T.mesh().objectRegistry::db(),
                 IOobject::MUST_READ,
                 IOobject::NO_WRITE
             )
@@ -75,7 +75,7 @@ autoPtr<radiationModel> radiationModel::New
             "radiationModel::New(const volScalarField&)"
         )   << "Unknown radiationModel type " << radiationModelTypeName
             << nl << nl
-            << "Valid radiationModel types are :" << nl
+            << "Valid radiationModel types are:" << nl
             << dictionaryConstructorTablePtr_->toc()
             << exit(FatalError);
     }
@@ -86,7 +86,7 @@ autoPtr<radiationModel> radiationModel::New
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace radiation
+} // End radiation
 } // End namespace Foam
 
 // ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.C

@@ -43,6 +43,30 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::radiation::radiationModel::radiationModel(const volScalarField& T)
+:
+    IOdictionary
+    (
+        IOobject
+        (
+            "radiationProperties",
+            T.time().constant(),
+            T.mesh().objectRegistry::db(),
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    mesh_(T.mesh()),
+    time_(T.time()),
+    T_(T),
+    radiation_(false),
+    coeffs_(dictionary::null),
+    solverFreq_(0),
+    absorptionEmission_(NULL),
+    scatter_(NULL)
+{}
+
+
 Foam::radiation::radiationModel::radiationModel
 (
     const word& type,
@@ -55,18 +79,22 @@ Foam::radiation::radiationModel::radiationModel
         (
             "radiationProperties",
             T.time().constant(),
-            T.db(),
+            T.mesh().objectRegistry::db(),
             IOobject::MUST_READ,
             IOobject::NO_WRITE
         )
     ),
-    T_(T),
     mesh_(T.mesh()),
+    time_(T.time()),
+    T_(T),
     radiation_(lookup("radiation")),
-    radiationModelCoeffs_(subDict(type + "Coeffs")),
+    coeffs_(subDict(type + "Coeffs")),
+    solverFreq_(readLabel(lookup("solverFreq"))),
     absorptionEmission_(absorptionEmissionModel::New(*this, mesh_)),
     scatter_(scatterModel::New(*this, mesh_))
-{}
+{
+    solverFreq_ = max(1, solverFreq_);
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor    * * * * * * * * * * * * * * //
@@ -82,7 +110,7 @@ bool Foam::radiation::radiationModel::read()
     if (regIOobject::read())
     {
         lookup("radiation") >> radiation_;
-        radiationModelCoeffs_ = subDict(type() + "Coeffs");
+        coeffs_ = subDict(type() + "Coeffs");
 
         return true;
     }
@@ -93,6 +121,20 @@ bool Foam::radiation::radiationModel::read()
 }
 
 
+void Foam::radiation::radiationModel::correct()
+{
+    if (!radiation_)
+    {
+        return;
+    }
+
+    if ((time_.timeIndex() == 0) || (time_.timeIndex() % solverFreq_ == 0))
+    {
+        calculate();
+    }
+}
+
+
 Foam::tmp<Foam::fvScalarMatrix> Foam::radiation::radiationModel::Sh
 (
     basicThermo& thermo

src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.H

@@ -28,7 +28,6 @@ Namespace
 Description
     Namespace for radiation modelling
 
-
 Class
     Foam::radiation::radiationModel
 
@@ -49,6 +48,7 @@ SourceFiles
 #include "volFields.H"
 #include "basicThermo.H"
 #include "fvMatrices.H"
+#include "blackBodyEmission.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -57,33 +57,40 @@ namespace Foam
 namespace radiation
 {
 
+// Forward declaration of classes
 class absorptionEmissionModel;
 class scatterModel;
 
 /*---------------------------------------------------------------------------*\
-                           Class radiationModel Declaration
+                       Class radiationModel Declaration
 \*---------------------------------------------------------------------------*/
 
 class radiationModel
 :
     public IOdictionary
 {
-
 protected:
 
     // Protected data
 
+        //- Reference to the mesh database
+        const fvMesh& mesh_;
+
+        //- Reference to the time database
+        const Time& time_;
+
         //- Reference to the temperature field
         const volScalarField& T_;
 
-        //- Reference to the mesh
-        const fvMesh& mesh_;
-
         //- Model specific dictionary input parameters
         Switch radiation_;
 
         //- Radiation model dictionary
-        dictionary radiationModelCoeffs_;
+        dictionary coeffs_;
+
+        //- Radiation solver frequency - number flow solver iterations per
+        //  radiation solver iteration
+        label solverFreq_;
 
 
         // References to the radiation sub-models
@@ -128,34 +135,32 @@ public:
 
     // Constructors
 
+        //- Null constructor
+        radiationModel(const volScalarField& T);
+
         //- Construct from components
-        radiationModel
-        (
-            const word& type,
-            const volScalarField& T
-        );
+        radiationModel(const word& type, const volScalarField& T);
 
 
     // Selectors
 
          //- Return a reference to the selected radiation model
-         static autoPtr<radiationModel> New
-         (
-             const volScalarField& T
-         );
+         static autoPtr<radiationModel> New(const volScalarField& T);
 
 
-    // Destructor
-
-        virtual ~radiationModel();
+    //- Destructor
+    virtual ~radiationModel();
 
 
     // Member Functions
 
         // Edit
 
+            //- Main update/correction routine
+            virtual void correct();
+
             //- Solve radiation equation(s)
-            virtual void correct() = 0;
+            virtual void calculate() = 0;
 
             //- Read radiationProperties dictionary
             virtual bool read() = 0;
@@ -170,10 +175,7 @@ public:
             virtual tmp<DimensionedField<scalar, volMesh> > Ru() const = 0;
 
             //- Enthalpy source term
-            virtual tmp<fvScalarMatrix> Sh
-            (
-                basicThermo& thermo
-            ) const;
+            virtual tmp<fvScalarMatrix> Sh(basicThermo& thermo) const;
 };
 
 

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.C

@@ -59,14 +59,14 @@ Foam::radiation::absorptionEmissionModel::~absorptionEmissionModel()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::a() const
+Foam::radiation::absorptionEmissionModel::a(const label bandI) const
 {
-    return aDisp() + aCont();
+    return aDisp(bandI) + aCont(bandI);
 }
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::aCont() const
+Foam::radiation::absorptionEmissionModel::aCont(const label bandI) const
 {
     return tmp<volScalarField>
     (
@@ -89,7 +89,7 @@ Foam::radiation::absorptionEmissionModel::aCont() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::aDisp() const
+Foam::radiation::absorptionEmissionModel::aDisp(const label bandI) const
 {
     return tmp<volScalarField>
     (
@@ -112,14 +112,14 @@ Foam::radiation::absorptionEmissionModel::aDisp() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::e() const
+Foam::radiation::absorptionEmissionModel::e(const label bandI) const
 {
-    return eDisp() + eCont();
+    return eDisp(bandI) + eCont(bandI);
 }
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::eCont() const
+Foam::radiation::absorptionEmissionModel::eCont(const label bandI) const
 {
     return tmp<volScalarField>
     (
@@ -142,7 +142,7 @@ Foam::radiation::absorptionEmissionModel::eCont() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::eDisp() const
+Foam::radiation::absorptionEmissionModel::eDisp(const label bandI) const
 {
     return tmp<volScalarField>
     (
@@ -165,14 +165,14 @@ Foam::radiation::absorptionEmissionModel::eDisp() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::E() const
+Foam::radiation::absorptionEmissionModel::E(const label bandI) const
 {
-    return EDisp() + ECont();
+    return EDisp(bandI) + ECont(bandI);
 }
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::ECont() const
+Foam::radiation::absorptionEmissionModel::ECont(const label bandI) const
 {
     return tmp<volScalarField>
     (
@@ -195,7 +195,7 @@ Foam::radiation::absorptionEmissionModel::ECont() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::EDisp() const
+Foam::radiation::absorptionEmissionModel::EDisp(const label bandI) const
 {
     return tmp<volScalarField>
     (
@@ -217,4 +217,45 @@ Foam::radiation::absorptionEmissionModel::EDisp() const
 }
 
 
+Foam::label Foam::radiation::absorptionEmissionModel::nBands() const
+{
+    return pTraits<label>::one;
+}
+
+
+const Foam::Vector2D<Foam::scalar>&
+Foam::radiation::absorptionEmissionModel::bands(const label n) const
+{
+    return Vector2D<scalar>::one;
+}
+
+
+bool Foam::radiation::absorptionEmissionModel::isGrey() const
+{
+    return false;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::absorptionEmissionModel::addIntensity
+(
+    const label rayI,
+    const volScalarField& ILambda
+) const
+{
+    return ILambda;
+}
+
+
+void Foam::radiation::absorptionEmissionModel::correct
+(
+    volScalarField& a,
+    PtrList<volScalarField>& aj
+) const
+{
+    a.internalField() = this->a();
+    aj[0].internalField() =  a.internalField();
+}
+
+
 // ************************************************************************* //

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.H

@@ -38,6 +38,7 @@ Description
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
 #include "volFields.H"
+#include "Vector2D.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -95,57 +96,94 @@ public:
 
 
     //- Selector
-
-        static autoPtr<absorptionEmissionModel> New
-        (
-            const dictionary& dict,
-            const fvMesh& mesh
-        );
+    static autoPtr<absorptionEmissionModel> New
+    (
+        const dictionary& dict,
+        const fvMesh& mesh
+    );
 
 
     //- Destructor
-
-        virtual ~absorptionEmissionModel();
+    virtual ~absorptionEmissionModel();
 
 
     // Member Functions
 
         // Access
 
+            //- Reference to the mesh
+            inline const fvMesh& mesh() const
+            {
+                return mesh_;
+            }
+
+            //- Reference to the dictionary
+            inline const dictionary& dict() const
+            {
+                return dict_;
+            }
+
+
             // Absorption coefficient
 
                 //- Absorption coefficient (net)
-                virtual tmp<volScalarField> a() const;
+                virtual tmp<volScalarField> a(const label bandI = 0) const;
 
                 //- Absorption coefficient for continuous phase
-                virtual tmp<volScalarField> aCont() const;
+                virtual tmp<volScalarField> aCont(const label bandI = 0) const;
 
                 //- Absorption coefficient for dispersed phase
-                virtual tmp<volScalarField> aDisp() const;
+                virtual tmp<volScalarField> aDisp(const label bandI = 0) const;
 
 
             // Emission coefficient
 
                 //- Emission coefficient (net)
-                virtual tmp<volScalarField> e() const;
+                virtual tmp<volScalarField> e(const label bandI = 0) const;
 
                 //- Return emission coefficient for continuous phase
-                virtual tmp<volScalarField> eCont() const;
+                virtual tmp<volScalarField> eCont(const label bandI = 0) const;
 
                 //- Return emission coefficient for dispersed phase
-                virtual tmp<volScalarField> eDisp() const;
+                virtual tmp<volScalarField> eDisp(const label bandI = 0) const;
 
 
             // Emission contribution
 
                 //- Emission contribution (net)
-                virtual tmp<volScalarField> E() const;
+                virtual tmp<volScalarField> E(const label bandI = 0) const;
 
                 //- Emission contribution for continuous phase
-                virtual tmp<volScalarField> ECont() const;
+                virtual tmp<volScalarField> ECont(const label bandI = 0) const;
 
                 //- Emission contribution for dispersed phase
-                virtual tmp<volScalarField> EDisp() const;
+                virtual tmp<volScalarField> EDisp(const label bandI = 0) const;
+
+
+            //- Const access to the number of bands - defaults to 1 for grey
+            //  absorption/emission
+            virtual label nBands() const;
+
+            //- Const access to the bands - defaults to Vector2D::one for grey
+            //  absorption/emission
+            virtual const Vector2D<scalar>& bands(const label n) const;
+
+            //- Flag for whether the absorption/emission is for a grey gas
+            virtual bool isGrey() const;
+
+            //- Add radiative intensity for ray i
+            virtual tmp<volScalarField> addIntensity
+            (
+                const label rayI,
+                const volScalarField& ILambda
+            ) const;
+
+            //- Correct absorption coefficients
+            virtual void correct
+            (
+                volScalarField& a,
+                PtrList<volScalarField>& aj
+            ) const;
 };
 
 

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/binaryAbsorptionEmission/binaryAbsorptionEmission.H

@@ -53,7 +53,6 @@ class binaryAbsorptionEmission
 :
     public absorptionEmissionModel
 {
-
     // Private data
 
         //- Coefficients dictionary
@@ -83,8 +82,7 @@ public:
 
 
     // Destructor
-
-        ~binaryAbsorptionEmission();
+    virtual ~binaryAbsorptionEmission();
 
 
     // Member Operators

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/constantAbsorptionEmission/constantAbsorptionEmission.C

@@ -70,7 +70,7 @@ Foam::radiation::constantAbsorptionEmission::~constantAbsorptionEmission()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::constantAbsorptionEmission::aCont() const
+Foam::radiation::constantAbsorptionEmission::aCont(const label bandI) const
 {
     tmp<volScalarField> ta
     (
@@ -95,7 +95,7 @@ Foam::radiation::constantAbsorptionEmission::aCont() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::constantAbsorptionEmission::eCont() const
+Foam::radiation::constantAbsorptionEmission::eCont(const label bandI) const
 {
     tmp<volScalarField> te
     (
@@ -120,7 +120,7 @@ Foam::radiation::constantAbsorptionEmission::eCont() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::constantAbsorptionEmission::ECont() const
+Foam::radiation::constantAbsorptionEmission::ECont(const label bandI) const
 {
     tmp<volScalarField> tE
     (

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/constantAbsorptionEmission/constantAbsorptionEmission.H

@@ -54,7 +54,6 @@ class constantAbsorptionEmission
 :
     public absorptionEmissionModel
 {
-
     // Private data
 
         //- Absorption model dictionary
@@ -87,8 +86,7 @@ public:
 
 
     // Destructor
-
-        ~constantAbsorptionEmission();
+    virtual ~constantAbsorptionEmission();
 
 
     // Member Operators
@@ -98,19 +96,27 @@ public:
             // Absorption coefficient
 
                 //- Absorption coefficient for continuous phase
-                tmp<volScalarField> aCont() const;
+                tmp<volScalarField> aCont(const label bandI = 0) const;
 
 
             // Emission coefficient
 
                 //- Emission coefficient for continuous phase
-                tmp<volScalarField> eCont() const;
+                tmp<volScalarField> eCont(const label bandI = 0) const;
 
 
             // Emission contribution
 
                 //- Emission contribution for continuous phase
-                tmp<volScalarField> ECont() const;
+                tmp<volScalarField> ECont(const label bandI = 0) const;
+
+
+    // Member Functions
+
+        inline bool isGrey() const
+        {
+            return true;
+        }
 };
 
 

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

@@ -0,0 +1,272 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "greyMeanAbsorptionEmission.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace radiation
+    {
+        defineTypeNameAndDebug(greyMeanAbsorptionEmission, 0);
+
+        addToRunTimeSelectionTable
+        (
+            absorptionEmissionModel,
+            greyMeanAbsorptionEmission,
+            dictionary
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::greyMeanAbsorptionEmission::greyMeanAbsorptionEmission
+(
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    absorptionEmissionModel(dict, mesh),
+    coeffsDict_((dict.subDict(typeName + "Coeffs"))),
+    speciesNames_(0),
+    specieIndex_(0),
+    lookUpTable_
+    (
+        fileName(coeffsDict_.lookup("lookUpTableFileName")),
+        mesh.time().constant(),
+        mesh
+    ),
+    thermo_(mesh.lookupObject<basicThermo>("thermophysicalProperties")),
+    EhrrCoeff_(readScalar(coeffsDict_.lookup("EhrrCoeff"))),
+    Yj_(nSpecies_)
+{
+    label nFunc = 0;
+    const dictionary& functionDicts = dict.subDict(typeName + "Coeffs");
+
+    forAllConstIter(dictionary, functionDicts, iter)
+    {
+        // safety:
+        if (!iter().isDict())
+        {
+            continue;
+        }
+        const word& key = iter().keyword();
+        speciesNames_.insert(key, nFunc);
+        const dictionary& dict = iter().dict();
+        coeffs_[nFunc].initialise(dict);
+        nFunc++;
+    }
+
+    // Check that all the species on the dictionary are present in the
+    // look-up table and save the corresponding indices of the look-up table
+
+    label j = 0;
+    forAllConstIter(HashTable<label>, speciesNames_, iter)
+    {
+        if (mesh.foundObject<volScalarField>("ft"))
+        {
+            if (lookUpTable_.found(iter.key()))
+            {
+                label index = lookUpTable_.findFieldIndex(iter.key());
+
+                Info<< "specie: " << iter.key() << " found on look-up table "
+                    << " with index: " << index << endl;
+
+                specieIndex_[iter()] = index;
+            }
+            else if (mesh.foundObject<volScalarField>(iter.key()))
+            {
+                volScalarField& Y =
+                    const_cast<volScalarField&>
+                    (
+                        mesh.lookupObject<volScalarField>(iter.key())
+                    );
+                Yj_.set(j, &Y);
+                specieIndex_[iter()] = 0;
+                j++;
+                Info<< "specie: " << iter.key() << " is being solved" << endl;
+            }
+            else
+            {
+                FatalErrorIn
+                (
+                    "Foam::radiation::greyMeanAbsorptionEmission(const"
+                    "dictionary& dict, const fvMesh& mesh)"
+                )   << "specie: " << iter.key()
+                    << " is neither in look-up table: "
+                    << lookUpTable_.tableName()
+                    << " nor is being solved" << nl
+                    << exit(FatalError);
+            }
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "Foam::radiation::greyMeanAbsorptionEmission(const"
+                "dictionary& dict, const fvMesh& mesh)"
+            )   << "specie ft is not present " << nl
+                << exit(FatalError);
+
+        }
+    }
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::radiation::greyMeanAbsorptionEmission::~greyMeanAbsorptionEmission()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::greyMeanAbsorptionEmission::aCont(const label bandI) const
+{
+    const volScalarField& T = thermo_.T();
+    const volScalarField& p = thermo_.p();
+    const volScalarField& ft = mesh_.lookupObject<volScalarField>("ft");
+
+    label nSpecies = speciesNames_.size();
+
+    tmp<volScalarField> ta
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "a",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("a", dimless/dimLength, 0.0)
+        )
+    );
+
+    scalarField& a = ta().internalField();
+
+    forAll(a, i)
+    {
+        const List<scalar>& species = lookUpTable_.lookUp(ft[i]);
+
+        for (label n=0; n<nSpecies; n++)
+        {
+            label l = 0;
+            scalar Yipi = 0;
+            if (specieIndex_[n] != 0)
+            {
+                //moles x pressure [atm]
+                Yipi = species[specieIndex_[n]]*p[i]*9.869231e-6;
+            }
+            else
+            {
+                // mass fraction
+                Yipi = Yj_[l][i];
+                l++;
+            }
+
+            const absorptionCoeffs::coeffArray& b = coeffs_[n].coeffs(T[i]);
+
+            scalar Ti = T[i];
+            // negative temperature exponents
+            if (coeffs_[n].invTemp())
+            {
+                Ti = 1./T[i];
+            }
+            a[i]+=
+                Yipi
+               *(
+                    ((((b[5]*Ti + b[4])*Ti + b[3])*Ti + b[2])*Ti + b[1])*Ti
+                  + b[0]
+                );
+        }
+    }
+    return ta;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::greyMeanAbsorptionEmission::eCont(const label bandI) const
+{
+    tmp<volScalarField> e
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "e",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("e", dimless/dimLength, 0.0)
+        )
+    );
+
+    return e;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::greyMeanAbsorptionEmission::ECont(const label bandI) const
+{
+    tmp<volScalarField> E
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "E",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
+        )
+    );
+
+    if (mesh_.foundObject<volScalarField>("hrr"))
+    {
+        const volScalarField& hrr = mesh_.lookupObject<volScalarField>("hrr");
+        E().internalField() = EhrrCoeff_*hrr.internalField();
+    }
+
+    return E;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.H

@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::radiation::greyMeanAbsorptionEmission
+
+Description
+    greyMeanAbsorptionEmission radiation absorption and emission coefficients
+    for continuous phase
+
+    The coefficients for the species in the Look up table have to be specified
+    for use in moles x P [atm], i.e. (k[i] = species[i]*p*9.869231e-6).
+
+    The coefficients for CO and soot or any other added are multiplied by the
+    respective mass fraction being solved
+
+    All the species in the dictionary need either to be in the look-up table or
+    being solved. Conversely, all the species solved do not need to be included
+    in the calculation of the absorption coefficient
+
+    The names of the species in the absorption dictionary must match exactly the
+    name in the look-up table or the name of the field being solved
+
+    The look-up table ("speciesTable") file should be in constant
+
+    i.e. dictionary
+
+    LookUpTableFileName     "speciesTable";
+
+    EhrrCoeff                0.0;
+
+    CO2
+    {
+        Tcommon         300.;   // Common Temp
+        invTemp         true;   // Is the polynomial using inverse temperature?
+        Tlow            300.;   // Low Temp
+        Thigh           2500.;  // High Temp
+
+        loTcoeffs       // coeffs for T < Tcommon
+        (
+            0           //  a0            +
+            0           //  a1*T          +
+            0           //  a2*T^(+/-)2   +
+            0           //  a3*T^(+/-)3   +
+            0           //  a4*T^(+/-)4   +
+            0           //  a5*T^(+/-)5   +
+        );
+        hiTcoeffs       // coeffs for T > Tcommon
+        (
+            18.741
+            -121.31e3
+            273.5e6
+            -194.05e9
+            56.31e12
+            -5.8169e15
+        );
+
+    }
+
+SourceFiles
+    greyMeanAbsorptionEmission.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef greyMeanAbsorptionEmission_H
+#define greyMeanAbsorptionEmission_H
+
+#include "interpolationLookUpTable.H"
+#include "absorptionEmissionModel.H"
+#include "HashTable.H"
+#include "absorptionCoeffs.H"
+#include "basicThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class constantAbsorptionEmission Declaration
+\*---------------------------------------------------------------------------*/
+
+class greyMeanAbsorptionEmission
+:
+    public absorptionEmissionModel
+{
+public:
+
+    // Public data
+
+        // - Maximum number of species considered for absorptivity
+        static const int nSpecies_ = 5;
+
+        //  Absorption Coefficients
+        absorptionCoeffs coeffs_[nSpecies_];
+
+
+private:
+
+    // Private data
+
+        //- Absorption model dictionary
+        dictionary coeffsDict_;
+
+        //- Hash table of species names
+        HashTable<label> speciesNames_;
+
+        // Indices of species in the look-up table
+        FixedList<label, nSpecies_> specieIndex_;
+
+        // Look-up table of species related to ft
+        mutable interpolationLookUpTable<scalar> lookUpTable_;
+
+        // Thermo package
+        const basicThermo& thermo_;
+
+        //- Emission constant coefficient
+        const scalar EhrrCoeff_;
+
+        //- Pointer list of species in the registry involved in the absorption
+        UPtrList<volScalarField> Yj_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("greyMeanAbsorptionEmission");
+
+
+    // Constructors
+
+        //- Construct from components
+        greyMeanAbsorptionEmission
+        (
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    // Destructor
+    virtual ~greyMeanAbsorptionEmission();
+
+
+    // Member Operators
+
+        // Access
+
+            // Absorption coefficient
+
+                //- Absorption coefficient for continuous phase
+                tmp<volScalarField> aCont(const label bandI = 0) const;
+
+
+            // Emission coefficient
+
+                //- Emission coefficient for continuous phase
+                tmp<volScalarField> eCont(const label bandI = 0) const;
+
+
+            // Emission contribution
+
+                //- Emission contribution for continuous phase
+                tmp<volScalarField> ECont(const label bandI = 0) const;
+
+
+    // Member Functions
+
+        inline bool isGrey() const
+        {
+            return true;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace radiation
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/noAbsorptionEmission/noAbsorptionEmission.H

@@ -71,8 +71,7 @@ public:
 
 
     // Destructor
-
-        ~noAbsorptionEmission();
+    virtual ~noAbsorptionEmission();
 };
 
 

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

@@ -0,0 +1,320 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "wideBandAbsorptionEmission.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace radiation
+    {
+        defineTypeNameAndDebug(wideBandAbsorptionEmission, 0);
+
+        addToRunTimeSelectionTable
+        (
+            absorptionEmissionModel,
+            wideBandAbsorptionEmission,
+            dictionary
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::wideBandAbsorptionEmission::wideBandAbsorptionEmission
+(
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    absorptionEmissionModel(dict, mesh),
+    coeffsDict_((dict.subDict(typeName + "Coeffs"))),
+    speciesNames_(0),
+    specieIndex_(0),
+    lookUpTable_
+    (
+        fileName(coeffsDict_.lookup("lookUpTableFileName")),
+        mesh.time().constant(),
+        mesh
+    ),
+    thermo_(mesh.lookupObject<basicThermo>("thermophysicalProperties")),
+    Yj_(nSpecies_),
+    totalWaveLength_(0)
+{
+    label nBand = 0;
+    const dictionary& functionDicts = dict.subDict(typeName +"Coeffs");
+    forAllConstIter(dictionary, functionDicts, iter)
+    {
+        // safety:
+        if (!iter().isDict())
+        {
+            continue;
+        }
+
+        const dictionary& dict = iter().dict();
+        dict.lookup("bandLimits") >> iBands_[nBand];
+        dict.lookup("EhrrCoeff") >> iEhrrCoeffs_[nBand];
+        totalWaveLength_ += (iBands_[nBand][1] - iBands_[nBand][0]);
+
+        label nSpec = 0;
+        const dictionary& specDicts = dict.subDict("species");
+        forAllConstIter(dictionary, specDicts, iter)
+        {
+            const word& key = iter().keyword();
+            if (nBand == 0)
+            {
+                speciesNames_.insert(key, nSpec);
+            }
+            else
+            {
+                if (!speciesNames_.found(key))
+                {
+                    FatalErrorIn
+                    (
+                        "Foam::radiation::wideBandAbsorptionEmission(const"
+                        "dictionary& dict, const fvMesh& mesh)"
+                    )   << "specie: " << key << "is not in all the bands"
+                        << nl << exit(FatalError);
+                }
+            }
+            coeffs_[nSpec][nBand].initialise(specDicts.subDict(key));
+            nSpec++;
+        }
+        nBand++;
+    }
+    nBands_ = nBand;
+
+    // Check that all the species on the dictionary are present in the
+    // look-up table  and save the corresponding indexes of the look-up table
+
+    label j = 0;
+    forAllConstIter(HashTable<label>, speciesNames_, iter)
+    {
+        if (lookUpTable_.found(iter.key()))
+        {
+            label index = lookUpTable_.findFieldIndex(iter.key());
+            Info<< "specie: " << iter.key() << " found in look-up table "
+                << " with index: " << index << endl;
+            specieIndex_[iter()] = index;
+        }
+        else if (mesh.foundObject<volScalarField>(iter.key()))
+        {
+            volScalarField& Y = const_cast<volScalarField&>
+                (mesh.lookupObject<volScalarField>(iter.key()));
+
+            Yj_.set(j, &Y);
+
+            specieIndex_[iter()] = 0.0;
+            j++;
+            Info<< "species: " << iter.key() << " is being solved" << endl;
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "radiation::wideBandAbsorptionEmission(const"
+                "dictionary& dict, const fvMesh& mesh)"
+            )   << "specie: " << iter.key()
+                << " is neither in look-up table : "
+                << lookUpTable_.tableName() << " nor is being solved"
+                << exit(FatalError);
+        }
+    }
+}
+
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::radiation::wideBandAbsorptionEmission::~wideBandAbsorptionEmission()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::wideBandAbsorptionEmission::aCont(const label bandI) const
+{
+    const volScalarField& T = thermo_.T();
+    const volScalarField& p = thermo_.p();
+    const volScalarField& ft = mesh_.lookupObject<volScalarField>("ft");
+
+    label nSpecies = speciesNames_.size();
+
+    tmp<volScalarField> ta
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "a",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("a",dimless/dimLength, 0.0)
+        )
+    );
+
+    scalarField& a = ta().internalField();
+
+    forAll(a, i)
+    {
+        const List<scalar>& species = lookUpTable_.lookUp(ft[i]);
+
+        for (label n=0; n<nSpecies; n++)
+        {
+            label l = 0;
+            scalar Yipi = 0;
+            if (specieIndex_[n] != 0)
+            {
+                // moles x pressure [atm]
+                Yipi = species[specieIndex_[n]]*p[i]*9.869231e-6;
+            }
+            else
+            {
+                // mass fraction from species being solved
+                Yipi = Yj_[l][i];
+                l++;
+            }
+
+            scalar Ti = T[i];
+
+            const absorptionCoeffs::coeffArray& b =
+                coeffs_[n][bandI].coeffs(T[i]);
+
+            if (coeffs_[n][bandI].invTemp())
+            {
+                Ti = 1./T[i];
+            }
+
+            a[i]+=
+                Yipi
+               *(
+                    ((((b[5]*Ti + b[4])*Ti + b[3])*Ti + b[2])*Ti + b[1])*Ti
+                  + b[0]
+                );
+        }
+    }
+
+    return ta;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::wideBandAbsorptionEmission::eCont(const label bandI) const
+{
+    tmp<volScalarField> e
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "e",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("e",dimless/dimLength, 0.0)
+        )
+    );
+
+    return e;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::wideBandAbsorptionEmission::ECont(const label bandI) const
+{
+    tmp<volScalarField> E
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "E",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
+        )
+    );
+
+    if (mesh().foundObject<volScalarField>("hrr"))
+    {
+        const volScalarField& hrr =  mesh().lookupObject<volScalarField>("hrr");
+        E().internalField() =
+            iEhrrCoeffs_[bandI]
+           *hrr.internalField()
+           *(iBands_[bandI][1] - iBands_[bandI][0])
+           /totalWaveLength_;
+    }
+
+    return E;
+}
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::wideBandAbsorptionEmission::addIntensity
+(
+    const label i,
+    const volScalarField& ILambda
+) const
+{
+    return ILambda*(iBands_[i][1] - iBands_[i][0])/totalWaveLength_;
+}
+
+
+void Foam::radiation::wideBandAbsorptionEmission::correct
+(
+    volScalarField& a,
+    PtrList<volScalarField>& aLambda
+) const
+{
+    a = dimensionedScalar("zero", dimless/dimLength, 0.0);
+
+    for (label j = 0; j < nBands_; j++)
+    {
+        Info<< "Calculating absorption in band: " << j << endl;
+        aLambda[j].internalField() = this->a(j);
+        Info<< "Calculated absorption in band: " << j << endl;
+        a.internalField() +=
+            aLambda[j].internalField()
+           *(iBands_[j][1] - iBands_[j][0])
+           /totalWaveLength_;
+    }
+
+}
+// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.H

@@ -0,0 +1,264 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::radiation::greyMeanAbsorptionEmission
+
+Description
+
+    wideBandAbsorptionEmission radiation absorption and emission coefficients
+    for continuous phase.
+
+    All the bands should have the same number of species and have to be entered
+    in the same order.
+
+    There is no check of continuity of the bands. They should not ovelap or
+    have gaps.
+
+    The black body emission power table(constant/blackBodyEmissivePower) range
+    of lambda * T = [1000; 10000] x 10E-6 (90% of the total emission).
+
+    The emission constant proportionality is specified per band (EhrrCoeff).
+
+    The coefficients for the species in the LookUpTable have to be specified
+    for use in moles x P [atm].i.e. (k[i] = species[i]*p*9.869231e-6).
+
+    The coefficients for CO and soot or any other added are multiplied by the
+    respective mass fraction being solved.
+
+    The look Up table file should be in the constant directory.
+
+    band dictionary:
+
+    band0
+    {
+        bandLimits (1.0e-6 2.63e-6);
+        EhrrCoeff       0.0;
+        species
+        {
+            CH4
+            {
+                Tcommon         300.;
+                Tlow            300.;
+                Thigh           2500.;
+                invTemp         false;
+                loTcoeffs (0 0 0 0 0 0) ;
+                hiTcoeffs (.1 0 0 0 0 0);
+            }
+            CO2
+            {
+                Tcommon         300.;
+                Tlow            300.;
+                Thigh           2500.;
+                invTemp         false;
+                loTcoeffs (0 0 0 0 0 0) ;
+                hiTcoeffs (.1 0 0 0 0 0);
+            }
+
+            H2O
+            {
+                Tcommon         300.;
+                Tlow            300.;
+                Thigh           2500.;
+                invTemp         false;
+                loTcoeffs (0 0 0 0 0 0) ;
+                hiTcoeffs (.1 0 0 0 0 0);
+            }
+
+            Ysoot
+            {
+                Tcommon         300.;
+                Tlow            300.;
+                Thigh           2500.;
+                invTemp         false;
+                loTcoeffs (0 0 0 0 0 0) ;
+                hiTcoeffs (.1 0 0 0 0 0);
+            }
+        }
+    }
+
+
+SourceFiles
+    wideBandAbsorptionEmission.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wideBandAbsorptionEmission_H
+#define wideBandAbsorptionEmission_H
+
+#include "interpolationLookUpTable.H"
+#include "absorptionEmissionModel.H"
+#include "HashTable.H"
+#include "absorptionCoeffs.H"
+#include "basicThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class wideBandAbsorptionEmission Declaration
+\*---------------------------------------------------------------------------*/
+
+class wideBandAbsorptionEmission
+:
+    public absorptionEmissionModel
+{
+public:
+
+    // Public data
+
+        //- Maximum number of species considered for absorptivity
+        static const int nSpecies_ = 5;
+
+        //- Maximum number of bands
+        static const int maxBands_ = 10;
+
+        //-  Absorption coefficients
+        FixedList<FixedList<absorptionCoeffs, nSpecies_>, maxBands_> coeffs_;
+
+
+private:
+
+    // Private data
+
+        //- Absorption model dictionary
+        dictionary coeffsDict_;
+
+        //- Hash table with species names
+        HashTable<label> speciesNames_;
+
+        //- Indices of species in the look-up table
+        FixedList<label,nSpecies_> specieIndex_;
+
+        //- Bands
+        FixedList<Vector2D<scalar>, maxBands_ > iBands_;
+
+        //- Proportion of the heat released rate emitted
+        FixedList<scalar, maxBands_ > iEhrrCoeffs_;
+
+        //- Look-up table of species related to ft
+        mutable interpolationLookUpTable<scalar> lookUpTable_;
+
+        //- Thermo package
+        const basicThermo& thermo_;
+
+        //- Bands
+        label nBands_ ;
+
+        //- Pointer list of species being solved involved in the absorption
+        UPtrList<volScalarField> Yj_;
+
+        // Total wave length covered by the bands
+        scalar totalWaveLength_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("wideBandAbsorptionEmission");
+
+
+    // Constructors
+
+        //- Construct from components
+        wideBandAbsorptionEmission
+        (
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    // Destructor
+    virtual ~wideBandAbsorptionEmission();
+
+
+    // Member Functions
+
+        // Access
+
+            // Absorption coefficient
+
+                //- Absorption coefficient for continuous phase
+                tmp<volScalarField> aCont(const label bandI = 0) const;
+
+
+            // Emission coefficient
+
+                //- Emission coefficient for continuous phase
+                tmp<volScalarField> eCont(const label bandI = 0) const;
+
+
+            // Emission contribution
+
+                //- Emission contribution for continuous phase
+                tmp<volScalarField> ECont(const label bandI = 0) const;
+
+
+        inline bool isGrey() const
+        {
+            return false;
+        }
+
+        //- Number of bands
+        inline label nBands() const
+        {
+            return nBands_;
+        }
+
+        //- Lower and upper limit of band i
+        inline const Vector2D<scalar>& bands(const label i) const
+        {
+            return iBands_[i];
+        }
+
+        //- Add contribution of ILambda to the total radiative intensity in
+        //  direction i
+        tmp<volScalarField> addIntensity
+        (
+            const label i,
+            const volScalarField& ILambda
+        ) const;
+
+        void correct
+        (
+            volScalarField& a_,
+            PtrList<volScalarField>& aLambda
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace radiation
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/triSurface/faceTriangulation/faceTriangulation.C

@@ -88,7 +88,7 @@ void Foam::faceTriangulation::calcHalfAngle
 
     scalar sin = (e0 ^ e1) & normal;
 
-    if (sin < -SMALL)
+    if (sin < -ROOTVSMALL)
     {
         // 3rd or 4th quadrant
         cosHalfAngle = - Foam::sqrt(0.5*(1 + cos));
@@ -366,7 +366,7 @@ Foam::label Foam::faceTriangulation::findStart
         const vector& rightEdge = edges[right(size, fp)];
         const vector leftEdge = -edges[left(size, fp)];
 
-        if (((rightEdge ^ leftEdge) & normal) < SMALL)
+        if (((rightEdge ^ leftEdge) & normal) < ROOTVSMALL)
         {
             scalar cos = rightEdge & leftEdge;
             if (cos < minCos)