Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C

@@ -37,7 +37,7 @@ Description
 #include "fvCFD.H"
 #include "psiThermo.H"
 #include "turbulenceModel.H"
-#include "MRFZones.H"
+#include "IOMRFZoneList.H"
 #include "IOporosityModelList.H"
 #include "IObasicSourceList.H"
 #include "fvcSmooth.H"

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/createZones.H

@@ -1,4 +1,4 @@
-    MRFZones mrfZones(mesh);
+    IOMRFZoneList mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U);
 
     IOporosityModelList pZones(mesh);

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C

@@ -37,7 +37,7 @@ Description
 #include "psiThermo.H"
 #include "turbulenceModel.H"
 #include "bound.H"
-#include "MRFZones.H"
+#include "IOMRFZoneList.H"
 #include "IOporosityModelList.H"
 #include "IObasicSourceList.H"
 #include "pimpleControl.H"

applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createZones.H

@@ -1,4 +1,4 @@
-    MRFZones mrfZones(mesh);
+    IOMRFZoneList mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U);
 
     IOporosityModelList pZones(mesh);

applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C

@@ -34,7 +34,7 @@ Description
 #include "fvCFD.H"
 #include "rhoThermo.H"
 #include "RASModel.H"
-#include "MRFZones.H"
+#include "IOMRFZoneList.H"
 #include "IObasicSourceList.H"
 #include "IOporosityModelList.H"
 #include "simpleControl.H"

applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H

@@ -8,6 +8,8 @@
 
     UEqn().relax();
 
+    mrfZones.addCoriolis(rho, UEqn());
+
     if (simple.momentumPredictor())
     {
         solve

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C

@@ -34,6 +34,7 @@ Description
 #include "RASModel.H"
 #include "fixedGradientFvPatchFields.H"
 #include "simpleControl.H"
+#include "IOMRFZoneList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,7 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createZones.H"
     #include "initContinuityErrs.H"
 
     simpleControl simple(mesh);

applications/solvers/heatTransfer/buoyantSimpleFoam/createZones.H

@@ -0,0 +1,3 @@
+    IOMRFZoneList mrfZones(mesh);
+    mrfZones.correctBoundaryVelocity(U);
+

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -17,6 +17,8 @@
         fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
     );
 
+    mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
     bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);
 
     phiHbyA += phig;

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -60,14 +60,9 @@ int main(int argc, char *argv[])
 
     #include "createFluidMeshes.H"
     #include "createSolidMeshes.H"
-    #include "createPorousFluidRegions.H"
-    #include "createPorousSolidMeshes.H"
 
     #include "createFluidFields.H"
     #include "createSolidFields.H"
-    #include "createPorousFluidFields.H"
-    #include "createPorousSolidFields.H"
-
 
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
@@ -116,24 +111,6 @@ int main(int argc, char *argv[])
                 #include "solveFluid.H"
             }
 
-            forAll(porousFluidRegions, i)
-            {
-                Info<< "\nSolving for fluid porous region "
-                    << porousFluidRegions[i].name() << endl;
-                #include "setPorousFluidFields.H"
-                #include "readPorousFluidRegionPIMPLEControls.H"
-                #include "solvePorousFluid.H"
-            }
-
-            forAll(porousSolidRegions, i)
-            {
-                Info<< "\nSolving for porous solid region "
-                    << porousSolidRegions[i].name() << endl;
-                #include "setPorousRegionSolidFields.H"
-                #include "readPorousSolidMultiRegionPIMPLEControls.H"
-                #include "solvePorousSolid.H"
-            }
-
             forAll(solidRegions, i)
             {
                 Info<< "\nSolving for solid region "

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options

@@ -1,8 +1,7 @@
 EXE_INC = \
     -Ifluid \
     -Isolid \
-    -I./porousFluid \
+    -I../solid \
-    -I./porousSolid \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/cfdTools \

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C

@@ -50,13 +50,9 @@ int main(int argc, char *argv[])
 
     #include "createFluidMeshes.H"
     #include "createSolidMeshes.H"
-    #include "createPorousFluidRegions.H"
-    #include "createPorousSolidMeshes.H"
 
     #include "createFluidFields.H"
     #include "createSolidFields.H"
-    #include "createPorousFluidFields.H"
-    #include "createPorousSolidFields.H"
 
     #include "initContinuityErrs.H"
 
@@ -74,24 +70,6 @@ int main(int argc, char *argv[])
             #include "solveFluid.H"
         }
 
-        forAll(porousFluidRegions, i)
-        {
-            Info<< "\nSolving for fluid porous region "
-                << porousFluidRegions[i].name() << endl;
-            #include "setPorousFluidFields.H"
-            #include "readPorousFluidRegionSIMPLEControls.H"
-            #include "solvePorousFluid.H"
-        }
-
-        forAll(porousSolidRegions, i)
-        {
-            Info<< "\nSolving for porous solid region "
-                << porousSolidRegions[i].name() << endl;
-            #include "setPorousRegionSolidFields.H"
-            #include "readPorousSolidMultiRegionSIMPLEControls.H"
-            #include "solvePorousSolid.H"
-        }
-
         forAll(solidRegions, i)
         {
             Info<< "\nSolving for solid region "

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/UPorousFluidEqn.H

@@ -1,11 +0,0 @@
-    // Solve the Momentum equation
-    tmp<fvVectorMatrix> porousUEqn
-    (
-        fvm::div(porousPhi, porousU)
-      + turbPorous.divDevRhoReff(porousU)
-      + porousSources(porousRho, porousU)
-    );
-
-    porousUEqn().relax();
-
-    solve(porousUEqn() == -fvc::grad(porousP));

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/createPorousFluidFields.H

@@ -1,181 +0,0 @@
-    // Initialise porous field pointer lists
-    PtrList<rhoThermo> thermoPorous(porousFluidRegions.size());
-    PtrList<volScalarField> rhoPorous(porousFluidRegions.size());
-    PtrList<volScalarField> kappaPorous(porousFluidRegions.size());
-    PtrList<volVectorField> UPorous(porousFluidRegions.size());
-    PtrList<surfaceScalarField> phiPorous(porousFluidRegions.size());
-    PtrList<compressible::turbulenceModel> turbulencePorous
-    (
-        porousFluidRegions.size()
-    );
-    PtrList<volScalarField> pPorous(porousFluidRegions.size());
-
-    List<scalar> initialMassFluidPorous(porousFluidRegions.size());
-    List<label> pRefCellFluidPorous(porousFluidRegions.size(),0);
-    List<scalar> pRefValueFluidPorous(porousFluidRegions.size(),0.0);
-
-    PtrList<dimensionedScalar> rhoMaxPorous(fluidRegions.size());
-    PtrList<dimensionedScalar> rhoMinPorous(fluidRegions.size());
-
-    PtrList<IObasicSourceList> heatPorousSources
-    (
-        porousFluidRegions.size()
-    );
-
-    forAll(porousFluidRegions, i)
-    {
-        Info<< "Reading fluid mesh thermophysical properties for porous "
-            << porousFluidRegions[i].name() << nl << endl;
-
-        Info<< "    Adding to thermoFluid porous\n" << endl;
-
-        thermoPorous.set
-        (
-            i,
-            rhoThermo::New(porousFluidRegions[i]).ptr()
-        );
-
-        Info<< "    Adding to rhoPorous\n" << endl;
-        rhoPorous.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "rho",
-                    runTime.timeName(),
-                    porousFluidRegions[i],
-                    IOobject::NO_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                thermoPorous[i].rho()
-            )
-        );
-
-        Info<< "    Adding to UPorous\n" << endl;
-        UPorous.set
-        (
-            i,
-            new volVectorField
-            (
-                IOobject
-                (
-                    "U",
-                    runTime.timeName(),
-                    porousFluidRegions[i],
-                    IOobject::MUST_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                porousFluidRegions[i]
-            )
-        );
-
-        Info<< "    Adding to phiPorous\n" << endl;
-        phiPorous.set
-        (
-            i,
-            new surfaceScalarField
-            (
-                IOobject
-                (
-                    "phi",
-                    runTime.timeName(),
-                    porousFluidRegions[i],
-                    IOobject::READ_IF_PRESENT,
-                    IOobject::AUTO_WRITE
-                ),
-                linearInterpolate(rhoPorous[i]*UPorous[i])
-                    & porousFluidRegions[i].Sf()
-            )
-        );
-
-        Info<< "    Adding turbulence to porous\n" << endl;
-        turbulencePorous.set
-        (
-            i,
-            compressible::turbulenceModel::New
-            (
-                rhoPorous[i],
-                UPorous[i],
-                phiPorous[i],
-                thermoPorous[i]
-            ).ptr()
-        );
-
-        Info<< "    Adding to kappaFluid\n" << endl;
-        kappaPorous.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "kappaPorous",
-                    runTime.timeName(),
-                    porousFluidRegions[i],
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                thermoPorous[i].Cp()*thermoPorous[i].alpha()
-            )
-        );
-
-        pPorous.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "p",
-                    runTime.timeName(),
-                    porousFluidRegions[i],
-                    IOobject::MUST_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                porousFluidRegions[i]
-            )
-        );
-
-        setRefCell
-        (
-            thermoPorous[i].p(),
-            pPorous[i],
-            porousFluidRegions[i].solutionDict().subDict("SIMPLE"),
-            pRefCellFluidPorous[i],
-            pRefValueFluidPorous[i]
-        );
-
-        rhoMaxPorous.set
-        (
-            i,
-            new dimensionedScalar
-            (
-                porousFluidRegions[i].solutionDict().subDict("SIMPLE").lookup
-                (
-                    "rhoMax"
-                )
-            )
-        );
-
-        rhoMinPorous.set
-        (
-            i,
-            new dimensionedScalar
-            (
-                porousFluidRegions[i].solutionDict().subDict("SIMPLE").lookup
-                (
-                    "rhoMin"
-                )
-            )
-        );
-
-        heatPorousSources.set
-        (
-            i,
-            new IObasicSourceList(porousFluidRegions[i])
-        );
-    }
-
-

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/createPorousFluidRegions.H

@@ -1,25 +0,0 @@
-    const wordList porousFluidNames(rp["porousFluid"]);
-
-    PtrList<fvMesh> porousFluidRegions(porousFluidNames.size());
-
-    forAll (porousFluidNames, iPorous)
-    {
-        const word porousFluidName = porousFluidNames[iPorous];
-        Info<< "Create porous fluid region " << porousFluidName
-            << nl << endl;
-
-        porousFluidRegions.set
-        (
-            iPorous,
-            new fvMesh
-            (
-                IOobject
-                (
-                    porousFluidName,
-                    runTime.timeName(),
-                    runTime,
-                    IOobject::MUST_READ
-                )
-            )
-        );
-    }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/hPorousFluidEqn.H

@@ -1,19 +0,0 @@
-{
-    fvScalarMatrix hPorousEqn
-    (
-        fvm::div(porousPhi, porousH)
-      - fvm::laplacian(turbPorous.alphaEff(), porousH)
-     ==
-      - fvc::div(porousPhi, 0.5*magSqr(porousU), "div(phi,K)")
-      + porousSources(porousRho, porousH)
-    );
-
-    hPorousEqn.relax();
-    hPorousEqn.solve();
-
-    porousThermo.correct();
-
-    Info<< "Min/max in the porous T:"
-        << min(porousThermo.T()).value() << ' '
-        << max(porousThermo.T()).value() << endl;
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/pPorousFluidEqn.H

@@ -1,53 +0,0 @@
-porousRho = porousThermo.rho();
-porousRho = max(porousRho, rhoMin);
-porousRho = min(porousRho, rhoMax);
-porousRho.relax();
-
-volScalarField rAUPorous(1.0/porousUEqn().A());
-
-porousU = rAUPorous*porousUEqn().H();
-porousUEqn.clear();
-
-bool closedVolume = false;
-
-porousPhi =
-    fvc::interpolate(porousRho)
-    *(fvc::interpolate(porousU) & porousMesh.Sf());
-
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-{
-    fvScalarMatrix pEqn
-    (
-        fvm::laplacian(porousRho*rAUPorous, porousP) ==  fvc::div(porousPhi)
-    );
-
-    pEqn.setReference(pRefCell, pRefValue);
-
-    pEqn.solve();
-
-    if (nonOrth == nNonOrthCorr)
-    {
-        porousPhi -= pEqn.flux();
-
-    }
-}
-
-porousP.relax();
-
-porousU -= rAUPorous*fvc::grad(porousP);
-porousU.correctBoundaryConditions();
-
-if (closedVolume)
-{
-    porousP += (initialMass - fvc::domainIntegrate(porousPsi*porousP))
-        /fvc::domainIntegrate(porousPsi);
-}
-
-porousRho = porousThermo.rho();
-porousRho = max(porousRho, rhoMin);
-porousRho = min(porousRho, rhoMax);
-porousRho.relax();
-
-Info<< "rho max/min : "
-    << max(porousRho).value() << " "
-    << min(porousRho).value() << endl;

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/readPorousFluidRegionSIMPLEControls.H

@@ -1,4 +0,0 @@
-    const dictionary& simple = porousMesh.solutionDict().subDict("SIMPLE");
-
-    const int nNonOrthCorr =
-        simple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/setPorousFluidFields.H

@@ -1,28 +0,0 @@
-    const fvMesh& porousMesh = porousFluidRegions[i];
-
-    rhoThermo& porousThermo = thermoPorous[i];
-    volScalarField& porousRho = rhoPorous[i];
-
-    volVectorField& porousU = UPorous[i];
-    surfaceScalarField& porousPhi = phiPorous[i];
-
-    compressible::turbulenceModel& turbPorous = turbulencePorous[i];
-
-    volScalarField& porousP = porousThermo.p();
-    const volScalarField& porousPsi = porousThermo.psi();
-    volScalarField& porousH = porousThermo.he();
-
-    const dimensionedScalar initialMass
-    (
-        "initialMass",
-        dimMass,
-        initialMassFluidPorous[i]
-    );
-
-    IObasicSourceList& porousSources = heatPorousSources[i];
-
-    const label pRefCell = pRefCellFluidPorous[i];
-    const scalar pRefValue = pRefValueFluidPorous[i];
-
-    const scalar rhoMax = rhoMaxPorous[i].value();
-    const scalar rhoMin = rhoMinPorous[i].value();

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/solvePorousFluid.H

@@ -1,11 +0,0 @@
-//  Pressure-velocity SIMPLE corrector
-
-    porousP.storePrevIter();
-    porousRho.storePrevIter();
-    {
-        #include "UPorousFluidEqn.H"
-        #include "hPorousFluidEqn.H"
-        #include "pPorousFluidEqn.H"
-    }
-
-    turbPorous.correct();

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/createPorousSolidFields.H

@@ -1,43 +0,0 @@
-    // Initialise solid field pointer lists
-    PtrList<solidThermo> porousSolidThermos(porousSolidRegions.size());
-    PtrList<IObasicSourceList> solidHeatSources(porousSolidRegions.size());
-    PtrList<volScalarField> betavSolid(porousSolidRegions.size());
-
-    // Populate solid field pointer lists
-    forAll(porousSolidRegions, i)
-    {
-        Info<< "*** Reading porous solid mesh thermophysical "
-            << "properties for region "
-            << porousSolidRegions[i].name() << nl << endl;
-
-        Info<< "    Adding to thermos\n" << endl;
-        porousSolidThermos.set
-        (
-            i,
-            solidThermo::New(porousSolidRegions[i])
-        );
-        Info<< "    Adding sources\n" << endl;
-        solidHeatSources.set
-        (
-            i,
-            new IObasicSourceList(porousSolidRegions[i])
-        );
-
-        Info<< "    Adding to betavSolid\n" << endl;
-        betavSolid.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "betavSolid",
-                    runTime.timeName(),
-                    porousSolidRegions[i],
-                    IOobject::MUST_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                porousSolidRegions[i]
-            )
-        );
-    }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/createPorousSolidMeshes.H

@@ -1,24 +0,0 @@
-    const wordList porousSolidNames(rp["porousSolid"]);
-
-    PtrList<fvMesh> porousSolidRegions(porousSolidNames.size());
-
-    forAll(porousSolidNames, i)
-    {
-        Info<< "Create solid mesh for region " << porousSolidNames[i]
-            << " for time = " << runTime.timeName() << nl << endl;
-
-        porousSolidRegions.set
-        (
-            i,
-            new fvMesh
-            (
-                IOobject
-                (
-                    porousSolidNames[i],
-                    runTime.timeName(),
-                    runTime,
-                    IOobject::MUST_READ
-                )
-            )
-        );
-    }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/readPorousSolidMultiRegionSIMPLEControls.H

@@ -1,4 +0,0 @@
-    const dictionary& simple = mesh.solutionDict().subDict("SIMPLE");
-
-    const int nNonOrthCorr =
-        simple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H

@@ -1,23 +0,0 @@
-    const fvMesh& mesh = porousSolidRegions[i];
-
-    solidThermo& thermo = porousSolidThermos[i];
-    const volScalarField& betav = betavSolid[i];
-
-    tmp<volScalarField> trho = thermo.rho();
-    const volScalarField& rho = trho();
-
-    tmp<volScalarField> tcp = thermo.Cp();
-    const volScalarField& cp = tcp();
-
-    tmp<volScalarField> tkappa = thermo.kappa();
-    //tmp<volSymmTensorField> tkappa = thermo.directionalKappa()*betav;
-
-    const volScalarField& kappa = tkappa();
-    //const volSymmTensorField& K = tK();
-
-    tmp<volScalarField> talpha = thermo.alpha();
-    const volScalarField& alpha = talpha();
-
-    volScalarField& h = thermo.he();
-
-    IObasicSourceList& sources = solidHeatSources[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solidPorousRegionDiffusionNo.H

@@ -1,18 +0,0 @@
-    scalar DiNum = -GREAT;
-
-    forAll(solidRegions, i)
-    {
-#       include "setRegionSolidFields.H"
-
-        DiNum = max
-        (
-            solidRegionDiffNo
-            (
-                solidRegions[i],
-                runTime,
-                rho*cp,
-                K
-            ),
-            DiNum
-        );
-    }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H

@@ -1,17 +0,0 @@
-{
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-    {
-        tmp<fvScalarMatrix> hEqn
-        (
-          - fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
-          + sources(rho, h)
-        );
-
-        hEqn().relax();
-        hEqn().solve();
-    }
-}
-
-thermo.correct();
-
-Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T()) << endl;

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/createSolidFields.H

@@ -1,20 +0,0 @@
-    // Initialise solid field pointer lists
-    PtrList<solidThermo> thermos(solidRegions.size());
-    PtrList<radiation::radiationModel> radiations(solidRegions.size());
-
-    // Populate solid field pointer lists
-    forAll(solidRegions, i)
-    {
-        Info<< "*** Reading solid mesh thermophysical properties for region "
-            << solidRegions[i].name() << nl << endl;
-
-        Info<< "    Adding to thermos\n" << endl;
-        thermos.set
-        (
-            i,
-            solidThermo::New(solidRegions[i])
-        );
-
-        Info<< "    Adding to radiations\n" << endl;
-        radiations.set(i, radiation::radiationModel::New(thermos[i].T()));
-    }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/createSolidMeshes.H

