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

Conflicts:
	src/postProcessing/functionObjects/field/Make/files

applications/solvers/lagrangian/sprayFoam/pEqn.H

@@ -53,7 +53,7 @@ else
         )
     );
 
-    fvOptions.makeRelative(fvc::interpolate(psi), phiHbyA);
+    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
     {

applications/test/ODE/Test-ODE.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -142,7 +142,7 @@ int main(int argc, char *argv[])
         scalar dxNext = dxEst;
 
         odeSolver->relTol() = relTol;
-        odeSolver->solve(ode, x, y, dxNext);
+        odeSolver->solve(x, y, dxNext);
 
         Info<< scientific << setw(13) << relTol;
         Info<< fixed << setw(11) << dxEst;
@@ -165,7 +165,7 @@ int main(int argc, char *argv[])
     scalar dxEst = 0.5;
 
     odeSolver->relTol() = 1e-4;
-    odeSolver->solve(ode, x, xEnd, y, dxEst);
+    odeSolver->solve(x, xEnd, y, dxEst);
 
     Info<< nl << "Analytical: y(2.0) = " << yEnd << endl;
     Info      << "Numerical:  y(2.0) = " << y << ", dxEst = " << dxEst << endl;

etc/bashrc

@@ -62,7 +62,7 @@ foamInstall=$HOME/$WM_PROJECT
 foamCompiler=system
 
 #- Compiler:
-#    WM_COMPILER = Gcc | Gcc43 | Gcc44 | Gcc45 | Gcc46 | Clang | Icc (Intel icc)
+#    WM_COMPILER = Gcc | Gcc45 | Gcc46 | Gcc47 | Clang | Icc (Intel icc)
 export WM_COMPILER=Gcc
 unset WM_COMPILER_ARCH WM_COMPILER_LIB_ARCH
 

etc/cshrc

@@ -61,7 +61,7 @@ if ( ! $?FOAM_INST_DIR ) setenv FOAM_INST_DIR $foamInstall
 set foamCompiler=system
 
 #- Compiler:
-#    WM_COMPILER = Gcc | Gcc43 | Gcc44 | Gcc45 | Gcc46 | Clang | Icc (Intel icc)
+#    WM_COMPILER = Gcc | Gcc45 | Gcc46 | Gcc47 | Clang | Icc (Intel icc)
 setenv WM_COMPILER Gcc
 setenv WM_COMPILER_ARCH # defined but empty
 unsetenv WM_COMPILER_LIB_ARCH

src/ODE/Make/files

@@ -8,13 +8,14 @@ ODESolvers/Trapezoid/Trapezoid.C
 ODESolvers/RKF45/RKF45.C
 ODESolvers/RKCK45/RKCK45.C
 ODESolvers/RKDP45/RKDP45.C
-ODESolvers/Rosenbrock21/Rosenbrock21.C
+ODESolvers/Rosenbrock12/Rosenbrock12.C
-ODESolvers/Rosenbrock32/Rosenbrock32.C
+ODESolvers/Rosenbrock23/Rosenbrock23.C
-ODESolvers/Rosenbrock43/Rosenbrock43.C
+ODESolvers/Rosenbrock34/Rosenbrock34.C
-ODESolvers/rodas32/rodas32.C
+ODESolvers/rodas23/rodas23.C
-ODESolvers/rodas43/rodas43.C
+ODESolvers/rodas34/rodas34.C
 ODESolvers/SIBS/SIBS.C
 ODESolvers/SIBS/SIMPR.C
 ODESolvers/SIBS/polyExtrapolate.C
+ODESolvers/seulex/seulex.C
 
 LIB = $(FOAM_LIBBIN)/libODE

src/ODE/ODESolvers/Euler/Euler.C

@@ -49,7 +49,6 @@ Foam::Euler::Euler(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::Euler::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -75,13 +74,12 @@ Foam::scalar Foam::Euler::solve
 
 void Foam::Euler::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/Euler/Euler.H

@@ -84,7 +84,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -95,7 +94,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/EulerSI/EulerSI.C

@@ -54,7 +54,6 @@ Foam::EulerSI::EulerSI(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::EulerSI::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -62,7 +61,7 @@ Foam::scalar Foam::EulerSI::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -95,13 +94,12 @@ Foam::scalar Foam::EulerSI::solve
 
 void Foam::EulerSI::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/EulerSI/EulerSI.H

@@ -91,7 +91,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -102,7 +101,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/ODESolver/ODESolver.C

@@ -38,6 +38,7 @@ namespace Foam
 
 Foam::ODESolver::ODESolver(const ODESystem& ode, const dictionary& dict)
 :
+    odes_(ode),
     n_(ode.nEqns()),
     absTol_(n_, dict.lookupOrDefault<scalar>("absTol", SMALL)),
     relTol_(n_, dict.lookupOrDefault<scalar>("relTol", 1e-4)),
@@ -52,6 +53,7 @@ Foam::ODESolver::ODESolver
     const scalarField& relTol
 )
 :
+    odes_(ode),
     n_(ode.nEqns()),
     absTol_(absTol),
     relTol_(relTol),
@@ -82,48 +84,70 @@ Foam::scalar Foam::ODESolver::normalizeError
 
 void Foam::ODESolver::solve
 (
-    const ODESystem& ode,
+    scalar& x,
+    scalarField& y,
+    stepState& step
+) const
+{
+    solve(x, y, step.dxTry);
+}
+
+
+void Foam::ODESolver::solve
+(
     const scalar xStart,
     const scalar xEnd,
     scalarField& y,
-    scalar& dxEst
+    scalar& dxTry
 ) const
 {
+    stepState step(dxTry);
     scalar x = xStart;
-    bool truncated = false;
 
     for (label nStep=0; nStep<maxSteps_; nStep++)
     {
-        // Store previous iteration dxEst
+        // Store previous iteration dxTry
-        scalar dxEst0 = dxEst;
+        scalar dxTry0 = step.dxTry;
 
-        // Check if this is a truncated step and set dxEst to integrate to xEnd
+        step.reject = false;
-        if ((x + dxEst - xEnd)*(x + dxEst - xStart) > 0)
+
+        // Check if this is a truncated step and set dxTry to integrate to xEnd
+        if ((x + step.dxTry - xEnd)*(x + step.dxTry - xStart) > 0)
         {
-            truncated = true;
+            step.last = true;
-            dxEst = xEnd - x;
+            step.dxTry = xEnd - x;
         }
 
-        // Integrate as far as possible up to dxEst
+        // Integrate as far as possible up to step.dxTry
-        solve(ode, x, y, dxEst);
+        solve(x, y, step);
 
         // Check if reached xEnd
         if ((x - xEnd)*(xEnd - xStart) >= 0)
         {
-            if (nStep > 0 && truncated)
+            if (nStep > 0 && step.last)
             {
-                dxEst = dxEst0;
+                step.dxTry = dxTry0;
             }
 
+            dxTry = step.dxTry;
+
             return;
         }
+
+        step.first = false;
+
+        // If the step.dxTry was reject set step.prevReject
+        if (step.reject)
+        {
+            step.prevReject = true;
+        }
     }
 
     FatalErrorIn
     (
         "ODESolver::solve"
-        "(const ODESystem& ode, const scalar xStart, const scalar xEnd,"
+        "(const scalar xStart, const scalar xEnd,"
-        "scalarField& y, scalar& dxEst) const"
+        "scalarField& y, scalar& dxTry) const"
     )   << "Integration steps greater than maximum " << maxSteps_
         << exit(FatalError);
 }

src/ODE/ODESolvers/ODESolver/ODESolver.H

@@ -55,6 +55,9 @@ protected:
 
     // Protected data
 
+        //- Reference to ODESystem
+        const ODESystem& odes_;
+
         //- Size of the ODESystem
         label n_;
 
@@ -90,6 +93,30 @@ public:
     //- Runtime type information
     TypeName("ODESolver");
 
+    class stepState
+    {
+        public:
+
+        const bool forward;
+        scalar dxTry;
+        scalar dxDid;
+        bool first;
+        bool last;
+        bool reject;
+        bool prevReject;
+
+        stepState(const scalar dx)
+        :
+            forward(dx > 0 ? true : false),
+            dxTry(dx),
+            dxDid(0),
+            first(true),
+            last(false),
+            reject(false),
+            prevReject(false)
+        {}
+    };
+
 
     // Declare run-time constructor selection table
 
@@ -150,17 +177,31 @@ public:
         //  Update the state and return an estimate for the next step in dxTry
         virtual void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry
-        ) const = 0;
+        ) const //= 0;
+        {
+            stepState step(dxTry);
+            solve(x, y, step);
+            dxTry = step.dxTry;
+        }
+
+        //- Solve the ODE system as far as possible upto dxTry
+        //  adjusting the step as necessary to provide a solution within
+        //  the specified tolerance.
+        //  Update the state and return an estimate for the next step in dxTry
+        virtual void solve
+        (
+            scalar& x,
+            scalarField& y,
+            stepState& step
+        ) const;
 
         //- Solve the ODE system from xStart to xEnd, update the state
         //  and return an estimate for the next step in dxTry
         virtual void solve
         (
-            const ODESystem& ode,
             const scalar xStart,
             const scalar xEnd,
             scalarField& y,

src/ODE/ODESolvers/RKCK45/RKCK45.C

@@ -89,7 +89,6 @@ Foam::RKCK45::RKCK45(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::RKCK45::solve
 (
-    const ODESystem& odes,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -102,21 +101,21 @@ Foam::scalar Foam::RKCK45::solve
         yTemp_[i] = y0[i] + a21*dx*dydx0[i];
     }
 
-    odes.derivatives(x0 + c2*dx, yTemp_, k2_);
+    odes_.derivatives(x0 + c2*dx, yTemp_, k2_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]);
     }
 
-    odes.derivatives(x0 + c3*dx, yTemp_, k3_);
+    odes_.derivatives(x0 + c3*dx, yTemp_, k3_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]);
     }
 
-    odes.derivatives(x0 + c4*dx, yTemp_, k4_);
+    odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
 
     forAll(yTemp_, i)
     {
@@ -124,7 +123,7 @@ Foam::scalar Foam::RKCK45::solve
           + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
     }
 
-    odes.derivatives(x0 + c5*dx, yTemp_, k5_);
+    odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
 
     forAll(yTemp_, i)
     {
@@ -133,7 +132,7 @@ Foam::scalar Foam::RKCK45::solve
            *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
     }
 
-    odes.derivatives(x0 + c6*dx, yTemp_, k6_);
+    odes_.derivatives(x0 + c6*dx, yTemp_, k6_);
 
     forAll(y, i)
     {
@@ -154,13 +153,12 @@ Foam::scalar Foam::RKCK45::solve
 
 void Foam::RKCK45::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/RKCK45/RKCK45.H

@@ -34,10 +34,7 @@ Description
         Cash, J.R.,
         Karp, A.H.
         ACM Transactions on Mathematical Software, vol. 16, 1990, pp. 201–222.
-    \endverbatim
 
-    Based on code from:
-    \verbatim
         "Solving Ordinary Differential Equations I: Nonstiff Problems,
          second edition",
         Hairer, E.,
@@ -110,7 +107,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -121,7 +117,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/RKDP45/RKDP45.C

@@ -94,7 +94,6 @@ Foam::RKDP45::RKDP45(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::RKDP45::solve
 (
-    const ODESystem& odes,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -107,21 +106,21 @@ Foam::scalar Foam::RKDP45::solve
         yTemp_[i] = y0[i] + a21*dx*dydx0[i];
     }
 
-    odes.derivatives(x0 + c2*dx, yTemp_, k2_);
+    odes_.derivatives(x0 + c2*dx, yTemp_, k2_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]);
     }
 
-    odes.derivatives(x0 + c3*dx, yTemp_, k3_);
+    odes_.derivatives(x0 + c3*dx, yTemp_, k3_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]);
     }
 
-    odes.derivatives(x0 + c4*dx, yTemp_, k4_);
+    odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
 
     forAll(yTemp_, i)
     {
@@ -129,7 +128,7 @@ Foam::scalar Foam::RKDP45::solve
           + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
     }
 
-    odes.derivatives(x0 + c5*dx, yTemp_, k5_);
+    odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
 
     forAll(yTemp_, i)
     {
@@ -138,7 +137,7 @@ Foam::scalar Foam::RKDP45::solve
            *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
     }
 
-    odes.derivatives(x0 + dx, yTemp_, k6_);
+    odes_.derivatives(x0 + dx, yTemp_, k6_);
 
     forAll(y, i)
     {
@@ -147,7 +146,7 @@ Foam::scalar Foam::RKDP45::solve
     }
 
     // Reuse k2_ for the derivative of the new state
-    odes.derivatives(x0 + dx, y, k2_);
+    odes_.derivatives(x0 + dx, y, k2_);
 
     forAll(err_, i)
     {
@@ -163,13 +162,12 @@ Foam::scalar Foam::RKDP45::solve
 
 void Foam::RKDP45::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/RKDP45/RKDP45.H

@@ -34,10 +34,7 @@ Description
         Prince, P. J.,
         Journal of Computational and Applied Mathematics,
         6 (1), 1980: pp. 19-26.
-    \endverbatim
 
-    Based on code from:
-    \verbatim
         "Solving Ordinary Differential Equations I: Nonstiff Problems,
          second edition",
         Hairer, E.,
@@ -114,7 +111,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -125,7 +121,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/RKF45/RKF45.C

@@ -90,7 +90,6 @@ Foam::RKF45::RKF45(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::RKF45::solve
 (
-    const ODESystem& odes,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -103,21 +102,21 @@ Foam::scalar Foam::RKF45::solve
         yTemp_[i] = y0[i] + a21*dx*dydx0[i];
     }
 
-    odes.derivatives(x0 + c2*dx, yTemp_, k2_);
+    odes_.derivatives(x0 + c2*dx, yTemp_, k2_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]);
     }
 
-    odes.derivatives(x0 + c3*dx, yTemp_, k3_);
+    odes_.derivatives(x0 + c3*dx, yTemp_, k3_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]);
     }
 
-    odes.derivatives(x0 + c4*dx, yTemp_, k4_);
+    odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
 
     forAll(yTemp_, i)
     {
@@ -125,7 +124,7 @@ Foam::scalar Foam::RKF45::solve
           + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
     }
 
-    odes.derivatives(x0 + c5*dx, yTemp_, k5_);
+    odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
 
     forAll(yTemp_, i)
     {
@@ -134,7 +133,7 @@ Foam::scalar Foam::RKF45::solve
            *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
     }
 
-    odes.derivatives(x0 + c6*dx, yTemp_, k6_);
+    odes_.derivatives(x0 + c6*dx, yTemp_, k6_);
 
     // Calculate the 5th-order solution
     forAll(y, i)
@@ -159,13 +158,12 @@ Foam::scalar Foam::RKF45::solve
 
 void Foam::RKF45::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/RKF45/RKF45.H

@@ -34,8 +34,6 @@ Description
         Fehlberg, E.,
         NASA Technical Report 315, 1969.
 
-    Based on code from:
-    \verbatim
         "Solving Ordinary Differential Equations I: Nonstiff Problems,
          second edition",
         Hairer, E.,
@@ -48,31 +46,6 @@ Description
     and allows to perform an adapdive step-size control using these two order
     without the need of re-evaluation.
 
-    \f{align}{
-     k_1 &= h f(t_k, y_k) \\
-     k_2 &= h f(t_k + \frac{1}{4}h, y_k + \frac{1}{4}k_1) \\
-     k_3 &= h f(t_k + \frac{3}{8}h, y_k + \frac{3}{32}k_1 + \frac{9}{32}k_2) \\
-     k_4 &= h f(t_k + \frac{12}{13}h, y_k + \frac{1932}{2197}k_1
-              \frac{7200}{2197}k_2 + \frac{7296}{2197}k_3) - \\
-     k_5 &= h f(t_k + h, y_k + \frac{439}{216}k_1 - 8k_2 + \frac{3680}{513}k_3 -
-           \frac{845}{4104}k_4) \\
-     k_6 &= h f(t_k + \frac{1}{2}h, y_k - \frac{8}{27}k_1 + 2k_2 -
-           \frac{3544}{2565}k_3 + \frac{1859}{4104}k_4 - \frac{11}{40}k_5) \\
-    \f}
-
-    Which yields the following update-steps for the 4th and 5th order:
-
-    \f{align}{
-     \Delta_4 &= \frac{25}{216}k_1 + \frac{1408}{2565}k_3 +
-                 \frac{2197}{4101}k_4 - \frac{1}{5}k_5 \\
-     \Delta_5 &= \frac{16}{135}k_1 + \frac{6656}{12825}k_3 +
-                 \frac{9}{50}k_5 + \frac{2}{55}k_6
-    \f}
-
-    The difference between the 5th and 4th order (\f$\epsilon = \left|\Delta_5
-    - \Delta_4\right|\f$) can be used to determine if the stepsize \f$h\f$ needs
-    to be adjusted.
-
 SourceFiles
     RKF45.C
 
@@ -138,7 +111,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -149,7 +121,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C

@@ -23,33 +23,33 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "Rosenbrock21.H"
+#include "Rosenbrock12.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(Rosenbrock21, 0);
+    defineTypeNameAndDebug(Rosenbrock12, 0);
-    addToRunTimeSelectionTable(ODESolver, Rosenbrock21, dictionary);
+    addToRunTimeSelectionTable(ODESolver, Rosenbrock12, dictionary);
 
 const scalar
-    Rosenbrock21::gamma = 1 + 1.0/std::sqrt(2.0),
+    Rosenbrock12::gamma = 1 + 1.0/std::sqrt(2.0),
-    Rosenbrock21::a21 = 1.0/gamma,
+    Rosenbrock12::a21 = 1.0/gamma,
-    Rosenbrock21::c2 = 1.0,
+    Rosenbrock12::c2 = 1.0,
-    Rosenbrock21::c21 = -2.0/gamma,
+    Rosenbrock12::c21 = -2.0/gamma,
-    Rosenbrock21::b1 = (3.0/2.0)/gamma,
+    Rosenbrock12::b1 = (3.0/2.0)/gamma,
-    Rosenbrock21::b2 = (1.0/2.0)/gamma,
+    Rosenbrock12::b2 = (1.0/2.0)/gamma,
-    Rosenbrock21::e1 = b1 - 1.0/gamma,
+    Rosenbrock12::e1 = b1 - 1.0/gamma,
-    Rosenbrock21::e2 = b2,
+    Rosenbrock12::e2 = b2,
-    Rosenbrock21::d1 = gamma,
+    Rosenbrock12::d1 = gamma,
-    Rosenbrock21::d2 = -gamma;
+    Rosenbrock12::d2 = -gamma;
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::Rosenbrock21::Rosenbrock21(const ODESystem& ode, const dictionary& dict)
+Foam::Rosenbrock12::Rosenbrock12(const ODESystem& ode, const dictionary& dict)
 :
     ODESolver(ode, dict),
     adaptiveSolver(ode, dict),
@@ -66,9 +66,8 @@ Foam::Rosenbrock21::Rosenbrock21(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::Rosenbrock21::solve
+Foam::scalar Foam::Rosenbrock12::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -76,7 +75,7 @@ Foam::scalar Foam::Rosenbrock21::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -104,7 +103,7 @@ Foam::scalar Foam::Rosenbrock21::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    ode.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, dydx_);
 
     forAll(k2_, i)
     {
@@ -124,15 +123,14 @@ Foam::scalar Foam::Rosenbrock21::solve
 }
 
 
-void Foam::Rosenbrock21::solve
+void Foam::Rosenbrock12::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::Rosenbrock21
+    Foam::Rosenbrock12
 
 Description
     L-stable embedded Rosenbrock ODE solver of order (1)2.
