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

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files

@@ -1,7 +1,6 @@
 fluid/compressibleCourantNo.C
 solid/solidRegionDiffNo.C
-
 chtMultiRegionFoam.C
 
-EXE = $(FOAM_APPBIN)/chtMultiRegionFoam
+EXE = $(FOAM_USER_APPBIN)/chtMultiRegionFoam
 

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options

@@ -1,14 +1,23 @@
+
 EXE_INC = \
-    -Ifluid \
-    -Isolid \
-    -Iinclude \
+    -I./fluid \
+    -I./solid \
+    -I./porousFluid \
+    -I./porousSolid \
+    -I./include \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude
+
 
 EXE_LIBS = \
     -lbasicThermophysicalModels \
@@ -16,7 +25,9 @@ EXE_LIBS = \
     -lspecie \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
-    -lcompressibleLESModels \
     -lmeshTools \
     -lfiniteVolume \
-    -lradiationModels
+    -lradiationModels \
+    -lfieldSources \
+    -lregionModels \
+    -lsampling

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -26,7 +26,13 @@ Application
 
 Description
     Combination of heatConductionFoam and buoyantFoam for conjugate heat
-    transfer between a solid region and fluid region
+    transfer between a solid region and fluid region. It includes
+    porous media in the primary fluid region treated explicitly.
+
+    It handles secondary fluid or solid circuits which can be coupled
+    thermally with the main fluid region. i.e radiators, etc.
+
+    The secondary fluid region is
 
 \*---------------------------------------------------------------------------*/
 
@@ -39,6 +45,9 @@ Description
 #include "solidRegionDiffNo.H"
 #include "basicSolidThermo.H"
 #include "radiationModel.H"
+#include "porousZones.H"
+#include "IObasicSourceList.H"
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,12 +60,16 @@ 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"
     #include "readSolidTimeControls.H"
 
@@ -65,14 +78,12 @@ int main(int argc, char *argv[])
     #include "solidRegionDiffusionNo.H"
     #include "setInitialMultiRegionDeltaT.H"
 
-
     while (runTime.run())
     {
         #include "readTimeControls.H"
         #include "readSolidTimeControls.H"
         #include "readPIMPLEControls.H"
 
-
         #include "compressibleMultiRegionCourantNo.H"
         #include "solidRegionDiffusionNo.H"
         #include "setMultiRegionDeltaT.H"
@@ -105,6 +116,24 @@ 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 "
@@ -113,6 +142,7 @@ int main(int argc, char *argv[])
                 #include "readSolidMultiRegionPIMPLEControls.H"
                 #include "solveSolid.H"
             }
+
         }
 
         runTime.write();

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

@@ -1,16 +1,21 @@
 EXE_INC = \
-    -I.. \
     -Ifluid \
     -Isolid \
+    -I./porousFluid \
+    -I./porousSolid \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
@@ -20,4 +25,7 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
-    -lradiationModels
+    -lradiationModels \
+    -lfieldSources \
+    -lregionModels \
+    -lsampling

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

@@ -36,6 +36,8 @@ Description
 #include "regionProperties.H"
 #include "basicSolidThermo.H"
 #include "radiationModel.H"
+#include "porousZones.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,9 +50,13 @@ 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"
 
@@ -68,6 +74,24 @@ 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/fluid/UEqn.H

@@ -3,10 +3,14 @@
     (
         fvm::div(phi, U)
       + turb.divDevRhoReff(U)
+      + sources(rho, U)
     );
 
     UEqn().relax();
 
+    // Add porous resistance
+    pZones.addResistance(UEqn());
+
     solve
     (
         UEqn()

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H

@@ -18,6 +18,9 @@
     PtrList<dimensionedScalar> rhoMax(fluidRegions.size());
     PtrList<dimensionedScalar> rhoMin(fluidRegions.size());
 
+    PtrList<IObasicSourceList> heatSources(fluidRegions.size());
+    PtrList<porousZones> porousZonesFluid(fluidRegions.size());
+
     // Populate fluid field pointer lists
     forAll(fluidRegions, i)
     {
@@ -209,6 +212,20 @@
                 )
             )
         );
+
+        Info<< "    Adding sources\n" << endl;
+        heatSources.set
+        (
+            i,
+            new IObasicSourceList(fluidRegions[i])
+        );
+
+        Info<< "    Adding porous zones\n" << endl;
+        porousZonesFluid.set
+        (
+            i,
+            new porousZones(fluidRegions[i])
+        );
     }
 
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidMeshes.H

@@ -1,8 +1,10 @@
-    PtrList<fvMesh> fluidRegions(rp.fluidRegionNames().size());
+    const wordList fluidNames(rp["fluid"]);
 
-    forAll(rp.fluidRegionNames(), i)
+    PtrList<fvMesh> fluidRegions(fluidNames.size());
+
+    forAll(fluidNames, i)
     {
-        Info<< "Create fluid mesh for region " << rp.fluidRegionNames()[i]
+        Info<< "Create fluid mesh for region " << fluidNames[i]
             << " for time = " << runTime.timeName() << nl << endl;
 
         fluidRegions.set
@@ -12,7 +14,7 @@
             (
                 IOobject
                 (
-                    rp.fluidRegionNames()[i],
+                    fluidNames[i],
                     runTime.timeName(),
                     runTime,
                     IOobject::MUST_READ

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H

@@ -7,6 +7,7 @@
      ==
       - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
       + rad.Sh(thermo)
+      + sources(rho, h)
     );
 
     hEqn.relax();

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H

@@ -3,9 +3,3 @@
     const int nNonOrthCorr =
         simple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
 
-    const bool momentumPredictor =
-        simple.lookupOrDefault("momentumPredictor", true);
-
-    const bool transonic =
-        simple.lookupOrDefault("transonic", false);
-

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H

@@ -12,6 +12,10 @@
     const volScalarField& psi = thermo.psi();
     volScalarField& h = thermo.h();
 
+    IObasicSourceList& sources = heatSources[i];
+
+    const porousZones& pZones = porousZonesFluid[i];
+
     const dimensionedScalar initialMass
     (
         "initialMass",

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

@@ -0,0 +1,11 @@
+    // 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

@@ -0,0 +1,181 @@
+    // Initialise porous field pointer lists
+    PtrList<basicRhoThermo> 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,
+            basicRhoThermo::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

@@ -0,0 +1,25 @@
+    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

@@ -0,0 +1,20 @@
+{
+    fvScalarMatrix hPorousEqn
+    (
+        fvm::div(porousPhi, porousH)
+      - fvm::Sp(fvc::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

@@ -0,0 +1,54 @@
+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

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

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

@@ -0,0 +1,31 @@
+
+    const fvMesh& porousMesh = porousFluidRegions[i];
+
+    basicRhoThermo& 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.h();
+
+    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

@@ -0,0 +1,11 @@
+//  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

@@ -0,0 +1,43 @@
+    // Initialise solid field pointer lists
+    PtrList<basicSolidThermo> 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,
+            basicSolidThermo::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

@@ -0,0 +1,24 @@
+    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

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

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

@@ -0,0 +1,20 @@
+    const fvMesh& mesh = porousSolidRegions[i];
+
+    basicSolidThermo& 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.K();
+    //tmp<volSymmTensorField> tkappa = thermo.directionalK()*betav;
+
+    const volScalarField& kappa = tkappa();
+    //const volSymmTensorField& K = tK();
+
+    volScalarField& T = thermo.T();
+
+    IObasicSourceList& sources = solidHeatSources[i];

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

@@ -0,0 +1,18 @@
+    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

@@ -0,0 +1,15 @@
+{
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        tmp<fvScalarMatrix> TEqn
+        (
+          - fvm::laplacian(betav*kappa, T, "laplacian(K,T)")
+          + sources(rho, T)
+        );
+
+        TEqn().relax();
+        TEqn().solve();
+    }
+
+    Info<< "Min/max T:" << min(T).value() << ' ' << max(T).value() << endl;
+}

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

@@ -1,8 +1,10 @@
-    PtrList<fvMesh> solidRegions(rp.solidRegionNames().size());
+    const wordList solidsNames(rp["solid"]);
 
-    forAll(rp.solidRegionNames(), i)
+    PtrList<fvMesh> solidRegions(solidsNames.size());
+
+    forAll(solidsNames, i)
     {
-        Info<< "Create solid mesh for region " << rp.solidRegionNames()[i]
+        Info<< "Create solid mesh for region " << solidsNames[i]
             << " for time = " << runTime.timeName() << nl << endl;
 
         solidRegions.set
@@ -12,16 +14,11 @@
             (
                 IOobject
                 (
-                    rp.solidRegionNames()[i],
+                    solidsNames[i],
                     runTime.timeName(),
                     runTime,
                     IOobject::MUST_READ
                 )
             )
         );
-
-        // Force calculation of geometric properties to prevent it being done
-        // later in e.g. some boundary evaluation
-        //(void)solidRegions[i].weights();
-        //(void)solidRegions[i].deltaCoeffs();
     }

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

@@ -4,10 +4,14 @@
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
       + turb.divDevRhoReff(U)
+      + sources(rho, U)
     );
 
     UEqn().relax();
 
+    // Add porous resistance
+    pZones.addResistance(UEqn());
+
     if (momentumPredictor)
     {
         solve

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

@@ -14,3 +14,18 @@
             CoNum
         );
     }
+
+    forAll (porousFluidRegions, porousI)
+    {
+        CoNum = max
+        (
+            compressibleCourantNo
+            (
+                porousFluidRegions[porousI],
+                runTime,
+                rhoPorous[porousI],
+                phiPorous[porousI]
+            ),
+            CoNum
+        );
+    }

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

@@ -15,6 +15,9 @@
 
     List<scalar> initialMassFluid(fluidRegions.size());
 
+    PtrList<IObasicSourceList> heatSources(fluidRegions.size());
+    PtrList<porousZones> porousZonesFluid(fluidRegions.size());
+
     // Populate fluid field pointer lists
     forAll(fluidRegions, i)
     {
@@ -22,6 +25,7 @@
             << fluidRegions[i].name() << nl << endl;
 
         Info<< "    Adding to thermoFluid\n" << endl;
+
         thermoFluid.set
         (
             i,
@@ -118,19 +122,16 @@
             )
         );
 
-        Info<< "    Adding to turbulence\n" << endl;
+        Info<< "    Adding turbulence \n" << endl;
         turbulence.set
         (
             i,
-            autoPtr<compressible::turbulenceModel>
-            (
             compressible::turbulenceModel::New
             (
                 rhoFluid[i],
                 UFluid[i],
                 phiFluid[i],
                 thermoFluid[i]
-                )
             ).ptr()
         );
 
@@ -197,4 +198,20 @@
                 fvc::ddt(thermoFluid[i].p())
             )
         );