@@ -1,24 +0,0 @@
-    const wordList solidsNames(rp["solid"]);
-
-    PtrList<fvMesh> solidRegions(solidsNames.size());
-
-    forAll(solidsNames, i)
-    {
-        Info<< "Create solid mesh for region " << solidsNames[i]
-            << " for time = " << runTime.timeName() << nl << endl;
-
-        solidRegions.set
-        (
-            i,
-            new fvMesh
-            (
-                IOobject
-                (
-                    solidsNames[i],
-                    runTime.timeName(),
-                    runTime,
-                    IOobject::MUST_READ
-                )
-            )
-        );
-    }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/setRegionSolidFields.H

@@ -1,18 +0,0 @@
-    fvMesh& mesh = solidRegions[i];
-    solidThermo& thermo = thermos[i];
-    const radiation::radiationModel& radiation = radiations[i];
-
-    tmp<volScalarField> trho = thermo.rho();
-    const volScalarField& rho = trho();
-
-    tmp<volScalarField> tcp = thermo.Cp();
-    const volScalarField& cp = tcp();
-
-    tmp<volScalarField> tkappa = thermo.kappa();
-    //tmp<volSymmTensorField> tkappa = thermo.directionalkappa();
-    const volScalarField& kappa = tkappa();
-
-    tmp<volScalarField> talpha = thermo.alpha();
-    const volScalarField& alpha = talpha();
-
-    volScalarField& h = thermo.he();

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H

@@ -3,7 +3,8 @@
     {
         fvScalarMatrix hEqn
         (
-            -fvm::laplacian(alpha, h)
+           - fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
+           + sources(rho, h)
         );
         hEqn.relax();
         hEqn.solve();

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

@@ -14,7 +14,7 @@
             CoNum
         );
     }
-
+/*
     forAll (porousFluidRegions, porousI)
     {
         CoNum = max
@@ -29,3 +29,4 @@
             CoNum
         );
     }
+*/

applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/UPorousFluidEqn.H

@@ -1,24 +0,0 @@
-
-// Solve the Momentum equation
-    tmp<fvVectorMatrix> porousUEqn
-    (
-        fvm::ddt(porousRho, porousU)
-      + fvm::div(porousPhi, porousU)
-      + turbPorous.divDevRhoReff(porousU)
-      + porousSources(porousRho, porousU)
-    );
-
-    porousUEqn().relax();
-
-    volScalarField rAUPorous(1.0/porousUEqn().A());
-
-    if (momentumPredictor)
-    {
-        solve(porousUEqn() == -fvc::grad(porousP));
-    }
-    else
-    {
-        porousU = rAUPorous*(porousUEqn().H() - fvc::grad(porousP));
-        porousU.correctBoundaryConditions();
-    }
-

applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/createPorousFluidFields.H

@@ -1,147 +0,0 @@
-    // Initialise porous field pointer lists
-    PtrList<rhoThermo> thermoPorous(porousFluidRegions.size());
-    PtrList<volScalarField> rhoPorous(porousFluidRegions.size());
-    PtrList<volVectorField> UPorous(porousFluidRegions.size());
-    PtrList<surfaceScalarField> phiPorous(porousFluidRegions.size());
-    PtrList<volScalarField> KPorous(porousFluidRegions.size());
-    PtrList<volScalarField> dpdtPorous(fluidRegions.size());
-
-    PtrList<compressible::turbulenceModel> turbulencePorous
-    (
-        porousFluidRegions.size()
-    );
-    PtrList<volScalarField> pPorous(porousFluidRegions.size());
-
-    PtrList<IObasicSourceList> heatPorousSources
-    (
-        porousFluidRegions.size()
-    );
-
-    forAll(porousFluidRegions, i)
-    {
-        Info<< "Reading fluid mesh thermophysical properties for porous "
-            << porousFluidRegions[i].name() << nl << endl;
-
-        Info<< "    Adding to thermoFluid porous\n" << endl;
-
-        thermoPorous.set
-        (
-            i,
-            rhoThermo::New(porousFluidRegions[i]).ptr()
-        );
-
-        Info<< "    Adding to rhoPorous\n" << endl;
-        rhoPorous.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "rho",
-                    runTime.timeName(),
-                    porousFluidRegions[i],
-                    IOobject::NO_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                thermoPorous[i].rho()
-            )
-        );
-
-        Info<< "    Adding to UPorous\n" << endl;
-        UPorous.set
-        (
-            i,
-            new volVectorField
-            (
-                IOobject
-                (
-                    "U",
-                    runTime.timeName(),
-                    porousFluidRegions[i],
-                    IOobject::MUST_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                porousFluidRegions[i]
-            )
-        );
-
-        Info<< "    Adding to phiPorous\n" << endl;
-        phiPorous.set
-        (
-            i,
-            new surfaceScalarField
-            (
-                IOobject
-                (
-                    "phi",
-                    runTime.timeName(),
-                    porousFluidRegions[i],
-                    IOobject::READ_IF_PRESENT,
-                    IOobject::AUTO_WRITE
-                ),
-                linearInterpolate(rhoPorous[i]*UPorous[i])
-                    & porousFluidRegions[i].Sf()
-            )
-        );
-
-        Info<< "    Adding turbulence to porous\n" << endl;
-        turbulencePorous.set
-        (
-            i,
-            compressible::turbulenceModel::New
-            (
-                rhoPorous[i],
-                UPorous[i],
-                phiPorous[i],
-                thermoPorous[i]
-            ).ptr()
-        );
-
-        Info<< "    Adding to KPorous\n" << endl;
-        KPorous.set
-        (
-            i,
-            new volScalarField
-            (
-                "KPorous",
-                0.5*magSqr(UPorous[i])
-            )
-        );
-
-        Info<< "    Adding to dpdtPorous\n" << endl;
-        dpdtPorous.set
-        (
-            i,
-            new volScalarField
-            (
-                "dpdtPorous",
-                fvc::ddt(thermoPorous[i].p())
-            )
-        );
-
-        pPorous.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "p",
-                    runTime.timeName(),
-                    porousFluidRegions[i],
-                    IOobject::MUST_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                porousFluidRegions[i]
-            )
-        );
-
-        heatPorousSources.set
-        (
-            i,
-            new IObasicSourceList(porousFluidRegions[i])
-        );
-    }
-
-

applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/createPorousFluidRegions.H

@@ -1,25 +0,0 @@
-    const wordList porousFluidNames(rp["porousFluid"]);
-
-    PtrList<fvMesh> porousFluidRegions(porousFluidNames.size());
-
-    forAll (porousFluidNames, iPorous)
-    {
-        const word porousFluidName = porousFluidNames[iPorous];
-        Info<< "Create porous fluid region " << porousFluidName
-            << nl << endl;
-
-        porousFluidRegions.set
-        (
-            iPorous,
-            new fvMesh
-            (
-                IOobject
-                (
-                    porousFluidName,
-                    runTime.timeName(),
-                    runTime,
-                    IOobject::MUST_READ
-                )
-            )
-        );
-    }

applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/hPorousFluidEqn.H

@@ -1,21 +0,0 @@
-{
-    fvScalarMatrix hPorousEqn
-    (
-        fvm::ddt(porousRho, porousH)
-      + fvm::div(porousPhi, porousH)
-      - fvm::laplacian(turbPorous.alphaEff(), porousH)
-     ==
-        porousdpdt
-      - (fvc::ddt(porousRho, porousK) + fvc::div(porousPhi, porousK))
-      + porousSources(porousRho, porousH)
-    );
-
-    hPorousEqn.relax();
-    hPorousEqn.solve();
-
-    porousThermo.correct();
-
-    Info<< "Min/max in the porous T:"
-        << min(porousThermo.T()).value() << ' '
-        << max(porousThermo.T()).value() << endl;
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/pPorousFluidEqn.H

@@ -1,63 +0,0 @@
-porousRho = porousThermo.rho();
-
-porousU = rAUPorous*porousUEqn().H();
-
-if (nCorr <= 1)
-{
-    porousUEqn.clear();
-}
-
-porousPhi =
-    fvc::interpolate(porousRho)*
-    (
-        (fvc::interpolate(porousU) & porousMesh.Sf())
-        + fvc::ddtPhiCorr(rAUPorous, porousRho, porousU, porousPhi)
-    );
-
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-{
-    // Pressure corrector
-    fvScalarMatrix pEqn
-    (
-        fvm::ddt(porousPsi, porousP)
-      + fvc::div(porousPhi)
-      - fvm::laplacian(porousRho*rAUPorous, porousP)
-    );
-
-    pEqn.solve
-    (
-        porousMesh.solver
-        (
-            porousP.select
-            (
-                oCorr == nOuterCorr-1
-                && corr == nCorr-1
-                && nonOrth == nNonOrthCorr
-            )
-        )
-    );
-
-    if (nonOrth == nNonOrthCorr)
-    {
-        porousPhi += pEqn.flux();
-    }
-}
-
-solve(fvm::ddt(porousRho) + fvc::div(porousPhi));
-
-// Explicitly relax pressure for momentum corrector
-porousP.relax();
-
-// Recalculate density from the relaxed pressure
-porousRho = porousThermo.rho();
-
-
-porousU -= rAUPorous*fvc::grad(porousP);
-porousU.correctBoundaryConditions();
-porousK = 0.5*magSqr(porousU);
-
-// Update pressure time derivative if needed
-if (porousThermo.dpdt())
-{
-    porousdpdt = fvc::ddt(porousP);
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/readPorousFluidRegionPIMPLEControls.H

@@ -1,11 +0,0 @@
-    const dictionary& pimple = porousMesh.solutionDict().subDict("PIMPLE");
-
-    const int nCorr =
-        pimple.lookupOrDefault<int>("nCorrectors", 1);
-
-    const int nNonOrthCorr =
-        pimple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-
-    const bool momentumPredictor =
-        pimple.lookupOrDefault("momentumPredictor", true);
-

applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/setPorousFluidFields.H

@@ -1,17 +0,0 @@
-    fvMesh& porousMesh = porousFluidRegions[i];
-
-    rhoThermo& porousThermo = thermoPorous[i];
-    volScalarField& porousRho = rhoPorous[i];
-
-    volVectorField& porousU = UPorous[i];
-    surfaceScalarField& porousPhi = phiPorous[i];
-
-    compressible::turbulenceModel& turbPorous = turbulencePorous[i];
-    volScalarField& porousK = KPorous[i];
-    volScalarField& porousdpdt = dpdtPorous[i];
-
-    volScalarField& porousP = porousThermo.p();
-    const volScalarField& porousPsi = porousThermo.psi();
-    volScalarField& porousH = porousThermo.he();
-
-    IObasicSourceList& porousSources = heatPorousSources[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/solvePorousFluid.H

@@ -1,28 +0,0 @@
-if (finalIter)
-{
-    porousMesh.data::add("finalIteration", true);
-}
-
-if (oCorr == 0)
-{
-    solve(fvm::ddt(porousRho) + fvc::div(porousPhi));
-}
-
-#include "UPorousFluidEqn.H"
-
-#include "hPorousFluidEqn.H"
-
-// --- PISO loop
-for (int corr=0; corr<nCorr; corr++)
-{
-    #include "pPorousFluidEqn.H"
-}
-
-turbPorous.correct();
-
-porousRho = porousThermo.rho();
-
-if (finalIter)
-{
-    porousMesh.data::remove("finalIteration");
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/createPorousSolidFields.H

@@ -1,43 +0,0 @@
-    // Initialise solid field pointer lists
-    PtrList<solidThermo> porousSolidThermos(porousSolidRegions.size());
-    PtrList<IObasicSourceList> solidHeatSources(porousSolidRegions.size());
-    PtrList<volScalarField> betavSolid(porousSolidRegions.size());
-
-    // Populate solid field pointer lists
-    forAll(porousSolidRegions, i)
-    {
-        Info<< "*** Reading porous solid mesh thermophysical "
-            << "properties for region "
-            << porousSolidRegions[i].name() << nl << endl;
-
-        Info<< "    Adding to thermos\n" << endl;
-        porousSolidThermos.set
-        (
-            i,
-            solidThermo::New(porousSolidRegions[i])
-        );
-        Info<< "    Adding sources\n" << endl;
-        solidHeatSources.set
-        (
-            i,
-            new IObasicSourceList(porousSolidRegions[i])
-        );
-
-        Info<< "    Adding to betavSolid\n" << endl;
-        betavSolid.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "betavSolid",
-                    runTime.timeName(),
-                    porousSolidRegions[i],
-                    IOobject::MUST_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                porousSolidRegions[i]
-            )
-        );
-    }

applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/createPorousSolidMeshes.H

@@ -1,24 +0,0 @@
-    const wordList porousSolidNames(rp["porousSolid"]);
-
-    PtrList<fvMesh> porousSolidRegions(porousSolidNames.size());
-
-    forAll(porousSolidNames, i)
-    {
-        Info<< "Create solid mesh for region " << porousSolidNames[i]
-            << " for time = " << runTime.timeName() << nl << endl;
-
-        porousSolidRegions.set
-        (
-            i,
-            new fvMesh
-            (
-                IOobject
-                (
-                    porousSolidNames[i],
-                    runTime.timeName(),
-                    runTime,
-                    IOobject::MUST_READ
-                )
-            )
-        );
-    }

applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/readPorousSolidMultiRegionPIMPLEControls.H

@@ -1,4 +0,0 @@
-    const dictionary& pimple = mesh.solutionDict().subDict("PIMPLE");
-
-    int nNonOrthCorr =
-        pimple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);

applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H

@@ -1,26 +0,0 @@
-    fvMesh& mesh = porousSolidRegions[i];
-
-    solidThermo& thermo = porousSolidThermos[i];
-    const volScalarField& betav = betavSolid[i];
-
-    tmp<volScalarField> trho = thermo.rho();
-    const volScalarField& rho = trho();
-
-    tmp<volScalarField> tcp = thermo.Cp();
-    const volScalarField& cp = tcp();
-
-    tmp<volScalarField> tkappa = thermo.kappa();
-    //tmp<volSymmTensorField> tkappa = thermo.directionalKappa()*betav;
-
-    const volScalarField& kappa = tkappa();
-    //const volSymmTensorField& K = tK();
-
-    //tmp<volScalarField> trhoCp = cp*rho;
-    //const volScalarField& rhoCp = trhoCp();
-
-    tmp<volScalarField> talpha = thermo.alpha();
-    const volScalarField& alpha = talpha();
-
-    volScalarField& h = thermo.he();
-
-    IObasicSourceList& sources = solidHeatSources[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solidPorousRegionDiffusionNo.H

@@ -1,18 +0,0 @@
-    scalar DiNum = -GREAT;
-
-    forAll(solidRegions, i)
-    {
-#       include "setRegionSolidFields.H"
-
-        DiNum = max
-        (
-            solidRegionDiffNo
-            (
-                solidRegions[i],
-                runTime,
-                rho*cp,
-                K
-            ),
-            DiNum
-        );
-    }

applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solvePorousSolid.H

@@ -1,28 +0,0 @@
-if (finalIter)
-{
-    mesh.data::add("finalIteration", true);
-}
-
-{
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-    {
-        tmp<fvScalarMatrix> hEqn
-        (
-            fvm::ddt(betav*rho, h)
-          - fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
-          + sources(rho, h)
-        );
-
-        hEqn().relax();
-        hEqn().solve(mesh.solver(h.select(finalIter)));
-    }
-}
-
-thermo.correct();
-
-Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T()) << endl;
-
-if (finalIter)
-{
-    mesh.data::remove("finalIteration");
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H

@@ -1,6 +1,8 @@
     // Initialise solid field pointer lists
     PtrList<solidThermo> thermos(solidRegions.size());
     PtrList<radiation::radiationModel> radiations(solidRegions.size());
+    PtrList<IObasicSourceList> solidHeatSources(solidRegions.size());
+    PtrList<volScalarField> betavSolid(solidRegions.size());
 
     // Populate solid field pointer lists
     forAll(solidRegions, i)
@@ -13,4 +15,49 @@
 
         Info<< "    Adding to radiations\n" << endl;
         radiations.set(i, radiation::radiationModel::New(thermos[i].T()));
+
+        Info<< "    Adding sources\n" << endl;
+        solidHeatSources.set
+        (
+            i,
+            new IObasicSourceList(solidRegions[i])
+        );
+
+        IOobject betavSolidIO
+        (
+            "betavSolid",
+            runTime.timeName(),
+            solidRegions[i],
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        );
+
+        if (betavSolidIO.headerOk())
+        {
+            betavSolid.set
+            (
+                i,
+                new volScalarField(betavSolidIO, solidRegions[i])
+            );
+        }
+        else
+        {
+            betavSolid.set
+            (
+                i,
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        "betavSolid",
+                        runTime.timeName(),
+                        solidRegions[i],
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    solidRegions[i],
+                    dimensionedScalar("1", dimless, scalar(1.0))
+                )
+            );
+        }
     }

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H

@@ -14,3 +14,7 @@
     const volScalarField& kappa = tkappa();
 
     volScalarField& h = thermo.he();
+
+    const volScalarField& betav = betavSolid[i];
+
+    IObasicSourceList& sources = solidHeatSources[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H

@@ -8,8 +8,9 @@ if (finalIter)
     {
         tmp<fvScalarMatrix> hEqn
         (
-            fvm::ddt(rho, h)
+            fvm::ddt(betav*rho, h)
-          - fvm::laplacian(alpha, h)
+          - fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
+          + sources(rho, h)
         );
         hEqn().relax();
         hEqn().solve(mesh.solver(h.select(finalIter)));

applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/MRFSimpleFoam.C

@@ -33,7 +33,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
-#include "MRFZones.H"
+#include "IOMRFZoneList.H"
 #include "simpleControl.H"
 #include "IObasicSourceList.H"
 
@@ -48,7 +48,7 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
-    MRFZones mrfZones(mesh);
+    IOMRFZoneList mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U);
 
     simpleControl simple(mesh);

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C

@@ -46,7 +46,7 @@ Description
 
 #include "pimpleControl.H"
 
-#include "MRFZones.H"
+#include "IOMRFZoneList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createMRFZones.H

@@ -1,4 +1,4 @@
-    MRFZones mrfZones(mesh);
+    IOMRFZoneList mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U1);
     mrfZones.correctBoundaryVelocity(U2);
     mrfZones.correctBoundaryVelocity(U);

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

@@ -42,7 +42,7 @@ Description
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
-#include "MRFZones.H"
+#include "IOMRFZoneList.H"
 #include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/multiphase/interFoam/MRFInterFoam/createMRFZones.H

@@ -1,2 +1,2 @@
-    MRFZones mrfZones(mesh);
+    IOMRFZoneList mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U);

applications/solvers/multiphase/multiphaseEulerFoam/createMRFZones.H

@@ -1,4 +1,4 @@
-    MRFZones mrfZones(mesh);
+    IOMRFZoneList mrfZones(mesh);
 
     forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
     {

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,7 +40,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "LESModel.H"
 
-#include "MRFZones.H"
+#include "IOMRFZoneList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C

@@ -35,7 +35,7 @@ Description
 #include "fvCFD.H"
 #include "multiphaseMixture.H"
 #include "turbulenceModel.H"
-#include "MRFZones.H"
+#include "IOMRFZoneList.H"
 #include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/multiphase/twoPhaseEulerFoam/createMRFZones.H

@@ -1,4 +1,4 @@
-    MRFZones mrfZones(mesh);
+    IOMRFZoneList mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U1);
     mrfZones.correctBoundaryVelocity(U2);
     mrfZones.correctBoundaryVelocity(U);

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C

@@ -46,7 +46,7 @@ Description
 #include "kineticTheoryModel.H"
 
 #include "pimpleControl.H"