@@ -39,12 +39,12 @@ Description
     \endverbatim
 
 SourceFiles
-    Rosenbrock21.C
+    Rosenbrock12.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Rosenbrock21_H
+#ifndef Rosenbrock12_H
-#define Rosenbrock21_H
+#define Rosenbrock12_H
 
 #include "ODESolver.H"
 #include "adaptiveSolver.H"
@@ -55,10 +55,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class Rosenbrock21 Declaration
+                           Class Rosenbrock12 Declaration
 \*---------------------------------------------------------------------------*/
 
-class Rosenbrock21
+class Rosenbrock12
 :
     public ODESolver,
     public adaptiveSolver
@@ -87,13 +87,13 @@ class Rosenbrock21
 public:
 
     //- Runtime type information
-    TypeName("Rosenbrock21");
+    TypeName("Rosenbrock12");
 
 
     // Constructors
 
         //- Construct from ODE
-        Rosenbrock21(const ODESystem& ode, const dictionary& dict);
+        Rosenbrock12(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -101,7 +101,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -112,7 +111,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C

@@ -23,45 +23,45 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "Rosenbrock32.H"
+#include "Rosenbrock23.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(Rosenbrock32, 0);
+    defineTypeNameAndDebug(Rosenbrock23, 0);
-    addToRunTimeSelectionTable(ODESolver, Rosenbrock32, dictionary);
+    addToRunTimeSelectionTable(ODESolver, Rosenbrock23, dictionary);
 
 const scalar
-    Rosenbrock32::a21 = 1,
+    Rosenbrock23::a21 = 1,
-    Rosenbrock32::a31 = 1,
+    Rosenbrock23::a31 = 1,
-    Rosenbrock32::a32 = 0,
+    Rosenbrock23::a32 = 0,
 
-    Rosenbrock32::c21 = -1.0156171083877702091975600115545,
+    Rosenbrock23::c21 = -1.0156171083877702091975600115545,
-    Rosenbrock32::c31 = 4.0759956452537699824805835358067,
+    Rosenbrock23::c31 = 4.0759956452537699824805835358067,
-    Rosenbrock32::c32 = 9.2076794298330791242156818474003,
+    Rosenbrock23::c32 = 9.2076794298330791242156818474003,
 
-    Rosenbrock32::b1 = 1,
+    Rosenbrock23::b1 = 1,
-    Rosenbrock32::b2 = 6.1697947043828245592553615689730,
+    Rosenbrock23::b2 = 6.1697947043828245592553615689730,
-    Rosenbrock32::b3 = -0.4277225654321857332623837380651,
+    Rosenbrock23::b3 = -0.4277225654321857332623837380651,
 
-    Rosenbrock32::e1 = 0.5,
+    Rosenbrock23::e1 = 0.5,
-    Rosenbrock32::e2 = -2.9079558716805469821718236208017,
+    Rosenbrock23::e2 = -2.9079558716805469821718236208017,
-    Rosenbrock32::e3 = 0.2235406989781156962736090927619,
+    Rosenbrock23::e3 = 0.2235406989781156962736090927619,
 
-    Rosenbrock32::gamma = 0.43586652150845899941601945119356,
+    Rosenbrock23::gamma = 0.43586652150845899941601945119356,
-    Rosenbrock32::c2 = 0.43586652150845899941601945119356,
+    Rosenbrock23::c2 = 0.43586652150845899941601945119356,
 
-    Rosenbrock32::d1 = 0.43586652150845899941601945119356,
+    Rosenbrock23::d1 = 0.43586652150845899941601945119356,
-    Rosenbrock32::d2 = 0.24291996454816804366592249683314,
+    Rosenbrock23::d2 = 0.24291996454816804366592249683314,
-    Rosenbrock32::d3 = 2.1851380027664058511513169485832;
+    Rosenbrock23::d3 = 2.1851380027664058511513169485832;
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::Rosenbrock32::Rosenbrock32(const ODESystem& ode, const dictionary& dict)
+Foam::Rosenbrock23::Rosenbrock23(const ODESystem& ode, const dictionary& dict)
 :
     ODESolver(ode, dict),
     adaptiveSolver(ode, dict),
@@ -79,9 +79,8 @@ Foam::Rosenbrock32::Rosenbrock32(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::Rosenbrock32::solve
+Foam::scalar Foam::Rosenbrock23::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -89,7 +88,7 @@ Foam::scalar Foam::Rosenbrock32::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -117,7 +116,7 @@ Foam::scalar Foam::Rosenbrock32::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    ode.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, dydx_);
 
     forAll(k2_, i)
     {
@@ -146,15 +145,14 @@ Foam::scalar Foam::Rosenbrock32::solve
 }
 
 
-void Foam::Rosenbrock32::solve
+void Foam::Rosenbrock23::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::Rosenbrock32
+    Foam::Rosenbrock23
 
 Description
     L-stable embedded Rosenbrock ODE solver of order (3)4.
@@ -42,12 +42,12 @@ Description
     \endverbatim
 
 SourceFiles
-    Rosenbrock32.C
+    Rosenbrock23.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Rosenbrock32_H
+#ifndef Rosenbrock23_H
-#define Rosenbrock32_H
+#define Rosenbrock23_H
 
 #include "ODESolver.H"
 #include "adaptiveSolver.H"
@@ -58,10 +58,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class Rosenbrock32 Declaration
+                           Class Rosenbrock23 Declaration
 \*---------------------------------------------------------------------------*/
 
-class Rosenbrock32
+class Rosenbrock23
 :
     public ODESolver,
     public adaptiveSolver
@@ -91,13 +91,13 @@ class Rosenbrock32
 public:
 
     //- Runtime type information
-    TypeName("Rosenbrock32");
+    TypeName("Rosenbrock23");
 
 
     // Constructors
 
         //- Construct from ODE
-        Rosenbrock32(const ODESystem& ode, const dictionary& dict);
+        Rosenbrock23(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -105,7 +105,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -116,7 +115,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C

@@ -23,89 +23,89 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "Rosenbrock43.H"
+#include "Rosenbrock34.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(Rosenbrock43, 0);
+    defineTypeNameAndDebug(Rosenbrock34, 0);
-    addToRunTimeSelectionTable(ODESolver, Rosenbrock43, dictionary);
+    addToRunTimeSelectionTable(ODESolver, Rosenbrock34, dictionary);
 
 const scalar
-    // L-Stable constants from Hairer et. al.
-    Rosenbrock43::a21 = 2,
-    Rosenbrock43::a31 = 1.867943637803922,
-    Rosenbrock43::a32 = 0.2344449711399156,
-
-    Rosenbrock43::c21 = -7.137615036412310,
-    Rosenbrock43::c31 = 2.580708087951457,
-    Rosenbrock43::c32 = 0.6515950076447975,
-    Rosenbrock43::c41 = -2.137148994382534,
-    Rosenbrock43::c42 = -0.3214669691237626,
-    Rosenbrock43::c43 = -0.6949742501781779,
-
-    Rosenbrock43::b1 = 2.255570073418735,
-    Rosenbrock43::b2 = 0.2870493262186792,
-    Rosenbrock43::b3 = 0.435317943184018,
-    Rosenbrock43::b4 = 1.093502252409163,
-
-    Rosenbrock43::e1 = -0.2815431932141155,
-    Rosenbrock43::e2 = -0.0727619912493892,
-    Rosenbrock43::e3 = -0.1082196201495311,
-    Rosenbrock43::e4 = -1.093502252409163,
-
-    Rosenbrock43::gamma = 0.57282,
-    Rosenbrock43::c2 = 1.14564,
-    Rosenbrock43::c3 = 0.65521686381559,
-
-    Rosenbrock43::d1 = 0.57282,
-    Rosenbrock43::d2 = -1.769193891319233,
-    Rosenbrock43::d3 = 0.7592633437920482,
-    Rosenbrock43::d4 = -0.1049021087100450;
-
     // Constants by Shampine
     // More accurate than the L-Stable coefficients for small step-size
     // but less stable for large step-size
+    Rosenbrock34::a21 = 2,
+    Rosenbrock34::a31 = 48.0/25.0,
+    Rosenbrock34::a32 = 6.0/25.0,
+
+    Rosenbrock34::c21 = -8,
+    Rosenbrock34::c31 = 372.0/25.0,
+    Rosenbrock34::c32 = 12.0/5.0,
+
+    Rosenbrock34::c41 = -112.0/125.0,
+    Rosenbrock34::c42 = -54.0/125.0,
+    Rosenbrock34::c43 = -2.0/5.0,
+
+    Rosenbrock34::b1 = 19.0/9.0,
+    Rosenbrock34::b2 = 1.0/2.0,
+    Rosenbrock34::b3 = 25.0/108.0,
+    Rosenbrock34::b4 = 125.0/108.0,
+
+    Rosenbrock34::e1 = 34.0/108.0,
+    Rosenbrock34::e2 = 7.0/36.0,
+    Rosenbrock34::e3 = 0,
+    Rosenbrock34::e4 = 125.0/108.0,
+
+    Rosenbrock34::gamma = 1.0/2.0,
+    Rosenbrock34::c2 = 1,
+    Rosenbrock34::c3  = 3.0/5.0,
+
+    Rosenbrock34::d1 = 1.0/2.0,
+    Rosenbrock34::d2 = -3.0/2.0,
+    Rosenbrock34::d3 = 605.0/250.0,
+    Rosenbrock34::d4 = 29.0/250.0;
+
     /*
-    Rosenbrock43::a21 = 2,
+    // L-Stable constants from Hairer et. al.
-    Rosenbrock43::a31 = 48.0/25.0,
+    Rosenbrock34::a21 = 2,
-    Rosenbrock43::a32 = 6.0/25.0,
+    Rosenbrock34::a31 = 1.867943637803922,
+    Rosenbrock34::a32 = 0.2344449711399156,
 
-    Rosenbrock43::c21 = -8,
+    Rosenbrock34::c21 = -7.137615036412310,
-    Rosenbrock43::c31 = 372.0/25.0,
+    Rosenbrock34::c31 = 2.580708087951457,
-    Rosenbrock43::c32 = 12.0/5.0,
+    Rosenbrock34::c32 = 0.6515950076447975,
+    Rosenbrock34::c41 = -2.137148994382534,
+    Rosenbrock34::c42 = -0.3214669691237626,
+    Rosenbrock34::c43 = -0.6949742501781779,
 
-    Rosenbrock43::c41 = -112.0/125.0,
+    Rosenbrock34::b1 = 2.255570073418735,
-    Rosenbrock43::c42 = -54.0/125.0,
+    Rosenbrock34::b2 = 0.2870493262186792,
-    Rosenbrock43::c43 = -2.0/5.0,
+    Rosenbrock34::b3 = 0.435317943184018,
+    Rosenbrock34::b4 = 1.093502252409163,
 
-    Rosenbrock43::b1 = 19.0/9.0,
+    Rosenbrock34::e1 = -0.2815431932141155,
-    Rosenbrock43::b2 = 1.0/2.0,
+    Rosenbrock34::e2 = -0.0727619912493892,
-    Rosenbrock43::b3 = 25.0/108.0,
+    Rosenbrock34::e3 = -0.1082196201495311,
-    Rosenbrock43::b4 = 125.0/108.0,
+    Rosenbrock34::e4 = -1.093502252409163,
 
-    Rosenbrock43::e1 = 34.0/108.0,
+    Rosenbrock34::gamma = 0.57282,
-    Rosenbrock43::e2 = 7.0/36.0,
+    Rosenbrock34::c2 = 1.14564,
-    Rosenbrock43::e3 = 0,
+    Rosenbrock34::c3 = 0.65521686381559,
-    Rosenbrock43::e4 = 125.0/108.0,
 
-    Rosenbrock43::gamma = 1.0/2.0,
+    Rosenbrock34::d1 = 0.57282,
-    Rosenbrock43::c2 = 1,
+    Rosenbrock34::d2 = -1.769193891319233,
-    Rosenbrock43::c3  = 3.0/5.0,
+    Rosenbrock34::d3 = 0.7592633437920482,
-
+    Rosenbrock34::d4 = -0.1049021087100450;
-    Rosenbrock43::d1 = 1.0/2.0,
-    Rosenbrock43::d2 = -3.0/2.0,
-    Rosenbrock43::d3 = 605.0/250.0,
-    Rosenbrock43::d4 = 29.0/250.0;
     */
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::Rosenbrock43::Rosenbrock43(const ODESystem& ode, const dictionary& dict)
+Foam::Rosenbrock34::Rosenbrock34(const ODESystem& ode, const dictionary& dict)
 :
     ODESolver(ode, dict),
     adaptiveSolver(ode, dict),
@@ -124,9 +124,8 @@ Foam::Rosenbrock43::Rosenbrock43(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::Rosenbrock43::solve
+Foam::scalar Foam::Rosenbrock34::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -134,7 +133,7 @@ Foam::scalar Foam::Rosenbrock43::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -162,7 +161,7 @@ Foam::scalar Foam::Rosenbrock43::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    ode.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, dydx_);
 
     forAll(k2_, i)
     {
@@ -177,7 +176,7 @@ Foam::scalar Foam::Rosenbrock43::solve
         y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
     }
 
-    ode.derivatives(x0 + c3*dx, y, dydx_);
+    odes_.derivatives(x0 + c3*dx, y, dydx_);
 
     forAll(k3_, i)
     {
@@ -206,15 +205,14 @@ Foam::scalar Foam::Rosenbrock43::solve
 }
 
 
-void Foam::Rosenbrock43::solve
+void Foam::Rosenbrock34::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.H

@@ -22,12 +22,11 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::Rosenbrock43
+    Foam::Rosenbrock34
 
 Description
     L-stable embedded Rosenbrock ODE solver of order (3)4.
 
-    Based on code from:
     \verbatim
         "Solving Ordinary Differential Equations II: Stiff
          and Differential-Algebraic Problems, second edition",
@@ -37,7 +36,7 @@ Description
         Springer-Verlag, Berlin. 1996.
     \endverbatim
 
-    Also provided are the constants from:
+    The default constants are from:
     \verbatim
         "Implementation of Rosenbrock Methods"
         Shampine, L. F.,
@@ -46,13 +45,15 @@ Description
     with which the scheme is more accurate than with the L-Stable coefficients
     for small step-size but less stable for large step-size.
 
+    The L-Stable scheme constants are provided commented-out in Rosenbrock34.C
+
 SourceFiles
-    Rosenbrock43.C
+    Rosenbrock34.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Rosenbrock43_H
+#ifndef Rosenbrock34_H
-#define Rosenbrock43_H
+#define Rosenbrock34_H
 
 #include "ODESolver.H"
 #include "adaptiveSolver.H"
@@ -63,10 +64,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class Rosenbrock43 Declaration
+                           Class Rosenbrock34 Declaration
 \*---------------------------------------------------------------------------*/
 
-class Rosenbrock43
+class Rosenbrock34
 :
     public ODESolver,
     public adaptiveSolver
@@ -98,13 +99,13 @@ class Rosenbrock43
 public:
 
     //- Runtime type information
-    TypeName("Rosenbrock43");
+    TypeName("Rosenbrock34");
 
 
     // Constructors
 
         //- Construct from ODE
-        Rosenbrock43(const ODESystem& ode, const dictionary& dict);
+        Rosenbrock34(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -112,7 +113,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -123,7 +123,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/SIBS/SIBS.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -70,13 +70,12 @@ Foam::SIBS::SIBS(const ODESystem& ode, const dictionary& dict)
 
 void Foam::SIBS::solve
 (
-    const ODESystem& ode,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    ode.derivatives(x, y, dydx0_);
+    odes_.derivatives(x, y, dydx0_);
 
     scalar h = dxTry;
     bool exitflag = false;
@@ -124,7 +123,7 @@ void Foam::SIBS::solve
     label k = 0;
     yTemp_ = y;
 
-    ode.jacobian(x, y, dfdx_, dfdy_);
+    odes_.jacobian(x, y, dfdx_, dfdy_);
 
     if (x != xNew_ || h != dxTry)
     {
@@ -151,7 +150,7 @@ void Foam::SIBS::solve
                     << exit(FatalError);
             }
 
-            SIMPR(ode, x, yTemp_, dydx0_, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
+            SIMPR(x, yTemp_, dydx0_, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
             scalar xest = sqr(h/nSeq_[k]);
 
             polyExtrapolate(k, xest, ySeq_, y, yErr_, x_p_, d_p_);

src/ODE/ODESolvers/SIBS/SIBS.H

@@ -84,7 +84,6 @@ class SIBS
 
         void SIMPR
         (
-            const ODESystem& ode,
             const scalar xStart,
             const scalarField& y,
             const scalarField& dydx,
@@ -123,7 +122,6 @@ public:
 
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/SIBS/SIMPR.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,7 +29,6 @@ License
 
 void Foam::SIBS::SIMPR
 (
-    const ODESystem& ode,
     const scalar xStart,
     const scalarField& y,
     const scalarField& dydx,
@@ -72,7 +71,7 @@ void Foam::SIBS::SIMPR
 
     scalar x = xStart + h;
 
-    ode.derivatives(x, ytemp, yEnd);
+    odes_.derivatives(x, ytemp, yEnd);
 
     for (register label nn=2; nn<=nSteps; nn++)
     {
@@ -90,7 +89,7 @@ void Foam::SIBS::SIMPR
 
         x += h;
 
-        ode.derivatives(x, ytemp, yEnd);
+        odes_.derivatives(x, ytemp, yEnd);
     }
     for (register label i=0; i<n_; i++)
     {

src/ODE/ODESolvers/Trapezoid/Trapezoid.C

@@ -49,7 +49,6 @@ Foam::Trapezoid::Trapezoid(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::Trapezoid::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -64,7 +63,7 @@ Foam::scalar Foam::Trapezoid::solve
     }
 
     // Evaluate the system for the predicted state
-    ode.derivatives(x0 + dx, y, err_);
+    odes_.derivatives(x0 + dx, y, err_);
 
     // Update the state as the average between the prediction and the correction
     // and estimate the error from the difference
@@ -80,13 +79,12 @@ Foam::scalar Foam::Trapezoid::solve
 
 void Foam::Trapezoid::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/Trapezoid/Trapezoid.H

@@ -74,7 +74,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -85,7 +84,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C

@@ -64,7 +64,7 @@ void Foam::adaptiveSolver::solve
     do
     {
         // Solve step and provide error estimate
-        err = solve(odes, x, y, dydx0_, dx, yTemp_);
+        err = solve(x, y, dydx0_, dx, yTemp_);
 
         // If error is large reduce dx
         if (err > 1)

src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H

@@ -77,7 +77,6 @@ public:
         //- Solve a single step dx and return the error
         virtual scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,

src/ODE/ODESolvers/rodas23/rodas23.C

@@ -23,35 +23,35 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "rodas32.H"
+#include "rodas23.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(rodas32, 0);
+    defineTypeNameAndDebug(rodas23, 0);
-    addToRunTimeSelectionTable(ODESolver, rodas32, dictionary);
+    addToRunTimeSelectionTable(ODESolver, rodas23, dictionary);
 
 const scalar
-    rodas32::c3 = 1,
+    rodas23::c3 = 1,
-    rodas32::d1 = 1.0/2.0,
+    rodas23::d1 = 1.0/2.0,
-    rodas32::d2 = 3.0/2.0,
+    rodas23::d2 = 3.0/2.0,
-    rodas32::a31 = 2,
+    rodas23::a31 = 2,
-    rodas32::a41 = 2,
+    rodas23::a41 = 2,
-    rodas32::c21 = 4,
+    rodas23::c21 = 4,
-    rodas32::c31 = 1,
+    rodas23::c31 = 1,
-    rodas32::c32 = -1,
+    rodas23::c32 = -1,
-    rodas32::c41 = 1,
+    rodas23::c41 = 1,
-    rodas32::c42 = -1,
+    rodas23::c42 = -1,
-    rodas32::c43 = -8.0/3.0,
+    rodas23::c43 = -8.0/3.0,
-    rodas32::gamma = 1.0/2.0;
+    rodas23::gamma = 1.0/2.0;
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::rodas32::rodas32(const ODESystem& ode, const dictionary& dict)
+Foam::rodas23::rodas23(const ODESystem& ode, const dictionary& dict)
 :
     ODESolver(ode, dict),
     adaptiveSolver(ode, dict),
@@ -70,9 +70,8 @@ Foam::rodas32::rodas32(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::rodas32::solve
+Foam::scalar Foam::rodas23::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -80,7 +79,7 @@ Foam::scalar Foam::rodas32::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -117,7 +116,7 @@ Foam::scalar Foam::rodas32::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    ode.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, dydx_);
 
     forAll(k3_, i)
     {
@@ -133,7 +132,7 @@ Foam::scalar Foam::rodas32::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    ode.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, dydx_);
 
     forAll(err_, i)
     {
@@ -151,15 +150,14 @@ Foam::scalar Foam::rodas32::solve
 }
 
 
-void Foam::rodas32::solve
+void Foam::rodas23::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/rodas23/rodas23.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::rodas32
+    Foam::rodas23
 
 Description
     L-stable, stiffly-accurate embedded Rosenbrock ODE solver of order (2)3.
@@ -42,12 +42,12 @@ Description
     \endverbatim
 
 SourceFiles
-    rodas32.C
+    rodas23.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef rodas32_H
+#ifndef rodas23_H
-#define rodas32_H
+#define rodas23_H
 
 #include "ODESolver.H"
 #include "adaptiveSolver.H"
@@ -58,10 +58,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class rodas32 Declaration
+                           Class rodas23 Declaration
 \*---------------------------------------------------------------------------*/
 
-class rodas32
+class rodas23
 :
     public ODESolver,
     public adaptiveSolver
@@ -91,13 +91,13 @@ class rodas32
 public:
 
     //- Runtime type information
-    TypeName("rodas32");
+    TypeName("rodas23");
 
 
     // Constructors
 
         //- Construct from ODE
-        rodas32(const ODESystem& ode, const dictionary& dict);
+        rodas23(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -105,7 +105,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -116,7 +115,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/rodas34/rodas34.C

@@ -23,56 +23,56 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "rodas43.H"
+#include "rodas34.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(rodas43, 0);
+    defineTypeNameAndDebug(rodas34, 0);
-    addToRunTimeSelectionTable(ODESolver, rodas43, dictionary);
+    addToRunTimeSelectionTable(ODESolver, rodas34, dictionary);
 
 const scalar
-    rodas43::c2 = 0.386,
+    rodas34::c2 = 0.386,
-    rodas43::c3 = 0.21,
+    rodas34::c3 = 0.21,
-    rodas43::c4 = 0.63,
+    rodas34::c4 = 0.63,
-    rodas43::d1 =  0.25,
+    rodas34::d1 =  0.25,
-    rodas43::d2 = -0.1043,
+    rodas34::d2 = -0.1043,
-    rodas43::d3 =  0.1035,
+    rodas34::d3 =  0.1035,
-    rodas43::d4 = -0.3620000000000023e-01,
+    rodas34::d4 = -0.3620000000000023e-01,
-    rodas43::a21 =  0.1544e1,
+    rodas34::a21 =  0.1544e1,
-    rodas43::a31 =  0.9466785280815826,
+    rodas34::a31 =  0.9466785280815826,
-    rodas43::a32 =  0.2557011698983284,
+    rodas34::a32 =  0.2557011698983284,
-    rodas43::a41 =  0.3314825187068521e1,
+    rodas34::a41 =  0.3314825187068521e1,
-    rodas43::a42 =  0.2896124015972201e1,
+    rodas34::a42 =  0.2896124015972201e1,
-    rodas43::a43 =  0.9986419139977817,
+    rodas34::a43 =  0.9986419139977817,
-    rodas43::a51 =  0.1221224509226641e1,
+    rodas34::a51 =  0.1221224509226641e1,
-    rodas43::a52 =  0.6019134481288629e1,
+    rodas34::a52 =  0.6019134481288629e1,
-    rodas43::a53 =  0.1253708332932087e2,
+    rodas34::a53 =  0.1253708332932087e2,
-    rodas43::a54 = -0.6878860361058950,
+    rodas34::a54 = -0.6878860361058950,
-    rodas43::c21 = -0.56688e1,
+    rodas34::c21 = -0.56688e1,
-    rodas43::c31 = -0.2430093356833875e1,
+    rodas34::c31 = -0.2430093356833875e1,
-    rodas43::c32 = -0.2063599157091915,
+    rodas34::c32 = -0.2063599157091915,
-    rodas43::c41 = -0.1073529058151375,
+    rodas34::c41 = -0.1073529058151375,
-    rodas43::c42 = -0.9594562251023355e1,
+    rodas34::c42 = -0.9594562251023355e1,
-    rodas43::c43 = -0.2047028614809616e2,
+    rodas34::c43 = -0.2047028614809616e2,
-    rodas43::c51 =  0.7496443313967647e1,
+    rodas34::c51 =  0.7496443313967647e1,
-    rodas43::c52 = -0.1024680431464352e2,
+    rodas34::c52 = -0.1024680431464352e2,
-    rodas43::c53 = -0.3399990352819905e2,
+    rodas34::c53 = -0.3399990352819905e2,
-    rodas43::c54 =  0.1170890893206160e2,
+    rodas34::c54 =  0.1170890893206160e2,
-    rodas43::c61 =  0.8083246795921522e1,
+    rodas34::c61 =  0.8083246795921522e1,
-    rodas43::c62 = -0.7981132988064893e1,
+    rodas34::c62 = -0.7981132988064893e1,
-    rodas43::c63 = -0.3152159432874371e2,
+    rodas34::c63 = -0.3152159432874371e2,
-    rodas43::c64 =  0.1631930543123136e2,
+    rodas34::c64 =  0.1631930543123136e2,
-    rodas43::c65 = -0.6058818238834054e1,
+    rodas34::c65 = -0.6058818238834054e1,
-    rodas43::gamma = 0.25;
+    rodas34::gamma = 0.25;
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::rodas43::rodas43(const ODESystem& ode, const dictionary& dict)
+Foam::rodas34::rodas34(const ODESystem& ode, const dictionary& dict)
 :
     ODESolver(ode, dict),
     adaptiveSolver(ode, dict),
@@ -93,9 +93,8 @@ Foam::rodas43::rodas43(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::rodas43::solve
+Foam::scalar Foam::rodas34::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -103,7 +102,7 @@ Foam::scalar Foam::rodas43::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -131,7 +130,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    ode.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, dydx_);
 
     forAll(k2_, i)
     {
@@ -146,7 +145,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
     }
 
-    ode.derivatives(x0 + c3*dx, y, dydx_);
+    odes_.derivatives(x0 + c3*dx, y, dydx_);
 
     forAll(k3_, i)
     {
@@ -161,7 +160,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + a41*k1_[i] + a42*k2_[i] + a43*k3_[i];
     }
 
-    ode.derivatives(x0 + c4*dx, y, dydx_);
+    odes_.derivatives(x0 + c4*dx, y, dydx_);
 
     forAll(k4_, i)
     {
@@ -178,7 +177,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    ode.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, dydx_);
 
     forAll(k5_, i)
     {
@@ -195,7 +194,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    ode.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, dydx_);
 
     forAll(err_, i)
     {
@@ -214,15 +213,14 @@ Foam::scalar Foam::rodas43::solve
 }
 
 
-void Foam::rodas43::solve
+void Foam::rodas34::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/rodas34/rodas34.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::rodas43
+    Foam::rodas34
 
 Description
     L-stable, stiffly-accurate embedded Rosenbrock ODE solver of order (3)4.
@@ -38,12 +38,12 @@ Description
     \endverbatim
 
 SourceFiles
-    rodas43.C
+    rodas34.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef rodas43_H
+#ifndef rodas34_H
-#define rodas43_H
+#define rodas34_H
 
 #include "ODESolver.H"
 #include "adaptiveSolver.H"
@@ -54,10 +54,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class rodas43 Declaration
+                           Class rodas34 Declaration
 \*---------------------------------------------------------------------------*/
 
-class rodas43
+class rodas34
 :
     public ODESolver,
     public adaptiveSolver
@@ -92,13 +92,13 @@ class rodas43
 public:
 
     //- Runtime type information
-    TypeName("rodas43");
+    TypeName("rodas34");
 
 
     // Constructors
 
         //- Construct from ODE
-        rodas43(const ODESystem& ode, const dictionary& dict);
+        rodas34(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -106,7 +106,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -117,7 +116,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/seulex/seulex.C

@@ -0,0 +1,463 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "seulex.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(seulex, 0);
+
+    addToRunTimeSelectionTable(ODESolver, seulex, dictionary);
+
+    const scalar
+        seulex::stepFactor1_ = 0.6,
+        seulex::stepFactor2_ = 0.93,
+        seulex::stepFactor3_ = 0.1,
+        seulex::stepFactor4_ = 4.0,
+        seulex::stepFactor5_ = 0.5,
+        seulex::kFactor1_ = 0.7,
+        seulex::kFactor2_ = 0.9;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict)
+:
+    ODESolver(ode, dict),
+    jacRedo_(min(1.0e-4, min(relTol_))),
+    nSeq_(iMaxx_),
+    cpu_(iMaxx_),
+    coeff_(iMaxx_, iMaxx_),
+    theta_(2.0*jacRedo_),
+    table_(kMaxx_, n_),
+    dfdx_(n_),
+    dfdy_(n_),
+    a_(n_),
+    pivotIndices_(n_),
+    dxOpt_(iMaxx_),
+    temp_(iMaxx_),
+    y0_(n_),
+    ySequence_(n_),
+    scale_(n_),
+    dy_(n_),
+    yTemp_(n_),
+    dydx_(n_)
+{
+    // The CPU time factors for the major parts of the algorithm
+    const scalar cpuFunc = 1.0, cpuJac = 5.0, cpuLU = 1.0, cpuSolve = 1.0;
+
+    nSeq_[0] = 2;
+    nSeq_[1] = 3;
+
+    for (int i=2; i<iMaxx_; i++)
+    {
+        nSeq_[i] = 2*nSeq_[i-2];
+    }
+    cpu_[0] = cpuJac + cpuLU + nSeq_[0]*(cpuFunc + cpuSolve);
+
+    for (int k=0; k<kMaxx_; k++)
+    {
+        cpu_[k+1] = cpu_[k] + (nSeq_[k+1]-1)*(cpuFunc + cpuSolve) + cpuLU;
+    }
+
+    // Set the extrapolation coefficients array
+    for (int k=0; k<iMaxx_; k++)
+    {
+        for (int l=0; l<k; l++)
+        {
+            scalar ratio = scalar(nSeq_[k])/nSeq_[l];
+            coeff_[k][l] = 1.0/(ratio - 1.0);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+bool Foam::seulex::seul
+(
+    const scalar x0,
+    const scalarField& y0,
+    const scalar dxTot,
+    const label k,
+    scalarField& y,
+    const scalarField& scale
+) const
+{
+    label nSteps = nSeq_[k];
+    scalar dx = dxTot/nSteps;
+
+    for (label i=0; i<n_; i++)
+    {
+        for (label j=0; j<n_; j++)
+        {
+            a_[i][j] = -dfdy_[i][j];
+        }
+
+        a_[i][i] += 1.0/dx;
+    }
+
+    LUDecompose(a_, pivotIndices_);
+
+    scalar xnew = x0 + dx;
+    odes_.derivatives(xnew, y0, dy_);
+    LUBacksubstitute(a_, pivotIndices_, dy_);
+
+    yTemp_ = y0;
+
+    for (label nn=1; nn<nSteps; nn++)
+    {
+        yTemp_ += dy_;
+        xnew += dx;
+
+        if (nn == 1 && k<=1)
+        {
+            scalar dy1 = 0.0;
+            for (label i=0; i<n_; i++)
+            {
+                dy1 += sqr(dy_[i]/scale[i]);
+            }
+            dy1 = sqrt(dy1);
+
+            odes_.derivatives(x0 + dx, yTemp_, dydx_);
+            for (label i=0; i<n_; i++)
+            {
+                dy_[i] = dydx_[i] - dy_[i]/dx;
+            }
+
+            LUBacksubstitute(a_, pivotIndices_, dy_);
+
+            scalar dy2 = 0.0;
+            for (label i=0; i<n_; i++)
+            {
+                dy2 += sqr(dy_[i]/scale[i]);
+            }
+            dy2 = sqrt(dy2);
+            theta_ = dy2/min(1.0, dy1 + SMALL);
+
+            if (theta_ > 1.0)
+            {
+                return false;
+            }
+        }
+
+        odes_.derivatives(xnew, yTemp_, dy_);
+        LUBacksubstitute(a_, pivotIndices_, dy_);
+    }
+
+    for (label i=0; i<n_; i++)
+    {
+        y[i] = yTemp_[i] + dy_[i];
+    }
+
+    return true;
+}
+
+
+void Foam::seulex::extrapolate
+(
+    const label k,
+    scalarRectangularMatrix& table,
+    scalarField& y
+) const
+{
+    for (int j=k-1; j>0; j--)
+    {
+        for (label i=0; i<n_; i++)
+        {
+            table[j-1][i] =
+                table[j][i] + coeff_[k][j]*(table[j][i] - table[j-1][i]);
+        }
+    }
+
+    for (int i=0; i<n_; i++)
+    {
+        y[i] = table[0][i] + coeff_[k][0]*(table[0][i] - y[i]);
+    }
+}
+
+
+void Foam::seulex::solve
+(
+    scalar& x,
+    scalarField& y,
+    stepState& step
+) const
+{
+    temp_[0] = GREAT;
+    scalar dx = step.dxTry;
+    y0_ = y;
+    dxOpt_[0] = mag(0.1*dx);
+
+    if (step.first || step.prevReject)
+    {
+        theta_ = 2.0*jacRedo_;
+    }
+
+    if (step.first)
+    {
+        // NOTE: the first element of relTol_ and absTol_ are used here.
+        scalar logTol = -log10(relTol_[0] + absTol_[0])*0.6 + 0.5;
+        kTarg_ = max(1, min(kMaxx_ - 1, int(logTol)));
+    }
+
+    forAll(scale_, i)
+    {
+        scale_[i] = absTol_[i] + relTol_[i]*mag(y[i]);
+    }
+
+    bool jacUpdated = false;
+
+    if (theta_ > jacRedo_)
+    {
+        odes_.jacobian(x, y, dfdx_, dfdy_);
+        jacUpdated = true;
+    }
+
+    int k;
+    scalar dxNew = mag(dx);
+    bool firstk = true;
+
+    while (firstk || step.reject)
+    {
+        dx = step.forward ? dxNew : -dxNew;
+        firstk = false;
+        step.reject = false;
+
+        if (mag(dx) <= mag(x)*SMALL)
+        {
+             WarningIn("seulex::solve(scalar& x, scalarField& y, stepState&")
+                    << "step size underflow :"  << dx << endl;
+        }
+
+        scalar errOld;
+
+        for (k=0; k<=kTarg_+1; k++)
+        {
+            bool success = seul(x, y0_, dx, k, ySequence_, scale_);
+
+            if (!success)
+            {
+                step.reject = true;
+                dxNew = mag(dx)*stepFactor5_;
+                break;
+            }
+
+            if (k == 0)
+            {
+                 y = ySequence_;
+            }
+            else
+            {
+                forAll(ySequence_, i)
+                {
+                    table_[k-1][i] = ySequence_[i];
+                }
+            }
+
+            if (k != 0)
+            {
+                extrapolate(k, table_, y);
+                scalar err = 0.0;
+                forAll(scale_, i)
+                {
+                    scale_[i] = absTol_[i] + relTol_[i]*mag(y0_[i]);
+                    err += sqr((y[i] - table_[0][i])/scale_[i]);
+                }
+                err = sqrt(err/n_);
+                if (err > 1.0/SMALL || (k > 1 && err >= errOld))
+                {
+                    step.reject = true;
+                    dxNew = mag(dx)*stepFactor5_;
+                    break;
+                }
+                errOld = min(4*err, 1);
+                scalar expo = 1.0/(k + 1);
+                scalar facmin = pow(stepFactor3_, expo);
+                scalar fac;
+                if (err == 0.0)
+                {
+                    fac = 1.0/facmin;
+                }
+                else
+                {
+                    fac = stepFactor2_/pow(err/stepFactor1_, expo);
+                    fac = max(facmin/stepFactor4_, min(1.0/facmin, fac));
+                }
+                dxOpt_[k] = mag(dx*fac);
+                temp_[k] = cpu_[k]/dxOpt_[k];
+
+                if ((step.first || step.last) && err <= 1.0)
+                {
+                    break;
+                }
+
+                if
+                (
+                    k == kTarg_ - 1
+                 && !step.prevReject
+                 && !step.first && !step.last
+                )
+                {
+                    if (err <= 1.0)
+                    {
+                        break;
+                    }
+                    else if (err > nSeq_[kTarg_]*nSeq_[kTarg_ + 1]*4.0)
+                    {
+                        step.reject = true;
+                        kTarg_ = k;
+                        if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
+                        {
+                            kTarg_--;
+                        }
+                        dxNew = dxOpt_[kTarg_];
+                        break;
+                    }
+                }
+
+                if (k == kTarg_)
+                {
+                    if (err <= 1.0)
+                    {
+                        break;
+                    }
+                    else if (err > nSeq_[k + 1]*2.0)
+                    {
+                        step.reject = true;
+                        if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
+                        {
+                            kTarg_--;
+                        }
+                        dxNew = dxOpt_[kTarg_];
+                        break;
+                    }
+                }
+
+                if (k == kTarg_+1)
+                {
+                    if (err > 1.0)
+                    {
+                        step.reject = true;
+                        if
+                        (
+                            kTarg_ > 1
+                         && temp_[kTarg_-1] < kFactor1_*temp_[kTarg_]
+                        )
+                        {
+                            kTarg_--;
+                        }
+                        dxNew = dxOpt_[kTarg_];
+                    }
+                    break;
+                }
+            }
+        }
+        if (step.reject)
+        {
+            step.prevReject = true;
+            if (!jacUpdated)
+            {
+                theta_ = 2.0*jacRedo_;
+
+                if (theta_ > jacRedo_ && !jacUpdated)
+                {
+                    odes_.jacobian(x, y, dfdx_, dfdy_);
+                    jacUpdated = true;
+                }
+            }
+        }
+    }
+
+    jacUpdated = false;
+
+    step.dxDid = dx;
+    x += dx;
+
+    label kopt;
+    if (k == 1)
+    {
+        kopt = 2;
+    }
+    else if (k <= kTarg_)
+    {
+        kopt=k;
+        if (temp_[k-1] < kFactor1_*temp_[k])
+        {
+            kopt = k - 1;
+        }
+        else if (temp_[k] < kFactor2_*temp_[k - 1])
+        {
+            kopt = min(k + 1, kMaxx_ - 1);
+        }
+    }
+    else
+    {
+        kopt = k - 1;
+        if (k > 2 && temp_[k-2] < kFactor1_*temp_[k - 1])
+        {
+            kopt = k - 2;
+        }
+        if (temp_[k] < kFactor2_*temp_[kopt])
+        {
+            kopt = min(k, kMaxx_ - 1);
+        }
+    }
+
+    if (step.prevReject)
+    {
+        kTarg_ = min(kopt, k);
+        dxNew = min(mag(dx), dxOpt_[kTarg_]);
+        step.prevReject = false;
+    }
+    else
+    {
+        if (kopt <= k)
+        {
+            dxNew = dxOpt_[kopt];
+        }
+        else
+        {
+            if (k < kTarg_ && temp_[k] < kFactor2_*temp_[k - 1])
+            {
+                dxNew = dxOpt_[k]*cpu_[kopt + 1]/cpu_[k];
+            }
+            else
+            {
+                dxNew = dxOpt_[k]*cpu_[kopt]/cpu_[k];
+            }
+        }
+        kTarg_ = kopt;
+    }
+
+    step.dxTry = step.forward ? dxNew : -dxNew;
+}
+
+
+// ************************************************************************* //

src/ODE/ODESolvers/seulex/seulex.H

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::seulex
+
+Description
+    An extrapolation-algorithm, based on the linearly implicit Euler method
+    with step size control and order selection.
+
+    The implementation is based on the SEULEX algorithm in
+    \verbatim
+        "Solving Ordinary Differential Equations II: Stiff
+         and Differential-Algebraic Problems, second edition",
+        Hairer, E.,
+        Nørsett, S.,
+        Wanner, G.,
+        Springer-Verlag, Berlin. 1996.
+    \endverbatim
+
+SourceFiles
+    seulex.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef seulex_H
+#define seulex_H
+
+#include "ODESolver.H"
+#include "scalarMatrices.H"
+#include "labelField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class seulex Declaration
+\*---------------------------------------------------------------------------*/
+
+class seulex
+:
+    public ODESolver
+{
+    // Private data
+
+        // Static constants
+
+            static const label kMaxx_ = 12;
+            static const label iMaxx_ = kMaxx_ + 1;
+
+            static const scalar
+                stepFactor1_, stepFactor2_, stepFactor3_,
+                stepFactor4_, stepFactor5_,
+                kFactor1_, kFactor2_;
+
+        // Evaluated constants
+
+            scalar jacRedo_;
+            labelField nSeq_;
+            scalarField cpu_;
+            scalarSquareMatrix coeff_;
+
+        // Temporary storage
+        // held to avoid dynamic memory allocation between calls
+        // and to transfer internal values between functions
+
+            mutable scalar theta_;
+            mutable label kTarg_;
+            mutable scalarRectangularMatrix table_;
+
+            mutable scalarField dfdx_;
+            mutable scalarSquareMatrix dfdy_;
+            mutable scalarSquareMatrix a_;
+            mutable labelList pivotIndices_;
+
+            // Fields space for "solve" function
+            mutable scalarField dxOpt_, temp_;
+            mutable scalarField y0_, ySequence_, scale_;
+
+            // Fields used in "seul" function
+            mutable scalarField dy_, yTemp_, dydx_;
+
+
+    // Private Member Functions
+
+        //- Computes the j-th line of the extrapolation table
+        bool seul
+        (
+            const scalar x0,
+            const scalarField& y0,
+            const scalar dxTot,
+            const label k,
+            scalarField& y,
+            const scalarField& scale
+        ) const;
+
+        //- Polynomial extrpolation
+        void extrapolate
+        (
+            const label k,
+            scalarRectangularMatrix& table,
+            scalarField& y
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("seulex");
+
+
+    // Constructors
+
+        //- Construct from ODE
+        seulex(const ODESystem& ode, const dictionary& dict);
+
+
+    // Member Functions
+
+        //- Solve the ODE system and the update the state
+        void solve
+        (
+            scalar& x,
+            scalarField& y,
+            stepState& step
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/db/Time/TimeIO.C

@@ -579,8 +579,8 @@ void Foam::Time::readModifiedObjects()
 
         if (controlDict_.readIfModified())
         {
-           readDict();
+            readDict();
-           functionObjects_.read();
+            functionObjects_.read();
         }
 
         bool registryModified = objectRegistry::modified();

src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,10 +59,10 @@ Foam::basicSymmetryPointPatchField<Type>::basicSymmetryPointPatchField
     const basicSymmetryPointPatchField<Type>& ptf,
     const pointPatch& p,
     const DimensionedField<Type, pointMesh>& iF,
-    const pointPatchFieldMapper&
+    const pointPatchFieldMapper& mapper
 )
 :
-    pointPatchField<Type>(p, iF)
+    pointPatchField<Type>(ptf, p, iF, mapper)
 {}
 
 

src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,10 +68,10 @@ calculatedPointPatchField<Type>::calculatedPointPatchField
     const calculatedPointPatchField<Type>& ptf,
     const pointPatch& p,
     const DimensionedField<Type, pointMesh>& iF,
-    const pointPatchFieldMapper&
+    const pointPatchFieldMapper& mapper
 )
 :
-    pointPatchField<Type>(p, iF)
+    pointPatchField<Type>(ptf, p, iF, mapper)
 {}
 
 

src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -61,10 +61,10 @@ coupledPointPatchField<Type>::coupledPointPatchField
     const coupledPointPatchField<Type>& ptf,
     const pointPatch& p,
     const DimensionedField<Type, pointMesh>& iF,
-    const pointPatchFieldMapper&
+    const pointPatchFieldMapper& mapper
 )
 :
-    pointPatchField<Type>(p, iF)
+    pointPatchField<Type>(ptf, p, iF, mapper)
 {}
 
 

src/OpenFOAM/fields/pointPatchFields/basic/value/valuePointPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -108,7 +108,7 @@ Foam::valuePointPatchField<Type>::valuePointPatchField
     const pointPatchFieldMapper& mapper
 )
 :
-    pointPatchField<Type>(p, iF),
+    pointPatchField<Type>(ptf, p, iF, mapper),
     Field<Type>(ptf, mapper)
 {}
 

src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -61,10 +61,10 @@ zeroGradientPointPatchField<Type>::zeroGradientPointPatchField
     const zeroGradientPointPatchField<Type>& ptf,
     const pointPatch& p,
     const DimensionedField<Type, pointMesh>& iF,
-    const pointPatchFieldMapper&
+    const pointPatchFieldMapper& mapper
 )
 :
-    pointPatchField<Type>(p, iF)
+    pointPatchField<Type>(ptf, p, iF, mapper)
 {}
 
 

src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,10 +77,10 @@ emptyPointPatchField<Type>::emptyPointPatchField
     const emptyPointPatchField<Type>& ptf,
     const pointPatch& p,
     const DimensionedField<Type, pointMesh>& iF,
-    const pointPatchFieldMapper&
+    const pointPatchFieldMapper& mapper
 )
 :
-    pointPatchField<Type>(p, iF)
+    pointPatchField<Type>(ptf, p, iF, mapper)
 {
     if (!isType<emptyPointPatch>(this->patch()))
     {

src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,10 +74,10 @@ Foam::wedgePointPatchField<Type>::wedgePointPatchField
     const wedgePointPatchField<Type>& ptf,
     const pointPatch& p,
     const DimensionedField<Type, pointMesh>& iF,
-    const pointPatchFieldMapper&
+    const pointPatchFieldMapper& mapper
 )
 :
-    pointPatchField<Type>(p, iF)
+    pointPatchField<Type>(ptf, p, iF, mapper)
 {
     if (!isType<wedgePointPatch>(this->patch()))
     {

src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,6 +63,22 @@ pointPatchField<Type>::pointPatchField
 {}
 
 
+template<class Type>
+Foam::pointPatchField<Type>::pointPatchField
+(
+    const pointPatchField<Type>& ptf,
+    const pointPatch& p,
+    const DimensionedField<Type, pointMesh>& iF,
+    const pointPatchFieldMapper&
+)
+:
+    patch_(p),
+    internalField_(iF),
+    updated_(false),
+    patchType_(ptf.patchType_)
+{}
+
+
 template<class Type>
 pointPatchField<Type>::pointPatchField
 (

src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H

@@ -167,6 +167,15 @@ public:
             const dictionary&
         );
 
+        //- Construct by mapping given patchField<Type> onto a new patch
+        pointPatchField
+        (
+            const pointPatchField<Type>&,
+            const pointPatch&,
+            const DimensionedField<Type, pointMesh>&,
+            const pointPatchFieldMapper&
+        );
+
         //- Construct as copy
         pointPatchField(const pointPatchField<Type>&);
 

src/OpenFOAM/primitives/Tensor/tensor/tensor.C

@@ -26,69 +26,68 @@ License
 #include "tensor.H"
 #include "mathematicalConstants.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+using namespace Foam::constant::mathematical;
-
-namespace Foam
-{
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-template<>
+namespace Foam
-const char* const tensor::typeName = "tensor";
-
-template<>
-const char* tensor::componentNames[] =
 {
-    "xx", "xy", "xz",
+    template<>
-    "yx", "yy", "yz",
+    const char* const tensor::typeName = "tensor";
-    "zx", "zy", "zz"
-};
 
-template<>
+    template<>
-const tensor tensor::zero
+    const char* tensor::componentNames[] =
-(
+    {
-    0, 0, 0,
+        "xx", "xy", "xz",
-    0, 0, 0,
+        "yx", "yy", "yz",
-    0, 0, 0
+        "zx", "zy", "zz"
-);
+    };
 
-template<>
+    template<>
-const tensor tensor::one
+    const tensor tensor::zero
-(
+    (
-    1, 1, 1,
+        0, 0, 0,
-    1, 1, 1,
+        0, 0, 0,
-    1, 1, 1
+        0, 0, 0
-);
+    );
 
-template<>
+    template<>
-const tensor tensor::max
+    const tensor tensor::one
-(
+    (
-    VGREAT, VGREAT, VGREAT,
+        1, 1, 1,
-    VGREAT, VGREAT, VGREAT,
+        1, 1, 1,
-    VGREAT, VGREAT, VGREAT
+        1, 1, 1
-);
+    );
 
-template<>
+    template<>
-const tensor tensor::min
+    const tensor tensor::max
-(
+    (
-    -VGREAT, -VGREAT, -VGREAT,
+        VGREAT, VGREAT, VGREAT,
-    -VGREAT, -VGREAT, -VGREAT,
+        VGREAT, VGREAT, VGREAT,
-    -VGREAT, -VGREAT, -VGREAT
+        VGREAT, VGREAT, VGREAT
-);
+    );
 
-template<>
+    template<>
-const tensor tensor::I
+    const tensor tensor::min
-(
+    (
-    1, 0, 0,
+        -VGREAT, -VGREAT, -VGREAT,
-    0, 1, 0,
+        -VGREAT, -VGREAT, -VGREAT,
-    0, 0, 1
+        -VGREAT, -VGREAT, -VGREAT
-);
+    );
+
+    template<>
+    const tensor tensor::I
+    (
+        1, 0, 0,
+        0, 1, 0,
+        0, 0, 1
+    );
+}
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Return eigenvalues in ascending order of absolute values
+Foam::vector Foam::eigenValues(const tensor& t)
-vector eigenValues(const tensor& t)
 {
     scalar i = 0;
     scalar ii = 0;
@@ -120,7 +119,7 @@ vector eigenValues(const tensor& t)
             + t.xy()*t.yx()*t.zz() + t.xx()*t.yz()*t.zy();
 
         // If there is a zero root
-        if (mag(c) < 1.0e-100)
+        if (mag(c) < ROOTVSMALL)
         {
             scalar disc = sqr(a) - 4*b;
 
@@ -157,14 +156,8 @@ vector eigenValues(const tensor& t)
                 scalar aBy3 = a/3;
 
                 i = m2SqrtQ*cos(theta/3) - aBy3;
-                ii =
+                ii = m2SqrtQ*cos((theta + twoPi)/3) - aBy3;
-                    m2SqrtQ
+                iii = m2SqrtQ*cos((theta - twoPi)/3) - aBy3;
-                   *cos((theta + constant::mathematical::twoPi)/3)
-                  - aBy3;
-                iii =
-                    m2SqrtQ
-                   *cos((theta - constant::mathematical::twoPi)/3)
-                  - aBy3;
             }
             else
             {
@@ -210,7 +203,7 @@ vector eigenValues(const tensor& t)
 }
 
 
-vector eigenVector(const tensor& t, const scalar lambda)
+Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda)
 {
     if (mag(lambda) < SMALL)
     {
@@ -273,7 +266,7 @@ vector eigenVector(const tensor& t, const scalar lambda)
 }
 
 
-tensor eigenVectors(const tensor& t)
+Foam::tensor Foam::eigenVectors(const tensor& t)
 {
     vector evals(eigenValues(t));
 
@@ -289,7 +282,7 @@ tensor eigenVectors(const tensor& t)
 
 
 // Return eigenvalues in ascending order of absolute values
-vector eigenValues(const symmTensor& t)
+Foam::vector Foam::eigenValues(const symmTensor& t)
 {
     scalar i = 0;
     scalar ii = 0;
@@ -321,7 +314,7 @@ vector eigenValues(const symmTensor& t)
             + t.xy()*t.xy()*t.zz() + t.xx()*t.yz()*t.yz();
 
         // If there is a zero root
-        if (mag(c) < 1.0e-100)
+        if (mag(c) < ROOTVSMALL)
         {
             scalar disc = sqr(a) - 4*b;
 
@@ -358,14 +351,8 @@ vector eigenValues(const symmTensor& t)
                 scalar aBy3 = a/3;
 
                 i = m2SqrtQ*cos(theta/3) - aBy3;
-                ii =
+                ii = m2SqrtQ*cos((theta + twoPi)/3) - aBy3;
-                    m2SqrtQ
+                iii = m2SqrtQ*cos((theta - twoPi)/3) - aBy3;
-                   *cos((theta + constant::mathematical::twoPi)/3)
-                  - aBy3;
-                iii =
-                    m2SqrtQ
-                   *cos((theta - constant::mathematical::twoPi)/3)
-                  - aBy3;
             }
             else
             {
@@ -411,7 +398,7 @@ vector eigenValues(const symmTensor& t)
 }
 
 
-vector eigenVector(const symmTensor& t, const scalar lambda)
+Foam::vector Foam::eigenVector(const symmTensor& t, const scalar lambda)
 {
     if (mag(lambda) < SMALL)
     {
@@ -474,7 +461,7 @@ vector eigenVector(const symmTensor& t, const scalar lambda)
 }
 
 
-tensor eigenVectors(const symmTensor& t)
+Foam::tensor Foam::eigenVectors(const symmTensor& t)
 {
     vector evals(eigenValues(t));
 
@@ -489,8 +476,4 @@ tensor eigenVectors(const symmTensor& t)
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/OpenFOAM/primitives/demandDrivenEntry/demandDrivenEntry.C

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "demandDrivenEntry.H"
+
+// * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::demandDrivenEntry<Type>::demandDrivenEntry
+(
+    const dictionary& dict,
+    const Type& value
+)
+:
+    dict_(dict),
+    keyword_("unknown-keyword"),
+    value_(value),
+    stored_(true)
+{}
+
+
+template<class Type>
+Foam::demandDrivenEntry<Type>::demandDrivenEntry
+(
+    const dictionary& dict,
+    const word& keyword
+)
+:
+    dict_(dict),
+    keyword_(keyword),
+    value_(pTraits<Type>::zero),
+    stored_(false)
+{}
+
+
+template<class Type>
+Foam::demandDrivenEntry<Type>::demandDrivenEntry
+(
+    const dictionary& dict,
+    const word& keyword,
+    const Type& defaultValue,
+    const bool readIfPresent
+)
+:
+    dict_(dict),
+    keyword_(keyword),
+    value_(defaultValue),
+    stored_(true)
+{
+    if (readIfPresent)
+    {
+        dict_.readIfPresent<Type>(keyword, value_);
+    }
+}
+
+
+template<class Type>
+Foam::demandDrivenEntry<Type>::demandDrivenEntry(const demandDrivenEntry& dde)
+:
+    dict_(dde.dict_),
+    keyword_(dde.keyword_),
+    value_(dde.value_),
+    stored_(dde.stored_)
+{}
+
+
+// ************************************************************************* //

src/OpenFOAM/primitives/demandDrivenEntry/demandDrivenEntry.H

@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::demandDrivenEntry
+
+Description
+    Class for demand-driven dictionary entries
+
+    Holds a reference to a dictionary, which is then queried if the value
+    is requested and has not already been cached
+
+SourceFiles
+    demandDrivenEntry.C
+    demandDrivenEntryI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef demandDrivenEntry_H
+#define demandDrivenEntry_H
+
+#include "dictionary.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class demandDrivenEntry Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class demandDrivenEntry
+{
+private:
+
+    // Private data
+
+        //- Reference to the dictionary
+        const dictionary& dict_;
+
+        //- Keyword to look up
+        const word keyword_;
+
+        //- Value
+        mutable Type value_;
+
+        //- Flag to say that the value has been stored
+        mutable bool stored_;
+
+
+public:
+
+    //- Constructors
+
+        //- Construct from dictionary and value - cannot be re-read
+        demandDrivenEntry
+        (
+            const dictionary& dict,
+            const Type& value
+        );
+
+
+        //- Construct from dictionary and keyword
+        demandDrivenEntry
+        (
+            const dictionary& dict,
+            const word& keyword
+        );
+
+
+        //- Construct from dictionary, keyword and default value
+        demandDrivenEntry
+        (
+            const dictionary& dict,
+            const word& keyword,
+            const Type& defaultValue,
+            const bool readIfPresent = true
+        );
+
+        //- Copy constructor
+        demandDrivenEntry(const demandDrivenEntry& dde);
+
+
+    // Public Member Functions
+
+        //- Initialise
+        inline void initialise() const;
+
+        //- Return the value
+        inline const Type& value() const;
+
+        //- Set the value
+        inline void setValue(const Type& value);
+
+        //- Reset the demand-driven entry
+        inline void reset();
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "demandDrivenEntryI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "demandDrivenEntry.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/primitives/demandDrivenEntry/demandDrivenEntryI.H

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "demandDrivenEntry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+inline void Foam::demandDrivenEntry<Type>::initialise() const
+{
+    if (!stored_)
+    {
+        dict_.lookup(keyword_) >> value_;
+        stored_ = true;
+    }
+}
+
+
+template<class Type>
+inline const Type& Foam::demandDrivenEntry<Type>::value() const
+{
+    initialise();
+
+    return value_;
+}
+
+
+template<class Type>
+inline void Foam::demandDrivenEntry<Type>::setValue(const Type& value)
+{
+//    dict_.set<Type>(keyword_, value);
+    value_ = value;
+    stored_ = true;
+}
+
+
+template<class Type>
+inline void Foam::demandDrivenEntry<Type>::reset()
+{
+    stored_ = false;
+}
+
+
+// ************************************************************************* //

src/combustionModels/PaSR/PaSR.C

@@ -42,7 +42,7 @@ Foam::combustionModels::PaSR<Type>::PaSR
     (
         IOobject
         (
-            typeName + ":kappa",
+            "PaSR:kappa",
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
@@ -117,7 +117,11 @@ template<class Type>
 Foam::tmp<Foam::volScalarField>
 Foam::combustionModels::PaSR<Type>::dQ() const
 {
-    return kappa_*laminar<Type>::dQ();
+    return
+        tmp<volScalarField>
+        (
+            new volScalarField("PaSR:dQ", kappa_*laminar<Type>::dQ())
+        );
 }
 
 
@@ -125,7 +129,11 @@ template<class Type>
 Foam::tmp<Foam::volScalarField>
 Foam::combustionModels::PaSR<Type>::Sh() const
 {
-    return kappa_*laminar<Type>::Sh();
+    return
+        tmp<volScalarField>
+        (
+            new volScalarField("PaSR:Sh", kappa_*laminar<Type>::Sh())
+        );
 }
 
 

src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.C

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

src/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChanger.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,7 +51,7 @@ void Foam::polyTopoChanger::readModifiers()
         {
             WarningIn("polyTopoChanger::readModifiers()")
                 << "Specified IOobject::MUST_READ_IF_MODIFIED but class"
-                << " does not support automatic rereading."
+                << " does not support automatic re-reading."
                 << endl;
         }
 
@@ -80,18 +80,13 @@ void Foam::polyTopoChanger::readModifiers()
         }
 
         // Check state of IOstream
-        is.check
+        is.check("polyTopoChanger::readModifiers()");
-        (
-            "polyTopoChanger::polyTopoChanger"
-            "(const IOobject&, const polyMesh&)"
-        );
 
         close();
     }
 }
 
 
-// Read constructor given IOobject and a polyMesh reference
 Foam::polyTopoChanger::polyTopoChanger
 (
     const IOobject& io,
@@ -106,7 +101,6 @@ Foam::polyTopoChanger::polyTopoChanger
 }
 
 
-// Read constructor given IOobject and a polyMesh reference
 Foam::polyTopoChanger::polyTopoChanger(polyMesh& mesh)
 :
     PtrList<polyMeshModifier>(),
@@ -133,7 +127,6 @@ Foam::polyTopoChanger::polyTopoChanger(polyMesh& mesh)
 }
 
 
-// Return a list of modifier types
 Foam::wordList Foam::polyTopoChanger::types() const
 {
     const PtrList<polyMeshModifier>& modifiers = *this;
@@ -149,7 +142,6 @@ Foam::wordList Foam::polyTopoChanger::types() const
 }
 
 
-// Return a list of modifier names
 Foam::wordList Foam::polyTopoChanger::names() const
 {
     const PtrList<polyMeshModifier>& modifiers = *this;
@@ -165,7 +157,6 @@ Foam::wordList Foam::polyTopoChanger::names() const
 }
 
 
-// Is topology change required
 bool Foam::polyTopoChanger::changeTopology() const
 {
     // Go through all mesh modifiers and accumulate the morphing information
@@ -211,7 +202,6 @@ bool Foam::polyTopoChanger::changeTopology() const
 }
 
 
-// Return topology change request
 Foam::autoPtr<Foam::polyTopoChange>
 Foam::polyTopoChanger::topoChangeRequest() const
 {
@@ -233,7 +223,6 @@ Foam::polyTopoChanger::topoChangeRequest() const
 }
 
 
-// Correct polyTopoChanger after moving points
 void Foam::polyTopoChanger::modifyMotionPoints(pointField& p) const
 {
     const PtrList<polyMeshModifier>& topoChanges = *this;
@@ -248,7 +237,6 @@ void Foam::polyTopoChanger::modifyMotionPoints(pointField& p) const
 }
 
 
-// Force recalculation of locally stored data on topological change
 void Foam::polyTopoChanger::update(const mapPolyMesh& m)
 {
     // Go through all mesh modifiers and accumulate the morphing information
@@ -299,7 +287,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChanger::changeMesh
 }
 
 
-// Add mesh modifiers to the morph engine
 void Foam::polyTopoChanger::addTopologyModifiers
 (
     const List<polyMeshModifier*>& tm
@@ -314,8 +301,10 @@ void Foam::polyTopoChanger::addTopologyModifiers
         {
             FatalErrorIn
             (
-                "void polyTopoChanger::addTopologyModifiers("
+                "void polyTopoChanger::addTopologyModifiers"
-                "const List<polyMeshModifier*>& tm)"
+                "("
+                    "const List<polyMeshModifier*>&"
+                ")"
             )   << "Mesh modifier created with different mesh reference."
                 << abort(FatalError);
         }
@@ -344,8 +333,7 @@ Foam::label Foam::polyTopoChanger::findModifierID
     // Modifier not found
     if (debug)
     {
-        Info<< "label polyTopoChanger::::findModifierID(const word& "
+        WarningIn("label polyTopoChanger::findModifierID(const word&) const")
-            << "modName) const"
             << "Modifier named " << modName << " not found.  "
             << "List of available modifier names: " << names() << endl;
     }
@@ -355,7 +343,6 @@ Foam::label Foam::polyTopoChanger::findModifierID
 }
 
 
-// writeData member function required by regIOobject
 bool Foam::polyTopoChanger::writeData(Ostream& os) const
 {
     os << *this;

src/fvOptions/sources/derived/MRFSource/MRFSource.C

@@ -175,6 +175,8 @@ bool Foam::fv::MRFSource::read(const dictionary& dict)
         coeffs_.readIfPresent("UName", UName_);
         coeffs_.readIfPresent("rhoName", rhoName_);
 
+        initialise();
+
         return true;
     }
     else

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

@@ -92,6 +92,7 @@ Foam::CollidingCloud<CloudType>::CollidingCloud
 )
 :
     CloudType(cloudName, rho, U, mu, g, false),
+    constProps_(this->particleProperties()),
     collisionModel_(NULL)
 {
     if (this->solution().steadyState())
@@ -246,4 +247,5 @@ void Foam::CollidingCloud<CloudType>::info()
         << rotationalKineticEnergy << nl;
 }
 
+
 // ************************************************************************* //

src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.H

@@ -97,6 +97,10 @@ protected:
 
     // Protected data
 
+        //- Thermo parcel constant properties
+        typename parcelType::constantProperties constProps_;
+
+
         // References to the cloud sub-models
 
             //- Collision model
@@ -180,6 +184,11 @@ public:
             //- Return a reference to the cloud copy
             inline const CollidingCloud& cloudCopy() const;
 
+            //- Return the constant properties
+            inline const typename parcelType::constantProperties&
+                constProps() const;
+
+
             //- If the collision model controls the wall interaction,
             //  then the wall impact distance should be zero.
             //  Otherwise, it should be allowed to be the value from

src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloudI.H

@@ -33,6 +33,14 @@ Foam::CollidingCloud<CloudType>::cloudCopy() const
 }
 
 
+template<class CloudType>
+inline const typename CloudType::particleType::constantProperties&
+Foam::CollidingCloud<CloudType>::constProps() const
+{
+    return constProps_;
+}
+
+
 template<class CloudType>
 inline const Foam::CollisionModel<Foam::CollidingCloud<CloudType> >&
 Foam::CollidingCloud<CloudType>::collision() const

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

@@ -308,7 +308,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
         )
     ),
     solution_(mesh_, particleProperties_.subDict("solution")),
-    constProps_(particleProperties_, solution_.active()),
+    constProps_(particleProperties_),
     subModelProperties_
     (
         particleProperties_.subOrEmptyDict("subModels", solution_.active())

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

@@ -106,7 +106,7 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
     CloudType(cloudName, rho, U, g, thermo, false),
     reactingCloud(),
     cloudCopyPtr_(NULL),
-    constProps_(this->particleProperties(), this->solution().active()),
+    constProps_(this->particleProperties()),
     compositionModel_(NULL),
     phaseChangeModel_(NULL),
     rhoTrans_(thermo.carrier().species().size())

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,7 +85,7 @@ Foam::ReactingMultiphaseCloud<CloudType>::ReactingMultiphaseCloud
     CloudType(cloudName, rho, U, g, thermo, false),
     reactingMultiphaseCloud(),
     cloudCopyPtr_(NULL),
-    constProps_(this->particleProperties(), this->solution().active()),
+    constProps_(this->particleProperties()),
     devolatilisationModel_(NULL),
     surfaceReactionModel_(NULL),
     dMassDevolatilisation_(0.0),

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

@@ -153,7 +153,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
     ),
     thermoCloud(),
     cloudCopyPtr_(NULL),
-    constProps_(this->particleProperties(), this->solution().active()),
+    constProps_(this->particleProperties()),
     thermo_(thermo),
     T_(thermo.thermo().T()),
     p_(thermo.thermo().p()),

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

@@ -74,6 +74,47 @@ class CollidingParcel
 :
     public ParcelType
 {
+public:
+
+    //- Class to hold thermo particle constant properties
+    class constantProperties
+    :
+        public ParcelType::constantProperties
+    {
+
+        // Private data
+
+            //- Young's modulus [N/m2]
+            demandDrivenEntry<scalar> youngsModulus_;
+
+            //- Poisson's ratio
+            demandDrivenEntry<scalar> poissonsRatio_;
+
+
+    public:
+
+        // Constructors
+
+            //- Null constructor
+            constantProperties();
+
+            //- Copy constructor
+            constantProperties(const constantProperties& cp);
+
+            //- Construct from dictionary
+            constantProperties(const dictionary& parentDict);
+
+
+        // Member functions
+
+            //- Return const access to Young's Modulus
+            inline scalar youngsModulus() const;
+
+            //- Return const access to Poisson's ratio
+            inline scalar poissonsRatio() const;
+    };
+
+
 protected:
 
     // Protected data

src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcelI.H

@@ -25,6 +25,40 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+template<class ParcelType>
+inline Foam::CollidingParcel<ParcelType>::constantProperties::
+constantProperties()
+:
+    ParcelType::constantProperties(),
+    youngsModulus_(this->dict_, 0.0),
+    poissonsRatio_(this->dict_, 0.0)
+{}
+
+
+template<class ParcelType>
+inline Foam::CollidingParcel<ParcelType>::constantProperties::constantProperties
+(
+    const constantProperties& cp
+)
+:
+    ParcelType::constantProperties(cp),
+    youngsModulus_(cp.youngsModulus_),
+    poissonsRatio_(cp.poissonsRatio_)
+{}
+
+
+template<class ParcelType>
+inline Foam::CollidingParcel<ParcelType>::constantProperties::constantProperties
+(
+    const dictionary& parentDict
+)
+:
+    ParcelType::constantProperties(parentDict),
+    youngsModulus_(this->dict_, "youngsModulus"),
+    poissonsRatio_(this->dict_, "poissonsRatio")
+{}
+
+
 template<class ParcelType>
 inline Foam::CollidingParcel<ParcelType>::CollidingParcel
 (
@@ -83,7 +117,25 @@ inline Foam::CollidingParcel<ParcelType>::CollidingParcel
 {}
 
 
-// * * * * * * * CollidingParcel Member Functions  * * * * * * * //
+// * * * * * * * * * constantProperties Member Functions * * * * * * * * * * //
+
+template<class ParcelType>
+inline Foam::scalar
+Foam::CollidingParcel<ParcelType>::constantProperties::youngsModulus() const
+{
+    return youngsModulus_.value();
+}
+
+
+template<class ParcelType>
+inline Foam::scalar
+Foam::CollidingParcel<ParcelType>::constantProperties::poissonsRatio() const
+{
+    return poissonsRatio_.value();
+}
+
+
+// * * * * * * * * * * CollidingParcel Member Functions  * * * * * * * * * * //
 
 template<class ParcelType>
 inline const Foam::vector& Foam::CollidingParcel<ParcelType>::f() const

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

@@ -48,6 +48,7 @@ SourceFiles
 #include "IOstream.H"
 #include "autoPtr.H"
 #include "interpolation.H"
+#include "demandDrivenEntry.H"
 
 // #include "ParticleForceList.H" // TODO
 
@@ -82,29 +83,30 @@ public:
     //- Class to hold kinematic particle constant properties
     class constantProperties
     {
-        // Private data
+    protected:
+
+        // Protected data
 
             //- Constant properties dictionary
             const dictionary dict_;
 
+
+    private:
+
+        // Private data
+
             //- Parcel type id - used for post-processing to flag the type
             //  of parcels issued by this cloud
-            label parcelTypeId_;
+            demandDrivenEntry<label> parcelTypeId_;
 
             //- Minimum density [kg/m3]
-            scalar rhoMin_;
+            demandDrivenEntry<scalar> rhoMin_;
 
             //- Particle density [kg/m3] (constant)
-            scalar rho0_;
+            demandDrivenEntry<scalar> rho0_;
 
             //- Minimum particle mass [kg]
-            scalar minParticleMass_;
+            demandDrivenEntry<scalar> minParticleMass_;
-
-            //- Young's modulus [N/m2]
-            scalar youngsModulus_;
-
-            //- Poisson's ratio
-            scalar poissonsRatio_;
 
 
     public:
@@ -117,23 +119,8 @@ public:
             //- Copy constructor
             constantProperties(const constantProperties& cp);
 
-            //- Constructor from dictionary
+            //- Construct from dictionary
-            constantProperties
+            constantProperties(const dictionary& parentDict);
-            (
-                const dictionary& parentDict,
-                const bool readFields = true
-            );
-
-            //- Construct from components
-            constantProperties
-            (
-                const label parcelTypeId,
-                const scalar rhoMin,
-                const scalar rho0,
-                const scalar minParticleMass,
-                const scalar youngsModulus,
-                const scalar poissonsRatio
-            );
 
 
         // Member functions
@@ -152,12 +139,6 @@ public:
 
             //- Return const access to the minimum particle mass
             inline scalar minParticleMass() const;
-
-            //- Return const access to Young's Modulus
-            inline scalar youngsModulus() const;
-
-            //- Return const access to Poisson's ratio
-            inline scalar poissonsRatio() const;
     };
 
 

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H

@@ -34,12 +34,10 @@ inline
 Foam::KinematicParcel<ParcelType>::constantProperties::constantProperties()
 :
     dict_(dictionary::null),
-    parcelTypeId_(-1),
+    parcelTypeId_(dict_, -1),
-    rhoMin_(0.0),
+    rhoMin_(dict_, 0.0),
-    rho0_(0.0),
+    rho0_(dict_, 0.0),
-    minParticleMass_(0.0),
+    minParticleMass_(dict_, 0.0)
-    youngsModulus_(0.0),
-    poissonsRatio_(0.0)
 {}
 
 
@@ -53,68 +51,21 @@ inline Foam::KinematicParcel<ParcelType>::constantProperties::constantProperties
     parcelTypeId_(cp.parcelTypeId_),
     rhoMin_(cp.rhoMin_),
     rho0_(cp.rho0_),
-    minParticleMass_(cp.minParticleMass_),
+    minParticleMass_(cp.minParticleMass_)
-    youngsModulus_(cp.youngsModulus_),
-    poissonsRatio_(cp.poissonsRatio_)
 {}
 
 
 template<class ParcelType>
 inline Foam::KinematicParcel<ParcelType>::constantProperties::constantProperties
 (
-    const dictionary& parentDict,
+    const dictionary& parentDict
-    const bool readFields
 )
 :
     dict_(parentDict.subOrEmptyDict("constantProperties")),
-    parcelTypeId_(1),
+    parcelTypeId_(dict_, "parcelTypeId", 1),
-    rhoMin_(1e-15),
+    rhoMin_(dict_, "rhoMin", 1e-15),
-    rho0_(0.0),
+    rho0_(dict_, "rho0"),
-    minParticleMass_(1e-15),
+    minParticleMass_(dict_, "minParticleMass", 1e-15)
-    youngsModulus_(0.0),
-    poissonsRatio_(0.0)
-{
-    if (readFields)
-    {
-        if (dict_.readIfPresent("parcelTypeId", parcelTypeId_))
-        {
-            Info<< "    employing parcel parcelTypeId of " << parcelTypeId_
-                << endl;
-        }
-        if (dict_.readIfPresent("rhoMin", rhoMin_))
-        {
-            Info<< "    employing parcel rhoMin of " << rhoMin_ << endl;
-        }
-        if (dict_.readIfPresent("minParticleMass", minParticleMass_))
-        {
-            Info<< "    employing parcel minParticleMass of "
-                << minParticleMass_ << endl;
-        }
-
-        dict_.lookup("rho0") >> rho0_;
-        dict_.lookup("youngsModulus") >> youngsModulus_;
-        dict_.lookup("poissonsRatio") >> poissonsRatio_;
-    }
-}
-
-template<class ParcelType>
-inline Foam::KinematicParcel<ParcelType>::constantProperties::constantProperties
-(
-    const label parcelTypeId,
-    const scalar rhoMin,
-    const scalar rho0,
-    const scalar minParticleMass,
-    const scalar youngsModulus,
-    const scalar poissonsRatio
-)
-:
-    dict_(dictionary::null),
-    parcelTypeId_(parcelTypeId),
-    rhoMin_(rhoMin),
-    rho0_(rho0),
-    minParticleMass_(minParticleMass),
-    youngsModulus_(youngsModulus),
-    poissonsRatio_(poissonsRatio)
 {}
 
 
@@ -192,7 +143,7 @@ template<class ParcelType>
 inline Foam::label
 Foam::KinematicParcel<ParcelType>::constantProperties::parcelTypeId() const
 {
-    return parcelTypeId_;
+    return parcelTypeId_.value();
 }
 
 
@@ -200,7 +151,7 @@ template<class ParcelType>
 inline Foam::scalar
 Foam::KinematicParcel<ParcelType>::constantProperties::rhoMin() const
 {
-    return rhoMin_;
+    return rhoMin_.value();
 }
 
 
@@ -208,7 +159,7 @@ template<class ParcelType>
 inline Foam::scalar
 Foam::KinematicParcel<ParcelType>::constantProperties::rho0() const
 {
-    return rho0_;
+    return rho0_.value();
 }
 
 
@@ -216,23 +167,7 @@ template<class ParcelType>
 inline Foam::scalar
 Foam::KinematicParcel<ParcelType>::constantProperties::minParticleMass() const
 {
-    return minParticleMass_;
+    return minParticleMass_.value();
-}
-
-
-template<class ParcelType>
-inline Foam::scalar
-Foam::KinematicParcel<ParcelType>::constantProperties::youngsModulus() const
-{
-    return youngsModulus_;
-}
-
-
-template<class ParcelType>
-inline Foam::scalar
-Foam::KinematicParcel<ParcelType>::constantProperties::poissonsRatio() const
-{
-    return poissonsRatio_;
 }
 
 

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

@@ -506,14 +506,19 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     scalarField& Cs
 ) const
 {
-    // Check that model is active, and that the parcel temperature is
+    // Check that model is active
-    // within necessary limits for devolatilisation to occur
+    if (!td.cloud().devolatilisation().active())
-    if
+    {
-    (
+        return;
-        !td.cloud().devolatilisation().active()
+    }
-     || T < td.cloud().constProps().TDevol()
+
-     || canCombust == -1
+    // Initialise demand-driven constants
-    )
+    (void)td.cloud().constProps().TDevol();
+    (void)td.cloud().constProps().LDevol();
+
+    // Check that the parcel temperature is within necessary limits for
+    // devolatilisation to occur
+    if (T < td.cloud().constProps().TDevol() || canCombust == -1)
     {
         return;
     }
@@ -604,11 +609,22 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
 ) const
 {
     // Check that model is active
-    if (!td.cloud().surfaceReaction().active() || (canCombust != 1))
+    if (!td.cloud().surfaceReaction().active())
     {
         return;
     }
 
+    // Initialise demand-driven constants
+    (void)td.cloud().constProps().hRetentionCoeff();
+    (void)td.cloud().constProps().TMax();
+
+    // Check that model is active
+    if (canCombust != 1)
+    {
+        return;
+    }
+
+
     // Update surface reactions
     const scalar hReaction = td.cloud().surfaceReaction().calculate
     (

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

@@ -40,6 +40,7 @@ SourceFiles
 
 #include "particle.H"
 #include "SLGThermo.H"
+#include "demandDrivenEntry.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -82,14 +83,14 @@ public:
         // Private data
 
             //- Devolatilisation activation temperature [K]
-            scalar TDevol_;
+            demandDrivenEntry<scalar> TDevol_;
 
             //- Latent heat of devolatilisation [J/kg]
-            scalar LDevol_;
+            demandDrivenEntry<scalar> LDevol_;
 
             //- Fraction of enthalpy retained by parcel due to surface
             //  reactions
-            scalar hRetentionCoeff_;
+            demandDrivenEntry<scalar> hRetentionCoeff_;
 
 
     public:
@@ -102,32 +103,9 @@ public:
             //- Copy constructor
             constantProperties(const constantProperties& cp);
 
-            //- Constructor from dictionary
+            //- Construct from dictionary
-            constantProperties
+            constantProperties(const dictionary& parentDict);
-            (
-                const dictionary& parentDict,
-                const bool readFields = true
-            );
 
-            //- Construct from components
-            constantProperties
-            (
-                const label parcelTypeId,
-                const scalar rhoMin,
-                const scalar rho0,
-                const scalar minParticleMass,
-                const scalar youngsModulus,
-                const scalar poissonsRatio,
-                const scalar T0,
-                const scalar TMin,
-                const scalar TMax,
-                const scalar Cp0,
-                const scalar epsilon0,
-                const scalar f0,
-                const scalar pMin,
-                const Switch& constantVolume,
-                const scalar TDevol
-            );
 
         // Access
 

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H

@@ -30,9 +30,9 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
 constantProperties()
 :
     ParcelType::constantProperties(),
-    TDevol_(0.0),
+    TDevol_(this->dict_, 0.0),
-    LDevol_(0.0),
+    LDevol_(this->dict_, 0.0),
-    hRetentionCoeff_(0.0)
+    hRetentionCoeff_(this->dict_, 0.0)
 {}
 
 
@@ -54,75 +54,13 @@ template<class ParcelType>
 inline Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
 constantProperties
 (
-    const dictionary& parentDict,
+    const dictionary& parentDict
-    const bool readFields
 )
 :
-    ParcelType::constantProperties(parentDict, readFields),
+    ParcelType::constantProperties(parentDict),
-    TDevol_(0.0),
+    TDevol_(this->dict_, "TDevol"),
-    LDevol_(0.0),
+    LDevol_(this->dict_, "LDevol"),
-    hRetentionCoeff_(0.0)
+    hRetentionCoeff_(this->dict_, "hRetentionCoeff")
-{
-    if (readFields)
-    {
-        this->dict().lookup("TDevol") >> TDevol_;
-        this->dict().lookup("LDevol") >> LDevol_;
-        this->dict().lookup("hRetentionCoeff") >> hRetentionCoeff_;
-
-        if ((hRetentionCoeff_ < 0) || (hRetentionCoeff_ > 1))
-        {
-            FatalErrorIn
-            (
-                "ReactingMultiphaseParcel<ParcelType>::constantProperties::"
-                "constantProperties"
-            )   << "hRetentionCoeff must be in the range 0 to 1" << nl
-                << exit(FatalError) << endl;
-        }
-
-        hRetentionCoeff_ = max(1e-06, hRetentionCoeff_);
-    }
-}
-
-
-template<class ParcelType>
-inline Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
-constantProperties
-(
-    const label parcelTypeId,
-    const scalar rhoMin,
-    const scalar rho0,
-    const scalar minParticleMass,
-    const scalar youngsModulus,
-    const scalar poissonsRatio,
-    const scalar T0,
-    const scalar TMin,
-    const scalar TMax,
-    const scalar Cp0,
-    const scalar epsilon0,
-    const scalar f0,
-    const scalar pMin,
-    const Switch& constantVolume,
-    const scalar TDevol
-)
-:
-    ParcelType::constantProperties
-    (
-        parcelTypeId,
-        rhoMin,
-        rho0,
-        minParticleMass,
-        youngsModulus,
-        poissonsRatio,
-        T0,
-        TMin,
-        TMax,
-        Cp0,
-        epsilon0,
-        f0,
-        pMin,
-        constantVolume
-    ),
-    TDevol_(TDevol)
 {}
 
 
@@ -198,7 +136,7 @@ template<class ParcelType>
 inline Foam::scalar
 Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::TDevol() const
 {
-    return TDevol_;
+    return TDevol_.value();
 }
 
 
@@ -206,7 +144,7 @@ template<class ParcelType>
 inline Foam::scalar
 Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::LDevol() const
 {
-    return LDevol_;
+    return LDevol_.value();
 }
 
 
@@ -215,7 +153,19 @@ inline Foam::scalar
 Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
 hRetentionCoeff() const
 {
-    return hRetentionCoeff_;
+    scalar value = hRetentionCoeff_.value();
+
+    if ((value < 0) || (value > 1))
+    {
+        FatalErrorIn
+        (
+            "ReactingMultiphaseParcel<ParcelType>::constantProperties::"
+            "constantProperties"
+        )   << "hRetentionCoeff must be in the range 0 to 1" << nl
+            << exit(FatalError) << endl;
+    }
+
+    return value;
 }
 
 

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

@@ -60,9 +60,14 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
     PhaseChangeModel<reactingCloudType>& phaseChange = td.cloud().phaseChange();
 
+    if (!phaseChange.active())
+    {
+        return;
+    }
+
     scalar Tvap = phaseChange.Tvap(YComponents);
 
-    if (!phaseChange.active() || T < Tvap || YPhase < SMALL)
+    if (T < Tvap || YPhase < SMALL)
     {
         return;
     }

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

@@ -40,6 +40,7 @@ SourceFiles
 
 #include "particle.H"
 #include "SLGThermo.H"
+#include "demandDrivenEntry.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -75,10 +76,10 @@ public:
         // Private data
 
             //- Minimum pressure [Pa]
-            scalar pMin_;
+            demandDrivenEntry<scalar> pMin_;
 
             //- Constant volume flag - e.g. during mass transfer
-            Switch constantVolume_;
+            demandDrivenEntry<bool> constantVolume_;
 
 
     public:
@@ -91,31 +92,8 @@ public:
             //- Copy constructor
             constantProperties(const constantProperties& cp);
 
-            //- Constructor from dictionary
+            //- Construct from dictionary
-            constantProperties
+            constantProperties(const dictionary& parentDict);
-            (
-                const dictionary& parentDict,
-                const bool readFields = true
-            );
-
-            //- Construct from components
-            constantProperties
-            (
-                const label parcelTypeId,
-                const scalar rhoMin,
-                const scalar rho0,
-                const scalar minParticleMass,
-                const scalar youngsModulus,
-                const scalar poissonsRatio,
-                const scalar T0,
-                const scalar TMin,
-                const scalar TMax,
-                const scalar Cp0,
-                const scalar epsilon0,
-                const scalar f0,
-                const scalar pMin,
-                const Switch& constantVolume
-            );
 
 
         // Access
@@ -124,7 +102,7 @@ public:
             inline scalar pMin() const;
 
             //- Return const access to the constant volume flag
-            inline Switch constantVolume() const;
+            inline bool constantVolume() const;
     };
 
 

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H

@@ -30,8 +30,8 @@ inline
 Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties()
 :
     ParcelType::constantProperties(),
-    pMin_(0.0),
+    pMin_(this->dict_, 0.0),
-    constantVolume_(false)
+    constantVolume_(this->dict_, false)
 {}
 
 
@@ -50,62 +50,12 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
 template<class ParcelType>
 inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
 (
-    const dictionary& parentDict,
+    const dictionary& parentDict
-    const bool readFields
 )
 :
-    ParcelType::constantProperties(parentDict, readFields),
+    ParcelType::constantProperties(parentDict),
-    pMin_(1000.0),
+    pMin_(this->dict_, "pMin", 1000.0),
-    constantVolume_(false)
+    constantVolume_(this->dict_, "constantVolume")
-{
-    if (readFields)
-    {
-        if (this->dict().readIfPresent("pMin", pMin_))
-        {
-            Info<< "    employing parcel pMin of " << pMin_ << endl;
-        }
-
-        this->dict().lookup("constantVolume") >> constantVolume_;
-    }
-}
-
-
-template<class ParcelType>
-inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
-(
-    const label parcelTypeId,
-    const scalar rhoMin,
-    const scalar rho0,
-    const scalar minParticleMass,
-    const scalar youngsModulus,
-    const scalar poissonsRatio,
-    const scalar T0,
-    const scalar TMin,
-    const scalar TMax,
-    const scalar Cp0,
-    const scalar epsilon0,
-    const scalar f0,
-    const scalar pMin,
-    const Switch& constantVolume
-)
-:
-    ParcelType::constantProperties
-    (
-        parcelTypeId,
-        rhoMin,
-        rho0,
-        minParticleMass,
-        youngsModulus,
-        poissonsRatio,
-        T0,
-        TMin,
-        TMax,
-        Cp0,
-        epsilon0,
-        f0
-    ),
-    pMin_(pMin),
-    constantVolume_(constantVolume)
 {}
 
 
@@ -178,15 +128,15 @@ template<class ParcelType>
 inline Foam::scalar
 Foam::ReactingParcel<ParcelType>::constantProperties::pMin() const
 {
-    return pMin_;
+    return pMin_.value();
 }
 
 
 template<class ParcelType>
-inline Foam::Switch
+inline bool
 Foam::ReactingParcel<ParcelType>::constantProperties::constantVolume() const
 {
-    return constantVolume_;
+    return constantVolume_.value();
 }
 
 

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

@@ -41,6 +41,7 @@ SourceFiles
 
 #include "particle.H"
 #include "SLGThermo.H"
+#include "demandDrivenEntry.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -77,22 +78,22 @@ public:
         // Private data
 
             //- Particle initial temperature [K]
-            scalar T0_;
+            demandDrivenEntry<scalar> T0_;
 
             //- Minimum temperature [K]
-            scalar TMin_;
+            demandDrivenEntry<scalar> TMin_;
 
             //- Maximum temperature [K]
-            scalar TMax_;
+            demandDrivenEntry<scalar> TMax_;
 
             //- Particle specific heat capacity [J/(kg.K)]
-            scalar Cp0_;
+            demandDrivenEntry<scalar> Cp0_;
 
             //- Particle emissivity [] (radiation)
-            scalar epsilon0_;
+            demandDrivenEntry<scalar> epsilon0_;
 
             //- Particle scattering factor [] (radiation)
-            scalar f0_;
+            demandDrivenEntry<scalar> f0_;
 
 
     public:
@@ -105,29 +106,8 @@ public:
             //- Copy constructor
             constantProperties(const constantProperties& cp);
 
-            //- Constructor from dictionary
+            //- Construct from dictionary
-            constantProperties
+            constantProperties(const dictionary& parentDict);
-            (
-                const dictionary& parentDict,
-                const bool readFields = true
-            );
-
-            //- Construct from components
-            constantProperties
-            (
-                const label parcelTypeId,
-                const scalar rhoMin,
-                const scalar rho0,
-                const scalar minParticleMass,
-                const scalar youngsModulus,
-                const scalar poissonsRatio,
-                const scalar T0,
-                const scalar TMin,
-                const scalar TMax,
-                const scalar Cp0,
-                const scalar epsilon0,
-                const scalar f0
-            );
 
 
         // Member functions
@@ -143,8 +123,8 @@ public:
                 //- Return const access to maximum temperature [K]
                 inline scalar TMax() const;
 
-                //- Return non-const access to maximum temperature [K]
+                //- Set the maximum temperature [K]
-                inline scalar& TMax();
+                inline void setTMax(const scalar TMax);
 
                 //- Return const access to the particle specific heat capacity
                 //  [J/(kg.K)]

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

@@ -29,12 +29,12 @@ template<class ParcelType>
 inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties()
 :
     ParcelType::constantProperties(),
-    T0_(0.0),
+    T0_(this->dict_, 0.0),
-    TMin_(0.0),
+    TMin_(this->dict_, 0.0),
-    TMax_(VGREAT),
+    TMax_(this->dict_, VGREAT),
-    Cp0_(0.0),
+    Cp0_(this->dict_, 0.0),
-    epsilon0_(0.0),
+    epsilon0_(this->dict_, 0.0),
-    f0_(0.0)
+    f0_(this->dict_, 0.0)
 {}
 
 
@@ -57,69 +57,16 @@ inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties
 template<class ParcelType>
 inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties
 (
-    const dictionary& parentDict,
+    const dictionary& parentDict
-    const bool readFields
 )
 :
-    ParcelType::constantProperties(parentDict, readFields),
+    ParcelType::constantProperties(parentDict),
-    T0_(0.0),
+    T0_(this->dict_, "T0"),
-    TMin_(200),
+    TMin_(this->dict_, "TMin", 200.0),
-    TMax_(5000),
+    TMax_(this->dict_, "TMax", 5000.0),
-    Cp0_(0.0),
+    Cp0_(this->dict_, "Cp0"),
-    epsilon0_(0.0),
+    epsilon0_(this->dict_, "epsilon0"),
-    f0_(0.0)
+    f0_(this->dict_, "f0")
-{
-    if (readFields)
-    {
-        if (this->dict().readIfPresent("TMin", TMin_))
-        {
-            Info<< "    employing parcel TMin of " << TMin_ << endl;
-        }
-        if (this->dict().readIfPresent("TMax", TMax_))
-        {
-            Info<< "    employing parcel TMax of " << TMax_ << endl;
-        }
-
-        this->dict().lookup("T0") >> T0_;
-        this->dict().lookup("Cp0") >> Cp0_;
-        this->dict().lookup("epsilon0") >> epsilon0_;
-        this->dict().lookup("f0") >> f0_;
-    }
-}
-
-
-template<class ParcelType>
-inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties
-(
-    const label parcelTypeId,
-    const scalar rhoMin,
-    const scalar rho0,
-    const scalar minParticleMass,
-    const scalar youngsModulus,
-    const scalar poissonsRatio,
-    const scalar T0,
-    const scalar TMin,
-    const scalar TMax,
-    const scalar Cp0,
-    const scalar epsilon0,
-    const scalar f0
-)
-:
-    ParcelType::constantProperties
-    (
-        parcelTypeId,
-        rhoMin,
-        rho0,
-        minParticleMass,
-        youngsModulus,
-        poissonsRatio
-    ),
-    T0_(T0),
-    TMin_(TMin),
-    TMax_(TMax),
-    Cp0_(Cp0),
-    epsilon0_(epsilon0),
-    f0_(f0)
 {}
 
 
@@ -190,7 +137,7 @@ template<class ParcelType>
 inline Foam::scalar
 Foam::ThermoParcel<ParcelType>::constantProperties::T0() const
 {
-    return T0_;
+    return T0_.value();
 }
 
 
@@ -198,7 +145,7 @@ template<class ParcelType>
 inline Foam::scalar
 Foam::ThermoParcel<ParcelType>::constantProperties::TMin() const
 {
-    return TMin_;
+    return TMin_.value();
 }
 
 
@@ -206,15 +153,15 @@ template<class ParcelType>
 inline Foam::scalar
 Foam::ThermoParcel<ParcelType>::constantProperties::TMax() const
 {
-    return TMax_;
+    return TMax_.value();
 }
 
 
 template<class ParcelType>
-inline Foam::scalar&
+inline void
-Foam::ThermoParcel<ParcelType>::constantProperties::TMax()
+Foam::ThermoParcel<ParcelType>::constantProperties::setTMax(const scalar TMax)
 {
-    return TMax_;
+    TMax_.setValue(TMax);
 }
 
 
@@ -222,7 +169,7 @@ template<class ParcelType>
 inline Foam::scalar
 Foam::ThermoParcel<ParcelType>::constantProperties::Cp0() const
 {
-    return Cp0_;
+    return Cp0_.value();
 }
 
 
@@ -230,7 +177,7 @@ template<class ParcelType>
 inline Foam::scalar
 Foam::ThermoParcel<ParcelType>::constantProperties::epsilon0() const
 {
-    return epsilon0_;
+    return epsilon0_.value();
 }
 
 
@@ -238,7 +185,7 @@ template<class ParcelType>
 inline Foam::scalar
 Foam::ThermoParcel<ParcelType>::constantProperties::f0() const
 {
-    return f0_;
+    return f0_.value();
 }
 
 

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

@@ -90,7 +90,7 @@ void Foam::SprayParcel<ParcelType>::calc
     }
 
     // set the maximum temperature limit
-    td.cloud().constProps().TMax() = TMax;
+    td.cloud().constProps().setTMax(TMax);
 
     // store the parcel properties
     const scalarField& Y(this->Y());

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

@@ -33,6 +33,7 @@ Description
 #define SprayParcel_H
 
 #include "particle.H"
+#include "demandDrivenEntry.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -69,10 +70,10 @@ public:
         // Private data
 
             //- Particle initial surface tension [N/m]
-            scalar sigma0_;
+            demandDrivenEntry<scalar> sigma0_;
 
             //- Particle initial dynamic viscosity [Pa.s]
-            scalar mu0_;
+            demandDrivenEntry<scalar> mu0_;
 
 
     public:
@@ -85,12 +86,8 @@ public:
             //- Copy constructor
             constantProperties(const constantProperties& cp);
 
-            //- Constructor from dictionary
+            //- Construct from dictionary
-            constantProperties
+            constantProperties(const dictionary& parentDict);
-            (
-                const dictionary& parentDict,
-                const bool readFields = true
-            );
 
             //- Construct from components
             constantProperties

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

@@ -29,8 +29,8 @@ template<class ParcelType>
 inline Foam::SprayParcel<ParcelType>::constantProperties::constantProperties()
 :
     ParcelType::constantProperties(),
-    sigma0_(0.0),
+    sigma0_(this->dict_, 0.0),
-    mu0_(0.0)
+    mu0_(this->dict_, 0.0)
 {}
 
 
@@ -49,19 +49,13 @@ inline Foam::SprayParcel<ParcelType>::constantProperties::constantProperties
 template<class ParcelType>
 inline Foam::SprayParcel<ParcelType>::constantProperties::constantProperties
 (
-    const dictionary& parentDict,
+    const dictionary& parentDict
-    const bool readFields
 )
 :
-    ParcelType::constantProperties(parentDict, readFields),
+    ParcelType::constantProperties(parentDict),
-    sigma0_(0.0),
+    sigma0_(this->dict_, "sigma0"),
-    mu0_(0.0)
+    mu0_(this->dict_, "mu0")
-{
+{}
-    if (readFields)
-    {
-        this->dict().lookup("sigma0") >> sigma0_;
-    }
-}
 
 
 template<class ParcelType>
@@ -104,8 +98,8 @@ inline Foam::SprayParcel<ParcelType>::constantProperties::constantProperties
         pMin,
         constantVolume
     ),
-    sigma0_(sigma0),
+    sigma0_(this->dict_, sigma0),
-    mu0_(mu0)
+    mu0_(this->dict_, mu0)
 {}
 
 
@@ -199,6 +193,24 @@ inline Foam::SprayParcel<ParcelType>::SprayParcel
 {}
 
 
+// * * * * * * * * * constantProperties Member Functions * * * * * * * * * * //
+
+template<class ParcelType>
+inline Foam::scalar
+Foam::SprayParcel<ParcelType>::constantProperties::sigma0() const
+{
+    return sigma0_.value();
+}
+
+
+template<class ParcelType>
+inline Foam::scalar
+Foam::SprayParcel<ParcelType>::constantProperties::mu0() const
+{
+    return mu0_.value();
+}
+
+
 // * * * * * * * * * * SprayParcel Member Functions * * * * * * * * * * * * //
 
 template<class ParcelType>

src/postProcessing/functionObjects/field/Make/options

@@ -7,11 +7,13 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/turbulenceModels \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
+    -lsurfMesh \
     -llagrangian \
     -lfileFormats \
     -lsurfMesh \
@@ -19,4 +21,3 @@ LIB_LIBS = \
     -lincompressibleTransportModels \
     -lcompressibleTurbulenceModel \
     -lincompressibleTurbulenceModel
-

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

@@ -48,11 +48,13 @@ namespace Foam
 
 
     template<>
-    const char* NamedEnum<fieldValues::faceSource::operationType, 12>::names[] =
+    const char* NamedEnum<fieldValues::faceSource::operationType, 14>::names[] =
     {
         "none",
         "sum",
+        "sumMag",
         "sumDirection",
+        "sumDirectionBalance",
         "average",
         "weightedAverage",
         "areaAverage",
@@ -75,7 +77,7 @@ namespace Foam
 const Foam::NamedEnum<Foam::fieldValues::faceSource::sourceType, 3>
     Foam::fieldValues::faceSource::sourceTypeNames_;
 
-const Foam::NamedEnum<Foam::fieldValues::faceSource::operationType, 12>
+const Foam::NamedEnum<Foam::fieldValues::faceSource::operationType, 14>
     Foam::fieldValues::faceSource::operationTypeNames_;
 
 
@@ -450,6 +452,11 @@ void Foam::fieldValues::faceSource::initialise(const dictionary& dict)
         Info<< "    weight field = " << weightFieldName_;
     }
 
+    if (dict.readIfPresent("scaleFactor", scaleFactor_))
+    {
+        Info<< "    scale factor = " << scaleFactor_;
+    }
+
     Info<< nl << endl;
 
     if (valueOutput_)
@@ -499,24 +506,15 @@ Foam::scalar Foam::fieldValues::faceSource::processValues
     {
         case opSumDirection:
         {
-            const vector direction(dict_.lookup("direction"));
+            vector n(dict_.lookup("direction"));
+            return sum(pos(values*(Sf & n))*mag(values));
+        }
+        case opSumDirectionBalance:
+        {
+            vector n(dict_.lookup("direction"));
+            const scalarField nv(values*(Sf & n));
 
-            scalar v = 0.0;
+            return sum(pos(nv)*mag(values) - neg(nv)*mag(values));
-
-            forAll(Sf, i)
-            {
-                scalar d = Sf[i] & direction;
-                if (d > 0)
-                {
-                    v += pos(values[i])*values[i];
-                }
-                else
-                {
-                    v += neg(values[i])*values[i];
-                }
-            }
-
-            return v;
         }
         default:
         {
@@ -539,8 +537,19 @@ Foam::vector Foam::fieldValues::faceSource::processValues
     {
         case opSumDirection:
         {
-            const vector direction(dict_.lookup("direction"));
+            vector n(dict_.lookup("direction"));
-            return sum(pos(values & direction)*values);
+            n /= mag(n) + ROOTVSMALL;
+            const scalarField nv(n & values);
+
+            return sum(pos(nv)*n*(nv));
+        }
+        case opSumDirectionBalance:
+        {
+            vector n(dict_.lookup("direction"));
+            n /= mag(n) + ROOTVSMALL;
+            const scalarField nv(n & values);
+
+            return sum(pos(nv)*n*(nv));
         }
         case opAreaNormalAverage:
         {
@@ -576,6 +585,7 @@ Foam::fieldValues::faceSource::faceSource
     source_(sourceTypeNames_.read(dict.lookup("source"))),
     operation_(operationTypeNames_.read(dict.lookup("operation"))),
     weightFieldName_("none"),
+    scaleFactor_(1.0),
     nFaces_(0),
     faceId_(),
     facePatchId_(),

src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H

@@ -71,6 +71,7 @@ Description
         sourceName   | name of face source if required  | no |
         operation    | operation to perform    | yes         |
         weightField  | name of field to apply weighting | no |
+        scaleFactor  | scale factor            | no          | 1
         fields       | list of fields to operate on | yes    |
     \endtable
 
@@ -87,7 +88,9 @@ Description
     \plaintable
        none          | no operation
        sum           | sum
+       sumMag        | sum of component magnitudes
        sumDirection  | sum values which are positive in given direction
+       sumDirectionBalance | sum of balance of values in given direction
        average       | ensemble average
        weightedAverage | weighted average
        areaAverage   | area weighted average
@@ -177,7 +180,9 @@ public:
         {
             opNone,
             opSum,
+            opSumMag,
             opSumDirection,
+            opSumDirectionBalance,
             opAverage,
             opWeightedAverage,
             opAreaAverage,
@@ -190,7 +195,7 @@ public:
         };
 
         //- Operation type names
-        static const NamedEnum<operationType, 12> operationTypeNames_;
+        static const NamedEnum<operationType, 14> operationTypeNames_;
 
 
 private:
@@ -237,6 +242,9 @@ protected:
         //- Weight field name - optional
         word weightFieldName_;
 
+        //- Scale factor - optional
+        scalar scaleFactor_;
+
         //- Global number of faces
         label nFaces_;
 

src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C

@@ -141,6 +141,11 @@ Type Foam::fieldValues::faceSource::processSameTypeValues
             result = sum(values);
             break;
         }
+        case opSumMag:
+        {
+            result = sum(cmptMag(values));
+            break;
+        }
         case opSumDirection:
         {
             FatalErrorIn
@@ -161,6 +166,26 @@ Type Foam::fieldValues::faceSource::processSameTypeValues
             result = pTraits<Type>::zero;
             break;
         }
+        case opSumDirectionBalance:
+        {
+            FatalErrorIn
+            (
+                "template<class Type>"
+                "Type Foam::fieldValues::faceSource::processSameTypeValues"
+                "("
+                    "const Field<Type>&, "
+                    "const vectorField&, "
+                    "const scalarField&"
+                ") const"
+            )
+                << "Operation " << operationTypeNames_[operation_]
+                << " not available for values of type "
+                << pTraits<Type>::typeName
+                << exit(FatalError);
+
+            result = pTraits<Type>::zero;
+            break;
+        }
         case opAverage:
         {
             result = sum(values)/values.size();
@@ -304,8 +329,8 @@ bool Foam::fieldValues::faceSource::writeValues(const word& fieldName)
             );
         }
 
-        // apply weight field
+        // apply scale factor and weight field
-        values *= weightField;
+        values *= scaleFactor_*weightField;
 
         if (Pstream::master())
         {

src/postProcessing/functionObjects/field/nearWallFields/findCellParticle.H

@@ -61,6 +61,7 @@ class findCellParticle
         //- passive data
         label data_;
 
+
 public:
 
     friend class Cloud<findCellParticle>;
@@ -75,7 +76,6 @@ public:
 
     public:
 
-
         // Constructors
 
             trackingData
@@ -97,6 +97,7 @@ public:
             {
                 return cellToData_;
             }
+
             List<List<point> >& cellToEnd()
             {
                 return cellToEnd_;
@@ -104,7 +105,6 @@ public:
     };
 
 
-
     // Constructors
 
         //- Construct from components
@@ -171,7 +171,6 @@ public:
         }
 
 
-
         // Tracking
 
             //- Track all particles to their end point

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

@@ -33,7 +33,7 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(nearWallFields, 0);
+    defineTypeNameAndDebug(nearWallFields, 0);
 }
 
 
@@ -72,8 +72,8 @@ void Foam::nearWallFields::calcAddressing()
         label patchI = iter.key();
         const fvPatch& patch = mesh.boundary()[patchI];
 
-        vectorField nf = patch.nf();
+        vectorField nf(patch.nf());
-        vectorField faceCellCentres = patch.patch().faceCellCentres();
+        vectorField faceCellCentres(patch.patch().faceCellCentres());
 
         forAll(patch, patchFaceI)
         {

src/postProcessing/functionObjects/field/surfaceInterpolateFields/surfaceInterpolateFields.C

@@ -92,6 +92,8 @@ void Foam::surfaceInterpolateFields::execute()
 {
     if (active_)
     {
+        Info<< type() << " " << name_ << " output:" << nl;
+
         // Clear out any previously loaded fields
         ssf_.clear();
         svf_.clear();
@@ -104,6 +106,8 @@ void Foam::surfaceInterpolateFields::execute()
         interpolateFields<sphericalTensor>(sSpheretf_);
         interpolateFields<symmTensor>(sSymmtf_);
         interpolateFields<tensor>(stf_);
+
+        Info<< endl;
     }
 }
 

src/regionModels/pyrolysisModels/noPyrolysis/noPyrolysis.C

@@ -56,6 +56,7 @@ void noPyrolysis::constructThermoChemistry()
     radiation_.reset(radiation::radiationModel::New(solidThermo_->T()).ptr());
 }
 
+
 bool noPyrolysis::read()
 {
     if (pyrolysisModel::read())
@@ -86,9 +87,14 @@ bool noPyrolysis::read(const dictionary& dict)
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-noPyrolysis::noPyrolysis(const word& modelType, const fvMesh& mesh)
+noPyrolysis::noPyrolysis
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& regionType
+)
 :
-    pyrolysisModel(mesh),
+    pyrolysisModel(mesh, regionType),
     solidChemistry_(NULL),
     solidThermo_(NULL),
     radiation_(NULL)
@@ -104,9 +110,11 @@ noPyrolysis::noPyrolysis
 (
     const word& modelType,
     const fvMesh& mesh,
-    const dictionary& dict
+    const dictionary& dict,
-):
+    const word& regionType
-    pyrolysisModel(mesh),
+)
+:
+    pyrolysisModel(mesh, regionType),
     solidChemistry_(NULL),
     solidThermo_(NULL),
     radiation_(NULL)
@@ -117,6 +125,7 @@ noPyrolysis::noPyrolysis
     }
 }
 
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 noPyrolysis::~noPyrolysis()

src/regionModels/pyrolysisModels/noPyrolysis/noPyrolysis.H

@@ -100,14 +100,20 @@ public:
     // Constructors
 
         //- Construct from type name and mesh
-        noPyrolysis(const word& modelType, const fvMesh& mesh);
+        noPyrolysis
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& regionType
+        );
 
         //- Construct from type name and mesh and dict
         noPyrolysis
         (
             const word& modelType,
             const fvMesh& mesh,
-            const dictionary& dict
+            const dictionary& dict,
+            const word& regionType
         );
 
 

src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,15 +80,20 @@ bool pyrolysisModel::read(const dictionary& dict)
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-pyrolysisModel::pyrolysisModel(const fvMesh& mesh)
+pyrolysisModel::pyrolysisModel(const fvMesh& mesh, const word& regionType)
 :
-    regionModel1D(mesh)
+    regionModel1D(mesh, regionType)
 {}
 
 
-pyrolysisModel::pyrolysisModel(const word& modelType, const fvMesh& mesh)
+pyrolysisModel::pyrolysisModel
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& regionType
+)
 :
-    regionModel1D(mesh, "pyrolysis", modelType)
+    regionModel1D(mesh, regionType, modelType)
 {
     if (active_)
     {
@@ -101,10 +106,11 @@ pyrolysisModel::pyrolysisModel
 (
     const word& modelType,
     const fvMesh& mesh,
-    const dictionary& dict
+    const dictionary& dict,
+    const word& regionType
 )
 :
-    regionModel1D(mesh, "pyrolysis", modelType, dict)
+    regionModel1D(mesh, regionType, modelType, dict)
 {
     if (active_)
     {

src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,9 +22,10 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::pyrolysisModel
+    Foam::regionModels::pyrolysisModels::pyrolysisModel
 
 Description
+    Base class for pyrolysis models
 
 SourceFiles
     pyrolysisModelI.H
@@ -104,9 +105,10 @@ public:
             mesh,
             (
                 const word& modelType,
-                const fvMesh& mesh
+                const fvMesh& mesh,
+                const word& regionType
             ),
-            (modelType, mesh)
+            (modelType, mesh, regionType)
         );
 
         declareRunTimeSelectionTable
@@ -117,26 +119,37 @@ public:
             (
                 const word& modelType,
                 const fvMesh& mesh,
-                const dictionary& dict
+                const dictionary& dict,
+                const word& regionType
             ),
-            (modelType, mesh, dict)
+            (modelType, mesh, dict, regionType)
         );
 
 
     // Constructors
 
         //- Construct null from mesh
-        pyrolysisModel(const fvMesh& mesh);
+        pyrolysisModel
+        (
+            const fvMesh& mesh,
+            const word& regionType
+        );
 
         //- Construct from type name and mesh
-        pyrolysisModel(const word& modelType, const fvMesh& mesh);
+        pyrolysisModel
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& regionType
+        );
 
         //- Construct from type name and mesh and dictionary
         pyrolysisModel
         (
             const word& modelType,
             const fvMesh& mesh,
-            const dictionary& dict
+            const dictionary& dict,
+            const word& regionType
         );
 
         //- Return clone
@@ -150,39 +163,20 @@ public:
     // Selectors
 
         //- Return a reference to the selected pyrolysis model
-        static autoPtr<pyrolysisModel> New(const fvMesh& mesh);
+        static autoPtr<pyrolysisModel> New
+        (
+            const fvMesh& mesh,
+            const word& regionType = "pyrolysis"
+        );
 
         //- Return a reference to a named selected pyrolysis model
         static autoPtr<pyrolysisModel> New
         (
             const fvMesh& mesh,
-            const dictionary& dict
+            const dictionary& dict,
+            const word& regionType = "pyrolysis"
         );
 
-         //- Return pointer to new pyrolysis created on freestore from Istream
-        class iNew
-        {
-            const fvMesh& mesh_;
-
-        public:
-
-            iNew(const fvMesh& mesh)
-            :
-                mesh_(mesh)
-            {}
-
-            autoPtr<pyrolysisModel> operator()(Istream& is) const
-            {
-                keyType key(is);
-                dictionary dict(is);
-
-                return autoPtr<pyrolysisModel>
-                (
-                    pyrolysisModel::New(mesh_, dict)
-                );
-            }
-        };
-
 
     //- Destructor
     virtual ~pyrolysisModel();
@@ -212,6 +206,7 @@ public:
                 //- Return the total gas mass flux to primary region [kg/m2/s]
                 virtual const surfaceScalarField& phiGas() const = 0;
 
+
             // Sources
 
                 //- External hook to add mass to the primary region

src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModelCollection.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,11 +30,13 @@ License
 
 namespace Foam
 {
-    defineTemplateTypeNameAndDebug
+    namespace regionModels
-    (
+    {
-        IOPtrList<regionModels::pyrolysisModels::pyrolysisModel>,
+        namespace pyrolysisModels
-        0
+        {
-    );
+            defineTypeNameAndDebug(pyrolysisModelCollection, 0);
+        }
+    }
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -48,12 +50,12 @@ namespace pyrolysisModels
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-pyrolysisModelCollection::pyrolysisModelCollection
+pyrolysisModelCollection::pyrolysisModelCollection(const fvMesh& mesh)
-(
-    const fvMesh& mesh
-)
 :
-    IOPtrList<pyrolysisModel>
+    PtrList<pyrolysisModel>()
+
+{
+    IOdictionary pyrolysisZonesDict
     (
         IOobject
         (
@@ -62,10 +64,32 @@ pyrolysisModelCollection::pyrolysisModelCollection
             mesh,
             IOobject::MUST_READ,
             IOobject::NO_WRITE
-        ),
+        )
-        pyrolysisModel::iNew(mesh)
+    );
-    ),
+
-    mesh_(mesh)
+    const wordList regions(pyrolysisZonesDict.toc());
+
+    setSize(regions.size());
+
+    for (label i = 0; i < regions.size(); i++)
+    {
+        set
+        (
+            i,
+            pyrolysisModel::New
+            (
+                mesh,
+                pyrolysisZonesDict.subDict(regions[i]),
+                regions[i]
+            )
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+
+pyrolysisModelCollection::~pyrolysisModelCollection()
 {}
 
 

src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModelCollection.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,16 +38,15 @@ SourceFiles
 #ifndef pyrolysisModelCollection_H
 #define pyrolysisModelCollection_H
 
-#include "IOPtrList.H"
+#include "PtrList.H"
 #include "pyrolysisModel.H"
-#include "fvc.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
+// Forward class declarations
 class fvMesh;
 
 namespace regionModels
@@ -56,19 +55,13 @@ namespace pyrolysisModels
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class pyrolysisModelCollection Declaration
+                  Class pyrolysisModelCollection Declaration
 \*---------------------------------------------------------------------------*/
 
 class pyrolysisModelCollection
 :
-    public IOPtrList<pyrolysisModel>
+    public PtrList<pyrolysisModel>
 {
-    // Private data
-
-        //- Reference to the finite volume mesh this zone is part of
-        const fvMesh& mesh_;
-
-
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -80,6 +73,9 @@ class pyrolysisModelCollection
 
 public:
 
+    // Runtime type information
+    TypeName("pyrolysisModelCollection");
+
 
     // Constructors
 
@@ -87,26 +83,29 @@ public:
         pyrolysisModelCollection(const fvMesh&);
 
 
+    //- Destructor
+    virtual ~pyrolysisModelCollection();
+
+
     // Member Functions
 
         //- Pre-evolve regions
-        void preEvolveRegion();
+        virtual void preEvolveRegion();
 
-         //- Evolve the pyrolysis equation regions
+        //- Evolve the pyrolysis equation regions
-        void evolveRegion();
+        virtual void evolveRegion();
 
         //- Evolve regions
-        void evolve();
+        virtual void evolve();
 
         //- Provide some feedback from pyrolysis regions
-        void info() const;
+        virtual void info() const;
 
         //- Return max diffusivity allowed in the solid
-        scalar maxDiff() const;
+        virtual scalar maxDiff() const;
 
         //- Mean diffusion number of the solid regions
-        scalar solidRegionDiffNo() const;
+        virtual scalar solidRegionDiffNo() const;
-
 };
 
 

src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModelNew.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,7 +38,11 @@ namespace pyrolysisModels
 
 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
 
-autoPtr<pyrolysisModel> pyrolysisModel::New(const fvMesh& mesh)
+autoPtr<pyrolysisModel> pyrolysisModel::New
+(
+    const fvMesh& mesh,
+    const word& regionType
+)
 {
     // get model name, but do not register the dictionary
     const word modelType
@@ -47,7 +51,7 @@ autoPtr<pyrolysisModel> pyrolysisModel::New(const fvMesh& mesh)
         (
             IOobject
             (
-                "pyrolysisProperties",
+                regionType + "Properties",
                 mesh.time().constant(),
                 mesh,
                 IOobject::MUST_READ,
@@ -64,21 +68,22 @@ autoPtr<pyrolysisModel> pyrolysisModel::New(const fvMesh& mesh)
 
     if (cstrIter == meshConstructorTablePtr_->end())
     {
-        FatalErrorIn("pyrolysisModel::New(const fvMesh&)")
+        FatalErrorIn("pyrolysisModel::New(const fvMesh&, const word&)")
             << "Unknown pyrolysisModel type " << modelType
             << nl << nl << "Valid pyrolisisModel types are:" << nl
             << meshConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }
 
-    return autoPtr<pyrolysisModel>(cstrIter()(modelType, mesh));
+    return autoPtr<pyrolysisModel>(cstrIter()(modelType, mesh, regionType));
 }
 
 
 autoPtr<pyrolysisModel> pyrolysisModel::New
 (
     const fvMesh& mesh,
-    const dictionary& dict
+    const dictionary& dict,
+    const word& regionType
 )
 {
 
@@ -91,14 +96,31 @@ autoPtr<pyrolysisModel> pyrolysisModel::New
 
     if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
-        FatalErrorIn("pyrolysisModel::New(const fvMesh&, const dictionary&)")
+        FatalErrorIn
+        (
+            "pyrolysisModel::New"
+            "("
+                "const fvMesh&, "
+                "const dictionary&, "
+                "const word&"
+            ")"
+        )
             << "Unknown pyrolysisModel type " << modelType
             << nl << nl << "Valid pyrolisisModel types are:" << nl
             << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }
 
-    return autoPtr<pyrolysisModel>(cstrIter()(modelType, mesh, dict));
+    return autoPtr<pyrolysisModel>
+    (
+        cstrIter()
+        (
+            modelType,
+            mesh,
+            dict,
+            regionType
+        )
+    );
 }
 
 

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C

@@ -382,9 +382,14 @@ void reactingOneDim::calculateMassTransfer()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
+reactingOneDim::reactingOneDim
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& regionType
+)
 :
-    pyrolysisModel(modelType, mesh),
+    pyrolysisModel(modelType, mesh, regionType),
     solidChemistry_(basicSolidChemistryModel::New(regionMesh())),
     solidThermo_(solidChemistry_->solidThermo()),
     radiation_(radiation::radiationModel::New(solidThermo_.T())),
@@ -482,10 +487,11 @@ reactingOneDim::reactingOneDim
 (
     const word& modelType,
     const fvMesh& mesh,
-    const dictionary& dict
+    const dictionary& dict,
+    const word& regionType
 )
 :
-    pyrolysisModel(modelType, mesh, dict),
+    pyrolysisModel(modelType, mesh, dict, regionType),
     solidChemistry_(basicSolidChemistryModel::New(regionMesh())),
     solidThermo_(solidChemistry_->solidThermo()),
     radiation_(radiation::radiationModel::New(solidThermo_.T())),

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H

@@ -203,14 +203,20 @@ public:
     // Constructors
 
         //- Construct from type name and mesh
-        reactingOneDim(const word& modelType, const fvMesh& mesh);
+        reactingOneDim
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& regionType
+        );
 
         //- Construct from type name, mesh and dictionary
         reactingOneDim
         (
             const word& modelType,
             const fvMesh& mesh,
-            const dictionary& dict
+            const dictionary& dict,
+            const word& regionType
         );
 
 

src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisRadiativeCoupledMixed/filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.C

@@ -38,13 +38,35 @@ const filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::filmModelType&
 filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::
 filmModel() const
 {
-    const regionModels::regionModel& model =
+    HashTable<const filmModelType*> models
-        db().time().lookupObject<regionModels::regionModel>
+        = db().time().lookupClass<filmModelType>();
-        (
-            "surfaceFilmProperties"
-        );
 
-    return dynamic_cast<const filmModelType&>(model);
+    forAllConstIter(HashTable<const filmModelType*>, models, iter)
+    {
+        if (iter()->regionMesh().name() == filmRegionName_)
+        {
+            return *iter();
+        }
+    }
+
+    DynamicList<word> modelNames;
+    forAllConstIter(HashTable<const filmModelType*>, models, iter)
+    {
+        modelNames.append(iter()->regionMesh().name());
+    }
+
+    FatalErrorIn
+    (
+        "const filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::"
+        "filmModelType& "
+        "filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::"
+        "filmModel() const"
+    )
+        << "Unable to locate film region " << filmRegionName_
+        << ".  Available regions include: " << modelNames
+        << abort(FatalError);
+
+    return **models.begin();
 }
 
 
@@ -53,13 +75,36 @@ pyrolysisModelType&
 filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::
 pyrModel() const
 {
-    const regionModels::regionModel& model =
+    HashTable<const pyrolysisModelType*> models =
-        db().time().lookupObject<regionModels::regionModel>
+        db().time().lookupClass<pyrolysisModelType>();
-        (
-            "pyrolysisProperties"
-        );
 
-    return dynamic_cast<const pyrolysisModelType&>(model);
+    forAllConstIter(HashTable<const pyrolysisModelType*>, models, iter)
+    {
+        if (iter()->regionMesh().name() == pyrolysisRegionName_)
+        {
+            return *iter();
+        }
+    }
+
+    DynamicList<word> modelNames;
+    forAllConstIter(HashTable<const pyrolysisModelType*>, models, iter)
+    {
+        modelNames.append(iter()->regionMesh().name());
+    }
+
+
+    FatalErrorIn
+    (
+        "const filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::"
+        "pyrolysisModelType& "
+        "filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::"
+        "pyrModel() const"
+    )
+        << "Unable to locate pyrolysis region " << pyrolysisRegionName_
+        << ".  Available regions include: " << modelNames
+        << abort(FatalError);
+
+    return **models.begin();
 }
 
 
@@ -74,6 +119,8 @@ filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(p, iF),
     temperatureCoupledBase(patch(), "undefined", "undefined-K"),
+    filmRegionName_("surfaceFilmProperties"),
+    pyrolysisRegionName_("pyrolysisProperties"),
     TnbrName_("undefined-Tnbr"),
     QrNbrName_("undefined-QrNbr"),
     QrName_("undefined-Qr"),
@@ -98,6 +145,8 @@ filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(psf, p, iF, mapper),
     temperatureCoupledBase(patch(), psf.KMethod(), psf.kappaName()),
+    filmRegionName_(psf.filmRegionName_),
+    pyrolysisRegionName_(psf.pyrolysisRegionName_),
     TnbrName_(psf.TnbrName_),
     QrNbrName_(psf.QrNbrName_),
     QrName_(psf.QrName_),
@@ -117,6 +166,14 @@ filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(p, iF),
     temperatureCoupledBase(patch(), dict),
+    filmRegionName_
+    (
+        dict.lookupOrDefault<word>("filmRegion", "surfaceFilmProperties")
+    ),
+    pyrolysisRegionName_
+    (
+        dict.lookupOrDefault<word>("pyrolysisRegion", "pyrolysisProperties")
+    ),
     TnbrName_(dict.lookup("Tnbr")),
     QrNbrName_(dict.lookup("QrNbr")),
     QrName_(dict.lookup("Qr")),
@@ -171,6 +228,8 @@ filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(psf, iF),
     temperatureCoupledBase(patch(), psf.KMethod(), psf.kappaName()),
+    filmRegionName_(psf.filmRegionName_),
+    pyrolysisRegionName_(psf.pyrolysisRegionName_),
     TnbrName_(psf.TnbrName_),
     QrNbrName_(psf.QrNbrName_),
     QrName_(psf.QrName_),
@@ -196,7 +255,6 @@ updateCoeffs()
 
     const label patchI = patch().index();
     const label nbrPatchI = mpp.samplePolyPatch().index();
-
     const polyMesh& mesh = patch().boundaryMesh().mesh();
     const polyMesh& nbrMesh = mpp.sampleMesh();
     const fvPatch& nbrPatch =
@@ -240,13 +298,13 @@ updateCoeffs()
 
     // Obtain Rad heat (Qr)
     scalarField Qr(patch().size(), 0.0);
-    if (QrName_ != "none") //region0
+    if (QrName_ != "none") // primary region (region0)
     {
         Qr = patch().lookupPatchField<volScalarField, scalar>(QrName_);
         myPatchINrbPatchI = nbrPatch.index();
     }
 
-    if (QrNbrName_ != "none") //pyrolysis
+    if (QrNbrName_ != "none") // pyrolysis region
     {
         Qr = nbrPatch.lookupPatchField<volScalarField, scalar>(QrNbrName_);
         mpp.distribute(Qr);
@@ -326,13 +384,12 @@ updateCoeffs()
             << "     convective heat[W] : " << Qc << nl
             << "     radiative heat [W] : " << Qr << nl
             << "     total heat     [W] : " << Qt << nl
-            << "     walltemperature "
+            << "     wall temperature "
             << " min:" << gMin(*this)
             << " max:" << gMax(*this)
             << " avg:" << gAverage(*this)
             << endl;
     }
-
 }
 
 
@@ -342,6 +399,20 @@ void filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::write
 ) const
 {
     mixedFvPatchScalarField::write(os);
+    writeEntryIfDifferent<word>
+    (
+        os,
+        "filmRegion",
+        "surfaceFilmProperties",
+        filmRegionName_
+    );
+    writeEntryIfDifferent<word>
+    (
+        os,
+        "pyrolysisRegion",
+        "pyrolysisProperties",
+        pyrolysisRegionName_
+    );
     os.writeKeyword("Tnbr")<< TnbrName_ << token::END_STATEMENT << nl;
     os.writeKeyword("QrNbr")<< QrNbrName_ << token::END_STATEMENT << nl;
     os.writeKeyword("Qr")<< QrName_ << token::END_STATEMENT << nl;

src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisRadiativeCoupledMixed/filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.H

@@ -112,6 +112,12 @@ private:
 
     // Private data
 
+        //- Name of film region
+        const word filmRegionName_;
+
+        //- Name of pyrolysis region
+        const word pyrolysisRegionName_;
+
         //- Name of field on the neighbour region
         const word TnbrName_;
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,6 +39,8 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(p, iF),
+    filmRegionName_("surfaceFilmProperties"),
+    pyrolysisRegionName_("pyrolysisProperties"),
     phiName_("phi"),
     rhoName_("rho")
 {}
@@ -54,6 +56,8 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    filmRegionName_(ptf.filmRegionName_),
+    pyrolysisRegionName_(ptf.pyrolysisRegionName_),
     phiName_(ptf.phiName_),
     rhoName_(ptf.rhoName_)
 {}
@@ -68,6 +72,14 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(p, iF),
+    filmRegionName_
+    (
+        dict.lookupOrDefault<word>("filmRegion", "surfaceFilmProperties")
+    ),
+    pyrolysisRegionName_
+    (
+        dict.lookupOrDefault<word>("pyrolysisRegion", "pyrolysisProperties")
+    ),
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
     rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
 {
@@ -82,6 +94,8 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(fptpsf),
+    filmRegionName_(fptpsf.filmRegionName_),
+    pyrolysisRegionName_(fptpsf.pyrolysisRegionName_),
     phiName_(fptpsf.phiName_),
     rhoName_(fptpsf.rhoName_)
 {}
@@ -95,6 +109,8 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(fptpsf, iF),
+    filmRegionName_(fptpsf.filmRegionName_),
+    pyrolysisRegionName_(fptpsf.pyrolysisRegionName_),
     phiName_(fptpsf.phiName_),
     rhoName_(fptpsf.rhoName_)
 {}
@@ -117,11 +133,10 @@ void Foam::filmPyrolysisTemperatureCoupledFvPatchScalarField::updateCoeffs()
     int oldTag = UPstream::msgType();
     UPstream::msgType() = oldTag+1;
 
-    bool filmOk =
+    bool filmOk = db().time().foundObject<filmModelType>(filmRegionName_);
-        db().time().foundObject<filmModelType>("surfaceFilmProperties");
 
 
-    bool pyrOk = db().time().foundObject<pyrModelType>("pyrolysisProperties");
+    bool pyrOk = db().time().foundObject<pyrModelType>(pyrolysisRegionName_);
 
     if (!filmOk || !pyrOk)
     {
@@ -135,7 +150,7 @@ void Foam::filmPyrolysisTemperatureCoupledFvPatchScalarField::updateCoeffs()
 
     // Retrieve film model
     const filmModelType& filmModel =
-        db().time().lookupObject<filmModelType>("surfaceFilmProperties");
+        db().time().lookupObject<filmModelType>(filmRegionName_);
 
     const label filmPatchI = filmModel.regionPatchID(patchI);
 
@@ -147,7 +162,7 @@ void Foam::filmPyrolysisTemperatureCoupledFvPatchScalarField::updateCoeffs()
 
     // Retrieve pyrolysis model
     const pyrModelType& pyrModel =
-        db().time().lookupObject<pyrModelType>("pyrolysisProperties");
+        db().time().lookupObject<pyrModelType>(pyrolysisRegionName_);
 
     const label pyrPatchI = pyrModel.regionPatchID(patchI);
 
@@ -171,6 +186,20 @@ void Foam::filmPyrolysisTemperatureCoupledFvPatchScalarField::write
 ) const
 {
     fvPatchScalarField::write(os);
+    writeEntryIfDifferent<word>
+    (
+        os,
+        "filmRegion",
+        "surfaceFilmProperties",
+        filmRegionName_
+    );
+    writeEntryIfDifferent<word>
+    (
+        os,
+        "pyrolysisRegion",
+        "pyrolysisProperties",
+        pyrolysisRegionName_
+    );
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
     writeEntry("value", os);

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -70,6 +70,12 @@ class filmPyrolysisTemperatureCoupledFvPatchScalarField
 {
     // Private data
 
+        //- Name of film region
+        const word filmRegionName_;
+
+        //- Name of pyrolysis region
+        const word pyrolysisRegionName_;
+
         //- Name of flux field
         word phiName_;
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,6 +39,8 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 )
 :
     fixedValueFvPatchVectorField(p, iF),
+    filmRegionName_("surfaceFilmProperties"),
+    pyrolysisRegionName_("pyrolysisProperties"),
     phiName_("phi"),
     rhoName_("rho")
 {}
@@ -54,6 +56,8 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 )
 :
     fixedValueFvPatchVectorField(ptf, p, iF, mapper),
+    filmRegionName_(ptf.filmRegionName_),
+    pyrolysisRegionName_(ptf.pyrolysisRegionName_),
     phiName_(ptf.phiName_),
     rhoName_(ptf.rhoName_)
 {}
@@ -68,6 +72,14 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 )
 :
     fixedValueFvPatchVectorField(p, iF),
+    filmRegionName_
+    (
+        dict.lookupOrDefault<word>("filmRegion", "surfaceFilmProperties")
+    ),
+    pyrolysisRegionName_
+    (
+        dict.lookupOrDefault<word>("pyrolysisRegion", "pyrolysisProperties")
+    ),
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
     rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
 {
@@ -82,6 +94,8 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 )
 :
     fixedValueFvPatchVectorField(fpvpvf),
+    filmRegionName_(fpvpvf.filmRegionName_),
+    pyrolysisRegionName_(fpvpvf.pyrolysisRegionName_),
     phiName_(fpvpvf.phiName_),
     rhoName_(fpvpvf.rhoName_)
 {}
@@ -95,6 +109,8 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 )
 :
     fixedValueFvPatchVectorField(fpvpvf, iF),
+    filmRegionName_(fpvpvf.filmRegionName_),
+    pyrolysisRegionName_(fpvpvf.pyrolysisRegionName_),
     phiName_(fpvpvf.phiName_),
     rhoName_(fpvpvf.rhoName_)
 {}
@@ -117,11 +133,10 @@ void Foam::filmPyrolysisVelocityCoupledFvPatchVectorField::updateCoeffs()
     int oldTag = UPstream::msgType();
     UPstream::msgType() = oldTag+1;
 
-    bool foundFilm =
+    bool foundFilm = db().time().foundObject<filmModelType>(filmRegionName_);
-        db().time().foundObject<filmModelType>("surfaceFilmProperties");
 
     bool foundPyrolysis =
-        db().time().foundObject<pyrModelType>("pyrolysisProperties");
+        db().time().foundObject<pyrModelType>(pyrolysisRegionName_);
 
     if (!foundFilm || !foundPyrolysis)
     {
@@ -135,7 +150,7 @@ void Foam::filmPyrolysisVelocityCoupledFvPatchVectorField::updateCoeffs()
 
     // Retrieve film model
     const filmModelType& filmModel =
-        db().time().lookupObject<filmModelType>("surfaceFilmProperties");
+        db().time().lookupObject<filmModelType>(filmRegionName_);
 
     const label filmPatchI = filmModel.regionPatchID(patchI);
 
@@ -147,7 +162,7 @@ void Foam::filmPyrolysisVelocityCoupledFvPatchVectorField::updateCoeffs()
 
     // Retrieve pyrolysis model
     const pyrModelType& pyrModel =
-        db().time().lookupObject<pyrModelType>("pyrolysisProperties");
+        db().time().lookupObject<pyrModelType>(pyrolysisRegionName_);
 
     const label pyrPatchI = pyrModel.regionPatchID(patchI);
 
@@ -201,6 +216,20 @@ void Foam::filmPyrolysisVelocityCoupledFvPatchVectorField::write
 ) const
 {
     fvPatchVectorField::write(os);
+    writeEntryIfDifferent<word>
+    (
+        os,
+        "filmRegion",
+        "surfaceFilmProperties",
+        filmRegionName_
+    );
+    writeEntryIfDifferent<word>
+    (
+        os,
+        "pyrolysisRegion",
+        "pyrolysisProperties",
+        pyrolysisRegionName_
+    );
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
     writeEntry("value", os);

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -70,6 +70,12 @@ class filmPyrolysisVelocityCoupledFvPatchVectorField
 {
     // Private data
 
+        //- Name of film region
+        word filmRegionName_;
+
+        //- Name of pyrolysis region
+        word pyrolysisRegionName_;
+
         //- Name of flux field
         word phiName_;
 

src/regionModels/regionModel/regionModel/regionModel.C

@@ -288,11 +288,41 @@ Foam::label Foam::regionModels::regionModel::nbrCoupledPatchID
     // boundary mesh
     const polyBoundaryMesh& nbrPbm = nbrRegionMesh.boundaryMesh();
 
-    const mappedPatchBase& mpb =
+    const polyBoundaryMesh& pbm = regionMesh().boundaryMesh();
-        refCast<const mappedPatchBase>
+
+    if (regionPatchI > pbm.size() - 1)
+    {
+        FatalErrorIn
         (
-            regionMesh().boundaryMesh()[regionPatchI]
+            "Foam::label Foam::regionModels::regionModel::nbrCoupledPatchID"
-        );
+            "("
+                "const regionModel&, "
+                "const label"
+            ") const"
+        )
+            << "region patch index out of bounds: "
+            << "region patch index = " << regionPatchI
+            << ", maximum index = " << pbm.size() - 1
+            << abort(FatalError);
+    }
+
+    const polyPatch& pp = regionMesh().boundaryMesh()[regionPatchI];
+
+    if (!isA<mappedPatchBase>(pp))
+    {
+        FatalErrorIn
+        (
+            "Foam::label Foam::regionModels::regionModel::nbrCoupledPatchID"
+            "("
+                "const regionModel&, "
+                "const label"
+            ") const"
+        )
+            << "Expected a " << mappedPatchBase::typeName
+            << " patch, but found a " << pp.type() << abort(FatalError);
+    }
+
+    const mappedPatchBase& mpb = refCast<const mappedPatchBase>(pp);
 
     // sample patch name on the primary region
     const word& primaryPatchName = mpb.samplePatch();
@@ -335,13 +365,17 @@ Foam::label Foam::regionModels::regionModel::nbrCoupledPatchID
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::regionModels::regionModel::regionModel(const fvMesh& mesh)
+Foam::regionModels::regionModel::regionModel
+(
+    const fvMesh& mesh,
+    const word& regionType
+)
 :
     IOdictionary
     (
         IOobject
         (
-            "regionModelProperties",
+            regionType + "Properties",
             mesh.time().constant(),
             mesh.time(),
             IOobject::NO_READ,

src/regionModels/regionModel/regionModel/regionModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -140,7 +140,7 @@ protected:
         virtual bool read(const dictionary& dict);
 
         //- Create or return a new inter-region AMI object
-        virtual const Foam::AMIPatchToPatchInterpolation& interRegionAMI
+        virtual const AMIPatchToPatchInterpolation& interRegionAMI
         (
             const regionModel& nbrRegion,
             const label regionPatchI,
@@ -148,6 +148,7 @@ protected:
             const bool flip
         );
 
+
 public:
 
     //- Runtime type information
@@ -157,7 +158,7 @@ public:
     // Constructors
 
         //- Construct null
-        regionModel(const fvMesh& mesh);
+        regionModel(const fvMesh& mesh, const word& regionType);
 
         //- Construct from mesh, region type and name
         regionModel