+
+        Info<< "    Adding sources\n" << endl;
+        heatSources.set
+        (
+            i,
+            new IObasicSourceList(fluidRegions[i])
+        );
+
+        Info<< "    Adding porous zones\n" << endl;
+        porousZonesFluid.set
+        (
+            i,
+            new porousZones(fluidRegions[i])
+        );
     }
+
+

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

@@ -1,8 +1,10 @@
-    PtrList<fvMesh> fluidRegions(rp.fluidRegionNames().size());
+    const wordList fluidNames(rp["fluid"]);
 
-    forAll(rp.fluidRegionNames(), i)
+    PtrList<fvMesh> fluidRegions(fluidNames.size());
+
+    forAll(fluidNames, i)
     {
-        Info<< "Create fluid mesh for region " << rp.fluidRegionNames()[i]
+        Info<< "Create fluid mesh for region " << fluidNames[i]
             << " for time = " << runTime.timeName() << nl << endl;
 
         fluidRegions.set
@@ -12,7 +14,7 @@
             (
                 IOobject
                 (
-                    rp.fluidRegionNames()[i],
+                    fluidNames[i],
                     runTime.timeName(),
                     runTime,
                     IOobject::MUST_READ

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

@@ -8,6 +8,7 @@
         dpdt
       - (fvc::ddt(rho, K) + fvc::div(phi, K))
       + rad.Sh(thermo)
+      + sources(rho, h)
     );
 
     hEqn.relax();
@@ -19,4 +20,5 @@
 
     Info<< "Min/max T:" << min(thermo.T()).value() << ' '
         << max(thermo.T()).value() << endl;
+
 }

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

@@ -68,14 +68,11 @@
     // Correct velocity field
     U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
     U.correctBoundaryConditions();
-    K = 0.5*magSqr(U);
-
     p = p_rgh + rho*gh;
 
     // Update pressure time derivative
     dpdt = fvc::ddt(p);
 
-
     // Solve continuity
     #include "rhoEqn.H"
 

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

@@ -8,4 +8,3 @@
 
     const bool momentumPredictor =
         pimple.lookupOrDefault("momentumPredictor", true);
-

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

@@ -20,6 +20,10 @@
 
     radiation::radiationModel& rad = radiation[i];
 
+    IObasicSourceList& sources = heatSources[i];
+
+    const porousZones& pZones = porousZonesFluid[i];
+
     const dimensionedScalar initialMass
     (
         "initialMass",

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

@@ -0,0 +1,24 @@
+
+// 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

@@ -0,0 +1,147 @@
+    // Initialise porous field pointer lists
+    PtrList<basicRhoThermo> 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,
+            basicRhoThermo::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

@@ -0,0 +1,25 @@
+    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

@@ -0,0 +1,21 @@
+{
+    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

@@ -0,0 +1,61 @@
+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
+porousdpdt = fvc::ddt(porousP);
+

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

@@ -0,0 +1,11 @@
+    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

@@ -0,0 +1,18 @@
+
+    fvMesh& porousMesh = porousFluidRegions[i];
+
+    basicRhoThermo& 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.h();
+
+    IObasicSourceList& porousSources = heatPorousSources[i];

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

@@ -0,0 +1,28 @@
+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

@@ -0,0 +1,43 @@
+    // Initialise solid field pointer lists
+    PtrList<basicSolidThermo> 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,
+            basicSolidThermo::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

@@ -0,0 +1,24 @@
+    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

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

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

@@ -0,0 +1,20 @@
+    fvMesh& mesh = porousSolidRegions[i];
+
+    basicSolidThermo& 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.K();
+    //tmp<volSymmTensorField> tkappa = thermo.directionalK()*betav;
+
+    const volScalarField& kappa = tkappa();
+    //const volSymmTensorField& K = tK();
+
+    volScalarField& T = thermo.T();
+
+    IObasicSourceList& sources = solidHeatSources[i];

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

@@ -0,0 +1,18 @@
+    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

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

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

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

applications/solvers/incompressible/pimpleFoam/Make/options

@@ -3,7 +3,10 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
 
 EXE_LIBS = \
     -lincompressibleTransportModels \
@@ -11,4 +14,6 @@ EXE_LIBS = \
     -lincompressibleRASModels \
     -lincompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lfieldSources \
+    -lsampling

applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/Make/options

@@ -3,7 +3,9 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \
     -lincompressibleTransportModels \
@@ -11,4 +13,6 @@ EXE_LIBS = \
     -lincompressibleRASModels \
     -lincompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lfieldSources \
+    -lsampling

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/Make/options

@@ -6,7 +6,10 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
 
 EXE_LIBS = \
     -ldynamicFvMesh \
@@ -16,4 +19,6 @@ EXE_LIBS = \
     -lincompressibleTransportModels \
     -lincompressibleRASModels \
     -lincompressibleLESModels \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lfieldSources \
+    -lsampling

applications/solvers/incompressible/potentialFreeSurfaceFoam/Make/options

@@ -3,7 +3,10 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
 
 EXE_LIBS = \
     -lincompressibleTransportModels \
@@ -11,4 +14,6 @@ EXE_LIBS = \
     -lincompressibleRASModels \
     -lincompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lfieldSources \
+    -lsampling

applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/Make/options

@@ -5,10 +5,15 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
 
 EXE_LIBS = \
     -lincompressibleRASModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lfieldSources \
+    -lsampling

applications/solvers/incompressible/simpleFoam/Make/options

@@ -4,11 +4,16 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
 
 EXE_LIBS = \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lfieldSources \
+    -lsampling

applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/Make/options

@@ -5,10 +5,15 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
 
 EXE_LIBS = \
     -lincompressibleRASModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lfieldSources \
+    -lsampling

applications/solvers/lagrangian/LTSReactingParcelFoam/Make/options

@@ -20,7 +20,10 @@ EXE_INC = \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-    -I$(LIB_SRC)/combustionModels/lnInclude
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
 
 EXE_LIBS = \
     -lfiniteVolume \
@@ -44,4 +47,6 @@ EXE_LIBS = \
     -lODE \
     -lregionModels \
     -lsurfaceFilmModels \
-    -lcombustionModels
+    -lcombustionModels \
+    -lfieldSources \
+    -lsampling

applications/solvers/lagrangian/coalChemistryFoam/Make/options

@@ -21,7 +21,10 @@ EXE_INC = \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude \
-    -I$(FOAM_SOLVERS)/combustion/reactingFoam
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
 
 
 EXE_LIBS = \
@@ -47,4 +50,7 @@ EXE_LIBS = \
     -lregionModels \
     -lsurfaceFilmModels \
     -lODE \
-    -lcombustionModels
+    -lcombustionModels \
+    -lfieldSources \
+    -lsampling
+

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/options

@@ -21,6 +21,8 @@ EXE_INC = \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam
 
 
@@ -46,4 +48,7 @@ EXE_LIBS = \
     -lODE \
     -lregionModels \
     -lsurfaceFilmModels \
-    -lcombustionModels
+    -lcombustionModels \
+    -lfieldSources \
+    -lsampling
+

applications/test/DLList/Test-DLList.C

@@ -65,30 +65,36 @@ int main(int argc, char *argv[])
     }
 
     myList.append(500.3);
+    myList.append(200.3);
     myList.append(100.3);
 
 
     Info<< nl << "And again using STL const_iterator: " << nl << endl;
 
-    const DLList<scalar>& const_myList = myList;
 
-    forAllConstIter(DLList<scalar>, const_myList, iter)
+    forAllConstIter(DLList<scalar>, myList, iter)
     {
         Info<< "element:" << *iter << endl;
     }
 
+    Info<< nl << "Testing swapUp and swapDown: " << endl;
+
+    Info<< nl << "swapUp" << endl;
+
     myList.swapUp(myList.DLListBase::first());
     myList.swapUp(myList.DLListBase::last());
 
-    forAllConstIter(DLList<scalar>, const_myList, iter)
+    forAllIter(DLList<scalar>, myList, iter)
     {
         Info<< "element:" << *iter << endl;
     }
 
+    Info<< nl << "swapDown" << endl;
+
     myList.swapDown(myList.DLListBase::first());
     myList.swapDown(myList.DLListBase::last());
 
-    forAllConstIter(DLList<scalar>, const_myList, iter)
+    forAllIter(DLList<scalar>, myList, iter)
     {
         Info<< "element:" << *iter << endl;
     }
@@ -103,8 +109,8 @@ int main(int argc, char *argv[])
     Info<< nl << "source: " << myList << nl
         << nl << "target: " << newList << endl;
 
-
     Info<< nl << "Done." << endl;
+
     return 0;
 }
 

applications/test/List/Test-List.C

@@ -37,6 +37,8 @@ Description
 #include "vector.H"
 #include "ListOps.H"
 
+#include<list>
+
 using namespace Foam;
 
 
@@ -117,6 +119,35 @@ int main(int argc, char *argv[])
         << "-wordList: " << wLst << nl
         << "-stringList: " << sLst << endl;
 
+
+    Info<< nl
+        << "Test List Iterator Constuctor" << endl;
+
+    List<vector> initialList(IStringStream("((0 1 2) (3 4 5) (6 7 8))")());
+
+    Info<< "    Initial List: " << initialList << endl;
+
+    List<vector> iteratorList(initialList.begin(), initialList.end());
+
+    Info<< "    Foam::List constructed from Foam::List: "
+        << iteratorList << endl;
+
+    std::list<vector> stlList(initialList.begin(), initialList.end());
+
+    Info<< "    std::list constructed from Foam::List: ";
+
+    std::list<vector>::iterator it;
+    for (it=stlList.begin(); it != stlList.end(); it++)
+    {
+        Info<< *it << " ";
+    }
+    Info<< endl;
+
+    List<vector> iteratorSTLList(stlList.begin(), stlList.end());
+
+    Info<< "    Foam::List constructed from std::list: "
+        << iteratorList << endl;
+
     return 0;
 }
 