-#include "MRFZones.H"
+#include "IOMRFZoneList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -226,37 +226,37 @@ castellatedMeshControls
 // Settings for the snapping.
 snapControls
 {
-    //- Number of patch smoothing iterations before finding correspondence
+    // Number of patch smoothing iterations before finding correspondence
-    //  to surface
+    // to surface
     nSmoothPatch 3;
 
-    //- Maximum relative distance for points to be attracted by surface.
+    // Maximum relative distance for points to be attracted by surface.
-    //  True distance is this factor times local maximum edge length.
+    // True distance is this factor times local maximum edge length.
     // Note: changed(corrected) w.r.t 17x! (17x used 2* tolerance)
     tolerance 2.0;
 
-    //- Number of mesh displacement relaxation iterations.
+    // Number of mesh displacement relaxation iterations.
     nSolveIter 30;
 
-    //- Maximum number of snapping relaxation iterations. Should stop
+    // Maximum number of snapping relaxation iterations. Should stop
-    //  before upon reaching a correct mesh.
+    // before upon reaching a correct mesh.
     nRelaxIter 5;
 
     // Feature snapping
 
-        //- Number of feature edge snapping iterations.
+        // Number of feature edge snapping iterations.
-        //  Leave out altogether to disable.
+        // Leave out altogether to disable.
         nFeatureSnapIter 10;
 
-        //- Detect (geometric only) features by sampling the surface
+        // Detect (geometric only) features by sampling the surface
-        //  (default=false).
+        // (default=false).
         implicitFeatureSnap false;
 
-        //- Use castellatedMeshControls::features (default = true)
+        // Use castellatedMeshControls::features (default = true)
         explicitFeatureSnap true;
 
-        //- Detect features between multiple surfaces
+        // Detect features between multiple surfaces
-        //  (only for explicitFeatureSnap, default = false)
+        // (only for explicitFeatureSnap, default = false)
         multiRegionFeatureSnap false;
 }
 
@@ -267,9 +267,43 @@ addLayersControls
     // size of the refined cell outside layer (true) or absolute sizes (false).
     relativeSizes true;
 
+    // Layer thickness specification. This can be specified in one of four ways
+    // - expansionRatio and finalLayerThickness (cell nearest internal mesh)
+    // - expansionRatio and firstLayerThickness (cell on surface)
+    // - overall thickness and firstLayerThickness
+    // - overall thickness and finalLayerThickness
+
+        // Expansion factor for layer mesh
+        expansionRatio 1.0;
+
+        // Wanted thickness of the layer furthest away from the wall.
+        // If relativeSizes this is relative to undistorted size of cell
+        // outside layer.
+        finalLayerThickness 0.3;
+
+        // Wanted thickness of the layer next to the wall.
+        // If relativeSizes this is relative to undistorted size of cell
+        // outside layer.
+        //firstLayerThickness 0.3;
+
+        // Wanted overall thickness of layers.
+        // If relativeSizes this is relative to undistorted size of cell
+        // outside layer.
+        //thickness 0.5
+
+
+    // Minimum overall thickness of total layers. If for any reason layer
+    // cannot be above minThickness do not add layer.
+    // If relativeSizes this is relative to undistorted size of cell
+    // outside layer..
+    minThickness 0.25;
+
+
     // Per final patch (so not geometry!) the layer information
     // Note: This behaviour changed after 21x. Any non-mentioned patches
-    //       now slide unless nSurfaceLayers is explicitly mentioned to be 0.
+    //       now slide unless:
+    //          - nSurfaceLayers is explicitly mentioned to be 0.
+    //          - angle to nearest surface < slipFeatureAngle (see below)
     layers
     {
         sphere.stl_firstSolid
@@ -281,9 +315,9 @@ addLayersControls
         {
             nSurfaceLayers 1;
             // Per patch layer data
-            expansionRatio 1.3;
+            expansionRatio      1.3;
             finalLayerThickness 0.3;
-            minThickness 0.1;
+            minThickness        0.1;
         }
 
         // Disable any mesh shrinking and layer addition on any point of
@@ -294,41 +328,24 @@ addLayersControls
         }
     }
 
-    // Expansion factor for layer mesh
+    // If points get not extruded do nGrow layers of connected faces that are
-    expansionRatio 1.0;
+    // also not grown. This helps convergence of the layer addition process
-
+    // close to features.
-    //- Wanted thickness of final added cell layer. If multiple layers
-    //  is the
-    //  thickness of the layer furthest away from the wall.
-    //  Relative to undistorted size of cell outside layer.
-    //  is the thickness of the layer furthest away from the wall.
-    //  See relativeSizes parameter.
-    finalLayerThickness 0.3;
-
-    //- Minimum thickness of cell layer. If for any reason layer
-    //  cannot be above minThickness do not add layer.
-    //  Relative to undistorted size of cell outside layer.
-    //  See relativeSizes parameter.
-    minThickness 0.25;
-
-    //- If points get not extruded do nGrow layers of connected faces that are
-    //  also not grown. This helps convergence of the layer addition process
-    //  close to features.
     // Note: changed(corrected) w.r.t 17x! (didn't do anything in 17x)
     nGrow 0;
 
     // Advanced settings
 
-    //- When not to extrude surface. 0 is flat surface, 90 is when two faces
+    // When not to extrude surface. 0 is flat surface, 90 is when two faces
-    //  make straight angle.
+    // make straight angle.
     featureAngle 60;
 
-    //- At non-patched sides allow mesh to slip if extrusion direction makes
+    // At non-patched sides allow mesh to slip if extrusion direction makes
-    //  angle larger than slipFeatureAngle.
+    // angle larger than slipFeatureAngle.
     slipFeatureAngle 30;
 
-    //- Maximum number of snapping relaxation iterations. Should stop
+    // Maximum number of snapping relaxation iterations. Should stop
-    //  before upon reaching a correct mesh.
+    // before upon reaching a correct mesh.
     nRelaxIter 5;
 
     // Number of smoothing iterations of surface normals
@@ -374,51 +391,51 @@ addLayersControls
 // where to undo.
 meshQualityControls
 {
-    //- Maximum non-orthogonality allowed. Set to 180 to disable.
+    // Maximum non-orthogonality allowed. Set to 180 to disable.
     maxNonOrtho 65;
 
-    //- Max skewness allowed. Set to <0 to disable.
+    // Max skewness allowed. Set to <0 to disable.
     maxBoundarySkewness 20;
     maxInternalSkewness 4;
 
-    //- Max concaveness allowed. Is angle (in degrees) below which concavity
+    // Max concaveness allowed. Is angle (in degrees) below which concavity
-    //  is allowed. 0 is straight face, <0 would be convex face.
+    // is allowed. 0 is straight face, <0 would be convex face.
-    //  Set to 180 to disable.
+    // Set to 180 to disable.
     maxConcave 80;
 
-    //- Minimum pyramid volume. Is absolute volume of cell pyramid.
+    // Minimum pyramid volume. Is absolute volume of cell pyramid.
-    //  Set to a sensible fraction of the smallest cell volume expected.
+    // Set to a sensible fraction of the smallest cell volume expected.
-    //  Set to very negative number (e.g. -1E30) to disable.
+    // Set to very negative number (e.g. -1E30) to disable.
     minVol 1e-13;
 
-    //- Minimum quality of the tet formed by the face-centre
+    // Minimum quality of the tet formed by the face-centre
-    //  and variable base point minimum decomposition triangles and
+    // and variable base point minimum decomposition triangles and
-    //  the cell centre. This has to be a positive number for tracking
+    // the cell centre. This has to be a positive number for tracking
-    //  to work. Set to very negative number (e.g. -1E30) to
+    // to work. Set to very negative number (e.g. -1E30) to
-    //  disable.
+    // disable.
-    //     <0 = inside out tet,
+    //    <0 = inside out tet,
-    //      0 = flat tet
+    //     0 = flat tet
-    //      1 = regular tet
+    //     1 = regular tet
     minTetQuality 1e-9;
 
-    //- Minimum face area. Set to <0 to disable.
+    // Minimum face area. Set to <0 to disable.
     minArea -1;
 
-    //- Minimum face twist. Set to <-1 to disable. dot product of face normal
+    // Minimum face twist. Set to <-1 to disable. dot product of face normal
-    //- and face centre triangles normal
+    // and face centre triangles normal
     minTwist 0.05;
 
-    //- minimum normalised cell determinant
+    // minimum normalised cell determinant
-    //- 1 = hex, <= 0 = folded or flattened illegal cell
+    // 1 = hex, <= 0 = folded or flattened illegal cell
     minDeterminant 0.001;
 
-    //- minFaceWeight (0 -> 0.5)
+    // minFaceWeight (0 -> 0.5)
     minFaceWeight 0.05;
 
-    //- minVolRatio (0 -> 1)
+    // minVolRatio (0 -> 1)
     minVolRatio 0.01;
 
-    //must be >0 for Fluent compatibility
+    // must be >0 for Fluent compatibility
     minTriangleTwist -1;
 
     //- if >0 : preserve single cells with all points on the surface if the
@@ -429,9 +446,9 @@ meshQualityControls
 
     // Advanced
 
-    //- Number of error distribution iterations
+    // Number of error distribution iterations
     nSmoothScale 4;
-    //- amount to scale back displacement at error points
+    // amount to scale back displacement at error points
     errorReduction 0.75;
 
     // Optional : some meshing phases allow usage of relaxed rules.

applications/utilities/postProcessing/sampling/sample/sampleDict

@@ -90,6 +90,8 @@ fields
 //      uniform: extra number of sampling points
 //      polyLine, cloud: list of coordinates
 //      patchCloud: list of coordinates and set of patches to look for nearest
+//      patchSeed: random sampling on set of patches. Points slightly off
+//                 face centre.
 sets
 (
     lineX1
@@ -135,6 +137,15 @@ sets
         maxDistance 0.1;    // maximum distance to search
         patches     (".*Wall.*");
     }
+
+    patchSeed
+    {
+        patches     (".*Wall.*");
+        // Number of points to seed. Divided amongst all processors according
+        // to fraction of patches they hold.
+        maxPoints   100;
+    }
+
 );
 
 

src/fieldSources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H

@@ -27,7 +27,6 @@ Class
 Description
     Base class for inter region heat exchange. The derived classes must
     provide the heat transfer coefficient (htc)
-    NOTE: mapToMap does to work in paralell
 
 \*---------------------------------------------------------------------------*/
 
@@ -44,7 +43,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                       Class interRegionHeatTransferModel Declaration
+                Class interRegionHeatTransferModel Declaration
 \*---------------------------------------------------------------------------*/
 
 class interRegionHeatTransferModel

src/finiteVolume/Make/files

@@ -379,7 +379,8 @@ $(porosity)/fixedCoeff/fixedCoeff.C
 
 MRF = $(general)/MRF
 $(MRF)/MRFZone.C
-$(MRF)/MRFZones.C
+$(MRF)/MRFZoneList.C
+$(MRF)/IOMRFZoneList.C
 
 SRF = $(general)/SRF
 $(SRF)/SRFModel/SRFModel/SRFModel.C

src/finiteVolume/cfdTools/general/MRF/IOMRFZoneList.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,46 +23,68 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "pureSolidMixture.H"
+#include "IOMRFZoneList.H"
 #include "fvMesh.H"
+#include "Time.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-namespace Foam
+Foam::IOobject Foam::IOMRFZoneList::createIOobject
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class ThermoType>
-pureSolidMixture<ThermoType>::pureSolidMixture
 (
-    const dictionary& thermoDict,
     const fvMesh& mesh
-)
+) const
-:
-    basicMixture(thermoDict, mesh),
-    mixture_(thermoDict.subDict("mixture"))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class ThermoType>
-pureSolidMixture<ThermoType>::~pureSolidMixture()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class ThermoType>
-void pureSolidMixture<ThermoType>::read(const dictionary& thermoDict)
 {
-    mixture_ = ThermoType(thermoDict.subDict("mixture"));
+    IOobject io
+    (
+        "MRFProperties",
+        mesh.time().constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (io.headerOk())
+    {
+        Info<< "Creating MRF zone list from " << io.name() << endl;
+
+        io.readOpt() = IOobject::MUST_READ_IF_MODIFIED;
+        return io;
+    }
+    else
+    {
+        Info<< "No MRF models present" << nl << endl;
+
+        io.readOpt() = IOobject::NO_READ;
+        return io;
+    }
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::IOMRFZoneList::IOMRFZoneList
+(
+    const fvMesh& mesh
+)
+:
+    IOdictionary(createIOobject(mesh)),
+    MRFZoneList(mesh, *this)
+{}
+
+
+bool Foam::IOMRFZoneList::read()
+{
+    if (regIOobject::read())
+    {
+        MRFZoneList::read(*this);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
 
-} // End namespace Foam
 
 // ************************************************************************* //
+

src/finiteVolume/cfdTools/general/MRF/IOMRFZoneList.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,21 +22,38 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::reactingSolidMixture
+    Foam::IOMRFZoneList
 
 Description
-    Foam::reactingSolidMixture
+    List of MRF zones with IO functionality.  MRF zones are specified by a list
+    of dictionary entries, e.g.
+
+    \verbatim
+    zone1
+    {
+        cellZone    rotor1;
+        active      yes;
+        ...
+    }
+
+    zone2
+    {
+        cellZone    rotor2;
+        active      yes;
+        ...
+    }
+    \endverbatim
 
 SourceFiles
-    reactingSolidMixture.C
+    IOMRFZoneList.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef reactingSolidMixture_H
+#ifndef IOMRFZoneList_H
-#define reactingSolidMixture_H
+#define IOMRFZoneList_H
 
-#include "multiComponentSolidMixture.H"
+#include "IOdictionary.H"
-#include "solidReaction.H"
+#include "MRFZoneList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,60 +61,53 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                    Class reactingSolidMixture Declaration
+                        Class IOMRFZoneList Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class ThermoSolidType>
+class IOMRFZoneList
-class reactingSolidMixture
 :
-    public multiComponentSolidMixture<ThermoSolidType>,
+    public IOdictionary,
-    public PtrList<solidReaction>
+    public MRFZoneList
 {
+private:
+
     // Private Member Functions
 
+        //- Create IO object if dictionary is present
+        IOobject createIOobject(const fvMesh& mesh) const;
+
         //- Disallow default bitwise copy construct
-        reactingSolidMixture(const reactingSolidMixture&);
+        IOMRFZoneList(const IOMRFZoneList&);
 
         //- Disallow default bitwise assignment
-        void operator=(const reactingSolidMixture&);
+        void operator=(const IOMRFZoneList&);
 
 
 public:
 
-    //- The type of thermo package this mixture is instantiated for
-    typedef ThermoSolidType thermoType;
-
-
     // Constructors
 
-        //- Construct from dictionary and mesh
+        //- Construct from mesh
-        reactingSolidMixture(const dictionary&, const fvMesh&);
+        IOMRFZoneList(const fvMesh& mesh);
 
 
-    //- Destructor
+        //- Destructor
-    virtual ~reactingSolidMixture()
+        virtual ~IOMRFZoneList()
-    {}
+        {}
 
 
-    // Member functions
+    // Member Functions
 
         //- Read dictionary
-        void read(const dictionary&);
+        virtual bool read();
 };
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-#   include "reactingSolidMixture.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/MRF/MRFZone.C

@@ -28,11 +28,9 @@ License
 #include "volFields.H"
 #include "surfaceFields.H"
 #include "fvMatrices.H"
-#include "syncTools.H"
 #include "faceSet.H"
 #include "geometricOneField.H"
 
-
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 defineTypeNameAndDebug(Foam::MRFZone, 0);
@@ -144,7 +142,7 @@ void Foam::MRFZone::setMRFFaces()
 
         forAll(pp, patchFacei)
         {
-            label faceI = pp.start()+patchFacei;
+            label faceI = pp.start() + patchFacei;
 
             if (faceType[faceI] == 1)
             {
@@ -173,7 +171,7 @@ void Foam::MRFZone::setMRFFaces()
 
         forAll(pp, patchFacei)
         {
-            label faceI = pp.start()+patchFacei;
+            label faceI = pp.start() + patchFacei;
 
             if (faceType[faceI] == 1)
             {
@@ -228,64 +226,71 @@ void Foam::MRFZone::setMRFFaces()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
+Foam::MRFZone::MRFZone
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
 :
     mesh_(mesh),
-    name_(is),
+    name_(name),
-    dict_(is),
+    coeffs_(dict),
-    cellZoneID_(mesh_.cellZones().findZoneID(name_)),
+    active_(readBool(coeffs_.lookup("active"))),
+    cellZoneName_(coeffs_.lookup("cellZone")),
+    cellZoneID_(mesh_.cellZones().findZoneID(cellZoneName_)),
     excludedPatchNames_
     (
-        dict_.lookupOrDefault("nonRotatingPatches", wordList(0))
+        coeffs_.lookupOrDefault("nonRotatingPatches", wordList(0))
     ),
-    origin_(dict_.lookup("origin")),
+    origin_(coeffs_.lookup("origin")),
-    axis_(dict_.lookup("axis")),
+    axis_(coeffs_.lookup("axis")),
-    omega_(DataEntry<scalar>::New("omega", dict_))
+    omega_(DataEntry<scalar>::New("omega", coeffs_))
 {
-    if (dict_.found("patches"))
+    if (!active_)
     {
-        WarningIn("MRFZone(const fvMesh&, Istream&)")
+        cellZoneID_ = -1;
-            << "Ignoring entry 'patches'\n"
-            << "    By default all patches within the rotating region rotate.\n"
-            << "    Optionally supply excluded patches "
-            << "using 'nonRotatingPatches'."
-            << endl;
     }
-
+    else
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-    axis_ = axis_/mag(axis_);
-
-    excludedPatchLabels_.setSize(excludedPatchNames_.size());
-
-    forAll(excludedPatchNames_, i)
     {
-        excludedPatchLabels_[i] = patches.findPatchID(excludedPatchNames_[i]);
+        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-        if (excludedPatchLabels_[i] == -1)
+        axis_ = axis_/mag(axis_);
+
+        excludedPatchLabels_.setSize(excludedPatchNames_.size());
+
+        forAll(excludedPatchNames_, i)
+        {
+            excludedPatchLabels_[i] =
+                patches.findPatchID(excludedPatchNames_[i]);
+
+            if (excludedPatchLabels_[i] == -1)
+            {
+                FatalErrorIn
+                (
+                    "MRFZone(const word&, const fvMesh&, const dictionary&)"
+                )
+                    << "cannot find MRF patch " << excludedPatchNames_[i]
+                    << exit(FatalError);
+            }
+        }
+
+        bool cellZoneFound = (cellZoneID_ != -1);
+
+        reduce(cellZoneFound, orOp<bool>());
+
+        if (!cellZoneFound)
         {
             FatalErrorIn
             (
-                "Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
+                "MRFZone(const word&, const fvMesh&, const dictionary&)"
-            )   << "cannot find MRF patch " << excludedPatchNames_[i]
+            )
+                << "cannot find MRF cellZone " << cellZoneName_
                 << exit(FatalError);
         }
+
+        setMRFFaces();
     }
-
-    bool cellZoneFound = (cellZoneID_ != -1);
-
-    reduce(cellZoneFound, orOp<bool>());
-
-    if (!cellZoneFound)
-    {
-        FatalErrorIn
-        (
-            "Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
-        )   << "cannot find MRF cellZone " << name_
-            << exit(FatalError);
-    }
-
-    setMRFFaces();
 }
 
 
@@ -497,24 +502,37 @@ void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
 }
 
 
-Foam::Ostream& Foam::operator<<(Ostream& os, const MRFZone& MRF)
+void Foam::MRFZone::writeData(Ostream& os) const
 {
-    os  << indent << nl;
+    os  << nl;
-    os.write(MRF.name_) << nl;
+    os.write(name_) << nl;
     os  << token::BEGIN_BLOCK << incrIndent << nl;
-    os.writeKeyword("origin") << MRF.origin_ << token::END_STATEMENT << nl;
+    os.writeKeyword("active") << active_ << token::END_STATEMENT << nl;
-    os.writeKeyword("axis") << MRF.axis_ << token::END_STATEMENT << nl;
+    os.writeKeyword("cellZone") << cellZoneName_ << token::END_STATEMENT << nl;
-    MRF.omega_->writeData(os);
+    os.writeKeyword("origin") << origin_ << token::END_STATEMENT << nl;
+    os.writeKeyword("axis") << axis_ << token::END_STATEMENT << nl;
+    omega_->writeData(os);
 
-    if (MRF.excludedPatchNames_.size())
+    if (excludedPatchNames_.size())
     {
-        os.writeKeyword("nonRotatingPatches") << MRF.excludedPatchNames_
+        os.writeKeyword("nonRotatingPatches") << excludedPatchNames_
             << token::END_STATEMENT << nl;
     }
 
     os  << decrIndent << token::END_BLOCK << nl;
-
-    return os;
 }
 
+
+bool Foam::MRFZone::read(const dictionary& dict)
+{
+    coeffs_ = dict;
+
+    active_ = readBool(coeffs_.lookup("active"));
+    coeffs_.lookup("cellZone") >> cellZoneName_;
+    cellZoneID_ = mesh_.cellZones().findZoneID(cellZoneName_);
+
+    return true;
+}
+
+
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/MRF/MRFZone.H

@@ -67,12 +67,22 @@ class MRFZone
 {
     // Private data
 
+        //- Reference to the mesh database
         const fvMesh& mesh_;
 
+        //- Name of the MRF region
         const word name_;
 
-        const dictionary dict_;
+        //- Coefficients dictionary
+        dictionary coeffs_;
 
+        //- MRF region active flag
+        bool active_;
+
+        //- Name of cell zone
+        word cellZoneName_;
+
+        //- Cell zone ID
         label cellZoneID_;
 
         const wordList excludedPatchNames_;
@@ -133,8 +143,8 @@ public:
 
     // Constructors
 
-        //- Construct from fvMesh and Istream
+        //- Construct from fvMesh
-        MRFZone(const fvMesh& mesh, Istream& is);
+        MRFZone(const word& name, const fvMesh& mesh, const dictionary& dict);
 
         //- Return clone
         autoPtr<MRFZone> clone() const
@@ -143,79 +153,81 @@ public:
             return autoPtr<MRFZone>(NULL);
         }
 
-        //- Return a pointer to a new MRFZone created on freestore
-        //  from Istream
-        class iNew
-        {
-            const fvMesh& mesh_;
-
-        public:
-
-            iNew(const fvMesh& mesh)
-            :
-                mesh_(mesh)
-            {}
-
-            autoPtr<MRFZone> operator()(Istream& is) const
-            {
-                return autoPtr<MRFZone>(new MRFZone(mesh_, is));
-            }
-        };
-
 
     // Member Functions
 
-        //- Update the mesh corresponding to given map
+        // Access
-        void updateMesh(const mapPolyMesh& mpm)
-        {
-            // Only updates face addressing
-            setMRFFaces();
-        }
 
-        //- Add the Coriolis force contribution to the acceleration field
+            //- Return const access to the MRF region name
-        void addCoriolis(const volVectorField& U, volVectorField& ddtU) const;
+            inline const word& name() const;
 
-        //- Add the Coriolis force contribution to the momentum equation
+            //- Return const access to the MRF active flag
-        void addCoriolis(fvVectorMatrix& UEqn) const;
+            inline bool active() const;
-
-        //- Add the Coriolis force contribution to the momentum equation
-        void addCoriolis(const volScalarField& rho, fvVectorMatrix& UEqn) const;
-
-        //- Make the given absolute velocity relative within the MRF region
-        void relativeVelocity(volVectorField& U) const;
-
-        //- Make the given relative velocity absolute within the MRF region
-        void absoluteVelocity(volVectorField& U) const;
-
-        //- Make the given absolute flux relative within the MRF region
-        void relativeFlux(surfaceScalarField& phi) const;
-
-        //- Make the given absolute mass-flux relative within the MRF region
-        void relativeFlux
-        (
-            const surfaceScalarField& rho,
-            surfaceScalarField& phi
-        ) const;
-
-        //- Make the given relative flux absolute within the MRF region
-        void absoluteFlux(surfaceScalarField& phi) const;
-
-        //- Make the given relative mass-flux absolute within the MRF region
-        void absoluteFlux
-        (
-            const surfaceScalarField& rho,
-            surfaceScalarField& phi
-        ) const;
-
-        //- Correct the boundary velocity for the roation of the MRF region
-        void correctBoundaryVelocity(volVectorField& U) const;
 
 
-    // IOstream operator
+        // Evaluation
 
-        friend Ostream& operator<<(Ostream& os, const MRFZone& MRF);
+            //- Update the mesh corresponding to given map
+            void updateMesh(const mapPolyMesh& mpm)
+            {
+                // Only updates face addressing
+                setMRFFaces();
+            }
+
+            //- Add the Coriolis force contribution to the acceleration field
+            void addCoriolis
+            (
+                const volVectorField& U,
+                volVectorField& ddtU
+            ) const;
+
+            //- Add the Coriolis force contribution to the momentum equation
+            void addCoriolis(fvVectorMatrix& UEqn) const;
+
+            //- Add the Coriolis force contribution to the momentum equation
+            void addCoriolis
+            (
+                const volScalarField& rho,
+                fvVectorMatrix& UEqn
+            ) const;
+
+            //- Make the given absolute velocity relative within the MRF region
+            void relativeVelocity(volVectorField& U) const;
+
+            //- Make the given relative velocity absolute within the MRF region
+            void absoluteVelocity(volVectorField& U) const;
+
+            //- Make the given absolute flux relative within the MRF region
+            void relativeFlux(surfaceScalarField& phi) const;
+
+            //- Make the given absolute mass-flux relative within the MRF region
+            void relativeFlux
+            (
+                const surfaceScalarField& rho,
+                surfaceScalarField& phi
+            ) const;
+
+            //- Make the given relative flux absolute within the MRF region
+            void absoluteFlux(surfaceScalarField& phi) const;
+
+            //- Make the given relative mass-flux absolute within the MRF region
+            void absoluteFlux
+            (
+                const surfaceScalarField& rho,
+                surfaceScalarField& phi
+            ) const;
+
+            //- Correct the boundary velocity for the roation of the MRF region
+            void correctBoundaryVelocity(volVectorField& U) const;
 
 
+    // I-O
+
+        //- Write
+        void writeData(Ostream& os) const;
+
+        //- Read MRF dictionary
+        bool read(const dictionary& dict);
 };
 
 
@@ -231,6 +243,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "MRFZoneI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/MRF/MRFZoneI.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,53 +23,15 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-
+inline const Foam::word& Foam::MRFZone::name() const
-inline Foam::PtrList<Foam::volScalarField>& Foam::basicSolidMixture::Y()
 {
-    return Y_;
+    return name_;
 }
 
 
-inline const Foam::PtrList<Foam::volScalarField>& Foam::basicSolidMixture::Y()
+inline bool Foam::MRFZone::active() const
-const
 {
-    return Y_;
+    return active_;
-}
-
-
-inline Foam::volScalarField& Foam::basicSolidMixture::Y(const label i)
-{
-    return Y_[i];
-}
-
-
-inline const Foam::volScalarField& Foam::basicSolidMixture::Y
-(
-    const label i
-) const
-{
-    return Y_[i];
-}
-
-
-inline Foam::volScalarField& Foam::basicSolidMixture::Y(const word& specieName)
-{
-    return Y_[components_[specieName]];
-}
-
-
-inline const Foam::volScalarField& Foam::basicSolidMixture::Y
-(
-    const word& specieName
-) const
-{
-    return Y_[components_[specieName]];
-}
-
-
-inline bool Foam::basicSolidMixture::contains(const word& specieName) const
-{
-    return components_.contains(specieName);
 }
 
 

src/finiteVolume/cfdTools/general/MRF/MRFZoneList.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,59 +23,109 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "MRFZones.H"
+#include "MRFZoneList.H"
-#include "Time.H"
+#include "volFields.H"
-#include "fvMesh.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTemplateTypeNameAndDebug(IOPtrList<MRFZone>, 0);
-}
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::MRFZones::MRFZones(const fvMesh& mesh)
+Foam::MRFZoneList::MRFZoneList
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
 :
-    IOPtrList<MRFZone>
+    PtrList<MRFZone>(),
-    (
-        IOobject
-        (
-            "MRFZones",
-            mesh.time().constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        ),
-        MRFZone::iNew(mesh)
-    ),
     mesh_(mesh)
 {
-    if
+    reset(dict);
-    (
+
-        Pstream::parRun()
+    active();
-     &&
-     (
-        regIOobject::fileModificationChecking == timeStampMaster
-     || regIOobject::fileModificationChecking == inotifyMaster
-     )
-    )
-    {
-        WarningIn("MRFZones(const fvMesh&)")
-            << "The MRFZones are not run time modifiable\n"
-            << "    using 'timeStampMaster' or 'inotifyMaster'\n"
-            << "    for the entry fileModificationChecking\n"
-            << "    in the etc/controlDict.\n"
-            << "    Use 'timeStamp' instead."
-            << endl;
-    }
 }
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::MRFZoneList::~MRFZoneList()
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::MRFZones::addCoriolis
+bool Foam::MRFZoneList::active() const
+{
+    bool a = false;
+    forAll(*this, i)
+    {
+        a = a || this->operator[](i).active();
+    }
+
+    if (!a)
+    {
+        Info<< "    No MRF zones active" << endl;
+    }
+
+    return a;
+}
+
+
+void Foam::MRFZoneList::reset(const dictionary& dict)
+{
+    label count = 0;
+    forAllConstIter(dictionary, dict, iter)
+    {
+        if (iter().isDict())
+        {
+            count++;
+        }
+    }
+
+    this->setSize(count);
+    label i = 0;
+    forAllConstIter(dictionary, dict, iter)
+    {
+        if (iter().isDict())
+        {
+            const word& name = iter().keyword();
+            const dictionary& modelDict = iter().dict();
+
+            Info<< "    creating MRF zone: " << name << endl;
+
+            this->set
+            (
+                i++,
+                new MRFZone(name, mesh_, modelDict)
+            );
+        }
+    }
+}
+
+
+bool Foam::MRFZoneList::read(const dictionary& dict)
+{
+    bool allOk = true;
+    forAll(*this, i)
+    {
+        MRFZone& pm = this->operator[](i);
+        bool ok = pm.read(dict.subDict(pm.name()));
+        allOk = (allOk && ok);
+    }
+    return allOk;
+}
+
+
+bool Foam::MRFZoneList::writeData(Ostream& os) const
+{
+    forAll(*this, i)
+    {
+        os  << nl;
+        this->operator[](i).writeData(os);
+    }
+
+    return os.good();
+}
+
+
+void Foam::MRFZoneList::addCoriolis
 (
     const volVectorField& U,
     volVectorField& ddtU
@@ -88,7 +138,7 @@ void Foam::MRFZones::addCoriolis
 }
 
 
-void Foam::MRFZones::addCoriolis(fvVectorMatrix& UEqn) const
+void Foam::MRFZoneList::addCoriolis(fvVectorMatrix& UEqn) const
 {
     forAll(*this, i)
     {
@@ -97,7 +147,7 @@ void Foam::MRFZones::addCoriolis(fvVectorMatrix& UEqn) const
 }
 
 
-void Foam::MRFZones::addCoriolis
+void Foam::MRFZoneList::addCoriolis
 (
     const volScalarField& rho,
     fvVectorMatrix& UEqn
@@ -110,7 +160,7 @@ void Foam::MRFZones::addCoriolis
 }
 
 
-void Foam::MRFZones::relativeVelocity(volVectorField& U) const
+void Foam::MRFZoneList::relativeVelocity(volVectorField& U) const
 {
     forAll(*this, i)
     {
@@ -119,7 +169,7 @@ void Foam::MRFZones::relativeVelocity(volVectorField& U) const
 }
 
 
-void Foam::MRFZones::absoluteVelocity(volVectorField& U) const
+void Foam::MRFZoneList::absoluteVelocity(volVectorField& U) const
 {
     forAll(*this, i)
     {
@@ -128,7 +178,7 @@ void Foam::MRFZones::absoluteVelocity(volVectorField& U) const
 }
 
 
-void Foam::MRFZones::relativeFlux(surfaceScalarField& phi) const
+void Foam::MRFZoneList::relativeFlux(surfaceScalarField& phi) const
 {
     forAll(*this, i)
     {
@@ -137,7 +187,7 @@ void Foam::MRFZones::relativeFlux(surfaceScalarField& phi) const
 }
 
 
-void Foam::MRFZones::relativeFlux
+void Foam::MRFZoneList::relativeFlux
 (
     const surfaceScalarField& rho,
     surfaceScalarField& phi
@@ -150,7 +200,7 @@ void Foam::MRFZones::relativeFlux
 }
 
 
-void Foam::MRFZones::absoluteFlux(surfaceScalarField& phi) const
+void Foam::MRFZoneList::absoluteFlux(surfaceScalarField& phi) const
 {
     forAll(*this, i)
     {
@@ -159,7 +209,7 @@ void Foam::MRFZones::absoluteFlux(surfaceScalarField& phi) const
 }
 
 
-void Foam::MRFZones::absoluteFlux
+void Foam::MRFZoneList::absoluteFlux
 (
     const surfaceScalarField& rho,
     surfaceScalarField& phi
@@ -172,7 +222,7 @@ void Foam::MRFZones::absoluteFlux
 }
 
 
-void Foam::MRFZones::correctBoundaryVelocity(volVectorField& U) const
+void Foam::MRFZoneList::correctBoundaryVelocity(volVectorField& U) const
 {
     forAll(*this, i)
     {
@@ -181,10 +231,16 @@ void Foam::MRFZones::correctBoundaryVelocity(volVectorField& U) const
 }
 
 
-bool Foam::MRFZones::readData(Istream& is)
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const MRFZoneList& models
+)
 {
-    PtrList<MRFZone>::read(is, MRFZone::iNew(mesh_));
+    models.writeData(os);
-    return is.good();
+    return os;
 }
 
 

src/finiteVolume/cfdTools/general/MRF/MRFZoneList.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,61 +22,77 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::MRFZones
+    Foam::MRFZoneList
 
 Description
-    Container class for a set of MRFZones with the MRFZone member functions
+    List container for MRF zomes
-    implemented to loop over the functions for each MRFZone.
 
 SourceFiles
-    MRFZones.C
+    MRFZoneList.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef MRFZones_H
+#ifndef MRFZoneList_H
-#define MRFZones_H
+#define MRFZoneList_H
 
+#include "fvMesh.H"
+#include "dictionary.H"
+#include "fvMatricesFwd.H"
 #include "MRFZone.H"
-#include "IOPtrList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class MRFZoneList;
+Ostream& operator<<(Ostream& os, const MRFZoneList& models);
+
 /*---------------------------------------------------------------------------*\
-                           Class MRFZones Declaration
+                         Class MRFZoneList Declaration
 \*---------------------------------------------------------------------------*/
 
-class MRFZones
+class MRFZoneList
 :
-    public IOPtrList<MRFZone>
+    PtrList<MRFZone>
 {
-    // Private data
+private:
-
-        //- Reference to mesh
-        const fvMesh& mesh_;
-
 
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
-        MRFZones(const MRFZones&);
+        MRFZoneList(const MRFZoneList&);
 
         //- Disallow default bitwise assignment
-        void operator=(const MRFZones&);
+        void operator=(const MRFZoneList&);
+
+
+protected:
+
+    // Protected data
+
+        //- Reference to the mesh database
+        const fvMesh& mesh_;
 
 
 public:
 
-    // Constructors
+    //- Constructor
+    MRFZoneList(const fvMesh& mesh, const dictionary& dict);
 
-        //- Construct from fvMesh
+    //- Destructor
-        MRFZones(const fvMesh& mesh);
+    ~MRFZoneList();
 
 
     // Member Functions
 
+        //- Return active status
+        bool active() const;
+
+        //- Reset the source list
+        void reset(const dictionary& dict);
+
         //- Add the Coriolis force contribution to the acceleration field
         void addCoriolis(const volVectorField& U, volVectorField& ddtU) const;
 
@@ -116,14 +132,22 @@ public:
         void correctBoundaryVelocity(volVectorField& U) const;
 
 
-    // I-O
+        // I-O
 
-        //- Read from Istream
+            //- Read dictionary
-        virtual bool readData(Istream&);
+            bool read(const dictionary& dict);
 
+            //- Write data to Ostream
+            bool writeData(Ostream& os) const;
+
+            //- Ostream operator
+            friend Ostream& operator<<
+            (
+                Ostream& os,
+                const MRFZoneList& models
+            );
 };
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C

@@ -27,7 +27,6 @@ License
 #include "fvMesh.H"
 #include "volFields.H"
 #include "surfaceFields.H"
-#include "geometricOneField.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.C

@@ -26,12 +26,6 @@ License
 #include "porosityModelList.H"
 #include "volFields.H"
 
-// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
-/*
-void Foam::porosityModelList::XXX()
-{}
-*/
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::porosityModelList::porosityModelList

src/finiteVolume/fields/fvPatchFields/derived/movingWallVelocity/movingWallVelocityFvPatchVectorField.C

@@ -105,33 +105,39 @@ void Foam::movingWallVelocityFvPatchVectorField::updateCoeffs()
         return;
     }
 
-    const fvPatch& p = patch();
-    const polyPatch& pp = p.patch();
     const fvMesh& mesh = dimensionedInternalField().mesh();
-    const pointField& oldPoints = mesh.oldPoints();
 
-    vectorField oldFc(pp.size());
+    if (mesh.changing())
-
-    forAll(oldFc, i)
     {
-        oldFc[i] = pp[i].centre(oldPoints);
+        const fvPatch& p = patch();
+        const polyPatch& pp = p.patch();
+        const pointField& oldPoints = mesh.oldPoints();
+
+        vectorField oldFc(pp.size());
+
+        forAll(oldFc, i)
+        {
+            oldFc[i] = pp[i].centre(oldPoints);
+        }
+
+        const scalar deltaT = mesh.time().deltaTValue();
+
+        const vectorField Up((pp.faceCentres() - oldFc)/deltaT);
+
+        const volVectorField& U = db().lookupObject<volVectorField>(UName_);
+        scalarField phip
+        (
+            p.patchField<surfaceScalarField, scalar>(fvc::meshPhi(U))
+        );
+
+        const vectorField n(p.nf());
+        const scalarField& magSf = p.magSf();
+        tmp<scalarField> Un = phip/(magSf + VSMALL);
+
+
+        vectorField::operator=(Up + n*(Un - (n & Up)));
     }
 
-    const vectorField Up((pp.faceCentres() - oldFc)/mesh.time().deltaTValue());
-
-    const volVectorField& U = db().lookupObject<volVectorField>(UName_);
-    scalarField phip
-    (
-        p.patchField<surfaceScalarField, scalar>(fvc::meshPhi(U))
-    );
-
-    const vectorField n(p.nf());
-    const scalarField& magSf = p.magSf();
-    tmp<scalarField> Un = phip/(magSf + VSMALL);
-
-
-    vectorField::operator=(Up + n*(Un - (n & Up)));
-
     fixedValueFvPatchVectorField::updateCoeffs();
 }
 

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

@@ -52,7 +52,8 @@ timeVaryingMappedFixedValueFvPatchField
     startAverage_(pTraits<Type>::zero),
     endSampleTime_(-1),
     endSampledValues_(0),
-    endAverage_(pTraits<Type>::zero)
+    endAverage_(pTraits<Type>::zero),
+    offSet_()
 {}
 
 
@@ -77,7 +78,8 @@ timeVaryingMappedFixedValueFvPatchField
     startAverage_(pTraits<Type>::zero),
     endSampleTime_(-1),
     endSampledValues_(0),
-    endAverage_(pTraits<Type>::zero)
+    endAverage_(pTraits<Type>::zero),
+    offSet_()
 {}
 
 
@@ -101,7 +103,8 @@ timeVaryingMappedFixedValueFvPatchField
     startAverage_(pTraits<Type>::zero),
     endSampleTime_(-1),
     endSampledValues_(0),
-    endAverage_(pTraits<Type>::zero)
+    endAverage_(pTraits<Type>::zero),
+    offSet_(DataEntry<Type>::New("offSet", dict))
 {
     dict.readIfPresent("fieldTableName", fieldTableName_);
 
@@ -134,7 +137,8 @@ timeVaryingMappedFixedValueFvPatchField
     startAverage_(ptf.startAverage_),
     endSampleTime_(ptf.endSampleTime_),
     endSampledValues_(ptf.endSampledValues_),
-    endAverage_(ptf.endAverage_)
+    endAverage_(ptf.endAverage_),
+    offSet_(ptf.offSet_().clone().ptr())
 {}
 
 
@@ -158,7 +162,8 @@ timeVaryingMappedFixedValueFvPatchField
     startAverage_(ptf.startAverage_),
     endSampleTime_(ptf.endSampleTime_),
     endSampledValues_(ptf.endSampledValues_),
-    endAverage_(ptf.endAverage_)
+    endAverage_(ptf.endAverage_),
+    offSet_(ptf.offSet_().clone().ptr())
 {}
 
 
@@ -276,7 +281,7 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::checkTable()
     {
         FatalErrorIn
         (
-            "timeVaryingMappedFixedValueFvPatchField<Type>::checkTable"
+            "timeVaryingMappedFixedValueFvPatchField<Type>::checkTable()"
         )   << "Cannot find starting sampling values for current time "
             << this->db().time().value() << nl
             << "Have sampling values for times "
@@ -366,6 +371,7 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::checkTable()
                         /sampleTimes_[endSampleTime_].name()
                     << endl;
             }
+
             // Reread values and interpolate
             AverageIOField<Type> vals
             (
@@ -392,7 +398,6 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::checkTable()
 template<class Type>
 void timeVaryingMappedFixedValueFvPatchField<Type>::updateCoeffs()
 {
-
     if (this->updated())
     {
         return;
@@ -423,7 +428,7 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::updateCoeffs()
         scalar start = sampleTimes_[startSampleTime_].value();
         scalar end = sampleTimes_[endSampleTime_].value();
 
-        scalar s = (this->db().time().value()-start)/(end-start);
+        scalar s = (this->db().time().value() - start)/(end - start);
 
         if (debug)
         {
@@ -434,8 +439,8 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::updateCoeffs()
                 << " with weight:" << s << endl;
         }
 
-        this->operator==((1-s)*startSampledValues_ + s*endSampledValues_);
+        this->operator==((1 - s)*startSampledValues_ + s*endSampledValues_);
-        wantedAverage = (1-s)*startAverage_ + s*endAverage_;
+        wantedAverage = (1 - s)*startAverage_ + s*endAverage_;
     }
 
     // Enforce average. Either by scaling (if scaling factor > 0.5) or by
@@ -465,7 +470,7 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::updateCoeffs()
                 Pout<< "updateCoeffs :"
                     << " offsetting with:" << offset << endl;
             }
-            this->operator==(fld+offset);
+            this->operator==(fld + offset);
         }
         else
         {
@@ -480,6 +485,10 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::updateCoeffs()
         }
     }
 
+    // apply offset to mapped values
+    const scalar t = this->db().time().timeOutputValue();
+    this->operator==(*this + offSet_->value(t));
+
     if (debug)
     {
         Pout<< "updateCoeffs : set fixedValue to min:" << gMin(*this)
@@ -503,6 +512,8 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::write(Ostream& os) const
             << token::END_STATEMENT << nl;
     }
 
+    offSet_->writeData(os);
+
     this->writeEntry("value", os);
 }
 

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

@@ -80,6 +80,7 @@ SourceFiles
 #include "FixedList.H"
 #include "instantList.H"
 #include "pointToPointPlanarInterpolation.H"
+#include "DataEntry.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -130,6 +131,9 @@ class timeVaryingMappedFixedValueFvPatchField
         //- If setAverage: end average value
         Type endAverage_;
 
+        //- Time varying offset values to interpolated data
+        autoPtr<DataEntry<Type> > offSet_;
+
 
 public:
 

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

@@ -318,7 +318,6 @@ void Foam::autoLayerDriver::handleNonManifolds
         }
     }
 
-
     Info<< "Set displacement to zero for all " << nNonManif
         << " non-manifold points" << endl;
 }
@@ -443,86 +442,6 @@ void Foam::autoLayerDriver::handleFeatureAngle
 }
 
 
-//Foam::tmp<Foam::scalarField> Foam::autoLayerDriver::undistortedEdgeLength
-//(
-//    const indirectPrimitivePatch& pp,
-//    const bool relativeSizes,
-//    const bool faceSize
-//)
-//{
-//    const fvMesh& mesh = meshRefiner_.mesh();
-//
-//    tmp<scalarField> tfld(new scalarField());
-//    scalarField& fld = tfld();
-//
-//    if (faceSize)
-//    {
-//        fld.setSize(pp.size());
-//    }
-//    else
-//    {
-//        fld.setSize(pp.nPoints());
-//    }
-//
-//
-//    if (relativeSizes)
-//    {
-//        const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
-//        const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
-//
-//        if (faceSize)
-//        {
-//            forAll(pp, i)
-//            {
-//                label faceI = pp.addressing()[i];
-//                label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
-//                fld[i] = edge0Len/(1<<ownLevel);
-//            }
-//        }
-//        else
-//        {
-//            // Determine per point the max cell level of connected cells
-//            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-//
-//            labelList maxPointLevel(pp.nPoints(), labelMin);
-//
-//            forAll(pp, i)
-//            {
-//                label faceI = pp.addressing()[i];
-//                label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
-//
-//                const face& f = pp.localFaces()[i];
-//                forAll(f, fp)
-//                {
-//                    maxPointLevel[f[fp]] =
-//                        max(maxPointLevel[f[fp]], ownLevel);
-//                }
-//            }
-//
-//            syncTools::syncPointList
-//            (
-//                mesh,
-//                pp.meshPoints(),
-//                maxPointLevel,
-//                maxEqOp<label>(),
-//                labelMin            // null value
-//            );
-//
-//
-//            forAll(maxPointLevel, pointI)
-//            {
-//                // Find undistorted edge size for this level.
-//                fld[i] = edge0Len/(1<<maxPointLevel[pointI]);
-//            }
-//        }
-//    }
-//    else
-//    {
-//        // Use actual cell size
-//    }
-//}
-
-
 // No extrusion on cells with warped faces. Calculates the thickness of the
 // layer and compares it to the space the warped face takes up. Disables
 // extrusion if layer thickness is more than faceRatio of the thickness of
@@ -702,7 +621,6 @@ void Foam::autoLayerDriver::handleWarpedFaces
 //}
 
 
-// No extrusion on faces with differing number of layers for points
 void Foam::autoLayerDriver::setNumLayers
 (
     const labelList& patchToNLayers,
@@ -865,15 +783,6 @@ Foam::autoLayerDriver::makeLayerDisplacementField
     }
 
 
-//Pout<< "*** makeLayerDisplacementField : boundary conditions:" << endl;
-//forAll(patchFieldTypes, patchI)
-//{
-//    Pout<< "\t" << patchI << " name:" << pointPatches[patchI].name()
-//        << " type:" << patchFieldTypes[patchI]
-//        << " nLayers:" << numLayers[patchI]
-//        << endl;
-//}
-
     const polyMesh& mesh = pMesh();
 
     // Note: time().timeName() instead of meshRefinement::timeName() since
@@ -1053,11 +962,10 @@ void Foam::autoLayerDriver::determineSidePatches
             patchDict.add("nFaces", 0);
             patchDict.add("startFace", mesh.nFaces());
 
-            Pout<< "Adding patch " << patchI
+            //Pout<< "Adding patch " << patchI
-                << " name:" << name
+            //    << " name:" << name
-                << " between " << Pstream::myProcNo()
+            //    << " between " << Pstream::myProcNo()
-                << " and " << nbrProcI
+            //    << " and " << nbrProcI << endl;
-                << endl;
 
             label procPatchI = meshRefiner_.appendPatch
             (
@@ -1090,12 +998,7 @@ void Foam::autoLayerDriver::calculateLayerThickness
 (
     const indirectPrimitivePatch& pp,
     const labelList& patchIDs,
-    const scalarField& patchExpansionRatio,
+    const layerParameters& layerParams,
-
-    const bool relativeSizes,
-    const scalarField& patchFinalLayerThickness,
-    const scalarField& patchMinThickness,
-
     const labelList& cellLevel,
     const labelList& patchNLayers,
     const scalar edge0Len,
@@ -1111,12 +1014,13 @@ void Foam::autoLayerDriver::calculateLayerThickness
 
     // Rework patch-wise layer parameters into minimum per point
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Note: only layer parameters consistent with layer specification
+    // method (see layerParameters) will be correct.
+    scalarField firstLayerThickness(pp.nPoints(), GREAT);
+    scalarField finalLayerThickness(pp.nPoints(), GREAT);
+    scalarField totalThickness(pp.nPoints(), GREAT);
+    scalarField expRatio(pp.nPoints(), GREAT);
 
-    // Reuse input fields
-    expansionRatio.setSize(pp.nPoints());
-    expansionRatio = GREAT;
-    thickness.setSize(pp.nPoints());
-    thickness = GREAT;
     minThickness.setSize(pp.nPoints());
     minThickness = GREAT;
 
@@ -1130,20 +1034,30 @@ void Foam::autoLayerDriver::calculateLayerThickness
         {
             label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
 
-            expansionRatio[ppPointI] = min
+            firstLayerThickness[ppPointI] = min
             (
-                expansionRatio[ppPointI],
+                firstLayerThickness[ppPointI],
-                patchExpansionRatio[patchI]
+                layerParams.firstLayerThickness()[patchI]
             );
-            thickness[ppPointI] = min
+            finalLayerThickness[ppPointI] = min
             (
-                thickness[ppPointI],
+                finalLayerThickness[ppPointI],
-                patchFinalLayerThickness[patchI]
+                layerParams.finalLayerThickness()[patchI]
+            );
+            totalThickness[ppPointI] = min
+            (
+                totalThickness[ppPointI],
+                layerParams.thickness()[patchI]
+            );
+            expRatio[ppPointI] = min
+            (
+                expRatio[ppPointI],
+                layerParams.expansionRatio()[patchI]
             );
             minThickness[ppPointI] = min
             (
                 minThickness[ppPointI],
-                patchMinThickness[patchI]
+                layerParams.minThickness()[patchI]
             );
         }
     }
@@ -1152,7 +1066,7 @@ void Foam::autoLayerDriver::calculateLayerThickness
     (
         mesh,
         pp.meshPoints(),
-        expansionRatio,
+        firstLayerThickness,
         minEqOp<scalar>(),
         GREAT               // null value
     );
@@ -1160,7 +1074,23 @@ void Foam::autoLayerDriver::calculateLayerThickness
     (
         mesh,
         pp.meshPoints(),
-        thickness,
+        finalLayerThickness,
+        minEqOp<scalar>(),
+        GREAT               // null value
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        totalThickness,
+        minEqOp<scalar>(),
+        GREAT               // null value
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        expRatio,
         minEqOp<scalar>(),
         GREAT               // null value
     );
@@ -1182,14 +1112,18 @@ void Foam::autoLayerDriver::calculateLayerThickness
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // by multiplying with the internal cell size.
 
-    if (relativeSizes)
+    if (layerParams.relativeSizes())
     {
-        if (min(patchMinThickness) < 0 || max(patchMinThickness) > 2)
+        if
+        (
+            min(layerParams.minThickness()) < 0
+         || max(layerParams.minThickness()) > 2
+        )
         {
             FatalErrorIn("calculateLayerThickness(..)")
                 << "Thickness should be factor of local undistorted cell size."
                 << " Valid values are [0..2]." << nl
-                << " minThickness:" << patchMinThickness
+                << " minThickness:" << layerParams.minThickness()
                 << exit(FatalError);
         }
 
@@ -1225,38 +1159,114 @@ void Foam::autoLayerDriver::calculateLayerThickness
         {
             // Find undistorted edge size for this level.
             scalar edgeLen = edge0Len/(1<<maxPointLevel[pointI]);
-            thickness[pointI] *= edgeLen;
+            firstLayerThickness[pointI] *= edgeLen;
+            finalLayerThickness[pointI] *= edgeLen;
+            totalThickness[pointI] *= edgeLen;
             minThickness[pointI] *= edgeLen;
         }
     }
 
 
 
-    // Rework thickness (of final layer) into overall thickness of all layers
+    // Rework thickness parameters into overall thickness
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    forAll(thickness, pointI)
+    forAll(firstLayerThickness, pointI)
     {
-        // Calculate layer thickness based on expansion ratio
+        thickness[pointI] = layerParams.layerThickness
-        // and final layer height
+        (
-        if (expansionRatio[pointI] == 1)
+            patchNLayers[pointI],
-        {
+            firstLayerThickness[pointI],
-            thickness[pointI] *= patchNLayers[pointI];
+            finalLayerThickness[pointI],
-        }
+            totalThickness[pointI],
-        else
+            expRatio[pointI]
-        {
+        );
 
-            scalar invExpansion = 1.0 / expansionRatio[pointI];
+        expansionRatio[pointI] = layerParams.layerExpansionRatio
-            label nLay = patchNLayers[pointI];
+        (
-            thickness[pointI] *=
+            patchNLayers[pointI],
-                (1.0 - pow(invExpansion, nLay))
+            firstLayerThickness[pointI],
-              / (1.0 - invExpansion);
+            finalLayerThickness[pointI],
-        }
+            totalThickness[pointI],
+            expRatio[pointI]
+        );
     }
 
-
     //Info<< "calculateLayerThickness : min:" << gMin(thickness)
     //    << " max:" << gMax(thickness) << endl;
+
+    // Print a bit
+    {
+        const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+        // Find maximum length of a patch name, for a nicer output
+        label maxPatchNameLen = 0;
+        forAll(patchIDs, i)
+        {
+            label patchI = patchIDs[i];
+            word patchName = patches[patchI].name();
+            maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
+        }
+
+        Info<< nl
+            << setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
+            << setw(0) << " faces    layers avg thickness[m]" << nl
+            << setf(ios_base::left) << setw(maxPatchNameLen) << " "
+            << setw(0) << "                 near-wall overall" << nl
+            << setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
+            << setw(0) << " -----    ------ --------- -------" << endl;
+
+        forAll(patchIDs, i)
+        {
+            label patchI = patchIDs[i];
+
+            const labelList& meshPoints = patches[patchI].meshPoints();
+
+            scalar sumThickness = 0;
+            scalar sumNearWallThickness = 0;
+
+            forAll(meshPoints, patchPointI)
+            {
+                label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
+
+                sumThickness += thickness[ppPointI];
+                sumNearWallThickness += layerParams.firstLayerThickness
+                (
+                    patchNLayers[ppPointI],
+                    firstLayerThickness[ppPointI],
+                    finalLayerThickness[ppPointI],
+                    thickness[ppPointI],
+                    expansionRatio[ppPointI]
+                );
+            }
+
+            label totNPoints = returnReduce(meshPoints.size(), sumOp<label>());
+
+            // For empty patches, totNPoints is 0.
+            scalar avgThickness = 0;
+            scalar avgNearWallThickness = 0;
+
+            if (totNPoints > 0)
+            {
+                avgThickness =
+                    returnReduce(sumThickness, sumOp<scalar>())
+                  / totNPoints;
+                avgNearWallThickness =
+                    returnReduce(sumNearWallThickness, sumOp<scalar>())
+                  / totNPoints;
+            }
+
+            Info<< setf(ios_base::left) << setw(maxPatchNameLen)
+                << patches[patchI].name() << setprecision(3)
+                << " " << setw(8)
+                << returnReduce(patches[patchI].size(), sumOp<scalar>())
+                << " " << setw(6) << layerParams.numLayers()[patchI]
+                << " " << setw(8) << avgNearWallThickness
+                << "  " << setw(8) << avgThickness
+                << endl;
+        }
+        Info<< endl;
+    }
 }
 
 
@@ -2618,7 +2628,8 @@ void Foam::autoLayerDriver::addLayers
     const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
     const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
 
-    // Determine (wanted) point-wise layer thickness and expansion ratio
+    // Determine (wanted) point-wise overall layer thickness and expansion
+    // ratio
     scalarField thickness(pp().nPoints());
     scalarField minThickness(pp().nPoints());
     scalarField expansionRatio(pp().nPoints());
@@ -2626,12 +2637,7 @@ void Foam::autoLayerDriver::addLayers
     (
         pp,
         meshMover().adaptPatchIDs(),
-        layerParams.expansionRatio(),
+        layerParams,
-
-        layerParams.relativeSizes(),        // thickness relative to cellsize?
-        layerParams.finalLayerThickness(),  // wanted thicknes
-        layerParams.minThickness(),         // minimum thickness
-
         cellLevel,
         patchNLayers,
         edge0Len,
@@ -2642,87 +2648,6 @@ void Foam::autoLayerDriver::addLayers
     );
 
 
-    // Print a bit
-    {
-        const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-        // Find maximum length of a patch name, for a nicer output
-        label maxPatchNameLen = 0;
-        forAll(meshMover().adaptPatchIDs(), i)
-        {
-            label patchI = meshMover().adaptPatchIDs()[i];
-            word patchName = patches[patchI].name();
-            maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
-        }
-
-        Info<< nl
-            << setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
-            << setw(0) << " faces    layers avg thickness[m]" << nl
-            << setf(ios_base::left) << setw(maxPatchNameLen) << " "
-            << setw(0) << "                 near-wall overall" << nl
-            << setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
-            << setw(0) << " -----    ------ --------- -------" << endl;
-
-        forAll(meshMover().adaptPatchIDs(), i)
-        {
-            label patchI = meshMover().adaptPatchIDs()[i];
-
-            const labelList& meshPoints = patches[patchI].meshPoints();
-
-            scalar sumThickness = 0;
-            scalar sumNearWallThickness = 0;
-
-            forAll(meshPoints, patchPointI)
-            {
-                label ppPointI = pp().meshPointMap()[meshPoints[patchPointI]];
-
-                sumThickness += thickness[ppPointI];
-
-                label nLay = patchNLayers[ppPointI];
-                if (nLay > 0)
-                {
-                    if (expansionRatio[ppPointI] == 1)
-                    {
-                        sumNearWallThickness += thickness[ppPointI]/nLay;
-                    }
-                    else
-                    {
-                        scalar s =
-                            (1.0-pow(expansionRatio[ppPointI], nLay))
-                          / (1.0-expansionRatio[ppPointI]);
-                        sumNearWallThickness += thickness[ppPointI]/s;
-                    }
-                }
-            }
-
-            label totNPoints = returnReduce(meshPoints.size(), sumOp<label>());
-
-            // For empty patches, totNPoints is 0.
-            scalar avgThickness = 0;
-            scalar avgNearWallThickness = 0;
-
-            if (totNPoints > 0)
-            {
-                avgThickness =
-                    returnReduce(sumThickness, sumOp<scalar>())
-                  / totNPoints;
-                avgNearWallThickness =
-                    returnReduce(sumNearWallThickness, sumOp<scalar>())
-                  / totNPoints;
-            }
-
-            Info<< setf(ios_base::left) << setw(maxPatchNameLen)
-                << patches[patchI].name() << setprecision(3)
-                << " " << setw(8)
-                << returnReduce(patches[patchI].size(), sumOp<scalar>())
-                << " " << setw(6) << layerParams.numLayers()[patchI]
-                << " " << setw(8) << avgNearWallThickness
-                << "  " << setw(8) << avgThickness
-                << endl;
-        }
-        Info<< endl;
-    }
-
 
     // Calculate wall to medial axis distance for smoothing displacement
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -2958,8 +2883,8 @@ void Foam::autoLayerDriver::addLayers
         // Determine per point/per face number of layers to extrude. Also
         // handles the slow termination of layers when going switching layers
 
-        labelList nPatchPointLayers(pp().nPoints(),-1);
+        labelList nPatchPointLayers(pp().nPoints(), -1);
-        labelList nPatchFaceLayers(pp().localFaces().size(),-1);
+        labelList nPatchFaceLayers(pp().size(), -1);
         setupLayerInfoTruncation
         (
             meshMover(),
@@ -2970,31 +2895,22 @@ void Foam::autoLayerDriver::addLayers
             nPatchFaceLayers
         );
 
-        // Calculate displacement for first layer for addPatchLayer.
+        // Calculate displacement for final layer for addPatchLayer.
-        // (first layer = layer of cells next to the original mesh)
+        // (layer of cells next to the original mesh)
-        vectorField firstDisp(patchNLayers.size(), vector::zero);
+        vectorField finalDisp(patchNLayers.size(), vector::zero);
 
         forAll(nPatchPointLayers, i)
         {
-            if (nPatchPointLayers[i] > 0)
+            scalar ratio = layerParams.finalLayerThicknessRatio
-            {
+            (
-                if (expansionRatio[i] == 1.0)
+                nPatchPointLayers[i],
-                {
+                expansionRatio[i]
-                    firstDisp[i] = patchDisp[i]/nPatchPointLayers[i];
+            );
-                }
+            finalDisp[i] = ratio*patchDisp[i];
-                else
-                {
-                    label nLay = nPatchPointLayers[i];
-                    scalar h =
-                        pow(expansionRatio[i], nLay - 1)
-                      * (1.0 - expansionRatio[i])
-                      / (1.0 - pow(expansionRatio[i], nLay));
-                    firstDisp[i] = h*patchDisp[i];
-                }
-            }
         }
 
-        const scalarField invExpansionRatio(1.0 / expansionRatio);
+
+        const scalarField invExpansionRatio(1.0/expansionRatio);
 
         // Add topo regardless of whether extrudeStatus is extruderemove.
         // Not add layer if patchDisp is zero.
@@ -3009,7 +2925,7 @@ void Foam::autoLayerDriver::addLayers
             labelList(0),       // exposed patchIDs, not used for adding layers
             nPatchFaceLayers,   // layers per face
             nPatchPointLayers,  // layers per point
-            firstDisp,          // thickness of layer nearest internal mesh
+            finalDisp,          // thickness of layer nearest internal mesh
             meshMod
         );
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H