applications/utilities/mesh/manipulation/objToVTK/objToVTK.C

@@ -273,13 +273,17 @@ int main(int argc, char *argv[])
         }
     }
 
+    if (!pointNormals.empty())
+    {
         outFile << nl << "NORMALS pointNormals float\n";
+
         forAll(pointNormals, i)
         {
             const vector& n = pointNormals[i];
 
             outFile << n.x() << ' ' << n.y() << ' ' << n.z() << nl;
         }
+    }
 
     Info<< "End\n" << endl;
 

src/Allwmake

@@ -73,4 +73,6 @@ fvAgglomerationMethods/Allwmake $*
 wmake $makeType fvMotionSolver
 wmake $makeType engine
 
+wmake $makeType fieldSources
+
 # ----------------------------------------------------------------- end-of-file

src/OpenFOAM/containers/LinkedLists/linkTypes/DLListBase/DLListBase.C

@@ -94,23 +94,23 @@ bool Foam::DLListBase::swapUp(DLListBase::link* a)
         if (ap == first_)
         {
             first_ = a;
+            ap->prev_ = a;
+        }
+        else
+        {
+            ap->prev_->next_ = a;
         }
 
         if (a == last_)
         {
             last_ = ap;
+            a->next_ = ap;
         }
-
-        if (a->next_)
+        else
         {
             a->next_->prev_ = ap;
         }
 
-        if (ap->prev_)
-        {
-            ap->prev_->next_ = a;
-        }
-
         a->prev_ = ap->prev_;
         ap->prev_ = a;
 
@@ -135,19 +135,19 @@ bool Foam::DLListBase::swapDown(DLListBase::link* a)
         if (a == first_)
         {
             first_ = an;
+            a->prev_ = an;
+        }
+        else
+        {
+            a->prev_->next_ = an;
         }
 
         if (an == last_)
         {
             last_ = a;
+            an->next_ = a;
         }
-
-        if (a->prev_)
-        {
-            a->prev_->next_ = an;
-        }
-
-        if (an->next_)
+        else
         {
             an->next_->prev_ = a;
         }

src/OpenFOAM/containers/Lists/List/List.C

@@ -200,7 +200,7 @@ Foam::List<T>::List(InputIterator first, InputIterator last)
         ++iter
     )
     {
-        this->operator[](s++) = iter();
+        this->operator[](s++) = *iter;
     }
 }
 

src/OpenFOAM/interpolations/interpolation2DTable/interpolation2DTable.C