@@ -237,12 +237,7 @@ class autoLayerDriver
                 (
                     const indirectPrimitivePatch& pp,
                     const labelList& patchIDs,
-
+                    const layerParameters& layerParams,
-                    const scalarField& patchExpansionRatio,
-                    const bool relativeSizes,
-                    const scalarField& patchFinalLayerThickness,
-                    const scalarField& patchMinThickness,
-
                     const labelList& cellLevel,
                     const labelList& patchNLayers,
                     const scalar edge0Len,

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

@@ -34,6 +34,72 @@ License
 const Foam::scalar Foam::layerParameters::defaultConcaveAngle = 90;
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::layerParameters::layerExpansionRatio
+(
+    const label n,
+    const scalar totalOverFirst
+) const
+{
+    if (n <= 1)
+    {
+        return 1.0;
+    }
+
+    //scalar totalOverFirst = totalThickness/firstLayerThickess;
+
+    const label maxIters = 10;
+    const scalar tol = 1e-8;
+
+
+    if (mag(n-totalOverFirst) < tol)
+    {
+        return 1.0;
+    }
+
+    // Calculate the bounds of the solution
+    scalar minR;
+    scalar maxR;
+
+    if (totalOverFirst < n)
+    {
+        minR = 0.0;
+        maxR = pow(totalOverFirst/n, 1/(n-1));
+    }
+    else
+    {
+        minR = pow(totalOverFirst/n, 1/(n-1));
+        maxR = totalOverFirst/(n - 1);
+    }
+
+    //Info<< "Solution bounds = (" << minR << ", " << maxR << ")" << nl << endl;
+
+    // Starting guess
+    scalar r = 0.5*(minR + maxR);
+
+    for (label i = 0; i < maxIters; ++i)
+    {
+        const scalar prevr = r;
+
+        const scalar fx = pow(r, n) - totalOverFirst*r - (1 - totalOverFirst);
+        const scalar dfx = n*pow(r, n - 1) - totalOverFirst;
+
+        r -= fx/dfx;
+
+        const scalar error = mag(r - prevr);
+
+        //Info<< i << " " << r << " Error = " << error << endl;
+
+        if (error < tol)
+        {
+            break;
+        }
+    }
+    return r;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from dictionary
@@ -44,17 +110,12 @@ Foam::layerParameters::layerParameters
 )
 :
     numLayers_(boundaryMesh.size(), -1),
-    expansionRatio_
-    (
-        boundaryMesh.size(),
-        readScalar(dict.lookup("expansionRatio"))
-    ),
     relativeSizes_(dict.lookup("relativeSizes")),
-    finalLayerThickness_
+    layerSpec_(ILLEGAL),
-    (
+    firstLayerThickness_(boundaryMesh.size(), -123),
-        boundaryMesh.size(),
+    finalLayerThickness_(boundaryMesh.size(), -123),
-        readScalar(dict.lookup("finalLayerThickness"))
+    thickness_(boundaryMesh.size(), -123),
-    ),
+    expansionRatio_(boundaryMesh.size(), -123),
     minThickness_
     (
         boundaryMesh.size(),
@@ -103,24 +164,108 @@ Foam::layerParameters::layerParameters
     nRelaxedIter_(labelMax),
     additionalReporting_(dict.lookupOrDefault("additionalReporting", false))
 {
-    if (nGrow_ > 0)
+    // Detect layer specification mode
+
+    label nSpec = 0;
+
+    bool haveFirst = dict.found("firstLayerThickness");
+    if (haveFirst)
     {
-        WarningIn("layerParameters::layerParameters(..)")
+        firstLayerThickness_ = scalarField
-            << "The nGrow parameter effect has changed with respect to 1.6.x."
+        (
-            << endl
+            boundaryMesh.size(),
-            << "Please set nGrow=0 for 1.6.x behaviour."
+            readScalar(dict.lookup("firstLayerThickness"))
+        );
+        nSpec++;
+    }
+    bool haveFinal = dict.found("finalLayerThickness");
+    if (haveFinal)
+    {
+        finalLayerThickness_ = scalarField
+        (
+            boundaryMesh.size(),
+            readScalar(dict.lookup("finalLayerThickness"))
+        );
+        nSpec++;
+    }
+    bool haveTotal = dict.found("thickness");
+    if (haveTotal)
+    {
+        thickness_ = scalarField
+        (
+            boundaryMesh.size(),
+            readScalar(dict.lookup("thickness"))
+        );
+        nSpec++;
+    }
+    bool haveExp = dict.found("expansionRatio");
+    if (haveExp)
+    {
+        expansionRatio_ = scalarField
+        (
+            boundaryMesh.size(),
+            readScalar(dict.lookup("expansionRatio"))
+        );
+        nSpec++;
+    }
+
+
+    if (haveFirst && haveTotal)
+    {
+        layerSpec_ = FIRST_AND_TOTAL;
+        Info<< "Layer thickness specified as first layer and overall thickness."
             << endl;
     }
+    else if (haveFirst && haveExp)
+    {
+        layerSpec_ = FIRST_AND_EXPANSION;
+        Info<< "Layer thickness specified as first layer and expansion ratio."
+            << endl;
+    }
+    else if (haveFinal && haveTotal)
+    {
+        layerSpec_ = FINAL_AND_TOTAL;
+        Info<< "Layer thickness specified as final layer and overall thickness."
+            << endl;
+    }
+    else if (haveFinal && haveExp)
+    {
+        layerSpec_ = FINAL_AND_EXPANSION;
+        Info<< "Layer thickness specified as final layer and expansion ratio."
+            << endl;
+    }
+
+
+    if (layerSpec_ == ILLEGAL || nSpec != 2)
+    {
+        FatalIOErrorIn
+        (
+            "layerParameters::layerParameters"
+            "(const dictionary&, const polyBoundaryMesh&)",
+            dict
+        )   << "Over- or underspecified layer thickness."
+            << " Please specify" << nl
+            << "    first layer thickness ('firstLayerThickness')"
+            << " and overall thickness ('thickness') or" << nl
+            << "    first layer thickness ('firstLayerThickness')"
+            << " and expansion ratio ('expansionRatio') or" << nl
+            << "    final layer thickness ('finalLayerThickness')"
+            << " and expansion ratio ('expansionRatio') or" << nl
+            << "    final layer thickness ('finalLayerThickness')"
+            << " and overall thickness ('thickness')"
+            << exit(FatalIOError);
+    }
+
 
     dict.readIfPresent("nRelaxedIter", nRelaxedIter_);
 
     if (nLayerIter_ < 0 || nRelaxedIter_ < 0)
     {
-        FatalErrorIn("layerParameters::layerParameters(..)")
+        FatalIOErrorIn("layerParameters::layerParameters(..)", dict)
             << "Layer iterations should be >= 0." << endl
             << "nLayerIter:" << nLayerIter_
             << " nRelaxedIter:" << nRelaxedIter_
-            << exit(FatalError);
+            << exit(FatalIOError);
     }
 
 
@@ -130,7 +275,7 @@ Foam::layerParameters::layerParameters
     {
         if (iter().isDict())
         {
-            const word& key = iter().keyword();
+            const keyType& key = iter().keyword();
             const labelHashSet patchIDs
             (
                 boundaryMesh.patchSet(List<wordRe>(1, key))
@@ -154,16 +299,69 @@ Foam::layerParameters::layerParameters
                     numLayers_[patchI] =
                         readLabel(layerDict.lookup("nSurfaceLayers"));
 
-                    layerDict.readIfPresent
+                    switch (layerSpec_)
-                    (
+                    {
-                        "expansionRatio",
+                        case FIRST_AND_TOTAL:
-                        expansionRatio_[patchI]
+                            layerDict.readIfPresent
-                    );
+                            (
-                    layerDict.readIfPresent
+                                "firstLayerThickness",
-                    (
+                                firstLayerThickness_[patchI]
-                        "finalLayerThickness",
+                            );
-                        finalLayerThickness_[patchI]
+                            layerDict.readIfPresent
-                    );
+                            (
+                                "thickness",
+                                thickness_[patchI]
+                            );
+                        break;
+
+                        case FIRST_AND_EXPANSION:
+                            layerDict.readIfPresent
+                            (
+                                "firstLayerThickness",
+                                firstLayerThickness_[patchI]
+                            );
+                            layerDict.readIfPresent
+                            (
+                                "expansionRatio",
+                                expansionRatio_[patchI]
+                            );
+                        break;
+
+                        case FINAL_AND_TOTAL:
+                            layerDict.readIfPresent
+                            (
+                                "finalLayerThickness",
+                                finalLayerThickness_[patchI]
+                            );
+                            layerDict.readIfPresent
+                            (
+                                "thickness",
+                                thickness_[patchI]
+                            );
+                        break;
+
+                        case FINAL_AND_EXPANSION:
+                            layerDict.readIfPresent
+                            (
+                                "finalLayerThickness",
+                                finalLayerThickness_[patchI]
+                            );
+                            layerDict.readIfPresent
+                            (
+                                "expansionRatio",
+                                expansionRatio_[patchI]
+                            );
+                        break;
+
+                        default:
+                            FatalIOErrorIn
+                            (
+                                "layerParameters::layerParameters(..)",
+                                dict
+                            )   << "problem." << exit(FatalIOError);
+                        break;
+                    }
+
                     layerDict.readIfPresent
                     (
                         "minThickness",
@@ -176,4 +374,190 @@ Foam::layerParameters::layerParameters
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::layerParameters::layerThickness
+(
+    const label nLayers,
+    const scalar firstLayerThickess,
+    const scalar finalLayerThickess,
+    const scalar totalThickness,
+    const scalar expansionRatio
+) const
+{
+    switch (layerSpec_)
+    {
+        case FIRST_AND_TOTAL:
+        case FINAL_AND_TOTAL:
+        {
+            return totalThickness;
+        }
+        break;
+
+        case FIRST_AND_EXPANSION:
+        {
+            if (mag(expansionRatio-1) < SMALL)
+            {
+                return firstLayerThickess * nLayers;
+            }
+            else
+            {
+                return firstLayerThickess *
+                    (1.0 - pow(expansionRatio, nLayers))
+                  / (1.0 - expansionRatio);
+            }
+        }
+        break;
+
+        case FINAL_AND_EXPANSION:
+        {
+            if (mag(expansionRatio-1) < SMALL)
+            {
+                return finalLayerThickess * nLayers;
+            }
+            else
+            {
+                scalar invExpansion = 1.0 / expansionRatio;
+                return finalLayerThickess *
+                    (1.0 - pow(invExpansion, nLayers))
+                  / (1.0 - invExpansion);
+            }
+        }
+        break;
+
+        default:
+        {
+            FatalErrorIn("layerParameters::layerThickness(..)")
+                << "Illegal thickness specification " <<    layerSpec_
+                << exit(FatalError);
+            return -VGREAT;
+        }
+    }
+}
+
+
+Foam::scalar Foam::layerParameters::layerExpansionRatio
+(
+    const label nLayers,
+    const scalar firstLayerThickess,
+    const scalar finalLayerThickess,
+    const scalar totalThickness,
+    const scalar expansionRatio
+) const
+{
+    switch (layerSpec_)
+    {
+        case FIRST_AND_EXPANSION:
+        case FINAL_AND_EXPANSION:
+        {
+            return expansionRatio;
+        }
+        break;
+
+        case FIRST_AND_TOTAL:
+        {
+            return layerExpansionRatio
+            (
+                nLayers,
+                totalThickness/firstLayerThickess
+            );
+        }
+        break;
+
+        case FINAL_AND_TOTAL:
+        {
+            return
+                1.0
+              / layerExpansionRatio
+                (
+                    nLayers,
+                    totalThickness/finalLayerThickess
+                );
+        }
+        break;
+
+        default:
+        {
+            FatalErrorIn("layerParameters::layerThickness(..)")
+                << "Illegal thickness specification" << exit(FatalError);
+            return -VGREAT;
+        }
+    }
+}
+
+
+Foam::scalar Foam::layerParameters::firstLayerThickness
+(
+    const label nLayers,
+    const scalar firstLayerThickess,
+    const scalar finalLayerThickess,
+    const scalar totalThickness,
+    const scalar expansionRatio
+) const
+{
+    switch (layerSpec_)
+    {
+        case FIRST_AND_EXPANSION:
+        case FIRST_AND_TOTAL:
+        {
+            return firstLayerThickess;
+        }
+
+        case FINAL_AND_EXPANSION:
+        {
+            return finalLayerThickess*pow(1.0/expansionRatio, nLayers-1);
+        }
+        break;
+
+        case FINAL_AND_TOTAL:
+        {
+            scalar r = layerExpansionRatio
+            (
+                nLayers,
+                firstLayerThickess,
+                finalLayerThickess,
+                totalThickness,
+                expansionRatio
+            );
+            return finalLayerThickess/pow(r, nLayers-1);
+        }
+        break;
+
+        default:
+        {
+            FatalErrorIn("layerParameters::layerThickness(..)")
+                << "Illegal thickness specification" << exit(FatalError);
+            return -VGREAT;
+        }
+    }
+}
+
+
+Foam::scalar Foam::layerParameters::finalLayerThicknessRatio
+(
+    const label nLayers,
+    const scalar expansionRatio
+) const
+{
+    if (nLayers > 0)
+    {
+        if (mag(expansionRatio-1) < SMALL)
+        {
+            return 1.0/nLayers;
+        }
+        else
+        {
+            return
+                pow(expansionRatio, nLayers - 1)
+              * (1.0 - expansionRatio)
+              / (1.0 - pow(expansionRatio, nLayers));
+        }
+    }
+    else
+    {
+        return 0.0;
+    }
+}
+
+
 // ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H

@@ -55,6 +55,27 @@ class refinementSurfaces;
 
 class layerParameters
 {
+public:
+
+    // Public data types
+
+        //- Enumeration defining the layer specification:
+        //  - first and total thickness specified
+        //  - first and expansion ratio specified
+        //  - final and total thickness specified
+        //  - final and expansion ratio specified
+        enum layerSpecification
+        {
+            ILLEGAL,
+            FIRST_AND_TOTAL,
+            FIRST_AND_EXPANSION,
+            FINAL_AND_TOTAL,
+            FINAL_AND_EXPANSION
+        };
+
+
+private:
+
     // Static data members
 
         //- Default angle for faces to be convcave
@@ -68,12 +89,21 @@ class layerParameters
             //- How many layers to add.
             labelList numLayers_;
 
-            scalarField expansionRatio_;
+            //- Are sizes relative to local cell size
-
             Switch relativeSizes_;
 
-            scalarField finalLayerThickness_;
+            //- How thickness is specified.
+            layerSpecification layerSpec_;
 
+                scalarField firstLayerThickness_;
+
+                scalarField finalLayerThickness_;
+
+                scalarField thickness_;
+
+                scalarField expansionRatio_;
+
+            //- Minimum total thickness
             scalarField minThickness_;
 
 
@@ -111,6 +141,14 @@ class layerParameters
 
     // Private Member Functions
 
+        //- Calculate expansion ratio from overall size v.s. thickness of
+        //  first layer.
+        scalar layerExpansionRatio
+        (
+            const label n,
+            const scalar totalOverFirst
+        ) const;
+
         //- Disallow default bitwise copy construct
         layerParameters(const layerParameters&);
 
@@ -128,154 +166,209 @@ public:
 
     // Member Functions
 
-        // Access
+        // Per patch information
 
-            // Per patch information
+            //- How many layers to add.
-
+            //  -1 : no specification. Assume 0 layers but allow sliding
-                //- How many layers to add.
+            //       to make layers
-                //  -1 : no specification. Assume 0 layers but allow sliding
+            //   0 : specified to have 0 layers. No sliding allowed.
-                //       to make layers
+            //  >0 : number of layers
-                //   0 : specified to have 0 layers. No sliding allowed.
+            const labelList& numLayers() const
-                //  >0 : number of layers
-                const labelList& numLayers() const
-                {
-                    return numLayers_;
-                }
-
-                // Expansion factor for layer mesh
-                const scalarField& expansionRatio() const
-                {
-                    return expansionRatio_;
-                }
-
-                //- Are size parameters relative to inner cell size or
-                //  absolute distances.
-                bool relativeSizes() const
-                {
-                    return relativeSizes_;
-                }
-
-                //- Wanted thickness of final added cell layer. If multiple
-                //  layers is the thickness of the layer furthest away
-                //  from the wall (i.e. nearest the original mesh)
-                //  If relativeSize() this number is relative to undistorted
-                //  size of the cell outside layer.
-                const scalarField& finalLayerThickness() const
-                {
-                    return finalLayerThickness_;
-                }
-
-                //- Minimum thickness of cell layer. If for any reason layer
-                //  cannot be above minThickness do not add layer.
-                //  If relativeSize() this number is relative to undistorted
-                //  size of the cell outside layer.
-                const scalarField& minThickness() const
-                {
-                    return minThickness_;
-                }
-
-
-            scalar featureAngle() const
             {
-                return featureAngle_;
+                return numLayers_;
             }
 
-            //- At non-patched sides allow mesh to slip if extrusion
+            //- Are size parameters relative to inner cell size or
-            //  direction makes angle larger than slipFeatureAngle.
+            //  absolute distances.
-            scalar slipFeatureAngle() const
+            bool relativeSizes() const
             {
-                return slipFeatureAngle_;
+                return relativeSizes_;
             }
 
-            scalar concaveAngle() const
+            // Expansion factor for layer mesh
+            const scalarField& expansionRatio() const
             {
-                return concaveAngle_;
+                return expansionRatio_;
             }
 
-            //- If points get not extruded do nGrow layers of connected faces
+            //- Wanted thickness of the layer furthest away
-            //  that are not grown. Is used to not do layers at all close to
+            //  from the wall (i.e. nearest the original mesh).
-            //  features.
+            //  If relativeSize() this number is relative to undistorted
-            label nGrow() const
+            //  size of the cell outside layer.
+            const scalarField& finalLayerThickness() const
             {
-                return nGrow_;
+                return finalLayerThickness_;
             }
 
-            //- Number of smoothing iterations of surface normals
+            //- Wanted thickness of the layer nearest to the wall.
-            label nSmoothSurfaceNormals() const
+            //  If relativeSize() this number is relative to undistorted
+            //  size of the cell outside layer.
+            const scalarField& firstLayerThickness() const
             {
-                return nSmoothSurfaceNormals_;
+                return firstLayerThickness_;
             }
 
-            //- Number of smoothing iterations of interior mesh movement
+            //- Wanted overall thickness of all layers.
-            //  direction
+            //  If relativeSize() this number is relative to undistorted
-            label nSmoothNormals() const
+            //  size of the cell outside layer.
+            const scalarField& thickness() const
             {
-                return nSmoothNormals_;
+                return thickness_;
             }
 
-            //- Stop layer growth on highly warped cells
+            //- Minimum overall thickness of cell layer. If for any reason layer
-            scalar maxFaceThicknessRatio() const
+            //  cannot be above minThickness do not add layer.
+            //  If relativeSize() this number is relative to undistorted
+            //  size of the cell outside layer.
+            const scalarField& minThickness() const
             {
-                return maxFaceThicknessRatio_;
+                return minThickness_;
-            }
-
-            scalar layerTerminationCos() const
-            {
-                return layerTerminationCos_;
-            }
-
-            //- Smooth layer thickness over surface patches
-            label nSmoothThickness() const
-            {
-                return nSmoothThickness_;
-            }
-
-            //- Reduce layer growth where ratio thickness to medial
-            //  distance is large
-            scalar maxThicknessToMedialRatio() const
-            {
-                return maxThicknessToMedialRatio_;
-            }
-
-            //- Angle used to pick up medial axis points
-            scalar minMedianAxisAngleCos() const
-            {
-                return minMedianAxisAngleCos_;
-            }
-
-            //- Create buffer region for new layer terminations
-            label nBufferCellsNoExtrude() const
-            {
-                return nBufferCellsNoExtrude_;
-            }
-
-            label nSnap() const
-            {
-                return nSnap_;
-            }
-
-            const Switch& additionalReporting() const
-            {
-                return additionalReporting_;
             }
 
 
-            // Overall
+        scalar featureAngle() const
+        {
+            return featureAngle_;
+        }
 
-                //- Number of overall layer addition iterations
+        //- At non-patched sides allow mesh to slip if extrusion
-                label nLayerIter() const
+        //  direction makes angle larger than slipFeatureAngle.
-                {
+        scalar slipFeatureAngle() const
-                    return nLayerIter_;
+        {
-                }
+            return slipFeatureAngle_;
+        }
 
-                //- Number of iterations after which relaxed motion rules
+        scalar concaveAngle() const
-                //  are to be used.
+        {
-                label nRelaxedIter() const
+            return concaveAngle_;
-                {
+        }
-                    return nRelaxedIter_;
+
-                }
+        //- If points get not extruded do nGrow layers of connected faces
+        //  that are not grown. Is used to not do layers at all close to
+        //  features.
+        label nGrow() const
+        {
+            return nGrow_;
+        }
+
+        //- Number of smoothing iterations of surface normals
+        label nSmoothSurfaceNormals() const
+        {
+            return nSmoothSurfaceNormals_;
+        }
+
+        //- Number of smoothing iterations of interior mesh movement
+        //  direction
+        label nSmoothNormals() const
+        {
+            return nSmoothNormals_;
+        }
+
+        //- Stop layer growth on highly warped cells
+        scalar maxFaceThicknessRatio() const
+        {
+            return maxFaceThicknessRatio_;
+        }
+
+        scalar layerTerminationCos() const
+        {
+            return layerTerminationCos_;
+        }
+
+        //- Smooth layer thickness over surface patches
+        label nSmoothThickness() const
+        {
+            return nSmoothThickness_;
+        }
+
+        //- Reduce layer growth where ratio thickness to medial
+        //  distance is large
+        scalar maxThicknessToMedialRatio() const
+        {
+            return maxThicknessToMedialRatio_;
+        }
+
+        //- Angle used to pick up medial axis points
+        scalar minMedianAxisAngleCos() const
+        {
+            return minMedianAxisAngleCos_;
+        }
+
+        //- Create buffer region for new layer terminations
+        label nBufferCellsNoExtrude() const
+        {
+            return nBufferCellsNoExtrude_;
+        }
+
+        label nSnap() const
+        {
+            return nSnap_;
+        }
+
+        const Switch& additionalReporting() const
+        {
+            return additionalReporting_;
+        }
 
 
+        // Overall
+
+            //- Number of overall layer addition iterations
+            label nLayerIter() const
+            {
+                return nLayerIter_;
+            }
+
+            //- Number of iterations after which relaxed motion rules
+            //  are to be used.
+            label nRelaxedIter() const
+            {
+                return nRelaxedIter_;
+            }
+
+
+        // Helper
+
+            //- Determine overall thickness. Uses two of the four parameters
+            //  according to the layerSpecification
+            scalar layerThickness
+            (
+                const label nLayers,
+                const scalar firstLayerThickess,
+                const scalar finalLayerThickess,
+                const scalar totalThickness,
+                const scalar expansionRatio
+            ) const;
+
+            //- Determine expansion ratio. Uses two of the four parameters
+            //  according to the layerSpecification
+            scalar layerExpansionRatio
+            (
+                const label nLayers,
+                const scalar firstLayerThickess,
+                const scalar finalLayerThickess,
+                const scalar totalThickness,
+                const scalar expansionRatio
+            ) const;
+
+            //- Determine first layer (near-wall) thickness. Uses two of the
+            //  four parameters according to the layerSpecification
+            scalar firstLayerThickness
+            (
+                const label nLayers,
+                const scalar firstLayerThickess,
+                const scalar finalLayerThickess,
+                const scalar totalThickness,
+                const scalar expansionRatio
+            ) const;
+
+            //- Determine ratio of final layer thickness to
+            //  overall layer thickness
+            scalar finalLayerThicknessRatio
+            (
+                const label nLayers,
+                const scalar expansionRatio
+            ) const;
 
 };
 

src/regionModels/pyrolysisModels/Make/options

@@ -3,10 +3,10 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/solid/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidChemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
@@ -20,12 +20,13 @@ LIB_LIBS = \
     -lmeshTools \
     -lchemistryModel \
     -lspecie \
-    -lsolidSpecie \
     -lfluidThermophysicalModels \
     -lsolidChemistryModel \
+    -lsolidThermo \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lLESdeltas \
     -lregionModels \
-    -lradiationModels
+    -lradiationModels \
+    -lreactionThermophysicalModels

src/regionModels/thermoBaffleModels/Make/options

@@ -1,7 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/solid/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
@@ -18,5 +17,4 @@ LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -lOpenFOAM \
-    -lsolidSpecie \
     -lradiationModels

src/sampling/Make/files

@@ -18,6 +18,7 @@ sampledSet/sampledSets/sampledSetsGrouping.C
 sampledSet/sampledSetsFunctionObject/sampledSetsFunctionObject.C
 sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C
 sampledSet/uniform/uniformSet.C
+sampledSet/array/arraySet.C
 
 setWriters = sampledSet/writers
 

src/sampling/sampledSet/array/arraySet.C

@@ -0,0 +1,195 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "arraySet.H"
+#include "sampledSet.H"
+#include "meshSearch.H"
+#include "DynamicList.H"
+#include "polyMesh.H"
+#include "addToRunTimeSelectionTable.H"
+#include "word.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(arraySet, 0);
+    addToRunTimeSelectionTable(sampledSet, arraySet, word);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::arraySet::calcSamples
+(
+    DynamicList<point>& samplingPts,
+    DynamicList<label>& samplingCells,
+    DynamicList<label>& samplingFaces,
+    DynamicList<label>& samplingSegments,
+    DynamicList<scalar>& samplingCurveDist
+) const
+{
+    const meshSearch& queryMesh = searchEngine();
+
+    label nTotalSamples
+    (
+        pointsDensity_.x()
+       *pointsDensity_.y()
+       *pointsDensity_.z()
+    );
+
+    List<point> sampleCoords(nTotalSamples);
+
+    const scalar deltax = spanBox_.x()/(pointsDensity_.x() + 1);
+    const scalar deltay = spanBox_.y()/(pointsDensity_.y() + 1);
+    const scalar deltaz = spanBox_.z()/(pointsDensity_.z() + 1);
+
+    label p(0);
+    for (label k=1; k<=pointsDensity_.z(); k++)
+    {
+        for (label j=1; j<=pointsDensity_.y(); j++)
+        {
+            for (label i=1; i<=pointsDensity_.x(); i++)
+            {
+                vector t(deltax*i , deltay*j, deltaz*k);
+                sampleCoords[p] = coordSys_.origin() + t;
+                p++;
+            }
+        }
+    }
+
+    forAll(sampleCoords, i)
+    {
+        sampleCoords[i] = transform(coordSys_.R(), sampleCoords[i]);
+    }
+
+    forAll(sampleCoords, sampleI)
+    {
+        label cellI = queryMesh.findCell(sampleCoords[sampleI]);
+
+        if (cellI != -1)
+        {
+            samplingPts.append(sampleCoords[sampleI]);
+            samplingCells.append(cellI);
+            samplingFaces.append(-1);
+            samplingSegments.append(0);
+            samplingCurveDist.append(1.0 * sampleI);
+        }
+    }
+}
+
+
+void Foam::arraySet::genSamples()
+{
+    // Storage for sample points
+    DynamicList<point> samplingPts;
+    DynamicList<label> samplingCells;
+    DynamicList<label> samplingFaces;
+    DynamicList<label> samplingSegments;
+    DynamicList<scalar> samplingCurveDist;
+
+    calcSamples
+    (
+        samplingPts,
+        samplingCells,
+        samplingFaces,
+        samplingSegments,
+        samplingCurveDist
+    );
+
+    samplingPts.shrink();
+    samplingCells.shrink();
+    samplingFaces.shrink();
+    samplingSegments.shrink();
+    samplingCurveDist.shrink();
+
+    setSamples
+    (
+        samplingPts,
+        samplingCells,
+        samplingFaces,
+        samplingSegments,
+        samplingCurveDist
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::arraySet::arraySet
+(
+    const word& name,
+    const polyMesh& mesh,
+    const meshSearch& searchEngine,
+    const word& axis,
+    const coordinateSystem& origin,
+    const Vector<label>& pointsDensity,
+    const Vector<scalar>& spanBox
+)
+:
+    sampledSet(name, mesh, searchEngine, axis),
+    coordSys_(origin),
+    pointsDensity_(pointsDensity),
+    spanBox_(spanBox)
+{
+    genSamples();
+
+    if (debug)
+    {
+        write(Info);
+    }
+}
+
+
+Foam::arraySet::arraySet
+(
+    const word& name,
+    const polyMesh& mesh,
+    const meshSearch& searchEngine,
+    const dictionary& dict
+)
+:
+    sampledSet(name, mesh, searchEngine, dict),
+    coordSys_(dict),
+    pointsDensity_(dict.lookup("pointsDensity")),
+    spanBox_(dict.lookup("spanBox"))
+{
+    genSamples();
+
+    if (debug)
+    {
+        write(Info);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::arraySet::~arraySet()
+{}
+
+
+// ************************************************************************* //

src/sampling/sampledSet/array/arraySet.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 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,103 +22,99 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::IrreversibleSolidReaction
+    Foam::arraySet
 
 Description
-    Simple extension of Reaction to handle reversible reactions
 
 SourceFiles
-    IrreversibleSolidReaction.C
+    arraySet.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef IrreversibleSolidReaction_H
+#ifndef arraySet_H
-#define IrreversibleSolidReaction_H
+#define arraySet_H
 
-#include "solidReaction.H"
+#include "sampledSet.H"
+#include "DynamicList.H"
+#include "coordinateSystem.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+// Forward declaration of classes
+class passiveParticle;
+template<class Type> class particle;
+
 /*---------------------------------------------------------------------------*\
-                      Class IrreversibleSolidReaction Declaration
+                           Class arraySet Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class ReactionRate>
+class arraySet
-class IrreversibleSolidReaction
 :
-    public solidReaction
+    public sampledSet
 {
     // Private data
 
-        // Reaction rate
+        //- Coordinate syste
-        ReactionRate k_;
+        coordinateSystem coordSys_;
 
-        // Reaction order
+        //- Point density vector
-        scalar nReact_;
+        Vector<label> pointsDensity_;
+
+        //- Span box
+        Vector<scalar> spanBox_;
 
 
     // Private Member Functions
 
-        //- Disallow default bitwise assignment
+        //- Samples all points in sampleCoords.
-        void operator=
+        void calcSamples
         (
-            const IrreversibleSolidReaction<ReactionRate>&
+            DynamicList<point>& samplingPts,
-        );
+            DynamicList<label>& samplingCells,
+            DynamicList<label>& samplingFaces,
+            DynamicList<label>& samplingSegments,
+            DynamicList<scalar>& samplingCurveDist
+        ) const;
+
+        //- Uses calcSamples to obtain samples. Copies them into *this.
+        void genSamples();
 
 
 public:
 
     //- Runtime type information
-    TypeName("irreversible");
+    TypeName("array");
 
 
     // Constructors
 
         //- Construct from components
-        IrreversibleSolidReaction
+        arraySet
         (
-            const solidReaction& reaction,
+            const word& name,
-            const ReactionRate& reactionRate,
+            const polyMesh& mesh,
-            const scalar nReact
+            const meshSearch& searchEngine,
+            const word& axis,
+            const coordinateSystem& coordSys,
+            const Vector<label>& pointsDensity,
+            const Vector<scalar>& spanBox
         );
 
-
+        //- Construct from dictionary
-        //- Construct from Istream
+        arraySet
-        IrreversibleSolidReaction
         (
-            const speciesTable& components,
+            const word& name,
-            Istream& is,
+            const polyMesh& mesh,
-            const speciesTable& pyrolysisGases
+            const meshSearch& searchEngine,
+            const dictionary& dict
         );
 
 
     //- Destructor
-    virtual ~IrreversibleSolidReaction()
+    virtual ~arraySet();
-    {}
-
-
-    // Member Functions
-
-        // IrreversibleSolidReaction rate coefficients
-
-            //- Forward rate constant
-            virtual scalar kf
-            (
-                const scalar p,
-                const scalar T,
-                const scalarField& c
-            ) const;
-
-
-            //- Reaction order
-            virtual scalar nReact() const;
-
-
-        //- Write
-        virtual void write(Ostream&) const;
 };
 
 
@@ -128,12 +124,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-#   include "IrreversibleSolidReaction.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C

@@ -164,17 +164,19 @@ void Foam::sampledSurfaces::write()
             writeGeometry();
         }
 
-        sampleAndWrite<volScalarField>(mesh_);
+        const IOobjectList objects(mesh_, mesh_.time().timeName());
-        sampleAndWrite<volVectorField>(mesh_);
-        sampleAndWrite<volSphericalTensorField>(mesh_);
-        sampleAndWrite<volSymmTensorField>(mesh_);
-        sampleAndWrite<volTensorField>(mesh_);
 
-        sampleAndWrite<surfaceScalarField>(mesh_);
+        sampleAndWrite<volScalarField>(objects);
-        sampleAndWrite<surfaceVectorField>(mesh_);
+        sampleAndWrite<volVectorField>(objects);
-        sampleAndWrite<surfaceSphericalTensorField>(mesh_);
+        sampleAndWrite<volSphericalTensorField>(objects);
-        sampleAndWrite<surfaceSymmTensorField>(mesh_);
+        sampleAndWrite<volSymmTensorField>(objects);
-        sampleAndWrite<surfaceTensorField>(mesh_);
+        sampleAndWrite<volTensorField>(objects);
+
+        sampleAndWrite<surfaceScalarField>(objects);
+        sampleAndWrite<surfaceVectorField>(objects);
+        sampleAndWrite<surfaceSphericalTensorField>(objects);
+        sampleAndWrite<surfaceSymmTensorField>(objects);
+        sampleAndWrite<surfaceTensorField>(objects);
     }
 }
 

src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H

@@ -160,7 +160,7 @@ class sampledSurfaces
 
         //- Sample and write all sampled fields
         template<class Type>
-        void sampleAndWrite(const fvMesh&);
+        void sampleAndWrite(const IOobjectList&);
 
         //- Disallow default bitwise copy construct and assignment
         sampledSurfaces(const sampledSurfaces&);

src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C

@@ -166,13 +166,18 @@ void Foam::sampledSurfaces::sampleAndWrite
 
 
 template<class GeoField>
-void Foam::sampledSurfaces::sampleAndWrite(const fvMesh& mesh)
+void Foam::sampledSurfaces::sampleAndWrite(const IOobjectList& allObjects)
 {
     forAll (fieldSelection_, fieldI)
     {
         const wordRe field = fieldSelection_[fieldI];
+        IOobject* fieldIOPtr = allObjects.lookup(field);
 
-        if (mesh.thisDb().foundObject<GeoField>(field))
+        if
+        (
+            fieldIOPtr != NULL
+         && fieldIOPtr->headerClassName() == GeoField::typeName
+        )
         {
             if (Pstream::master() && verbose_)
             {
@@ -181,17 +186,25 @@ void Foam::sampledSurfaces::sampleAndWrite(const fvMesh& mesh)
 
             if (loadFromFiles_)
             {
-                const GeoField& geoField =
+                const GeoField geoField
-                    mesh.thisDb().lookupObject<GeoField>(field);
+                (
+                    IOobject
+                    (
+                        field,
+                        mesh_.time().timeName(),
+                        mesh_,
+                        IOobject::MUST_READ
+                    ),
+                    mesh_
+                );
 
-                const_cast<GeoField&>(geoField).readOpt() = IOobject::MUST_READ;
                 sampleAndWrite(geoField);
             }
             else
             {
                 sampleAndWrite
                 (
-                    mesh.thisDb().lookupObject<GeoField>(field)
+                    mesh_.thisDb().lookupObject<GeoField>(field)
                 );
             }
         }

src/thermophysicalModels/Allwmake

@@ -4,6 +4,7 @@ makeType=${1:-libso}
 set -x
 
 wmake $makeType specie
+wmake $makeType solidSpecie
 wmake $makeType thermophysicalFunctions
 ./properties/Allwmake $*
 
@@ -14,7 +15,7 @@ wmake $makeType chemistryModel
 wmake $makeType barotropicCompressibilityModel
 wmake $makeType SLGThermo
 
-wmake $makeType solidSpecie
+
 wmake $makeType solidThermo
 wmake $makeType solidChemistryModel
 

src/thermophysicalModels/basic/mixtures/pureMixture/pureMixture.H

@@ -87,6 +87,27 @@ public:
             return mixture_;
         }
 
+        const ThermoType& cellVolMixture
+        (
+            const scalar,
+            const scalar,
+            const label
+        ) const
+        {
+            return mixture_;
+        }
+
+        const ThermoType& patchFaceVolMixture
+        (
+            const scalar,
+            const scalar,
+            const label,
+            const label
+        ) const
+        {
+            return mixture_;
+        }
+
         //- Read dictionary
         void read(const dictionary&);
 };

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

@@ -27,7 +27,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "unitConversion.H"
 #include "zeroGradientFvPatchFields.H"
-#include "basicSolidMixture.H"
+
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -71,7 +71,7 @@ greyMeanSolidAbsorptionEmission::X(const word specie) const
         }
     }
     const scalarField& Yj = mixture_.Y(specie);
-    const label mySpecieI = mixture_.components()[specie];
+    const label mySpecieI = mixture_.species()[specie];
     forAll(Xj, iCell)
     {
         Xj[iCell] = Yj[iCell]/mixture_.rho(mySpecieI, p[iCell], T[iCell]);
@@ -93,10 +93,10 @@ greyMeanSolidAbsorptionEmission
     coeffsDict_((dict.subDict(typeName + "Coeffs"))),
     thermo_(mesh.lookupObject<solidThermo>("thermophysicalProperties")),
     speciesNames_(0),
-    mixture_(dynamic_cast<const basicSolidMixture&>(thermo_)),
+    mixture_(dynamic_cast<const basicMultiComponentMixture&>(thermo_)),
     solidData_(mixture_.Y().size())
 {
-    if (!isA<basicSolidMixture>(thermo_))
+    if (!isA<basicMultiComponentMixture>(thermo_))
     {
         FatalErrorIn
         (

src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.H

@@ -90,7 +90,7 @@ private:
         HashTable<label> speciesNames_;
 
         //- Basic multicomponent mixture
-        const basicSolidMixture& mixture_;
+        const basicMultiComponentMixture& mixture_;
 
         //- List of solid species data
         List<FixedList<scalar, 2> > solidData_;

src/thermophysicalModels/reactionThermo/Make/options

@@ -2,10 +2,13 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
+
 LIB_LIBS = \
     -lfiniteVolume \
     -lfluidThermophysicalModels \
     -lspecie \
+    -lsolidSpecie \
     -lmeshTools

src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/makeChemistryReaders.C

@@ -25,6 +25,7 @@ License
 
 #include "makeReactionThermo.H"
 #include "thermoPhysicsTypes.H"
+#include "solidThermoPhysicsTypes.H"
 
 #include "chemistryReader.H"
 #include "foamChemistryReader.H"
@@ -41,6 +42,8 @@ makeChemistryReader(gasThermoPhysics);
 makeChemistryReader(constIncompressibleGasThermoPhysics);
 makeChemistryReader(incompressibleGasThermoPhysics);
 makeChemistryReader(icoPoly8ThermoPhysics);
+makeChemistryReader(hConstSolidThermoPhysics);
+makeChemistryReader(hExponentialSolidThermoPhysics);
 
 makeChemistryReaderType(foamChemistryReader, constGasThermoPhysics);
 makeChemistryReaderType(foamChemistryReader, gasThermoPhysics);
@@ -51,6 +54,8 @@ makeChemistryReaderType
 );
 makeChemistryReaderType(foamChemistryReader, incompressibleGasThermoPhysics);
 makeChemistryReaderType(foamChemistryReader, icoPoly8ThermoPhysics);
+makeChemistryReaderType(foamChemistryReader, hConstSolidThermoPhysics);
+makeChemistryReaderType(foamChemistryReader, hExponentialSolidThermoPhysics);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -41,6 +41,7 @@ License
 #include "powerSeriesReactionRate.H"
 #include "addToRunTimeSelectionTable.H"
 
+
 /* * * * * * * * * * * * * * * * * Static data * * * * * * * * * * * * * * * */
 
 namespace Foam
@@ -178,7 +179,8 @@ void Foam::chemkinReader::addReactionType
         {
             reactions_.append
             (
-                new IrreversibleReaction<gasThermoPhysics, ReactionRateType>
+                new IrreversibleReaction
+                <Reaction, gasThermoPhysics, ReactionRateType>
                 (
                     Reaction<gasThermoPhysics>
                     (
@@ -197,7 +199,8 @@ void Foam::chemkinReader::addReactionType
         {
             reactions_.append
             (
-                new ReversibleReaction<gasThermoPhysics, ReactionRateType>
+                new ReversibleReaction
+                <Reaction, gasThermoPhysics, ReactionRateType>
                 (
                     Reaction<gasThermoPhysics>
                     (
@@ -496,7 +499,7 @@ void Foam::chemkinReader::addReaction
                 reactions_.append
                 (
                     new NonEquilibriumReversibleReaction
-                        <gasThermoPhysics, ArrheniusReactionRate>
+                        <Reaction, gasThermoPhysics, ArrheniusReactionRate>
                     (
                         Reaction<gasThermoPhysics>
                         (
@@ -549,7 +552,11 @@ void Foam::chemkinReader::addReaction
                 reactions_.append
                 (
                     new NonEquilibriumReversibleReaction
-                        <gasThermoPhysics, thirdBodyArrheniusReactionRate>
+                    <
+                        Reaction,
+                        gasThermoPhysics,
+                        thirdBodyArrheniusReactionRate
+                    >
                     (
                         Reaction<gasThermoPhysics>
                         (
@@ -654,7 +661,7 @@ void Foam::chemkinReader::addReaction
                 reactions_.append
                 (
                     new NonEquilibriumReversibleReaction
-                        <gasThermoPhysics, LandauTellerReactionRate>
+                        <Reaction, gasThermoPhysics, LandauTellerReactionRate>
                     (
                         Reaction<gasThermoPhysics>
                         (

src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C

@@ -77,7 +77,8 @@ Foam::multiComponentMixture<ThermoType>::multiComponentMixture
 :
     basicMultiComponentMixture(thermoDict, specieNames, mesh),
     speciesData_(species_.size()),
-    mixture_("mixture", *thermoData[specieNames[0]])
+    mixture_("mixture", *thermoData[specieNames[0]]),
+    mixtureVol_("volMixture", *thermoData[specieNames[0]])
 {
     forAll(species_, i)
     {
@@ -101,7 +102,8 @@ Foam::multiComponentMixture<ThermoType>::multiComponentMixture
 :
     basicMultiComponentMixture(thermoDict, thermoDict.lookup("species"), mesh),
     speciesData_(species_.size()),
-    mixture_("mixture", constructSpeciesData(thermoDict))
+    mixture_("mixture", constructSpeciesData(thermoDict)),
+    mixtureVol_("volMixture", speciesData_[0])
 {
     correctMassFractions();
 }
@@ -148,6 +150,65 @@ const ThermoType& Foam::multiComponentMixture<ThermoType>::patchFaceMixture
 }
 
 
+template<class ThermoType>
+const ThermoType& Foam::multiComponentMixture<ThermoType>::cellVolMixture
+(
+    const scalar p,
+    const scalar T,
+    const label celli
+) const
+{
+    scalar rhoInv = 0.0;
+    forAll(speciesData_, i)
+    {
+        rhoInv += Y_[i][celli]/speciesData_[i].rho(p, T);
+    }
+
+    mixtureVol_ =
+        Y_[0][celli]/speciesData_[0].rho(p, T)/rhoInv*speciesData_[0];
+
+    for (label n=1; n<Y_.size(); n++)
+    {
+        mixtureVol_ +=
+            Y_[n][celli]/speciesData_[n].rho(p, T)/rhoInv*speciesData_[n];
+    }
+
+    return mixtureVol_;
+}
+
+
+template<class ThermoType>
+const ThermoType& Foam::multiComponentMixture<ThermoType>::
+patchFaceVolMixture
+(
+    const scalar p,
+    const scalar T,
+    const label patchi,
+    const label facei
+) const
+{
+    scalar rhoInv = 0.0;
+    forAll(speciesData_, i)
+    {
+        rhoInv +=
+            Y_[i].boundaryField()[patchi][facei]/speciesData_[i].rho(p, T);
+    }
+
+    mixtureVol_ =
+        Y_[0].boundaryField()[patchi][facei]/speciesData_[0].rho(p, T)/rhoInv
+      * speciesData_[0];
+
+    for (label n=1; n<Y_.size(); n++)
+    {
+        mixtureVol_ +=
+            Y_[n].boundaryField()[patchi][facei]/speciesData_[n].rho(p,T)
+          / rhoInv*speciesData_[n];
+    }
+
+    return mixtureVol_;
+}
+
+
 template<class ThermoType>
 void Foam::multiComponentMixture<ThermoType>::read
 (

src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.H

@@ -60,6 +60,10 @@ class multiComponentMixture
         //- Temporary storage for the cell/face mixture thermo data
         mutable ThermoType mixture_;
 
+        //- Temporary storage for the volume weighted
+        //  cell/face mixture thermo data
+        mutable ThermoType mixtureVol_;
+
 
     // Private Member Functions
 
@@ -110,6 +114,21 @@ public:
             const label facei
         ) const;
 
+        const ThermoType& cellVolMixture
+        (
+            const scalar p,
+            const scalar T,
+            const label celli
+        ) const;
+
+        const ThermoType& patchFaceVolMixture
+        (
+            const scalar p,
+            const scalar T,
+            const label patchi,
+            const label facei
+        ) const;
+
         //- Return the raw specie thermodynamic data
         const PtrList<ThermoType>& speciesData() const
         {

src/thermophysicalModels/solidChemistryModel/Make/options

@@ -1,10 +1,11 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude  \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude
 
@@ -12,4 +13,5 @@ EXE_INC = \
 LIB_LIBS = \
     -lchemistryModel \
     -lfiniteVolume \
-    -lODE
+    -lODE\
+    -lreactionThermophysicalModels

src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C

@@ -24,7 +24,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "ODESolidChemistryModel.H"
-#include "reactingSolidMixture.H"
+#include "reactingMixture.H"
+#include "solidReaction.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -38,24 +39,20 @@ ODESolidChemistryModel
     CompType(mesh),
     ODE(),
     Ys_(this->solidThermo().composition().Y()),
-    pyrolisisGases_
-    (
-        mesh.lookupObject<dictionary>
-            ("thermophysicalProperties").lookup("gaseousSpecies")
-    ),
     reactions_
     (
-        dynamic_cast<const reactingSolidMixture<SolidThermo>& >
+        dynamic_cast<const reactingMixture<SolidThermo>& >
         (
             this->solidThermo()
         )
     ),
+    pyrolisisGases_(reactions_[0].gasSpecies()),
     solidThermo_
     (
-        dynamic_cast<const reactingSolidMixture<SolidThermo>& >
+        dynamic_cast<const reactingMixture<SolidThermo>& >
         (
             this->solidThermo()
-        ).solidData()
+        ).speciesData()
     ),
     gasThermo_(pyrolisisGases_.size()),
     nGases_(pyrolisisGases_.size()),
@@ -184,7 +181,7 @@ ODESolidChemistryModel
             mesh.lookupObject<dictionary>
             (
                 "thermophysicalProperties"
-            ).subDict(pyrolisisGases_[gasI] + "Coeffs");
+            ).subDict(pyrolisisGases_[gasI]);
 
         gasThermo_.set
         (
@@ -193,14 +190,13 @@ ODESolidChemistryModel
         );
     }
 
-    Info<< "ODESolidChemistryModel: Number of solids = " << nSolids_
+    Info<< "solidChemistryModel: Number of solids = " << nSolids_ << endl;
-        << " and reactions = " << nReaction_ << endl;
 
-    Info<< "Number of gases from pyrolysis = " << nGases_ << endl;
+    Info<< "solidChemistryModel: Number of gases = " << nGases_ << endl;
 
     forAll(reactions_, i)
     {
-        Info<< indent << "Reaction " << i << nl << reactions_[i] << nl;
+        Info<< indent << reactions_[i] << nl;
     }
 
 }
@@ -235,23 +231,23 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
 
     forAll(reactions_, i)
     {
-        const solidReaction& R = reactions_[i];
+        const Reaction<SolidThermo>& R = reactions_[i];
 
         scalar omegai = omega
         (
             R, c, T, p, pf, cf, lRef, pr, cr, rRef
         );
         scalar rhoL = 0.0;
-        forAll(R.slhs(), s)
+        forAll(R.lhs(), s)
         {
-            label si = R.slhs()[s];
+            label si = R.lhs()[s].index;
             om[si] -= omegai;
             rhoL = solidThermo_[si].rho(p, T);
         }
         scalar sr = 0.0;
-        forAll(R.srhs(), s)
+        forAll(R.rhs(), s)
         {
-            label si = R.srhs()[s];
+            label si = R.rhs()[s].index;
             scalar rhoR = solidThermo_[si].rho(p, T);
             sr = rhoR/rhoL;
             om[si] += sr*omegai;
@@ -263,7 +259,7 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
         }
         forAll(R.grhs(), g)
         {
-            label gi = R.grhs()[g];
+            label gi = R.grhs()[g].index;
             om[gi + nSolids_] += (1.0 - sr)*omegai;
         }
     }
@@ -276,7 +272,7 @@ template<class CompType, class SolidThermo, class GasThermo>
 Foam::scalar
 Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
 (
-    const solidReaction& R,
+    const Reaction<SolidThermo>& R,
     const scalarField& c,
     const scalar T,
     const scalar p,
@@ -299,16 +295,17 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
 
     scalar kf = R.kf(p, T, c1);
 
-    const scalar exponent = R.nReact();
+    //const scalar exponent = R.nReact();
 
-    const label Nl = R.slhs().size();
+    const label Nl = R.lhs().size();
 
     for (label s=0; s<Nl; s++)
     {
-        label si = R.slhs()[s];
+        label si = R.lhs()[s].index;
+        const scalar exp = R.lhs()[si].exponent;
 
         kf *=
-            pow(c1[si]/Ys0_[si][cellI], exponent)
+            pow(c1[si]/Ys0_[si][cellI], exp)
            *(Ys0_[si][cellI]);
     }
 
@@ -390,18 +387,18 @@ void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::jacobian
 
     for (label ri=0; ri<reactions_.size(); ri++)
     {
-        const solidReaction& R = reactions_[ri];
+        const Reaction<SolidThermo>& R = reactions_[ri];
 
         scalar kf0 = R.kf(p, T, c2);
 
-        forAll(R.slhs(), j)
+        forAll(R.lhs(), j)
         {
-            label sj = R.slhs()[j];
+            label sj = R.lhs()[j].index;
             scalar kf = kf0;
-            forAll(R.slhs(), i)
+            forAll(R.lhs(), i)
             {
-                label si = R.slhs()[i];
+                label si = R.lhs()[i].index;
-                scalar exp = R.nReact();
+                scalar exp = R.lhs()[i].exponent;
                 if (i == j)
                 {
                     if (exp < 1.0)
@@ -428,14 +425,14 @@ void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::jacobian
                 }
             }
 
-            forAll(R.slhs(), i)
+            forAll(R.lhs(), i)
             {
-                label si = R.slhs()[i];
+                label si = R.lhs()[i].index;
                 dfdc[si][sj] -= kf;
             }
-            forAll(R.srhs(), i)
+            forAll(R.rhs(), i)
             {
-                label si = R.srhs()[i];
+                label si = R.rhs()[i].index;
                 dfdc[si][sj] += kf;
             }
         }
@@ -551,7 +548,8 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::nEqns() const
 
 
 template<class CompType, class SolidThermo, class GasThermo>
-void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::calculate()
+void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
+calculate()
 {
     if (!this->chemistry_)
     {

src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.H

@@ -38,7 +38,7 @@ SourceFiles
 #ifndef ODESolidChemistryModel_H
 #define ODESolidChemistryModel_H
 
-#include "solidReaction.H"
+#include "Reaction.H"
 #include "ODE.H"
 #include "volFieldsFwd.H"
 #include "DimensionedField.H"
@@ -72,12 +72,12 @@ protected:
         //- Reference to solid mass fractions
         PtrList<volScalarField>& Ys_;
 
+        //- Reactions
+        const PtrList<Reaction<SolidThermo> >& reactions_;
+
         //- List of gas species present in reaction system
         speciesTable pyrolisisGases_;
 
-        //- Reactions
-        const PtrList<solidReaction>& reactions_;
-
         //- Thermodynamic data of solids
         const PtrList<SolidThermo>& solidThermo_;
 
@@ -149,7 +149,7 @@ public:
     // Member Functions
 
         //- The reactions
-        inline const PtrList<solidReaction>& reactions() const;
+        inline const PtrList<Reaction<SolidThermo> >& reactions() const;
 
         //- Thermodynamic data of gases
         inline const PtrList<GasThermo>& gasThermo() const;
@@ -180,7 +180,7 @@ public:
         //  species and charateristic times
         virtual scalar omega
         (
-            const solidReaction& r,
+            const Reaction<SolidThermo>& r,
             const scalarField& c,
             const scalar T,
             const scalar p,

src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModelI.H

@@ -45,7 +45,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg()
 
 
 template<class CompType, class SolidThermo, class GasThermo>
-inline const Foam::PtrList<Foam::solidReaction>&
+inline const Foam::PtrList<Foam::Reaction<SolidThermo> >&
 Foam::ODESolidChemistryModel<CompType, SolidThermo,GasThermo>::reactions() const
 {
     return reactions_;

src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolvers.C

@@ -51,6 +51,7 @@ namespace Foam
         hExponentialSolidThermoPhysics,
         gasThermoPhysics
     )
+
 }