@@ -0,0 +1,800 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "IFstream.H"
+#include "openFoamTableReader.H"
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::interpolation2DTable<Type>::readTable()
+{
+    fileName fName(fileName_);
+    fName.expand();
+
+    // Read data from file
+    reader_()(fName, *this);
+    //IFstream(fName)() >> *this;
+
+    if (this->empty())
+    {
+        FatalErrorIn
+        (
+            "Foam::interpolation2DTable<Type>::readTable()"
+        )   << "table read from " << fName << " is empty" << nl
+            << exit(FatalError);
+    }
+
+    // Check that the data are okay
+    check();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::interpolation2DTable<Type>::interpolation2DTable()
+:
+    List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >(),
+    boundsHandling_(interpolation2DTable::WARN),
+    fileName_("fileNameIsUndefined"),
+    reader_(NULL)
+{}
+
+
+template<class Type>
+Foam::interpolation2DTable<Type>::interpolation2DTable
+(
+    const List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >& values,
+    const boundsHandling bounds,
+    const fileName& fName
+)
+:
+    List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >(values),
+    boundsHandling_(bounds),
+    fileName_(fName),
+    reader_(NULL)
+{}
+
+
+template<class Type>
+Foam::interpolation2DTable<Type>::interpolation2DTable(const fileName& fName)
+:
+    List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >(),
+    boundsHandling_(interpolation2DTable::WARN),
+    fileName_(fName),
+    reader_(new openFoamTableReader<Type>(dictionary()))
+{
+    readTable();
+}
+
+
+template<class Type>
+Foam::interpolation2DTable<Type>::interpolation2DTable(const dictionary& dict)
+:
+    List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >(),
+    boundsHandling_(wordToBoundsHandling(dict.lookup("outOfBounds"))),
+    fileName_(dict.lookup("fileName")),
+    reader_(tableReader<Type>::New(dict))
+{
+    readTable();
+}
+
+
+template<class Type>
+Foam::interpolation2DTable<Type>::interpolation2DTable
+(
+     const interpolation2DTable& interpTable
+)
+:
+    List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >(interpTable),
+    boundsHandling_(interpTable.boundsHandling_),
+    fileName_(interpTable.fileName_),
+    reader_(interpTable.reader_)    // note: steals reader. Used in write().
+{}
+
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Type Foam::interpolation2DTable<Type>::interpolateValue
+(
+    const List<Tuple2<scalar, Type> >& data,
+    const scalar lookupValue
+) const
+{
+    label n = data.size();
+
+    scalar minLimit = data[0].first();
+    scalar maxLimit = data[n-1].first();
+
+    if (lookupValue < minLimit)
+    {
+        switch (boundsHandling_)
+        {
+            case interpolation2DTable::ERROR:
+            {
+                FatalErrorIn
+                (
+                    "Foam::interpolation2DTable<Type>::interpolateValue("
+                    "List<Tuple2<scalar, Type> > data,"
+                    "const scalar lookupValue)"
+                )   << "value (" << lookupValue << ") underflow" << nl
+                    << exit(FatalError);
+                break;
+            }
+            case interpolation2DTable::WARN:
+            {
+                WarningIn
+                (
+                    "Foam::interpolation2DTable<Type>::interpolateValue("
+                    "List<Tuple2<scalar, Type> > data,"
+                    "const scalar lookupValue)"
+                )   << "value (" << lookupValue << ") underflow" << nl
+                    << "    Continuing with the first entry"
+                    << endl;
+                // fall-through to 'CLAMP'
+            }
+            case interpolation2DTable::CLAMP:
+            {
+                return data[0].second();
+                break;
+            }
+        }
+    }
+    else if (lookupValue >= maxLimit)
+    {
+        switch (boundsHandling_)
+        {
+            case interpolation2DTable::ERROR:
+            {
+                FatalErrorIn
+                (
+                    "Foam::interpolation2DTable<Type>::interpolateValue("
+                    "List<Tuple2<scalar, Type> > data,"
+                    "const scalar lookupValue)"
+                )   << "value (" << lookupValue << ") overflow" << nl
+                    << exit(FatalError);
+                break;
+            }
+            case interpolation2DTable::WARN:
+            {
+                WarningIn
+                (
+                    "Foam::interpolation2DTable<Type>::interpolateValue("
+                    "List<Tuple2<scalar, Type> > data,"
+                    "const scalar lookupValue)"
+                )   << "value (" << lookupValue << ") overflow" << nl
+                    << "    Continuing with the last entry"
+                    << endl;
+                // fall-through to 'CLAMP'
+            }
+            case interpolation2DTable::CLAMP:
+            {
+                return data[n-1].second();
+                break;
+            }
+        }
+    }
+
+    // look for the correct range in X
+    label lo = 0;
+    label hi = 0;
+
+    for (label i = 0; i < n; ++i)
+    {
+        if (lookupValue >= data[i].first())
+        {
+            lo = hi = i;
+        }
+        else
+        {
+            hi = i;
+            break;
+        }
+    }
+
+    if (lo == hi)
+    {
+        return data[lo].second();
+    }
+    else
+    {
+        // normal interpolation
+        return
+        (
+            data[lo].second()
+          + (
+                data[hi].second()
+              - data[lo].second()
+            )
+           *(
+                lookupValue
+              - data[lo].first()
+            )
+           /(
+                data[hi].first()
+              - data[lo].first()
+            )
+        );
+    }
+
+}
+
+
+template<class Type>
+Foam::word Foam::interpolation2DTable<Type>::boundsHandlingToWord
+(
+     const boundsHandling& bound
+) const
+{
+    word enumName("warn");
+
+    switch (bound)
+    {
+        case interpolation2DTable::ERROR:
+        {
+            enumName = "error";
+            break;
+        }
+        case interpolation2DTable::WARN:
+        {
+            enumName = "warn";
+            break;
+        }
+        case interpolation2DTable::CLAMP:
+        {
+            enumName = "clamp";
+            break;
+        }
+    }
+
+    return enumName;
+}
+
+
+template<class Type>
+typename Foam::interpolation2DTable<Type>::boundsHandling
+Foam::interpolation2DTable<Type>::wordToBoundsHandling
+(
+    const word& bound
+) const
+{
+    if (bound == "error")
+    {
+        return interpolation2DTable::ERROR;
+    }
+    else if (bound == "warn")
+    {
+        return interpolation2DTable::WARN;
+    }
+    else if (bound == "clamp")
+    {
+        return interpolation2DTable::CLAMP;
+    }
+    else
+    {
+        WarningIn
+        (
+            "Foam::interpolation2DTable<Type>::wordToBoundsHandling"
+            "("
+            "    const word&"
+            ")"
+        )   << "bad outOfBounds specifier " << bound << " using 'warn'" << endl;
+
+        return interpolation2DTable::WARN;
+    }
+}
+
+
+template<class Type>
+typename Foam::interpolation2DTable<Type>::boundsHandling
+Foam::interpolation2DTable<Type>::outOfBounds
+(
+    const boundsHandling& bound
+)
+{
+    boundsHandling prev = boundsHandling_;
+    boundsHandling_ = bound;
+    return prev;
+}
+
+
+template<class Type>
+void Foam::interpolation2DTable<Type>::check() const
+{
+    label n = this->size();
+    typedef  List<Tuple2<scalar, List<Tuple2<scalar, Type> > > > matrix;
+    scalar prevValue = matrix::operator[](0).first();
+
+    for (label i=1; i<n; ++i)
+    {
+        const scalar currValue = matrix::operator[](i).first();
+
+        // avoid duplicate values (divide-by-zero error)
+        if (currValue <= prevValue)
+        {
+            FatalErrorIn
+            (
+                "Foam::interpolation2DTable<Type>::check() const"
+            )   << "out-of-order value: "
+                << currValue << " at index " << i << nl
+                << exit(FatalError);
+        }
+        prevValue = currValue;
+    }
+}
+
+
+template<class Type>
+void Foam::interpolation2DTable<Type>::write(Ostream& os) const
+{
+    os.writeKeyword("fileName")
+        << fileName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("outOfBounds")
+        << boundsHandlingToWord(boundsHandling_) << token::END_STATEMENT << nl;
+
+    if (reader_.valid())
+    {
+        reader_->write(os);
+    }
+    //*this >> os;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class Type>
+const Foam::List<Foam::Tuple2<Foam::scalar, Type> >&
+Foam::interpolation2DTable<Type>::operator[](const label i) const
+{
+    label ii = i;
+    label n  = this->size();
+
+    if (n <= 1)
+    {
+        ii = 0;
+    }
+    else if (ii < 0)
+    {
+        switch (boundsHandling_)
+        {
+            case interpolation2DTable::ERROR:
+            {
+                FatalErrorIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const label) const"
+                )   << "index (" << ii << ") underflow" << nl
+                    << exit(FatalError);
+                break;
+            }
+            case interpolation2DTable::WARN:
+            {
+                WarningIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const label) const"
+                )   << "index (" << ii << ") underflow" << nl
+                    << "    Continuing with the first entry"
+                    << endl;
+                // fall-through to 'CLAMP'
+            }
+            case interpolation2DTable::CLAMP:
+            {
+                ii = 0;
+                break;
+            }
+        }
+    }
+    else if (ii >= n)
+    {
+        switch (boundsHandling_)
+        {
+            case interpolation2DTable::ERROR:
+            {
+                FatalErrorIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const label) const"
+                )   << "index (" << ii << ") overflow" << nl
+                    << exit(FatalError);
+                break;
+            }
+            case interpolation2DTable::WARN:
+            {
+                WarningIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const label) const"
+                )   << "index (" << ii << ") overflow" << nl
+                    << "    Continuing with the last entry"
+                    << endl;
+                // fall-through to 'CLAMP'
+            }
+            case interpolation2DTable::CLAMP:
+            {
+                ii = n - 1;
+                break;
+            }
+        }
+    }
+
+    return List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >::operator[](ii);
+}
+
+
+template<class Type>
+Type Foam::interpolation2DTable<Type>::operator()
+(
+    const scalar valueX,
+    const scalar valueY
+) const
+{
+    typedef  List<Tuple2<scalar, List<Tuple2<scalar, Type> > > > matrix;
+    label nX = this->size();
+
+    if (nX <= 1)
+    {
+        const List<Tuple2<scalar, Type> >& dataY =
+            matrix::operator[](0).second();
+
+        return interpolateValue(dataY, valueY);
+    }
+
+
+    scalar minLimit = matrix::operator[](0).first();
+    scalar maxLimit = matrix::operator[](nX-1).first();
+    scalar lookupValue = valueX;
+
+    if (lookupValue < minLimit)
+    {
+        switch (boundsHandling_)
+        {
+            case interpolation2DTable::ERROR:
+            {
+                FatalErrorIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const scalar, const scalar) const"
+                )   << "value (" << lookupValue << ") underflow" << nl
+                    << exit(FatalError);
+                break;
+            }
+            case interpolation2DTable::WARN:
+            {
+                WarningIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const scalar, const scalar) const"
+                )   << "value (" << lookupValue << ") underflow" << nl
+                    << "    Continuing with the first entry"
+                    << endl;
+                // fall-through to 'CLAMP'
+            }
+            case interpolation2DTable::CLAMP:
+            {
+                return interpolateValue
+                (
+                    matrix::operator[](0).second(), valueY
+                );
+                break;
+            }
+        }
+    }
+    else if (lookupValue >= maxLimit)
+    {
+        switch (boundsHandling_)
+        {
+            case interpolation2DTable::ERROR:
+            {
+                FatalErrorIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const label, const scalar) const"
+                )   << "value (" << lookupValue << ") overflow" << nl
+                    << exit(FatalError);
+                break;
+            }
+            case interpolation2DTable::WARN:
+            {
+                WarningIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const label, const scalar) const"
+                )   << "value (" << lookupValue << ") overflow" << nl
+                    << "    Continuing with the last entry"
+                    << endl;
+                // fall-through to 'CLAMP'
+            }
+            case interpolation2DTable::CLAMP:
+            {
+                return interpolateValue
+                (
+                    matrix::operator[](nX-1).second(), valueY
+                );
+                break;
+            }
+        }
+    }
+
+    label loX = 0;
+    label hiX = 0;
+
+    // look for the correct range in X
+    for (label i = 0; i < nX; ++i)
+    {
+        if (lookupValue >= matrix::operator[](i).first())
+        {
+            loX = hiX = i;
+        }
+        else
+        {
+            hiX = i;
+            break;
+        }
+    }
+
+    // look for the correct range in y
+    lookupValue = valueY;
+    label loY1 = 0;
+    label hiY1 = 0;
+
+    label nY = matrix::operator[](loX).second().size();
+
+    minLimit = matrix::operator[](loX).second()[0].first();
+    maxLimit = matrix::operator[](loX).second()[nY-1].first();
+
+    if (lookupValue < minLimit)
+    {
+        switch (boundsHandling_)
+        {
+            case interpolation2DTable::ERROR:
+            {
+                FatalErrorIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const scalar, const scalar) const"
+                )   << "value (" << lookupValue << ") underflow" << nl
+                    << exit(FatalError);
+                break;
+            }
+            case interpolation2DTable::WARN:
+            {
+                WarningIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const scalar, const scalar) const"
+                )   << "value (" << lookupValue << ") underflow" << nl
+                    << "    Continuing with the first entry"
+                    << endl;
+                // fall-through to 'CLAMP'
+            }
+            case interpolation2DTable::CLAMP:
+            {
+                hiY1 = 0;
+                hiY1 = 1;
+                break;
+            }
+        }
+    }
+    else if (lookupValue >= maxLimit)
+    {
+        switch (boundsHandling_)
+        {
+            case interpolation2DTable::ERROR:
+            {
+                FatalErrorIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const scalar, const scalar) const"
+                )   << "value (" << lookupValue << ") overflow" << nl
+                    << exit(FatalError);
+                break;
+            }
+            case interpolation2DTable::WARN:
+            {
+                WarningIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const scalar, const scalar) const"
+                )   << "value (" << lookupValue << ") overflow" << nl
+                    << "    Continuing with the last entry"
+                    << endl;
+                // fall-through to 'CLAMP'
+            }
+            case interpolation2DTable::CLAMP:
+            {
+                hiY1 = nY-1;
+                hiY1 = nY;
+                break;
+            }
+        }
+    }
+    else
+    {
+        // Finds the lo and hi of Y on the lowest x
+        for (label i = 0; i < nY; ++i)
+        {
+            if
+            (
+                lookupValue >= matrix::operator[](loX).second()[i].first()
+            )
+            {
+                loY1 = hiY1 = i;
+            }
+            else
+            {
+                hiY1 = i;
+                break;
+            }
+        }
+    }
+
+    label loY2 = 0;
+    label hiY2 = 0;
+
+    nY = matrix::operator[](hiX).second().size();
+
+    minLimit = matrix::operator[](loX).second()[0].first();
+    maxLimit = matrix::operator[](loX).second()[nY-1].first();
+
+    if (lookupValue < minLimit)
+    {
+        switch (boundsHandling_)
+        {
+            case interpolation2DTable::ERROR:
+            {
+                FatalErrorIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const scalar, const scalar) const"
+                )   << "value (" << lookupValue << ") underflow" << nl
+                    << exit(FatalError);
+                break;
+            }
+            case interpolation2DTable::WARN:
+            {
+                WarningIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const scalar, const scalar) const"
+                )   << "value (" << lookupValue << ") underflow" << nl
+                    << "    Continuing with the first entry"
+                    << endl;
+                // fall-through to 'CLAMP'
+            }
+            case interpolation2DTable::CLAMP:
+            {
+                loY2 = 0;
+                loY2 = 1;
+                break;
+            }
+        }
+    }
+    else if (lookupValue >= maxLimit)
+    {
+        switch (boundsHandling_)
+        {
+            case interpolation2DTable::ERROR:
+            {
+                FatalErrorIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const scalar, const scalar) const"
+                )   << "value (" << lookupValue << ") overflow" << nl
+                    << exit(FatalError);
+                break;
+            }
+            case interpolation2DTable::WARN:
+            {
+                WarningIn
+                (
+                    "Foam::interpolation2DTable<Type>::operator[]"
+                    "(const scalar, const scalar) const"
+                )   << "value (" << lookupValue << ") overflow" << nl
+                    << "    Continuing with the last entry"
+                    << endl;
+                // fall-through to 'CLAMP'
+            }
+            case interpolation2DTable::CLAMP:
+            {
+                loY2 = nY-1;
+                loY2 = nY;
+                break;
+            }
+        }
+    }
+    else
+    {
+        // Finds the lo and hi of Y on the high x
+        for (label i = 0; i < nY; ++i)
+        {
+            if
+            (
+                lookupValue >= matrix::operator[](hiX).second()[i].first()
+            )
+            {
+                loY2 = hiY2 = i;
+            }
+            else
+            {
+                hiY2 = i;
+                break;
+            }
+        }
+    }
+
+    if (loX == hiX)
+    {
+        // we are at the end of the table - or there is only a single entry
+        return (interpolateValue(matrix::operator[](hiX).second(), valueY));
+    }
+    else
+    {
+        Type loXData = matrix::operator[](loX).second()[loY1].second();
+        Type hiXData = matrix::operator[](hiX).second()[loY1].second();
+
+        Type hiYData = matrix::operator[](loX).second()[hiY1].second();
+
+        Type refValue = matrix::operator[](loX).second()[loY1].second();
+
+        // normal interpolation on x
+        refValue +=
+            (
+                hiXData
+              - loXData
+            )
+            *(
+                valueX
+              - matrix::operator[](loX).first()
+             )
+            /(
+                matrix::operator[](hiX).first()
+              - matrix::operator[](loX).first()
+            );
+
+        // normal interpolation on y
+        refValue +=
+            (
+                hiYData
+              - loXData
+            )
+            *(
+                valueY
+              - matrix::operator[](loX).second()[loY1].first()
+            )
+            /(
+                matrix::operator[](loX).second()[hiY1].first()
+              - matrix::operator[](loX).second()[loY1].first()
+            );
+
+        return refValue;
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/interpolations/interpolation2DTable/interpolation2DTable.H

@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::interpolation2DTable
+
+Description
+    2D table interpolation. The data must be in ascending order in both
+    dimensions x and y.
+
+SourceFiles
+    interpolation2DTable.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interpolation2DTable_H
+#define interpolation2DTable_H
+
+#include "List.H"
+#include "Tuple2.H"
+#include "tableReader.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class interpolation2DTable Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class interpolation2DTable
+:
+    public List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >
+{
+public:
+
+    // Public data types
+
+        //- Enumeration for handling out-of-bound values
+        enum boundsHandling
+        {
+            ERROR,          /*!< Exit with a FatalError */
+            WARN,           /*!< Issue warning and clamp value (default) */
+            CLAMP           /*!< Clamp value to the start/end value */
+        };
+
+
+private:
+
+    // Private data
+
+        //- Enumeration for handling out-of-bound values
+        boundsHandling boundsHandling_;
+
+        //- File name
+        fileName fileName_;
+
+        //- the actual reader
+        autoPtr<tableReader<Type> > reader_;
+
+
+    // Private Member Functions
+
+        //- Read the table of data from file
+        void readTable();
+
+        //- Return interpolated value in List
+        Type interpolateValue
+        (
+            const List<Tuple2<scalar, Type> >& data,
+            const scalar
+        ) const;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        interpolation2DTable();
+
+        //- Construct from components
+        interpolation2DTable
+        (
+            const List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >& values,
+            const boundsHandling bounds,
+            const fileName& fName
+        );
+
+        //- Construct given the name of the file containing the table of data
+        interpolation2DTable(const fileName& fName);
+
+        //- Construct by reading the fileName and boundsHandling from dictionary
+        interpolation2DTable(const dictionary& dict);
+
+        //- Construct copy
+        interpolation2DTable(const interpolation2DTable& interpTable);
+
+
+    // Member Functions
+
+        //- Return the out-of-bounds handling as a word
+        word boundsHandlingToWord(const boundsHandling& bound) const;
+
+        //- Return the out-of-bounds handling as an enumeration
+        boundsHandling wordToBoundsHandling(const word& bound) const;
+
+        //- Set the out-of-bounds handling from enum, return previous setting
+        boundsHandling outOfBounds(const boundsHandling& bound);
+
+        //- Check that list is monotonically increasing
+        //  Exit with a FatalError if there is a problem
+        void check() const;
+
+        //- Write
+        void write(Ostream& os) const;
+
+
+    // Member Operators
+
+        //- Return an element of constant Tuple2<scalar, Type>
+        const List<Tuple2<scalar, Type> >& operator[](const label) const;
+
+        //- Return an interpolated value
+        Type operator()(const scalar, const scalar) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "interpolation2DTable.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/interpolations/interpolationTable/tableReaders/csv/csvTableReader.C

@@ -162,6 +162,17 @@ void Foam::csvTableReader<Type>::operator()
 }
 
 
+template<class Type>
+void Foam::csvTableReader<Type>::operator()
+(
+    const fileName& fName,
+    List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >& data
+)
+{
+    notImplemented("csvTableReader<Type>::operator()");
+}
+
+
 template<class Type>
 void Foam::csvTableReader<Type>::write(Ostream& os) const
 {

src/OpenFOAM/interpolations/interpolationTable/tableReaders/csv/csvTableReader.H

@@ -103,6 +103,13 @@ public:
         //- Read the table
         virtual void operator()(const fileName&, List<Tuple2<scalar, Type> >&);
 
+        //- Read 2D table
+        virtual void operator()
+        (
+            const fileName&,
+            List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >&
+        );
+
         //- write the remaining parameters
         virtual void write(Ostream& os) const;
 };

src/OpenFOAM/interpolations/interpolationTable/tableReaders/openFoam/openFoamTableReader.C

@@ -56,4 +56,16 @@ void Foam::openFoamTableReader<Type>::operator()
 }
 
 
+template<class Type>
+void Foam::openFoamTableReader<Type>::operator()
+(
+    const fileName& fName,
+    List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >& data
+)
+{
+    // Read data from file
+    IFstream(fName)() >> data;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/interpolations/interpolationTable/tableReaders/openFoam/openFoamTableReader.H

@@ -84,6 +84,13 @@ public:
 
         //- Read the table
         virtual void operator()(const fileName&, List<Tuple2<scalar, Type> > &);
+
+        //- Read 2D table
+        virtual void operator()
+        (
+            const fileName&,
+            List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >&
+        );
 };
 
 

src/OpenFOAM/interpolations/interpolationTable/tableReaders/tableReader.H

@@ -103,6 +103,13 @@ public:
             List<Tuple2<scalar, Type> >&
         ) = 0;
 
+        //- Read the 2D table
+        virtual void operator()
+        (
+            const fileName&,
+            List<Tuple2<scalar, List<Tuple2<scalar, Type> > > >&
+        ) = 0;
+
         //- Write additional information
         virtual void write(Ostream& os) const;
 };

src/fieldSources/Make/files

@@ -0,0 +1,28 @@
+basicSource/basicSource/basicSource.C
+basicSource/basicSource/basicSourceIO.C
+basicSource/basicSource/basicSourceList.C
+basicSource/basicSource/IObasicSourceList.C
+
+basicSource/pressureGradientExplicitSource/pressureGradientExplicitSource.C
+basicSource/pressureGradientExplicitSource/pressureGradientExplicitSourceIO.C
+
+basicSource/explicitSource/explicitSource.C
+basicSource/explicitSetValue/explicitSetValue.C
+
+basicSource/rotorDiskSource/rotorDiskSource.C
+basicSource/rotorDiskSource/bladeModel/bladeModel.C
+basicSource/rotorDiskSource/profileModel/profileModel.C
+basicSource/rotorDiskSource/profileModel/profileModelList.C
+basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C
+basicSource/rotorDiskSource/profileModel/series/seriesProfile.C
+
+basicSource/actuationDiskSource/actuationDiskSource.C
+basicSource/radialActuationDiskSource/radialActuationDiskSource.C
+
+interRegion = basicSource/interRegionHeatTransferModel
+$(interRegion)/interRegionHeatTransferModel/interRegionHeatTransferModel.C
+$(interRegion)/constantHeatTransfer/constantHeatTransfer.C
+$(interRegion)/tabulatedHeatTransfer/tabulatedHeatTransfer.C
+$(interRegion)/variableHeatTransfer/variableHeatTransfer.C
+
+LIB = $(FOAM_LIBBIN)/libfieldSources

src/fieldSources/Make/options

@@ -0,0 +1,17 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels
+
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    -lsampling \
+    -lmeshTools \
+    -lbasicSolidThermo \
+    -lcompressibleTurbulenceModel

src/fieldSources/basicSource/basicSource/basicSource.C

@@ -38,16 +38,17 @@ namespace Foam
     template<> const char* NamedEnum
     <
         basicSource::selectionModeType,
-        4
+        5
         >::names[] =
     {
         "points",
         "cellSet",
         "cellZone",
+        "mapRegion",
         "all"
     };
 
-    const NamedEnum<basicSource::selectionModeType, 4>
+    const NamedEnum<basicSource::selectionModeType, 5>
         basicSource::selectionModeTypeNames_;
 }
 
@@ -73,6 +74,13 @@ void Foam::basicSource::setSelection(const dictionary& dict)
             dict.lookup("cellZone") >> cellSetName_;
             break;
         }
+        case smMapRegion:
+        {
+            dict_.lookup("secondarySourceName") >> secondarySourceName_;
+            dict_.lookup("mapRegionName") >> mapRegionName_;
+            master_  = readBool(dict_.lookup("master"));
+            break;
+        }
         case smAll:
         {
             break;
@@ -151,6 +159,47 @@ void Foam::basicSource::setCellSet()
 
             break;
         }
+        case smMapRegion:
+        {
+            if(active_)
+            {
+                Info<< indent << "- selecting inter region mapping" << endl;
+                const fvMesh& secondaryMesh =
+                    mesh_.time().lookupObject<fvMesh>(mapRegionName_);
+                const boundBox primaryBB = mesh_.bounds();
+                const boundBox secondaryBB = secondaryMesh.bounds();
+                if (secondaryBB.overlaps(primaryBB))
+                {
+
+                    // Dummy patches
+                    wordList cuttingPatches;
+                    HashTable<word> patchMap;
+
+                    secondaryToPrimaryInterpPtr_.reset
+                    (
+                        new meshToMesh
+                        (
+                            secondaryMesh,
+                            mesh_,
+                            patchMap,
+                            cuttingPatches
+                        )
+                    );
+                }
+                else
+                {
+                    FatalErrorIn
+                    (
+                        "Foam::basicSource::setCellSet()"
+                    )   << "regions dont overlap "
+                        << secondaryMesh.name()
+                        << " in region " << mesh_.name()
+                        << nl
+                        << exit(FatalError);
+                }
+            }
+            break;
+        }
         case smAll:
         {
             Info<< indent << "- selecting all cells" << endl;
@@ -169,6 +218,8 @@ void Foam::basicSource::setCellSet()
     }
 
     // Set volume information
+    if(selectionMode_ != smMapRegion)
+    {
         V_ = 0.0;
         forAll(cells_, i)
         {
@@ -180,6 +231,7 @@ void Foam::basicSource::setCellSet()
             << returnReduce(cells_.size(), sumOp<label>())
             << " cell(s) with volume " << V_ << nl << decrIndent << endl;
     }
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -205,6 +257,11 @@ Foam::basicSource::basicSource
     ),
     cellSetName_("none"),
     V_(0.0),
+    secondaryToPrimaryInterpPtr_(),
+    secondarySourceName_("none"),
+    mapRegionName_("none"),
+    master_(false),
+
     fieldNames_(),
     applied_()
 {
@@ -246,6 +303,13 @@ Foam::autoPtr<Foam::basicSource> Foam::basicSource::New
     return autoPtr<basicSource>(cstrIter()(name, modelType, coeffs, mesh));
 }
 
+Foam::basicSource::~basicSource()
+{
+    if (!secondaryToPrimaryInterpPtr_.empty())
+    {
+        secondaryToPrimaryInterpPtr_.clear();
+    }
+}
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 

src/fieldSources/basicSource/basicSource/basicSource.H

@@ -32,6 +32,7 @@ Description
         timeStart       0.0;                    // start time
         duration        1000.0;                 // duration
         selectionMode   cellSet;                // cellSet // points //cellZone
+                                                // mapRegion
 
     On evaluation, source expects to be added to the rhs of the equation
 
@@ -48,6 +49,8 @@ SourceFiles
 #include "fvMatricesFwd.H"
 #include "cellSet.H"
 #include "autoPtr.H"
+#include "meshToMesh.H"
+
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -57,6 +60,7 @@ namespace Foam
 
 class fvMesh;
 
+
 /*---------------------------------------------------------------------------*\
                          Class basicSource Declaration
 \*---------------------------------------------------------------------------*/
@@ -73,11 +77,13 @@ public:
             smPoints,
             smCellSet,
             smCellZone,
+            smMapRegion,
             smAll
         };
 
         //- Word list of selection mode type names
-        static const NamedEnum<selectionModeType, 4> selectionModeTypeNames_;
+        static const NamedEnum<selectionModeType, 5>
+            selectionModeTypeNames_;
 
 
 protected:
@@ -120,6 +126,21 @@ protected:
         //- Sum of cell volumes
         scalar V_;
 
+        // Data for smMapRegion only
+
+            //- Mesh to mesh mapping for map optiom
+            autoPtr<meshToMesh> secondaryToPrimaryInterpPtr_;
+
+            //- Name of the source in the secondary mesh
+            word secondarySourceName_;
+
+            //- Name of the region to map
+            word mapRegionName_;
+
+            //- Master or slave region
+            bool master_;
+
+
         //- Field names to apply source to - populated by derived models
         wordList fieldNames_;
 
@@ -227,8 +248,8 @@ public:
 
 
     //- Destructor
-    virtual ~basicSource()
-    {}
+    virtual ~basicSource();
+
 
 
     // Member Functions
@@ -266,6 +287,16 @@ public:
             //- Return const access to the total cell volume
             inline scalar V() const;
 
+            //- Return const access to the secondarySourceName
+            inline const word secondarySourceName() const;
+
+            //- Return const access to the mapToMap Ptr
+            inline const autoPtr<meshToMesh>
+                secondaryToPrimaryInterpPtr() const;
+
+            //- Return const referenc to the mapRegion
+            inline const word mapRegionName() const;
+
             //- Return const access to the cell set
             inline const labelList& cells() const;
 

src/fieldSources/basicSource/basicSource/basicSourceI.H

@@ -118,4 +118,21 @@ inline Foam::scalar& Foam::basicSource::duration()
 }
 
 
+inline const Foam::word Foam::basicSource::secondarySourceName() const
+{
+    return secondarySourceName_;
+}
+
+
+inline const Foam::autoPtr<Foam::meshToMesh> Foam::basicSource::
+secondaryToPrimaryInterpPtr() const
+{
+    return secondaryToPrimaryInterpPtr_;
+}
+
+inline const Foam::word Foam::basicSource::mapRegionName() const
+{
+    return mapRegionName_;
+}
+
 // ************************************************************************* //

src/fieldSources/basicSource/basicSource/basicSourceIO.C

@@ -72,6 +72,10 @@ void Foam::basicSource::writeData(Ostream& os) const
         {
             break;
         }
+        case smMapRegion:
+        {
+            break;
+        }
         default:
         {
             FatalErrorIn("basicSource::writeData(Ostream&) const")

src/fieldSources/basicSource/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C

@@ -0,0 +1,141 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "constantHeatTransfer.H"
+#include "fvm.H"
+#include "IObasicSourceList.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvcVolumeIntegrate.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(constantHeatTransfer, 0);
+    addToRunTimeSelectionTable
+    (
+        basicSource,
+        constantHeatTransfer,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::constantHeatTransfer::constantHeatTransfer
+(
+    const word& name,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    interRegionHeatTransferModel(name, modelType, dict, mesh),
+    htCoeffs_(),
+    area_()
+{
+    if (master_)
+    {
+        htCoeffs_.reset
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "htCoeffs",
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                mesh_
+            )
+        );
+
+        area_.reset
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "area",
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                mesh_
+            )
+        );
+
+        htc_.internalField() = htCoeffs_()*area_()/mesh_.V();
+        htc_.correctBoundaryConditions();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::constantHeatTransfer::~constantHeatTransfer()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::tmp<Foam::volScalarField>Foam::constantHeatTransfer::
+calculateHtc()
+{
+    return htc_;
+}
+
+
+void Foam::constantHeatTransfer::writeData(Ostream& os) const
+{
+    os  << indent << token::BEGIN_BLOCK << incrIndent << nl;
+    interRegionHeatTransferModel::writeData(os);
+
+    os << indent << "constantHeatTransfer";
+
+    dict_.write(os);
+
+    os << decrIndent << indent << token::END_BLOCK << endl;
+}
+
+
+bool Foam::constantHeatTransfer::read(const dictionary& dict)
+{
+    if (basicSource::read(dict))
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/fieldSources/basicSource/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H

@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::constantHeatTransfer
+
+Description
+    Constant heat transfer model. htCoeffs [W/m2/K] and area [m2] must be
+    provided.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constantHeatTransfer_H
+#define constantHeatTransfer_H
+
+#include "interRegionHeatTransferModel.H"
+#include "autoPtr.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                        Class constantHeatTransfer Declaration
+\*---------------------------------------------------------------------------*/
+
+class constantHeatTransfer
+:
+    public interRegionHeatTransferModel
+{
+private:
+
+    // Private data
+
+        //- Heat transfer coefficient
+        autoPtr<volScalarField> htCoeffs_;
+
+        //- Area of heat exchange
+        autoPtr<volScalarField> area_;
+
+
+public:
+
+    //- Runtime type information
+        TypeName("constantHeatTransfer");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        constantHeatTransfer
+        (
+            const word& name,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    //- Destructor
+    virtual ~constantHeatTransfer();
+
+
+    // Public Functions
+
+        //- Calculate the heat transfer coefficient
+        virtual const tmp<volScalarField> calculateHtc();
+
+
+        // I-O
+
+            //- Write data
+            virtual void writeData(Ostream&) const;
+
+            //- Read dictionary
+            virtual bool read(const dictionary& dict) ;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#endif
+
+// ************************************************************************* //

src/fieldSources/basicSource/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C

@@ -0,0 +1,270 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "interRegionHeatTransferModel.H"
+#include "basicThermo.H"
+#include "fvm.H"
+#include "IObasicSourceList.H"
+#include "zeroGradientFvPatchFields.H"
+#include "fvcVolumeIntegrate.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(interRegionHeatTransferModel, 0);
+};
+
+
+// * * * * * * * * * * * *  Private member functions * * * * * * * * * * * //
+
+void Foam::interRegionHeatTransferModel::check()
+{
+    const fvMesh& secondaryMesh =
+            mesh_.time().lookupObject<fvMesh>(mapRegionName_);
+
+    const basicSourceList& IObsl =
+        secondaryMesh.lookupObject<basicSourceList>
+        (
+            "sourcesProperties"
+        );
+
+    const PtrList<basicSource>& bsl = IObsl;
+
+    bool secSourceFound(false);
+
+    forAll(bsl, i)
+    {
+        if (bsl[i].name() == secondarySourceName_)
+        {
+            secIrht_ = &const_cast<interRegionHeatTransferModel&>
+            (
+                refCast<const interRegionHeatTransferModel>(bsl[i])
+            );
+            secSourceFound = true;
+            break;
+        }
+    }
+
+    if(!secSourceFound)
+    {
+        FatalErrorIn
+        (
+            "constantHeatTransfer::interRegionHeatTransferModel"
+            "("
+            "   const word& name,"
+            "   const word& modelType,"
+            "   const dictionary& dict,"
+            "   const fvMesh& mesh"
+            ")"
+        )   << "Secondary source name not found" << secondarySourceName_
+            << " in region " << secondaryMesh.name()
+            << nl
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::interRegionHeatTransferModel::interRegionHeatTransferModel
+(
+    const word& name,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    basicSource(name, modelType, dict, mesh),
+    secIrht_(),
+    firstIter_(true),
+    htc_
+    (
+        IOobject
+        (
+            "htc",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar
+        (
+            "htc",
+            dimEnergy/dimTime/dimTemperature/dimVolume,
+            0.0
+        ),
+        zeroGradientFvPatchScalarField::typeName
+    )
+{
+    coeffs_.lookup("fieldNames") >> fieldNames_;
+    applied_.setSize(fieldNames_.size(), false);
+}
+
+
+Foam::interRegionHeatTransferModel::~interRegionHeatTransferModel()
+{}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+void Foam::interRegionHeatTransferModel::addSup
+(
+    fvMatrix<scalar>& eEqn,
+    const label fieldI
+)
+{
+    if (firstIter_)
+    {
+        check();
+        firstIter_ = false;
+    }
+
+    const volScalarField& h = eEqn.psi();
+
+    tmp<volScalarField> tTmapped
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Tmapped" + mesh_.name(),
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("T", dimTemperature, 0.0)
+        )
+    );
+
+    volScalarField& Tmapped = tTmapped();
+
+    const fvMesh& secondaryMesh =
+        mesh_.time().lookupObject<fvMesh>(mapRegionName_);
+
+    const volScalarField& Tsecondary =
+        secondaryMesh.lookupObject<volScalarField>("T");
+
+    secondaryToPrimaryInterpPtr_->interpolateInternalField
+    (
+        Tmapped,
+        Tsecondary,
+        meshToMesh::MAP,
+        eqOp<scalar>()
+    );
+
+    if (!master_)
+    {
+        secondaryToPrimaryInterpPtr_->interpolateInternalField
+        (
+            htc_,
+            secIrht_->calculateHtc(),
+            meshToMesh::CELL_VOLUME_WEIGHT,
+            eqOp<scalar>()
+        );
+    }
+
+    if (debug)
+    {
+        Info<< " Volumetric integral of htc : "
+            << fvc::domainIntegrate(htc_).value()
+            << endl;
+    }
+
+    //SAF: temporarily output
+    if (mesh_.time().outputTime())
+    {
+        Tmapped.write();
+        htc_.write();
+    }
+
+    if (h.dimensions() == dimEnergy/dimMass)
+    {
+        const basicThermo& primaryThermo =
+            mesh_.lookupObject<basicThermo>("thermophysicalProperties");
+
+        eEqn += htc_*Tmapped - fvm::Sp(htc_/primaryThermo.Cp(), h);
+
+        if (debug)
+        {
+            Info<< " Energy exchange from region " << secondaryMesh.name()
+                << " To " << mesh_.name() << " : "
+                <<  fvc::domainIntegrate
+                    (
+                        htc_*(h/primaryThermo.Cp() - Tmapped)
+                    ).value()
+                << endl;
+        }
+    }
+    else if(h.dimensions() == dimTemperature)
+    {
+        eEqn += htc_*Tmapped - fvm::Sp(htc_, h);
+
+        if (debug)
+        {
+            Info<< " Enegy exchange from region " << secondaryMesh.name()
+                << " To " << mesh_.name() << " : "
+                <<  fvc::domainIntegrate(htc_*(h - Tmapped)).value()
+                << endl;
+        }
+    }
+}
+
+
+void Foam::interRegionHeatTransferModel::writeData(Ostream& os) const
+{
+    os.writeKeyword("name") << this->name() << token::END_STATEMENT << nl;
+    os.writeKeyword("mapRegionName") << mapRegionName_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("secondarySourceName") << secondarySourceName_
+        << token::END_STATEMENT << nl;
+
+    os.writeKeyword("master") << master_ << token::END_STATEMENT << nl;
+
+    if (dict_.found("note"))
+    {
+        os.writeKeyword("note") << string(dict_.lookup("note"))
+            << token::END_STATEMENT << nl;
+    }
+
+    dict_.write(os);
+}
+
+
+bool Foam::interRegionHeatTransferModel::read(const dictionary& dict)
+{
+    if (basicSource::read(dict))
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+// ************************************************************************* //

src/fieldSources/basicSource/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H

@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::interRegionHeatTransferModel
+
+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
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interRegionHeatTransferModel_H
+#define interRegionHeatTransferModel_H
+
+#include "basicSource.H"
+#include "volFields.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class interRegionHeatTransferModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class interRegionHeatTransferModel
+:
+    public  basicSource
+{
+
+private:
+
+    // Private data
+
+        //- Pointer to secondary interRegionHeatTransferModel
+        interRegionHeatTransferModel* secIrht_;
+
+        //- First iteration
+        bool firstIter_;
+
+
+    // Private members
+
+        //- Check coupled interRegionHeatTransferModel
+        void check();
+
+
+protected:
+
+    // Protected data
+
+        //- Heat transfer coefficient [W/m2/k] by area/volume [1/m]
+        // registered on the master mesh
+        volScalarField htc_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("interRegionHeatTransferModel");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        interRegionHeatTransferModel
+        (
+            const word& name,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    //- Destructor
+    virtual ~interRegionHeatTransferModel();
+
+
+    // Member Functions
+
+
+        // Access
+
+            //- Return the heat transfer coefficient
+            const volScalarField& htc() const
+            {
+                return  htc_;
+            }
+
+            //-Source term to fvMatrix<scalar>
+            virtual void addSup(fvMatrix<scalar>& eEqn, const label fieldI);
+
+            //- Calculate heat transfer coefficient
+            virtual const tmp<volScalarField> calculateHtc() = 0;
+
+
+        // I-O
+
+            //- Write data
+            virtual void writeData(Ostream&) const;
+
+            //- Read dictionary
+            virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#endif
+
+// ************************************************************************* //