TUT: Add and update overset tutorials

COMP: Compilation of src with Allwmake
ENH: Update overset solvers
2021-08-23 19:31:22 -07:00 · 2021-08-23 19:30:43 -07:00 · 2021-08-23 19:29:54 -07:00 · 2021-08-23 12:06:18 -07:00 · 2021-08-23 12:05:36 -07:00 · 2021-08-23 12:05:35 -07:00
1080 changed files with 37704 additions and 7173 deletions
--- a/META-INFO/api-info
+++ b/META-INFO/api-info
@ -1,2 +1,2 @@
-api=2106
+api=2107
 patch=0
--- a/applications/solvers/basic/potentialFoam/overPotentialFoam/overPotentialFoam.C
+++ b/applications/solvers/basic/potentialFoam/overPotentialFoam/overPotentialFoam.C
@ -149,7 +149,7 @@ int main(int argc, char *argv[])

    mesh.update();

-    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
+    //surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));

    // Since solver contains no time loop it would never execute
    // function objects so do it ourselves
--- a/applications/solvers/combustion/fireFoam/Make/options
+++ b/applications/solvers/combustion/fireFoam/Make/options
@ -1,5 +1,6 @@
 EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I${LIB_SRC}/meshTools/lnInclude \
    -I${LIB_SRC}/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
@ -14,7 +15,6 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/solidChemistryModel/lnInclude \
    -I$(LIB_SRC)/combustionModels/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
@ -23,7 +23,9 @@ EXE_INC = \
    -I$(LIB_SRC)/regionModels/pyrolysisModels/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/ODE/lnInclude
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -50,4 +52,7 @@ EXE_LIBS = \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
-    -lODE
+    -lODE \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/compressible/rhoPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/Make/options
@ -8,7 +8,8 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/regionFaModels\lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -23,4 +24,5 @@ EXE_LIBS = \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
    -lsampling \
-    -latmosphericModels
+    -latmosphericModels \
+    -lregionFaModels
--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createControls.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createControls.H
@ -1,4 +0,0 @@
-bool ddtCorr
-(
-    pimple.dict().getOrDefault("ddtCorr", true)
-);
--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createFields.H
@ -69,6 +69,8 @@ mesh.setFluxRequired(p.name());

 // Mask field for zeroing out contributions on hole cells
 #include "createCellMask.H"
+// Create bool field with interpolated cells
+#include "createInterpolatedCells.H"

 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/overRhoPimpleDyMFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/overRhoPimpleDyMFoam.C
@ -43,7 +43,6 @@ Description
 #include "dynamicFvMesh.H"
 #include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
-#include "bound.H"
 #include "pimpleControl.H"
 #include "pressureControl.H"
 #include "CorrectPhi.H"
@ -89,10 +88,8 @@ int main(int argc, char *argv[])

    while (runTime.run())
    {
-        #include "readControls.H"
        #include "readDyMControls.H"

-
        // Store divrhoU from the previous mesh so that it can be mapped
        // and used in correctPhi to ensure the corrected phi has the
        // same divergence
@ -128,7 +125,6 @@ int main(int argc, char *argv[])
        {
            if (pimple.firstIter() || moveMeshOuterCorrectors)
            {
-
                // Do any mesh changes
                mesh.update();

@ -137,52 +133,22 @@ int main(int argc, char *argv[])
                    MRF.update();

                    #include "setCellMask.H"
-
-                    const surfaceScalarField faceMaskOld
-                    (
-                        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-                    );
-
-                    // Zero Uf on old faceMask (H-I)
-                    rhoUf() *= faceMaskOld;
-
-                    surfaceVectorField rhoUfint(fvc::interpolate(rho*U));
-
-                    // Update Uf and phi on new C-I faces
-                    rhoUf() += (1-faceMaskOld)*rhoUfint;
-
-                    // Update Uf boundary
-                    forAll(rhoUf().boundaryField(), patchI)
-                    {
-                        rhoUf().boundaryFieldRef()[patchI] =
-                            rhoUfint.boundaryField()[patchI];
-                    }
-
-                    // Calculate absolute flux from the mapped surface velocity
-                    phi = mesh.Sf() & rhoUf();
+                    #include "setInterpolatedCells.H"
+                    #include "correctRhoPhiFaceMask.H"

                    if (correctPhi)
                    {
+                        // Corrects flux on separated regions
                        #include "correctPhi.H"
                    }

-                    // Zero phi on current H-I
-                    const surfaceScalarField faceMask
-                    (
-                        localMin<scalar>(mesh).interpolate(cellMask)
-                    );
-
-                    phi *= faceMask;
-                    U   *= cellMask;
-
                     // Make the fluxes relative to the mesh-motion
                    fvc::makeRelative(phi, rho, U);

-                }
-
-                if (checkMeshCourantNo)
-                {
-                    #include "meshCourantNo.H"
+                    if (checkMeshCourantNo)
+                    {
+                        #include "meshCourantNo.H"
+                    }
                }
            }

--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H
@ -25,17 +25,6 @@ surfaceScalarField phiHbyA
    fvc::interpolate(rho)*fvc::flux(HbyA)
 );

-if (ddtCorr)
-{
-    surfaceScalarField faceMaskOld
-    (
-        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-    );
-
-    phiHbyA +=
-        faceMaskOld*MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf));
-}
-
 fvc::makeRelative(phiHbyA, rho, U);
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);

@ -134,8 +123,4 @@ if (thermo.dpdt())
    }
 }

-surfaceScalarField faceMask
-(
-    localMin<scalar>(mesh).interpolate(cellMask)
-);
 phi *= faceMask;
--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/readControls.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/readControls.H
@ -1,9 +0,0 @@
-#include "readTimeControls.H"
-
-correctPhi = pimple.dict().getOrDefault("correctPhi", false);
-
-checkMeshCourantNo =
-    pimple.dict().getOrDefault("checkMeshCourantNo", false);
-
-
-ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);
--- a/applications/solvers/compressible/rhoSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/Make/options
@ -7,6 +7,7 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/regionFaModels\lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -18,4 +19,5 @@ EXE_LIBS = \
    -lspecie \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -latmosphericModels
+    -latmosphericModels \
+    -lregionFaModels
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
@ -6,7 +6,8 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/regionFaModels\lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -19,4 +20,5 @@ EXE_LIBS = \
    -lspecie \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -latmosphericModels
+    -latmosphericModels \
+    -lregionFaModels
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/createFields.H
@ -124,3 +124,6 @@ dimensionedScalar initialMass("initialMass", fvc::domainIntegrate(rho));

 // Mask field for zeroing out contributions on hole cells
 #include "createCellMask.H"
+
+// Create bool field with interpolated cells
+#include "createInterpolatedCells.H"
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/overBuoyantPimpleDyMFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/overBuoyantPimpleDyMFoam.C
@ -86,9 +86,6 @@ int main(int argc, char *argv[])

    while (runTime.run())
    {
-        #include "readTimeControls.H"
-
-        #include "readControls.H"
        #include "readDyMControls.H"

        #include "compressibleCourantNo.H"
@ -128,45 +125,14 @@ int main(int argc, char *argv[])
                    MRF.update();

                    #include "setCellMask.H"
-
-                    const surfaceScalarField faceMaskOld
-                    (
-                        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-                    );
-
-                    // Zero Uf on old faceMask (H-I)
-                    rhoUf() *= faceMaskOld;
-
-                    //fvc::correctRhoUf(rhoUfint, rho, U, phi);
-                    surfaceVectorField rhoUfint(fvc::interpolate(rho*U));
-
-                    // Update Uf and phi on new C-I faces
-                    rhoUf() += (1-faceMaskOld)*rhoUfint;
-
-                    // Update Uf boundary
-                    forAll(rhoUf().boundaryField(), patchI)
-                    {
-                        rhoUf().boundaryFieldRef()[patchI] =
-                            rhoUfint.boundaryField()[patchI];
-                    }
-
-                    // Calculate absolute flux from the mapped surface velocity
-                    phi = mesh.Sf() & rhoUf();
+                    #include "setInterpolatedCells.H"
+                    #include "correctRhoPhiFaceMask.H"

                    if (correctPhi)
                    {
                        #include "correctPhi.H"
                    }

-                    // Zero phi on current H-I
-                    const surfaceScalarField faceMask
-                    (
-                        localMin<scalar>(mesh).interpolate(cellMask)
-                    );
-
-                    phi *= faceMask;
-                    U   *= cellMask;
-
                     // Make the fluxes relative to the mesh-motion
                    fvc::makeRelative(phi, rho, U);
                }
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/pEqn.H
@ -21,16 +21,6 @@ surfaceScalarField phiHbyA
    fvc::flux(rho*HbyA) + phig
 );

-if (ddtCorr)
-{
-    surfaceScalarField faceMaskOld
-    (
-        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-    );
-
-    phiHbyA +=
-        faceMaskOld*MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi));
-}

 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);

@ -122,8 +112,4 @@ if (thermo.dpdt())
    }
 }

-surfaceScalarField faceMask
-(
-    localMin<scalar>(mesh).interpolate(cellMask)
-);
 phi *= faceMask;
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/readControls.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/readControls.H
@ -1,9 +0,0 @@
-#include "readTimeControls.H"
-
-correctPhi = pimple.dict().getOrDefault("correctPhi", false);
-
-checkMeshCourantNo =
-    pimple.dict().getOrDefault("checkMeshCourantNo", false);
-
-
-ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options
@ -6,7 +6,8 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/regionFaModels\lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -19,4 +20,5 @@ EXE_LIBS = \
    -lradiationModels \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -latmosphericModels
+    -latmosphericModels \
+    -lregionFaModels
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
@ -18,7 +18,9 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-    -I$(LIB_SRC)/regionModels/regionModel/lnInclude
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionFaModels\lnInclude
+

 EXE_LIBS = \
    -lfiniteVolume \
@ -36,4 +38,5 @@ EXE_LIBS = \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
    -lradiationModels \
-    -lregionModels
+    -lregionModels \
+    -lregionFaModels
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -54,6 +54,7 @@ Description
 #include "loopControl.H"
 #include "pressureControl.H"

+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

 int main(int argc, char *argv[])
@ -80,6 +81,8 @@ int main(int argc, char *argv[])
    #include "solidRegionDiffusionNo.H"
    #include "setInitialMultiRegionDeltaT.H"

+    #include "createCoupledRegions.H"
+
    while (runTime.run())
    {
        #include "readTimeControls.H"
@ -109,8 +112,6 @@ int main(int argc, char *argv[])

            forAll(fluidRegions, i)
            {
-                Info<< "\nSolving for fluid region "
-                    << fluidRegions[i].name() << endl;
                #include "setRegionFluidFields.H"
                #include "readFluidMultiRegionPIMPLEControls.H"
                #include "solveFluid.H"
@ -118,13 +119,41 @@ int main(int argc, char *argv[])

            forAll(solidRegions, i)
            {
-                Info<< "\nSolving for solid region "
-                    << solidRegions[i].name() << endl;
                #include "setRegionSolidFields.H"
                #include "readSolidMultiRegionPIMPLEControls.H"
                #include "solveSolid.H"
            }

+            if (coupled)
+            {
+                Info<< "\nSolving energy coupled regions " << endl;
+                fvMatrixAssemblyPtr->solve();
+                #include "correctThermos.H"
+
+                forAll(fluidRegions, i)
+                {
+                    #include "setRegionFluidFields.H"
+                    #include "readFluidMultiRegionPIMPLEControls.H"
+                    if (!frozenFlow)
+                    {
+                        Info<< "\nSolving for fluid region "
+                            << fluidRegions[i].name() << endl;
+                        // --- PISO loop
+                        for (int corr=0; corr<nCorr; corr++)
+                        {
+                            #include "pEqn.H"
+                        }
+                        turbulence.correct();
+                    }
+
+                    rho = thermo.rho();
+                    Info<< "Min/max T:" << min(thermo.T()).value() << ' '
+                        << max(thermo.T()).value() << endl;
+                }
+
+                fvMatrixAssemblyPtr->clear();
+            }
+
            // Additional loops for energy solution only
            if (!oCorr && nOuterCorr > 1)
            {
@ -152,6 +181,21 @@ int main(int argc, char *argv[])
                        #include "readSolidMultiRegionPIMPLEControls.H"
                        #include "solveSolid.H"
                    }
+
+                    if (coupled)
+                    {
+                        Info<< "\nSolving energy coupled regions " << endl;
+                        fvMatrixAssemblyPtr->solve();
+                        #include "correctThermos.H"
+
+                        forAll(fluidRegions, i)
+                        {
+                            #include "setRegionFluidFields.H"
+                            rho = thermo.rho();
+                        }
+
+                        fvMatrixAssemblyPtr->clear();
+                    }
                }
            }
        }
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options
@ -3,6 +3,7 @@ EXE_INC = \
    -I./fluid \
    -I./solid \
    -I../solid \
+    -I./../include \
    -I$(LIB_SRC)/finiteVolume/cfdTools \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
@ -35,3 +36,4 @@ EXE_LIBS = \
    -lregionModels \
    -lsampling \
    -lregionFaModels
+
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -66,6 +66,7 @@ int main(int argc, char *argv[])
    #include "createTime.H"
    #include "createMeshes.H"
    #include "createFields.H"
+    #include "createCoupledRegions.H"
    #include "initContinuityErrs.H"

    while (runTime.loop())
@ -83,13 +84,32 @@ int main(int argc, char *argv[])

        forAll(solidRegions, i)
        {
-            Info<< "\nSolving for solid region "
-                << solidRegions[i].name() << endl;
            #include "setRegionSolidFields.H"
            #include "readSolidMultiRegionSIMPLEControls.H"
            #include "solveSolid.H"
        }

+
+        if (coupled)
+        {
+            Info<< "\nSolving energy coupled regions" << endl;
+            fvMatrixAssemblyPtr->solve();
+            #include "correctThermos.H"
+
+            forAll(fluidRegions, i)
+            {
+                #include "setRegionFluidFields.H"
+                #include "readSolidMultiRegionSIMPLEControls.H"
+                if (!frozenFlow)
+                {
+                    #include "pEqn.H"
+                    turb.correct();
+                }
+            }
+
+            fvMatrixAssemblyPtr->clear();
+        }
+
        // Additional loops for energy solution only
        {
            loopControl looping(runTime, "SIMPLE", "energyCoupling");
@ -116,6 +136,15 @@ int main(int argc, char *argv[])
                    #include "readSolidMultiRegionSIMPLEControls.H"
                    #include "solveSolid.H"
                }
+
+                if (coupled)
+                {
+                    Info<< "\nSolving energy coupled regions.. " << endl;
+                    fvMatrixAssemblyPtr->solve();
+                    #include "correctThermos.H"
+
+                    fvMatrixAssemblyPtr->clear();
+                }
            }
        }

--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/EEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/EEqn.H
@ -20,13 +20,20 @@

    fvOptions.constrain(EEqn);

-    EEqn.solve();
+    if (coupled)
+    {
+        fvMatrixAssemblyPtr->addFvMatrix(EEqn);
+    }
+    else
+    {
+        EEqn.solve();

-    fvOptions.correct(he);
+        fvOptions.correct(he);

-    thermo.correct();
-    rad.correct();
+        thermo.correct();
+        rad.correct();

-    Info<< "Min/max T:" << min(thermo.T()).value() << ' '
-        << max(thermo.T()).value() << endl;
+        Info<< "Min/max T:" << min(thermo.T()).value() << ' '
+            << max(thermo.T()).value() << endl;
+    }
 }
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H
@ -2,7 +2,7 @@

    MRF.correctBoundaryVelocity(U);

-    tmp<fvVectorMatrix> tUEqn
+    UEqn =
    (
        fvm::div(phi, U)
      + MRF.DDt(rho, U)
@ -10,7 +10,6 @@
     ==
        fvOptions(rho, U)
    );
-    fvVectorMatrix& UEqn = tUEqn.ref();

    UEqn.relax();

--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
@ -21,6 +21,8 @@ PtrList<dimensionedScalar> rhoMin(fluidRegions.size());
 PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
 PtrList<fv::options> fluidFvOptions(fluidRegions.size());

+PtrList<fvVectorMatrix> UEqFluid(fluidRegions.size());
+
 const uniformDimensionedVectorField& g = meshObjects::gravity::New(runTime);

 // Populate fluid field pointer lists
@ -222,5 +224,11 @@ forAll(fluidRegions, i)
        new fv::options(fluidRegions[i])
    );

+    UEqFluid.set
+    (
+        i,
+        new fvVectorMatrix(UFluid[i], dimForce)
+    );
+
    turbulence[i].validate();
 }
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@ -2,7 +2,7 @@
    volScalarField rAU("rAU", 1.0/UEqn.A());
    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
-    tUEqn.clear();
+    //tUEqn.clear();

    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());

--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
@ -22,6 +22,8 @@
    IOMRFZoneList& MRF = MRFfluid[i];
    fv::options& fvOptions = fluidFvOptions[i];

+    fvVectorMatrix& UEqn = UEqFluid[i];
+
    const dimensionedScalar initialMass
    (
        "initialMass",
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H
@ -12,8 +12,10 @@

            #include "UEqn.H"
            #include "EEqn.H"
-            #include "pEqn.H"
-
-            turb.correct();
+            if (!coupled)
+            {
+                #include "pEqn.H"
+                turb.correct();
+            }
        }
    }
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H
@ -16,13 +16,21 @@

        fvOptions.constrain(hEqn);

-        hEqn.solve();
+        if (coupled)
+        {
+            fvMatrixAssemblyPtr->addFvMatrix(hEqn);
+        }
+        else
+        {
+            Info<< "\nSolving for solid region "
+                << solidRegions[i].name() << endl;

-        fvOptions.correct(h);
+            hEqn.solve();
+            fvOptions.correct(h);
+            thermo.correct();
+
+            Info<< "Min/max T:" << min(thermo.T()).value() << ' '
+                << max(thermo.T()).value() << endl;
+        }
    }
-
-    thermo.correct();
-
-    Info<< "Min/max T:" << min(thermo.T()).value() << ' '
-        << max(thermo.T()).value() << endl;
 }
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/Make/options
@ -7,6 +7,7 @@ EXE_INC = \
    -I${phaseSystem}/twoPhaseCompressibleTurbulenceModels/lnInclude \
    -I${phaseSystem}/multiphaseSystem/lnInclude \
    -I./fluid \
+    -I./solid \
    -I../solid \
    -I../fluid \
    -I../include \
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/solid/solveSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/solid/solveSolid.H
@ -0,0 +1,39 @@
+if (finalIter)
+{
+    mesh.data::add("finalIteration", true);
+}
+
+{
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; ++nonOrth)
+    {
+        fvScalarMatrix hEqn
+        (
+            fvm::ddt(betav*rho, h)
+          - (
+               thermo.isotropic()
+             ? fvm::laplacian(betav*thermo.alpha(), h, "laplacian(alpha,h)")
+             : fvm::laplacian(betav*taniAlpha(), h, "laplacian(alpha,h)")
+            )
+          ==
+            fvOptions(rho, h)
+        );
+
+        hEqn.relax();
+
+        fvOptions.constrain(hEqn);
+
+        hEqn.solve(mesh.solver(h.select(finalIter)));
+
+        fvOptions.correct(h);
+    }
+
+    thermo.correct();
+
+    Info<< "Min/max T:" << min(thermo.T()).value() << ' '
+        << max(thermo.T()).value() << endl;
+}
+
+if (finalIter)
+{
+    mesh.data::remove("finalIteration");
+}
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H
@ -27,13 +27,19 @@

    fvOptions.constrain(EEqn);

-    EEqn.solve(mesh.solver(he.select(finalIter)));
+    if (coupled)
+    {
+        fvMatrixAssemblyPtr->addFvMatrix(EEqn);
+    }
+    else
+    {
+        EEqn.solve(mesh.solver(he.select(finalIter)));
+        fvOptions.correct(he);

-    fvOptions.correct(he);
+        thermo.correct();
+        rad.correct();

-    thermo.correct();
-    rad.correct();
-
-    Info<< "Min/max T:" << min(thermo.T()).value() << ' '
-        << max(thermo.T()).value() << endl;
+        Info<< "Min/max T:" << min(thermo.T()).value() << ' '
+            << max(thermo.T()).value() << endl;
+    }
 }
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
@ -2,7 +2,7 @@

    MRF.correctBoundaryVelocity(U);

-    tmp<fvVectorMatrix> tUEqn
+    UEqn =
    (
        fvm::ddt(rho, U) + fvm::div(phi, U)
      + MRF.DDt(rho, U)
@ -10,7 +10,6 @@
     ==
        fvOptions(rho, U)
    );
-    fvVectorMatrix& UEqn = tUEqn.ref();

    UEqn.relax();

--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@ -16,6 +16,8 @@ PtrList<multivariateSurfaceInterpolationScheme<scalar>::fieldTable>
    fieldsFluid(fluidRegions.size());
 PtrList<volScalarField> QdotFluid(fluidRegions.size());

+PtrList<fvVectorMatrix> UEqFluid(fluidRegions.size());
+
 List<scalar> initialMassFluid(fluidRegions.size());
 List<bool> frozenFlowFluid(fluidRegions.size(), false);

@ -293,6 +295,12 @@ forAll(fluidRegions, i)
        new fv::options(fluidRegions[i])
    );

+    UEqFluid.set
+    (
+        i,
+        new fvVectorMatrix(UFluid[i], dimForce)
+    );
+
    turbulenceFluid[i].validate();

    pRefCellFluid[i] = -1;
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@ -49,6 +49,8 @@
    IOMRFZoneList& MRF = MRFfluid[i];
    fv::options& fvOptions = fluidFvOptions[i];

+    fvVectorMatrix& UEqn = UEqFluid[i];
+
    const dimensionedScalar initialMass
    (
        "initialMass",
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
@ -18,15 +18,20 @@ else
    #include "YEqn.H"
    #include "EEqn.H"

-    // --- PISO loop
-    for (int corr=0; corr<nCorr; corr++)
+    if (!coupled)
    {
-        #include "pEqn.H"
+        Info<< "\nSolving for fluid region " << fluidRegions[i].name() << endl;
+
+        // --- PISO loop
+        for (int corr=0; corr<nCorr; corr++)
+        {
+            #include "pEqn.H"
+        }
+
+        turbulence.correct();
+
+        rho = thermo.rho();
    }
-
-    turbulence.correct();
-
-    rho = thermo.rho();
 }

 if (finalIter)
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/include/correctThermos.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/include/correctThermos.H
@ -0,0 +1,20 @@
+forAll(fluidRegions, i)
+{
+    rhoThermo& thermo = thermoFluid[i];
+    radiation::radiationModel& rad = radiation[i];
+    fv::options& fvOptions = fluidFvOptions[i];
+    volScalarField& he = thermo.he();
+    fvOptions.correct(he);
+    thermo.correct();
+    rad.correct();
+}
+
+forAll(solidRegions, i)
+{
+    solidThermo& thermo = thermos[i];
+    fv::options& fvOptions = solidHeatSources[i];
+    volScalarField& h = thermo.he();
+    fvOptions.correct(h);
+    thermo.correct();
+}
+
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/include/createCoupledRegions.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/include/createCoupledRegions.H
@ -0,0 +1,51 @@
+fvSolution solutionDict(runTime);
+
+bool coupled(solutionDict.getOrDefault("coupledEnergyField", false));
+
+autoPtr<fvMatrix<scalar>> fvMatrixAssemblyPtr;
+
+forAll(fluidRegions, i)
+{
+    const rhoThermo& thermo = refCast<const rhoThermo>(thermoFluid[i]);
+    const auto& bpsi = thermo.T().boundaryField();
+
+    forAll (bpsi, patchI)
+    {
+        if (bpsi[patchI].useImplicit())
+        {
+            coupled = true;
+        }
+    }
+}
+
+forAll(solidRegions, i)
+{
+    solidThermo& thermo = thermos[i];
+    const auto& bpsi = thermo.T().boundaryField();
+
+    forAll (bpsi, patchI)
+    {
+        if (bpsi[patchI].useImplicit())
+        {
+            coupled = true;
+        }
+    }
+}
+
+forAll(fluidRegions, i)
+{
+    const rhoThermo& thermo = refCast<const rhoThermo>(thermoFluid[i]);
+    if (coupled)
+    {
+        Info << "Create fvMatrixAssembly." << endl;
+        fvMatrixAssemblyPtr.reset
+        (
+            new fvMatrix<scalar>
+            (
+                thermo.he(),
+                dimEnergy/dimTime
+            )
+        );
+        break;
+    }
+}
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
@ -1,39 +1,45 @@
-if (finalIter)
 {
-    mesh.data::add("finalIteration", true);
-}
+    fvScalarMatrix hEqn
+    (
+        fvm::ddt(betav*rho, h)
+        - (
+            thermo.isotropic()
+            ? fvm::laplacian(betav*thermo.alpha(), h, "laplacian(alpha,h)")
+            : fvm::laplacian(betav*taniAlpha(), h, "laplacian(alpha,h)")
+        )
+        ==
+        fvOptions(rho, h)
+    );

-{
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; ++nonOrth)
+    hEqn.relax();
+
+    fvOptions.constrain(hEqn);
+
+    if (coupled)
    {
-        fvScalarMatrix hEqn
-        (
-            fvm::ddt(betav*rho, h)
-          - (
-               thermo.isotropic()
-             ? fvm::laplacian(betav*thermo.alpha(), h, "laplacian(alpha,h)")
-             : fvm::laplacian(betav*taniAlpha(), h, "laplacian(alpha,h)")
-            )
-          ==
-            fvOptions(rho, h)
-        );
+        fvMatrixAssemblyPtr->addFvMatrix(hEqn);
+    }
+    else
+    {
+        Info<< "\nSolving for solid region "<< solidRegions[i].name() << endl;

-        hEqn.relax();
-
-        fvOptions.constrain(hEqn);
+        if (finalIter)
+        {
+            mesh.data::add("finalIteration", true);
+        }

        hEqn.solve(mesh.solver(h.select(finalIter)));

        fvOptions.correct(h);
+
+        thermo.correct();
+
+        Info<< "Min/max T:" << min(thermo.T()).value() << ' '
+            << max(thermo.T()).value() << endl;
+
+        if (finalIter)
+        {
+            mesh.data::remove("finalIteration");
+        }
    }
-
-    thermo.correct();
-
-    Info<< "Min/max T:" << min(thermo.T()).value() << ' '
-        << max(thermo.T()).value() << endl;
-}
-
-if (finalIter)
-{
-    mesh.data::remove("finalIteration");
 }
--- a/applications/solvers/incompressible/pimpleFoam/Make/options
+++ b/applications/solvers/incompressible/pimpleFoam/Make/options
@ -7,7 +7,8 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/regionFaModels\lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -20,4 +21,5 @@ EXE_LIBS = \
    -ldynamicMesh \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
-    -latmosphericModels
+    -latmosphericModels \
+    -lregionFaModels
--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/CourantNo.H
+++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/CourantNo.H
@ -37,7 +37,10 @@ scalar meanCoNum = 0.0;

 if (mesh.nInternalFaces())
 {
-    surfaceScalarField phiMask(localMin<scalar>(mesh).interpolate(cellMask));
+    surfaceScalarField phiMask
+    (
+        localMin<scalar>(mesh).interpolate(cellMask + interpolatedCells)
+    );

    scalarField sumPhi(fvc::surfaceSum(mag(phiMask*phi))().internalField());

--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/UEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/UEqn.H
@ -1,5 +1,4 @@
 // Solve the Momentum equation
-
 MRF.correctBoundaryVelocity(U);

 tmp<fvVectorMatrix> tUEqn
--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createControls.H
+++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createControls.H
@ -1,26 +0,0 @@
-#include "createTimeControls.H"
-
-bool correctPhi
-(
-    pimple.dict().getOrDefault("correctPhi", false)
-);
-
-bool checkMeshCourantNo
-(
-    pimple.dict().getOrDefault("checkMeshCourantNo", false)
-);
-
-bool massFluxInterpolation
-(
-    pimple.dict().getOrDefault("massFluxInterpolation", false)
-);
-
-bool adjustFringe
-(
-    pimple.dict().getOrDefault("oversetAdjustPhi", false)
-);
-
-bool ddtCorr
-(
-    pimple.dict().getOrDefault("ddtCorr", true)
-);
--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H
+++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H
@ -1,273 +0,0 @@
-// Interpolation used
-interpolationCellPoint<vector> UInterpolator(HbyA);
-
-// Determine faces on outside of interpolated cells
-bitSet isOwnerInterpolatedFace(mesh.nInternalFaces());
-bitSet isNeiInterpolatedFace(mesh.nInternalFaces());
-
-// Determine donor cells
-labelListList donorCell(mesh.nInternalFaces());
-
-scalarListList weightCellCells(mesh.nInternalFaces());
-
-// Interpolated HbyA faces
-vectorField UIntFaces(mesh.nInternalFaces(), Zero);
-
-// Determine receptor neighbour cells
-labelList receptorNeigCell(mesh.nInternalFaces(), -1);
-
-{
-    const cellCellStencilObject& overlap = Stencil::New(mesh);
-    const labelList& cellTypes = overlap.cellTypes();
-    const labelIOList& zoneID = overlap.zoneID();
-
-    label nZones = gMax(zoneID)+1;
-    PtrList<fvMeshSubset> meshParts(nZones);
-    labelList nCellsPerZone(nZones, Zero);
-
-    // A mesh subset for each zone
-    forAll(meshParts, zonei)
-    {
-        meshParts.set
-        (
-            zonei,
-            // Select cells where the zoneID == zonei
-            new fvMeshSubset(mesh, zonei, zoneID)
-        );
-    }
-
-    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
-    {
-        label ownType = cellTypes[mesh.faceOwner()[faceI]];
-        label neiType = cellTypes[mesh.faceNeighbour()[faceI]];
-        if
-        (
-            ownType == cellCellStencil::INTERPOLATED
-            && neiType == cellCellStencil::CALCULATED
-        )
-        {
-            isOwnerInterpolatedFace.set(faceI);
-
-            const vector& fc = mesh.faceCentres()[faceI];
-
-            for (label zoneI = 0; zoneI < nZones; zoneI++)
-            {
-                if (zoneI != zoneID[mesh.faceOwner()[faceI]])
-                {
-                    const fvMesh& partMesh = meshParts[zoneI].subMesh();
-                    const labelList& cellMap = meshParts[zoneI].cellMap();
-                    label cellI = partMesh.findCell(fc);
-
-                    if (cellI != -1)
-                    {
-                         // Determine weights
-                        labelList stencil(partMesh.cellCells()[cellI]);
-
-                        stencil.append(cellI);
-
-                        label st = stencil.size();
-
-                        donorCell[faceI].setSize(st);
-
-                        weightCellCells[faceI].setSize(st);
-
-                        scalarField weights(st);
-
-                        forAll(stencil, i)
-                        {
-                            scalar d = mag
-                            (
-                                partMesh.cellCentres()[stencil[i]]
-                              - fc
-                            );
-                            weights[i] = 1.0/d;
-                            donorCell[faceI][i] = cellMap[stencil[i]];
-                        }
-                        weights /= sum(weights);
-
-                        weightCellCells[faceI] = weights;
-
-                        forAll(stencil, i)
-                        {
-                            UIntFaces[faceI] +=
-                                weightCellCells[faceI][i]
-                               *UInterpolator.interpolate
-                                (
-                                    fc,
-                                    donorCell[faceI][i]
-                                );
-                        }
-
-                        break;
-                    }
-                }
-            }
-
-            receptorNeigCell[faceI] = mesh.faceNeighbour()[faceI];
-        }
-        else if
-        (
-            ownType == cellCellStencil::CALCULATED
-            && neiType == cellCellStencil::INTERPOLATED
-        )
-        {
-            isNeiInterpolatedFace.set(faceI);
-
-            const vector& fc = mesh.faceCentres()[faceI];
-            for (label zoneI = 0; zoneI < nZones; zoneI++)
-            {
-                if (zoneI != zoneID[mesh.faceNeighbour()[faceI]])
-                {
-                    const fvMesh& partMesh = meshParts[zoneI].subMesh();
-                    const labelList& cellMap = meshParts[zoneI].cellMap();
-                    label cellI = partMesh.findCell(fc);
-
-                    if (cellI != -1)
-                    {
-                        // Determine weights
-                        labelList stencil(partMesh.cellCells()[cellI]);
-
-                        stencil.append(cellI);
-
-                        label st = stencil.size();
-
-                        donorCell[faceI].setSize(st);
-
-                        weightCellCells[faceI].setSize(st);
-
-                        scalarField weights(st);
-
-                        forAll(stencil, i)
-                        {
-                            scalar d = mag
-                            (
-                                partMesh.cellCentres()[stencil[i]]
-                              - fc
-                            );
-                            weights[i] = 1.0/d;
-                            donorCell[faceI][i] = cellMap[stencil[i]];
-                        }
-                        weights /= sum(weights);
-
-                        weightCellCells[faceI] = weights;
-
-                        forAll(stencil, i)
-                        {
-                            UIntFaces[faceI] +=
-                                weightCellCells[faceI][i]
-                               *UInterpolator.interpolate
-                                (
-                                    fc,
-                                    donorCell[faceI][i]
-                                );
-                        }
-
-                        break;
-                    }
-                }
-            }
-
-            receptorNeigCell[faceI] = mesh.faceOwner()[faceI];
-        }
-    }
-}
-
-// contravariant U
-vectorField U1Contrav(mesh.nInternalFaces(), Zero);
-
-surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf());
-
-forAll(isNeiInterpolatedFace, faceI)
-{
-    label cellId = -1;
-    if (isNeiInterpolatedFace.test(faceI))
-    {
-        cellId = mesh.faceNeighbour()[faceI];
-    }
-    else if (isOwnerInterpolatedFace.test(faceI))
-    {
-        cellId = mesh.faceOwner()[faceI];
-    }
-
-    if (cellId != -1)
-    {
-        const vector& n = faceNormals[faceI];
-        vector n1(Zero);
-
-        // 2-D cases
-        if (mesh.nSolutionD() == 2)
-        {
-            for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
-            {
-                if (mesh.geometricD()[cmpt] == -1)
-                {
-                    switch (cmpt)
-                    {
-                        case vector::X:
-                        {
-                            n1 = vector(0, n.z(), -n.y());
-                            break;
-                        }
-
-                        case vector::Y:
-                        {
-                            n1 = vector(n.z(), 0, -n.x());
-                            break;
-                        }
-
-                        case vector::Z:
-                        {
-                            n1 = vector(n.y(), -n.x(), 0);
-                            break;
-                        }
-                    }
-                }
-            }
-        }
-        else if (mesh.nSolutionD() == 3)
-        {
-            //Determine which is the primary direction
-            if (mag(n.x()) > mag(n.y()) && mag(n.x()) > mag(n.z()))
-            {
-                n1 = vector(n.y(), -n.x(), 0);
-            }
-            else if (mag(n.y()) > mag(n.z()))
-            {
-                n1 = vector(0, n.z(), -n.y());
-            }
-            else
-            {
-                n1 = vector(-n.z(), 0, n.x());
-            }
-        }
-        n1.normalise();
-
-        const vector n2 = normalised(n ^ n1);
-
-        tensor rot =
-            tensor
-            (
-               n.x() ,n.y(), n.z(),
-               n1.x() ,n1.y(), n1.z(),
-               n2.x() ,n2.y(), n2.z()
-            );
-
-//         tensor rot =
-//             tensor
-//             (
-//                n  & x ,n  & y, n  & z,
-//                n1 & x ,n1 & y, n1 & z,
-//                n2 & x ,n2 & y, n2 & z
-//             );
-
-        U1Contrav[faceI].x() =
-           2*transform(rot, UIntFaces[faceI]).x()
-           - transform(rot, HbyA[receptorNeigCell[faceI]]).x();
-
-        U1Contrav[faceI].y() = transform(rot, HbyA[cellId]).y();
-
-        U1Contrav[faceI].z() = transform(rot, HbyA[cellId]).z();
-
-        HbyA[cellId] = transform(inv(rot), U1Contrav[faceI]);
-    }
-}
--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/overPimpleDyMFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/overPimpleDyMFoam.C
@ -46,11 +46,7 @@ Description
 #include "fvOptions.H"

 #include "cellCellStencilObject.H"
-#include "zeroGradientFvPatchFields.H"
 #include "localMin.H"
-#include "interpolationCellPoint.H"
-#include "transform.H"
-#include "fvMeshSubset.H"
 #include "oversetAdjustPhi.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -68,10 +64,9 @@ int main(int argc, char *argv[])
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createDynamicFvMesh.H"
+    #include "createDyMControls.H"
    #include "initContinuityErrs.H"

-    pimpleControl pimple(mesh);
-
    #include "createFields.H"
    #include "createUf.H"
    #include "createMRF.H"
@ -88,7 +83,7 @@ int main(int argc, char *argv[])

    while (runTime.run())
    {
-        #include "readControls.H"
+        #include "readDyMControls.H"
        #include "CourantNo.H"

        #include "setDeltaT.H"
@ -97,45 +92,26 @@ int main(int argc, char *argv[])

        Info<< "Time = " << runTime.timeName() << nl << endl;

-        bool changed = mesh.update();
+        mesh.update();

-        if (changed)
+        if (mesh.changing())
        {
            #include "setCellMask.H"
            #include "setInterpolatedCells.H"
+            #include "correctPhiFaceMask.H"

-            surfaceScalarField faceMaskOld
-            (
-                localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-            );
+            if (correctPhi)
+            {
+                // Corrects flux on separated regions
+                #include "correctPhi.H"
+            }

-            // Zero Uf on old faceMask (H-I)
-            Uf *= faceMaskOld;
-            // Update Uf and phi on new C-I faces
-            Uf += (1-faceMaskOld)*fvc::interpolate(U);
-            phi = mesh.Sf() & Uf;
+            fvc::makeRelative(phi, U);

-            // Zero phi on current H-I
-            surfaceScalarField faceMask
-            (
-                localMin<scalar>(mesh).interpolate(cellMask)
-            );
-            phi *= faceMask;
-        }
-
-
-        if (mesh.changing() && correctPhi)
-        {
-            // Calculate absolute flux from the mapped surface velocity
-            #include "correctPhi.H"
-        }
-
-        // Make the flux relative to the mesh motion
-        fvc::makeRelative(phi, U);
-
-        if (mesh.changing() && checkMeshCourantNo)
-        {
-            #include "meshCourantNo.H"
+            if (checkMeshCourantNo)
+            {
+                #include "meshCourantNo.H"
+            }
        }

        // --- Pressure-velocity PIMPLE corrector loop
--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H
@ -1,29 +1,11 @@
-// Option 1: interpolate rAU, do not block out rAU on blocked cells
 volScalarField rAU("rAU", 1.0/UEqn.A());
-mesh.interpolate(rAU);

-// Option 2: do not interpolate rAU but block out rAU
-//surfaceScalarField rAUf("rAUf", fvc::interpolate(blockedCells*rAU));
-
-
-// Option 3: do not interpolate rAU but zero out rAUf on faces on holes
-// But what about:
-//
-//   H
-// H I C C C C
-//   H
-//
 surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 volVectorField H("H", UEqn.H());

 volVectorField HbyA("HbyA", U);
 HbyA = constrainHbyA(rAU*H, U, p);

-if (massFluxInterpolation)
-{
-    #include "interpolatedFaces.H"
-}
-
 if (runTime.outputTime())
 {
    H.write();
@ -38,18 +20,14 @@ if (pimple.nCorrPISO() <= 1)

 phiHbyA = fvc::flux(HbyA);

-if (ddtCorr)
+if (runTime.outputTime())
 {
-    surfaceScalarField faceMaskOld
-    (
-        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-    );
-    phiHbyA += rAUf*faceMaskOld*fvc::ddtCorr(U, Uf);
+    volScalarField divPhiHbyA("divPhiHbyA", fvc::div(phiHbyA));
+    divPhiHbyA.write();
 }

 MRF.makeRelative(phiHbyA);

-// WIP
 if (p.needReference())
 {
    fvc::makeRelative(phiHbyA, U);
@ -57,7 +35,8 @@ if (p.needReference())
    fvc::makeAbsolute(phiHbyA, U);
 }

-
+// WIP: To adjust phi on fringe faces to help mass
+// conservation
 if (adjustFringe)
 {
    fvc::makeRelative(phiHbyA, U);
@ -79,27 +58,29 @@ while (pimple.correctNonOrthogonal())
    if (pimple.finalNonOrthogonalIter())
    {
        phi = phiHbyA - pEqn.flux();
-        // option 2:
-        // rAUf*fvc::snGrad(p)*mesh.magSf();
+        pEqn.relax();
+        U =
+            cellMask*
+            (
+                HbyA -
+                rAU*fvc::reconstruct((pEqn.flux())/rAUf)
+            );
+        U.correctBoundaryConditions();
+        fvOptions.correct(U);
    }
 }

-
+// Excludes error in interpolated/hole cells
 #include "continuityErrs.H"

-// Explicitly relax pressure for momentum corrector
-p.relax();
-volVectorField gradP(fvc::grad(p));
+if (runTime.outputTime())
+{
+    volVectorField Ucorrected("Ucorrected", U);
+    Ucorrected.write();

-// Option 2: zero out velocity on blocked out cells
-//U = HbyA - rAU*cellMask*gradP;
-// Option 3: zero out velocity on blocked out cells
-// This is needed for the scalar Eq (k,epsilon, etc)
-// which can use U as source term
-U = cellMask*(HbyA - rAU*gradP);
-U.correctBoundaryConditions();
-
-fvOptions.correct(U);
+    volScalarField divPhiCon("divPhiCon", fvc::div(phi));
+    divPhiCon.write();
+}

 {
    Uf = fvc::interpolate(U);
@ -109,9 +90,4 @@ fvOptions.correct(U);

 // Make the fluxes relative to the mesh motion
 fvc::makeRelative(phi, U);
-
-surfaceScalarField faceMask
-(
-    localMin<scalar>(mesh).interpolate(cellMask)
-);
 phi *= faceMask;
--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/readControls.H
+++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/readControls.H
@ -1,10 +0,0 @@
-#include "readTimeControls.H"
-
-correctPhi = pimple.dict().getOrDefault("correctPhi", false);
-
-checkMeshCourantNo = pimple.dict().getOrDefault("checkMeshCourantNo", false);
-
-massFluxInterpolation =
-    pimple.dict().getOrDefault("massFluxInterpolation", false);
-
-ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);
--- a/applications/solvers/incompressible/pimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pEqn.H
@ -6,10 +6,7 @@ if (pimple.ddtCorr())
 {
    phiHbyA += MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi, Uf));
 }
-else
-{
-    phiHbyA += MRF.zeroFilter(fvc::interpolate(rAU));
-}
+

 MRF.makeRelative(phiHbyA);

--- a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createOversetFields.H
+++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createOversetFields.H
@ -24,7 +24,3 @@ bool adjustFringe
 (
    simple.dict().getOrDefault("oversetAdjustPhi", false)
 );
-bool massFluxInterpolation
-(
-    simple.dict().getOrDefault("massFluxInterpolation", false)
-);
--- a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/pEqn.H
@ -1,5 +1,5 @@
 {
-    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
+    //surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));

    volScalarField rAU(1.0/UEqn.A());
    surfaceScalarField rAUf("rAUf", faceMask*fvc::interpolate(rAU));
@ -7,12 +7,6 @@
    volVectorField HbyA("HbyA", U);
    HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);

-    //mesh.interpolate(HbyA);
-    if (massFluxInterpolation)
-    {
-        #include "interpolatedFaces.H"
-    }
-
    tUEqn.clear();

    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options
@ -3,6 +3,7 @@ EXE_INC = \
    -I../.. \
    -I../../DPMTurbulenceModels \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
@ -19,7 +20,9 @@ EXE_INC = \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -39,4 +42,7 @@ EXE_LIBS = \
    -ldynamicMesh \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
-    -lsampling
+    -lsampling \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/options
@ -2,6 +2,7 @@ EXE_INC = \
    -I.. \
    -I../DPMTurbulenceModels \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
@ -18,7 +19,9 @@ EXE_INC = \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -38,4 +41,7 @@ EXE_LIBS = \
    -ldynamicMesh \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
-    -lsampling
+    -lsampling \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/DPMIncompressibleTurbulenceModel.C
+++ b/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/DPMIncompressibleTurbulenceModel.C
@ -149,6 +149,17 @@ Foam::DPMIncompressibleTurbulenceModel<TransportModel>::devReff() const
 }


+template<class TransportModel>
+Foam::tmp<Foam::volSymmTensorField>
+Foam::DPMIncompressibleTurbulenceModel<TransportModel>::devReff
+(
+    const volVectorField& U
+) const
+{
+    return devRhoReff(U);
+}
+
+
 template<class TransportModel>
 Foam::tmp<Foam::fvVectorMatrix>
 Foam::DPMIncompressibleTurbulenceModel<TransportModel>::divDevReff
@ -170,6 +181,19 @@ Foam::DPMIncompressibleTurbulenceModel<TransportModel>::devRhoReff() const
 }


+template<class TransportModel>
+Foam::tmp<Foam::volSymmTensorField>
+Foam::DPMIncompressibleTurbulenceModel<TransportModel>::devRhoReff
+(
+    const volVectorField& U
+) const
+{
+    NotImplemented;
+
+    return nullptr;
+}
+
+
 template<class TransportModel>
 Foam::tmp<Foam::fvVectorMatrix>
 Foam::DPMIncompressibleTurbulenceModel<TransportModel>::divDevRhoReff
--- a/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/DPMIncompressibleTurbulenceModel.H
+++ b/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/DPMIncompressibleTurbulenceModel.H
@ -116,12 +116,24 @@ public:
        //- Return the effective stress tensor
        virtual tmp<volSymmTensorField> devReff() const;

+        //- Return the effective stress tensor based on a given velocity field
+        virtual tmp<volSymmTensorField> devReff
+        (
+            const volVectorField& U
+        ) const;
+
        //- Return the source term for the momentum equation
        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;

        //- Return the effective stress tensor
        virtual tmp<volSymmTensorField> devRhoReff() const;

+        //- Return the effective stress tensor based on a given velocity field
+        virtual tmp<volSymmTensorField> devRhoReff
+        (
+            const volVectorField& U
+        ) const;
+
        //- Return the source term for the momentum equation
        virtual tmp<fvVectorMatrix> divDevRhoReff(volVectorField& U) const;
 };
--- a/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options
@ -2,6 +2,7 @@ EXE_INC = \
    -I.. \
    -I../DPMTurbulenceModels \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
@ -17,6 +18,8 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -33,4 +36,7 @@ EXE_LIBS = \
    -lDPMTurbulenceModels \
    -lregionModels \
    -lsurfaceFilmModels \
-    -lsampling
+    -lsampling \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/DPMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/Make/options
@ -1,6 +1,7 @@
 EXE_INC = \
    -I./DPMTurbulenceModels \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
@ -14,7 +15,9 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -31,4 +34,7 @@ EXE_LIBS = \
    -lDPMTurbulenceModels \
    -lregionModels \
    -lsurfaceFilmModels \
-    -lsampling
+    -lsampling \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/coalChemistryFoam/Make/options
+++ b/applications/solvers/lagrangian/coalChemistryFoam/Make/options
@ -1,5 +1,6 @@
 EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I${LIB_SRC}/meshTools/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
@ -11,13 +12,14 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude \
    -I$(LIB_SRC)/ODE/lnInclude \
    -I$(LIB_SRC)/combustionModels/lnInclude \
    -I$(FOAM_SOLVERS)/combustion/reactingFoam \
@ -45,4 +47,7 @@ EXE_LIBS = \
    -lsurfaceFilmModels \
    -lODE \
    -lcombustionModels \
-    -lsampling
+    -lsampling \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Make/options
@ -1,5 +1,6 @@
 EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
@ -13,7 +14,9 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -30,4 +33,7 @@ EXE_LIBS = \
    -lincompressibleTurbulenceModels \
    -lincompressibleTransportModels \
    -lregionModels \
-    -lsurfaceFilmModels
+    -lsurfaceFilmModels \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options
@ -1,6 +1,7 @@
 EXE_INC = \
    -I.. \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
@ -16,7 +17,9 @@ EXE_INC = \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -36,4 +39,7 @@ EXE_LIBS = \
    -lsurfaceFilmModels \
    -ldynamicMesh \
    -ldynamicFvMesh \
-    -ltopoChangerFvMesh
+    -ltopoChangerFvMesh \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/kinematicParcelFoam/Make/files
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/Make/files
@ -0,0 +1,3 @@
+kinematicParcelFoam.C
+
+EXE = $(FOAM_APPBIN)/kinematicParcelFoam
--- a/applications/solvers/lagrangian/kinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/Make/options
@ -0,0 +1,42 @@
+EXE_INC = \
+    -I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude
+
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -ldynamicMesh \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -latmosphericModels \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lsurfaceFilmDerivedFvPatchFields \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -llagrangianTurbulence \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/kinematicParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/UEqn.H
@ -0,0 +1,22 @@
+    MRF.correctBoundaryVelocity(U);
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(U) + fvm::div(phi, U)
+      + MRF.DDt(U)
+      + turbulence->divDevReff(U)
+     ==
+        parcels.SU(U, true)
+      + fvOptions(U)
+    );
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve(UEqn == -fvc::grad(p));
+
+        fvOptions.correct(U);
+    }
--- a/applications/solvers/lagrangian/kinematicParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/createClouds.H
@ -0,0 +1,16 @@
+const word kinematicCloudName
+(
+    args.getOrDefault<word>("cloud", "kinematicCloud")
+);
+
+Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
+
+basicKinematicCloud parcels
+(
+    kinematicCloudName,
+    rhoInf,
+    U,
+    muc,
+    g
+);
+
--- a/applications/solvers/lagrangian/kinematicParcelFoam/createFieldRefs.H
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/createFieldRefs.H
@ -0,0 +1 @@
+regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
--- a/applications/solvers/lagrangian/kinematicParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/createFields.H
@ -0,0 +1,85 @@
+#include "readGravitationalAcceleration.H"
+
+Info<< "Reading field p\n" << endl;
+volScalarField p
+(
+    IOobject
+    (
+        "p",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "\nReading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "createPhi.H"
+
+singlePhaseTransportModel laminarTransport(U, phi);
+
+dimensionedScalar rhoInfValue
+(
+    "rhoInf",
+    dimDensity,
+    laminarTransport
+);
+
+volScalarField rhoInf
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    rhoInfValue
+);
+
+volScalarField muc
+(
+    IOobject
+    (
+        "muc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    rhoInf*laminarTransport.nu()
+);
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<incompressible::turbulenceModel> turbulence
+(
+    incompressible::turbulenceModel::New(U, phi, laminarTransport)
+);
+
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+setRefCell(p, pimple.dict(), pRefCell, pRefValue);
+mesh.setFluxRequired(p.name());
+
+#include "createMRF.H"
+#include "createClouds.H"
+#include "createSurfaceFilmModel.H"
+#include "createFvOptions.H"
+
+
--- a/applications/solvers/lagrangian/kinematicParcelFoam/kinematicParcelFoam.C
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/kinematicParcelFoam.C
@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Application
+    kinematicParcelFoam
+
+Group
+    grpLagrangianSolvers
+
+Description
+    Transient solver for incompressible, turbulent flow with kinematic,
+    particle cloud, and surface film modelling.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "surfaceFilmModel.H"
+#include "basicKinematicCloud.H"
+#include "fvOptions.H"
+#include "pimpleControl.H"
+#include "CorrectPhi.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Transient solver for incompressible, turbulent flow"
+        " with kinematic particle clouds"
+        " and surface film modelling."
+    );
+
+    #define CREATE_MESH createMeshesPostProcess.H
+    #include "postProcess.H"
+
+    #include "addCheckCaseOptions.H"
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "initContinuityErrs.H"
+    #include "createDyMControls.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createRegionControls.H"
+    #include "createUfIfPresent.H"
+
+    turbulence->validate();
+
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readDyMControls.H"
+        #include "CourantNo.H"
+        #include "setMultiRegionDeltaT.H"
+
+        ++runTime;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // Store the particle positions
+        parcels.storeGlobalPositions();
+
+        // Do any mesh changes
+        mesh.update();
+
+        if (solvePrimaryRegion && mesh.changing())
+        {
+            MRF.update();
+
+            if (correctPhi)
+            {
+                // Calculate absolute flux
+                // from the mapped surface velocity
+                phi = mesh.Sf() & Uf();
+
+                #include "../../incompressible/pimpleFoam/correctPhi.H"
+
+                // Make the fluxes relative to the mesh-motion
+                fvc::makeRelative(phi, U);
+            }
+
+            if (checkMeshCourantNo)
+            {
+                #include "meshCourantNo.H"
+            }
+        }
+
+        parcels.evolve();
+        surfaceFilm.evolve();
+
+        if (solvePrimaryRegion)
+        {
+            // --- PIMPLE loop
+            while (pimple.loop())
+            {
+                #include "UEqn.H"
+
+                // --- Pressure corrector loop
+                while (pimple.correct())
+                {
+                    #include "pEqn.H"
+                }
+
+                if (pimple.turbCorr())
+                {
+                    laminarTransport.correct();
+                    turbulence->correct();
+                }
+            }
+        }
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/lagrangian/kinematicParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/pEqn.H
@ -0,0 +1,58 @@
+
+volScalarField rAU(1.0/UEqn.A());
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
+
+if (pimple.ddtCorr())
+{
+    phiHbyA += MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi, Uf));
+}
+
+MRF.makeRelative(phiHbyA);
+
+if (p.needReference())
+{
+    fvc::makeRelative(phiHbyA, U);
+    adjustPhi(phiHbyA, U, p);
+    fvc::makeAbsolute(phiHbyA, U);
+}
+
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, U, phiHbyA, rAU, MRF);
+
+// Non-orthogonal pressure corrector loop
+while (pimple.correctNonOrthogonal())
+{
+    fvScalarMatrix pEqn
+    (
+        fvm::laplacian(rAU, p)
+     ==
+        fvc::div(phiHbyA)
+    );
+
+    pEqn.setReference(pRefCell, pRefValue);
+
+    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+    if (pimple.finalNonOrthogonalIter())
+    {
+        phi = phiHbyA - pEqn.flux();
+    }
+}
+
+#include "continuityErrs.H"
+
+p.relax();
+
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+
+// Correct rhoUf if the mesh is moving
+fvc::correctUf(Uf, U, phi);
+
+// Make the fluxes relative to the mesh motion
+fvc::makeRelative(phi, U);
+
--- a/applications/solvers/lagrangian/reactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/Make/options
@ -1,6 +1,7 @@
 EXE_INC = \
    -I../reactingParcelFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I${LIB_SRC}/meshTools/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
@ -12,13 +13,14 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
    -I$(LIB_SRC)/ODE/lnInclude \
@ -49,4 +51,7 @@ EXE_LIBS = \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
    -lODE \
-    -lcombustionModels
+    -lcombustionModels \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
@ -11,7 +11,11 @@ basicSpecieMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();

 const word inertSpecie(thermo.get<word>("inertSpecie"));
-if (!composition.species().found(inertSpecie))
+if
+(
+    !composition.species().found(inertSpecie)
+  && composition.species().size() > 0
+)
 {
    FatalIOErrorIn(args.executable().c_str(), thermo)
        << "Inert specie " << inertSpecie << " not found in available species "
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/options
@ -1,6 +1,7 @@
 EXE_INC = \
    -I.. \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I${LIB_SRC}/sampling/lnInclude \
    -I${LIB_SRC}/meshTools/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
@ -12,13 +13,14 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
    -I$(LIB_SRC)/ODE/lnInclude \
@ -48,4 +50,7 @@ EXE_LIBS = \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
    -lODE \
-    -lcombustionModels
+    -lcombustionModels \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
@ -1,5 +1,6 @@
 EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I${LIB_SRC}/meshTools/lnInclude \
    -I$(LIB_SRC)/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
@ -12,7 +13,6 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
@ -20,6 +20,8 @@ EXE_INC = \
    -I$(LIB_SRC)/ODE/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude \
    -I$(LIB_SRC)/combustionModels/lnInclude \
    -I$(FOAM_SOLVERS)/combustion/reactingFoam

@ -44,4 +46,7 @@ EXE_LIBS = \
    -lregionModels \
    -lsurfaceFilmModels \
    -lcombustionModels \
-    -lsampling
+    -lsampling \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/simpleCoalParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/simpleCoalParcelFoam/Make/options
@ -1,5 +1,6 @@
 EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I$(LIB_SRC)/fvOptions/lnInclude \
    -I${LIB_SRC}/meshTools/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
@ -12,7 +13,6 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
@ -21,7 +21,9 @@ EXE_INC = \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
    -I$(LIB_SRC)/combustionModels/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -44,4 +46,7 @@ EXE_LIBS = \
    -lsurfaceFilmModels \
    -lcombustionModels \
    -lsampling \
-    -lcoalCombustion
+    -lcoalCombustion \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/sprayFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/Make/options
@ -1,6 +1,7 @@
 EXE_INC = \
    -I../reactingParcelFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I${LIB_SRC}/meshTools/lnInclude \
    -I${LIB_SRC}/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
@ -13,7 +14,6 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
@ -21,6 +21,8 @@ EXE_INC = \
    -I$(LIB_SRC)/ODE/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude \
    -I$(LIB_SRC)/combustionModels/lnInclude

 EXE_LIBS = \
@ -45,4 +47,7 @@ EXE_LIBS = \
    -lregionModels \
    -lsurfaceFilmModels \
    -lcombustionModels \
-    -lsampling
+    -lsampling \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options
@ -3,6 +3,7 @@ EXE_INC = \
    -I.. \
    -I../../reactingParcelFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I${LIB_SRC}/meshTools/lnInclude \
    -I${LIB_SRC}/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
@ -15,7 +16,6 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
@ -24,7 +24,9 @@ EXE_INC = \
    -I$(LIB_SRC)/engine/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)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -49,4 +51,7 @@ EXE_LIBS = \
    -lengine \
    -lregionModels \
    -lsurfaceFilmModels \
-    -lcombustionModels
+    -lcombustionModels \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/Make/options
@ -1,6 +1,7 @@
 EXE_INC = \
    -I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam/simpleReactingParcelFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I${LIB_SRC}/meshTools/lnInclude \
    -I$(LIB_SRC)/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
@ -13,7 +14,6 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
@ -21,7 +21,9 @@ 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)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -45,4 +47,7 @@ EXE_LIBS = \
    -lregionModels \
    -lsurfaceFilmModels \
    -lcombustionModels \
-    -lsampling
+    -lsampling \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/Make/options
@ -3,6 +3,7 @@ EXE_INC = \
    -I../../reactingParcelFoam \
    -I../../../compressible/rhoPimpleFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I${LIB_SRC}/meshTools/lnInclude \
    -I${LIB_SRC}/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
@ -15,7 +16,6 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
@ -26,7 +26,9 @@ EXE_INC = \
    -I$(LIB_SRC)/combustionModels/lnInclude \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -52,4 +54,7 @@ EXE_LIBS = \
    -lcombustionModels \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
-    -ldynamicMesh
+    -ldynamicMesh \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/Make/options
@ -1,5 +1,6 @@
 EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
@ -11,7 +12,9 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude \

 EXE_LIBS = \
    -lfiniteVolume \
@ -27,4 +30,7 @@ EXE_LIBS = \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
    -lregionModels \
-    -lsurfaceFilmModels
+    -lsurfaceFilmModels \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
@ -1,6 +1,7 @@
 EXE_INC = \
    -I.. \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
@ -14,7 +15,9 @@ EXE_INC = \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -33,4 +36,7 @@ EXE_LIBS = \
    -lsurfaceFilmModels \
    -ldynamicMesh \
    -ldynamicFvMesh \
-    -ltopoChangerFvMesh
+    -ltopoChangerFvMesh \
+    -lregionFaModels \
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/multiphase/MPPICInterFoam/Make/options
+++ b/applications/solvers/multiphase/MPPICInterFoam/Make/options
@ -2,6 +2,7 @@ EXE_INC = \
    -I../VoF \
    -I$(FOAM_SOLVERS)/multiphase/interFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I$(LIB_SRC)/fvOptions/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/sampling/lnInclude \
@ -20,7 +21,9 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -38,3 +41,4 @@ EXE_LIBS = \
    -lregionModels \
    -lsurfaceFilmModels \
    -lcompressibleTwoPhaseMixtureTurbulenceModels
+
--- a/applications/solvers/multiphase/VoF/alphaEqn.H
+++ b/applications/solvers/multiphase/VoF/alphaEqn.H
@ -161,12 +161,11 @@

        surfaceScalarField phir(phic*mixture.nHatf());

-        tmp<surfaceScalarField> talphaPhi1Un
-        (
+        alphaPhiUn =
            fvc::flux
            (
-                phiCN(),
-                cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(),
+                phi,
+                alpha1,
                alphaScheme
            )
          + fvc::flux
@ -174,19 +173,18 @@
               -fvc::flux(-phir, alpha2, alpharScheme),
                alpha1,
                alpharScheme
-            )
-        );
+            );

        if (MULESCorr)
        {
-            tmp<surfaceScalarField> talphaPhi1Corr(talphaPhi1Un() - alphaPhi10);
+            tmp<surfaceScalarField> talphaPhi1Corr(alphaPhiUn - alphaPhi10);
            volScalarField alpha10("alpha10", alpha1);

            MULES::correct
            (
                geometricOneField(),
                alpha1,
-                talphaPhi1Un(),
+                alphaPhiUn,
                talphaPhi1Corr.ref(),
                Sp,
                (-Sp*alpha1)(),
@ -207,7 +205,7 @@
        }
        else
        {
-            alphaPhi10 = talphaPhi1Un;
+            alphaPhi10 = alphaPhiUn;

            MULES::explicitSolve
            (
--- a/applications/solvers/multiphase/VoF/createAlphaFluxes.H
+++ b/applications/solvers/multiphase/VoF/createAlphaFluxes.H
@ -24,3 +24,18 @@ surfaceScalarField alphaPhi10

 // MULES Correction
 tmp<surfaceScalarField> talphaPhi1Corr0;
+
+// MULES compressed flux is registered in case scalarTransport FO needs it.
+surfaceScalarField alphaPhiUn
+(
+    IOobject
+    (
+        "alphaPhiUn",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar(phi.dimensions(), Zero)
+);
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/options
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/options
@ -14,7 +14,6 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2017 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -122,7 +122,7 @@ Foam::fv::VoFSolidificationMeltingSource::VoFSolidificationMeltingSource
    const fvMesh& mesh
 )
 :
-    cellSetOption(sourceName, modelType, dict, mesh),
+    fv::cellSetOption(sourceName, modelType, dict, mesh),
    alphaSolidT_(Function1<scalar>::New("alphaSolidT", coeffs_)),
    L_("L", dimEnergy/dimMass, coeffs_),
    relax_(coeffs_.getOrDefault("relax", 0.9)),
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H
@ -99,7 +99,7 @@ namespace fv

 class VoFSolidificationMeltingSource
 :
-    public cellSetOption
+    public fv::cellSetOption
 {
    // Private data

--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceIO.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceIO.C
@ -31,7 +31,7 @@ License

 bool Foam::fv::VoFSolidificationMeltingSource::read(const dictionary& dict)
 {
-    if (cellSetOption::read(dict))
+    if (fv::cellSetOption::read(dict))
    {
        alphaSolidT_ = Function1<scalar>::New("alphaSolidT", coeffs_);
        coeffs_.readEntry("L", L_);
--- a/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/createFields.H
@ -129,10 +129,5 @@ compressibleInterPhaseTransportModel turbulence
 // Mask field for zeroing out contributions on hole cells
 #include "createCellMask.H"

-surfaceScalarField faceMask
-(
-    localMin<scalar>(mesh).interpolate(cellMask)
-);
-
 // Create bool field with interpolated cells
 #include "createInterpolatedCells.H"
--- a/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/overCompressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/overCompressibleInterDyMFoam.C
@ -80,9 +80,6 @@ int main(int argc, char *argv[])
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createDynamicFvMesh.H"
-    pimpleControl pimple(mesh);
-
-    #include "createTimeControls.H"
    #include "createDyMControls.H"
    #include "createFields.H"

@ -109,7 +106,7 @@ int main(int argc, char *argv[])

    while (runTime.run())
    {
-        #include "readControls.H"
+        #include "readDyMControls.H"

        if (LTS)
        {
@ -154,25 +151,7 @@ int main(int argc, char *argv[])
                    // Update cellMask field for blocking out hole cells
                    #include "setCellMask.H"
                    #include "setInterpolatedCells.H"
-
-                    faceMask =
-                        localMin<scalar>(mesh).interpolate(cellMask.oldTime());
-
-                    // Zero Uf on old faceMask (H-I)
-                    Uf *= faceMask;
-
-                    const surfaceVectorField Uint(fvc::interpolate(U));
-                    // Update Uf and phi on new C-I faces
-                    Uf += (1-faceMask)*Uint;
-
-                    // Update Uf boundary
-                    forAll(Uf.boundaryField(), patchI)
-                    {
-                        Uf.boundaryFieldRef()[patchI] =
-                            Uint.boundaryField()[patchI];
-                    }
-
-                    phi = mesh.Sf() & Uf;
+                    #include "correctPhiFaceMask.H"

                    // Correct phi on individual regions
                    if (correctPhi)
@ -182,12 +161,6 @@ int main(int argc, char *argv[])

                    mixture.correct();

-                    // Zero phi on current H-I
-                    faceMask = localMin<scalar>(mesh).interpolate(cellMask);
-
-                    phi *= faceMask;
-                    U *= cellMask;
-
                    // Make the flux relative to the mesh motion
                    fvc::makeRelative(phi, U);

@ -202,10 +175,6 @@ int main(int argc, char *argv[])
            #include "alphaControls.H"
            #include "compressibleAlphaEqnSubCycle.H"

-            const surfaceScalarField faceMask
-            (
-                localMin<scalar>(mesh).interpolate(cellMask)
-            );
            rhoPhi *= faceMask;

            turbulence.correctPhasePhi();
--- a/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/pEqn.H
@ -10,19 +10,6 @@
        fvc::flux(HbyA)
    );

-    if (ddtCorr)
-    {
-        surfaceScalarField faceMaskOld
-        (
-            localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-        );
-        phiHbyA +=
-            MRF.zeroFilter
-            (
-                fvc::interpolate(rho*rAU)*faceMaskOld*fvc::ddtCorr(U, Uf)
-            );
-    }
-
    MRF.makeRelative(phiHbyA);

    surfaceScalarField phig
--- a/applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/BinghamPlastic/BinghamPlastic.C
+++ b/applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/BinghamPlastic/BinghamPlastic.C
@ -6,6 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2014-2015 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -130,9 +131,9 @@ bool Foam::mixtureViscosityModels::BinghamPlastic::read
 {
    plastic::read(viscosityProperties);

-    plasticCoeffs_.readEntry("yieldStressCoeff", yieldStressCoeff_);
-    plasticCoeffs_.readEntry("yieldStressExponent", yieldStressExponent_);
-    plasticCoeffs_.readEntry("yieldStressOffset", yieldStressOffset_);
+    plasticCoeffs_.readEntry("BinghamCoeff", yieldStressCoeff_);
+    plasticCoeffs_.readEntry("BinghamExponent", yieldStressExponent_);
+    plasticCoeffs_.readEntry("BinghamOffset", yieldStressOffset_);

    return true;
 }
--- a/applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/mixtureViscosityModel/mixtureViscosityModel.C
+++ b/applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/mixtureViscosityModel/mixtureViscosityModel.C
@ -27,6 +27,7 @@ License

 #include "mixtureViscosityModel.H"
 #include "volFields.H"
+#include "surfaceMesh.H"

 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H
@ -120,21 +120,6 @@ if (p_rgh.needReference())
 mesh.setFluxRequired(p_rgh.name());
 mesh.setFluxRequired(alpha1.name());

-// MULES compressed flux is registered in case scalarTransport FO needs it.
-surfaceScalarField alphaPhiUn
-(
-    IOobject
-    (
-        "alphaPhiUn",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    mesh,
-    dimensionedScalar(phi.dimensions(), Zero)
-);
-
 #include "createMRF.H"
 #include "createFvOptions.H"

--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/Make/options
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/Make/options
@ -1,5 +1,6 @@
 EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
@ -9,6 +10,7 @@ EXE_INC = \

 LIB_LIBS = \
    -lfiniteVolume \
+    -lmeshTools \
    -lgeometricVoF \
    -ltwoPhaseMixture \
    -linterfaceProperties \
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@ -126,20 +126,5 @@ if (p_rgh.needReference())
 mesh.setFluxRequired(p_rgh.name());
 mesh.setFluxRequired(alpha1.name());

-// MULES compressed flux is registered in case scalarTransport FO needs it.
-surfaceScalarField alphaPhiUn
-(
-    IOobject
-    (
-        "alphaPhiUn",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    mesh,
-    dimensionedScalar(phi.dimensions(), Zero)
-);
-
 #include "createMRF.H"
 #include "createFvOptions.H"
--- a/applications/solvers/multiphase/interFoam/overInterDyMFoam/CourantNo.H
+++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/CourantNo.H
@ -36,7 +36,10 @@ scalar meanCoNum = 0.0;

 if (mesh.nInternalFaces())
 {
-    surfaceScalarField phiMask(localMin<scalar>(mesh).interpolate(cellMask));
+    surfaceScalarField phiMask
+    (
+        localMin<scalar>(mesh).interpolate(cellMask + interpolatedCells)
+    );

    scalarField sumPhi
    (
--- a/applications/solvers/multiphase/interFoam/overInterDyMFoam/createDyMControls.H
+++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/createDyMControls.H
@ -1,30 +0,0 @@
-bool correctPhi
-(
-    pimple.dict().getOrDefault("correctPhi", true)
-);
-
-bool checkMeshCourantNo
-(
-    pimple.dict().getOrDefault("checkMeshCourantNo", false)
-);
-
-bool moveMeshOuterCorrectors
-(
-    pimple.dict().getOrDefault("moveMeshOuterCorrectors", false)
-);
-
-
-bool massFluxInterpolation
-(
-    pimple.dict().getOrDefault("massFluxInterpolation", false)
-);
-
-bool adjustFringe
-(
-    pimple.dict().getOrDefault("oversetAdjustPhi", false)
-);
-
-bool ddtCorr
-(
-    pimple.dict().getOrDefault("ddtCorr", true)
-);
--- a/applications/solvers/multiphase/interFoam/overInterDyMFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/createFields.H
@ -147,19 +147,4 @@ if (p_rgh.needReference())
 mesh.setFluxRequired(p_rgh.name());
 mesh.setFluxRequired(alpha1.name());

-// MULES compressed flux is registered in case scalarTransport FO needs it.
-surfaceScalarField alphaPhiUn
-(
-    IOobject
-    (
-        "alphaPhiUn",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    mesh,
-    dimensionedScalar(phi.dimensions(), Zero)
-);
-
 #include "createMRF.H"
--- a/applications/solvers/multiphase/interFoam/overInterDyMFoam/overInterDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/overInterDyMFoam.C
@ -49,13 +49,9 @@ Description
 #include "turbulentTransportModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
-#include "CorrectPhi.H"
 #include "fvcSmooth.H"
 #include "cellCellStencilObject.H"
 #include "localMin.H"
-#include "interpolationCellPoint.H"
-#include "transform.H"
-#include "fvMeshSubset.H"
 #include "oversetAdjustPhi.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -76,8 +72,7 @@ int main(int argc, char *argv[])
    #include "createTime.H"
    #include "createDynamicFvMesh.H"
    #include "initContinuityErrs.H"
-    pimpleControl pimple(mesh);
-    #include "createTimeControls.H"
+
    #include "createDyMControls.H"
    #include "createFields.H"
    #include "createAlphaFluxes.H"
@ -102,6 +97,7 @@ int main(int argc, char *argv[])
        #include "correctPhi.H"
    }
    #include "createUf.H"
+    #include "createControls.H"

    #include "setCellMask.H"
    #include "setInterpolatedCells.H"
@ -119,7 +115,7 @@ int main(int argc, char *argv[])

    while (runTime.run())
    {
-        #include "readControls.H"
+        #include "readDyMControls.H"

        if (LTS)
        {
@ -164,27 +160,7 @@ int main(int argc, char *argv[])
                    // Update cellMask field for blocking out hole cells
                    #include "setCellMask.H"
                    #include "setInterpolatedCells.H"
-
-                    const surfaceScalarField faceMaskOld
-                    (
-                        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-                    );
-
-                    // Zero Uf on old faceMask (H-I)
-                    Uf *= faceMaskOld;
-
-                    const surfaceVectorField Uint(fvc::interpolate(U));
-                    // Update Uf and phi on new C-I faces
-                    Uf += (1-faceMaskOld)*Uint;
-
-                    // Update Uf boundary
-                    forAll(Uf.boundaryField(), patchI)
-                    {
-                        Uf.boundaryFieldRef()[patchI] =
-                            Uint.boundaryField()[patchI];
-                    }
-
-                    phi = mesh.Sf() & Uf;
+                    #include "correctPhiFaceMask.H"

                    // Correct phi on individual regions
                    if (correctPhi)
@ -194,17 +170,8 @@ int main(int argc, char *argv[])

                    mixture.correct();

-                    // Zero phi on current H-I
-                    const surfaceScalarField faceMask
-                    (
-                        localMin<scalar>(mesh).interpolate(cellMask)
-                    );
-                    phi *= faceMask;
-                    U   *= cellMask;
-
                    // Make the flux relative to the mesh motion
                    fvc::makeRelative(phi, U);
-
                }

                if (mesh.changing() && checkMeshCourantNo)
@ -213,14 +180,9 @@ int main(int argc, char *argv[])
                }
            }

-
            #include "alphaControls.H"
            #include "alphaEqnSubCycle.H"

-            const surfaceScalarField faceMask
-            (
-                localMin<scalar>(mesh).interpolate(cellMask)
-            );
            rhoPhi *= faceMask;

            mixture.correct();
--- a/applications/solvers/multiphase/interFoam/overInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/pEqn.H
@ -1,8 +1,5 @@
 {
    rAU = 1.0/UEqn.A();
-    //mesh.interpolate(rAU);
-
-    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));

    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

@ -12,22 +9,8 @@
    //HbyA = rAU*UEqn.H();
    HbyA = constrainHbyA(rAU*H, U, p_rgh);

-    if (massFluxInterpolation)
-    {
-        #include "interpolatedFaces.H"
-    }
-
    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));

-    if (ddtCorr)
-    {
-        surfaceScalarField faceMaskOld
-        (
-            localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-        );
-        phiHbyA +=
-            fvc::interpolate(rho*rAU)*faceMaskOld*fvc::ddtCorr(U, Uf);
-    }
    MRF.makeRelative(phiHbyA);

    if (p_rgh.needReference())
@ -37,13 +20,6 @@
        fvc::makeAbsolute(phiHbyA, U);
    }

-    if (adjustFringe)
-    {
-        fvc::makeRelative(phiHbyA, U);
-        oversetAdjustPhi(phiHbyA, U);
-        fvc::makeAbsolute(phiHbyA, U);
-    }
-
    surfaceScalarField phig
    (
        (
--- a/applications/solvers/multiphase/interFoam/overInterDyMFoam/readControls.H
+++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/readControls.H
@ -1,16 +0,0 @@
-#include "readTimeControls.H"
-
-correctPhi = pimple.dict().getOrDefault("correctPhi", false);
-
-checkMeshCourantNo =
-    pimple.dict().getOrDefault("checkMeshCourantNo", false);
-
-moveMeshOuterCorrectors =
-    pimple.dict().getOrDefault("moveMeshOuterCorrectors", false);
-
-massFluxInterpolation =
-    pimple.dict().getOrDefault("massFluxInterpolation", false);
-
-ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);
-
-adjustFringe =  pimple.dict().getOrDefault("oversetAdjustPhi", false);
--- a/applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/createFields.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/createFields.H
@ -149,10 +149,5 @@ surfaceScalarField alphaPhi10
 // Mask field for zeroing out contributions on hole cells
 #include "createCellMask.H"

-surfaceScalarField faceMask
-(
-    localMin<scalar>(mesh).interpolate(cellMask)
-);
-
 // Create bool field with interpolated cells
 #include "createInterpolatedCells.H"
--- a/applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/overInterPhaseChangeDyMFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/overInterPhaseChangeDyMFoam.C
@ -82,9 +82,7 @@ int main(int argc, char *argv[])
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createDynamicFvMesh.H"
-    pimpleControl pimple(mesh);

-    #include "createTimeControls.H"
    #include "createDyMControls.H"
    #include "initContinuityErrs.H"
    #include "createFields.H"
@ -118,7 +116,7 @@ int main(int argc, char *argv[])

    while (runTime.run())
    {
-        #include "readControls.H"
+        #include "readDyMControls.H"

        // Store divU from the previous mesh so that it can be mapped
        // and used in correctPhi to ensure the corrected phi has the
@ -153,26 +151,7 @@ int main(int argc, char *argv[])
                     // Update cellMask field for blocking out hole cells
                    #include "setCellMask.H"
                    #include "setInterpolatedCells.H"
-
-                    faceMask =
-                        localMin<scalar>(mesh).interpolate(cellMask.oldTime());
-
-
-                    // Zero Uf on old faceMask (H-I)
-                    Uf *= faceMask;
-
-                    const surfaceVectorField Uint(fvc::interpolate(U));
-                    // Update Uf and phi on new C-I faces
-                    Uf += (1-faceMask)*Uint;
-
-                    // Update Uf boundary
-                    forAll(Uf.boundaryField(), patchI)
-                    {
-                        Uf.boundaryFieldRef()[patchI] =
-                            Uint.boundaryField()[patchI];
-                    }
-
-                    phi = mesh.Sf() & Uf;
+                    #include "correctPhiFaceMask.H"

                    if (correctPhi)
                    {
@ -181,12 +160,6 @@ int main(int argc, char *argv[])

                    mixture->correct();

-                    // Zero phi on current H-I
-                    faceMask = localMin<scalar>(mesh).interpolate(cellMask);
-
-                    phi *= faceMask;
-                    U   *= cellMask;
-
                    // Make the flux relative to the mesh motion
                    fvc::makeRelative(phi, U);
                }
--- a/applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/pEqn.H
@ -8,16 +8,6 @@
        fvc::flux(HbyA)
    );

-
-    if (ddtCorr)
-    {
-        surfaceScalarField faceMaskOld
-        (
-            localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-        );
-        phiHbyA += faceMaskOld*fvc::ddtCorr(U, Uf);
-    }
-
    if (p_rgh.needReference())
    {
        fvc::makeRelative(phiHbyA, U);
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
sergio	95fe43b079	TUT: Add and update overset tutorials	2021-08-23 19:31:22 -07:00
sergio	10b7479dca	COMP: Compilation of src with Allwmake	2021-08-23 19:30:43 -07:00
sergio	054e3ebc8b	ENH: Update overset solvers	2021-08-23 19:29:54 -07:00
sergio	3dc63dc67c	ENH: Updating tutorials	2021-08-23 12:06:18 -07:00
sergio	141b8376f3	ENh: Updating overset solvers 1) Deleting ddtCorr option 2) Deleting fluxInterpolation flux option 3) Reducing control readers files	2021-08-23 12:05:36 -07:00
sergio	72b9a002a2	ENH: adding cellTypes for use in overset to cellSetOption	2021-08-23 12:05:35 -07:00
sergio	808565a673	ENH: Delete unused ot updated files	2021-08-23 12:03:33 -07:00
sergio	5da36a3386	ENH: overset lib modifications 1) Creating static/dynamic oversetFvMesh 2) Adding helper functions for overset solvers in include 3) Changes to allow overlapping patches from the inset mesh on the background mesh for cellVolumeWeightCell, trackingInverseDistanc and inverseDistance stencils.	2021-08-23 12:03:33 -07:00
Andrew Heather	8be3f3d910	BUG: Corrected ddtCorr switch - see #2183	2021-08-17 15:36:34 +01:00
mattijs	9678b1dc22	ENH: splitMeshRegions: combine cellZones. See #2046 addZones: specifies the name instead of synthesising	2021-08-11 16:14:23 +01:00
sergio	9bf030d86a	BUG: Fixing frozenFlow option with implicit fvMatrix in CHT solvers	2021-08-06 12:54:44 -07:00
sergio	b826e23d2f	TUT: Rename tutorial scripts to Allrun	2021-08-06 11:15:31 -07:00
sergio	93b0a57791	BUG: Fixes #2175 and #2176	2021-08-06 11:07:48 -07:00
Andrew Heather	b78dccdd51	ENH: noiseModel - protected log10 calls	2021-08-05 16:06:33 +01:00
Andrew Heather	d8eecc5e2a	Merge branch 'issue-2093-custom-U-in-devRhoReff-and-devReff' into 'develop' ENH: enable custom input of U for devRhoReff and devReff See merge request Development/openfoam!479	2021-08-05 12:41:24 +00:00
Kutalmis Bercin	a870be6511	BUG: forces: ensure UName entry can be used (fixes #2093 ) Previously, for basic incompressible and compressible simulations, the "force" function object has not been using the user-specified "UName" for the "devRhoReff" computation (affecting the tangential component), but using the "U" of the latest available step. In contrast, the user-specified "pName" has always been being used correctly. This has been causing issues for users when they wish to use a specific "UMean" field in various force and forceCoeff function object computations.	2021-08-05 12:40:47 +00:00
Kutalmis Bercin	3207e315b2	ENH: enable custom input of U for devRhoReff and devReff	2021-08-05 12:40:47 +00:00
Andrew Heather	c8a52fad96	Merge branch 'issue-2132-turbulenceFields' into 'develop' ENH: turbulenceFields: various improvements See merge request Development/openfoam!476	2021-08-04 11:37:16 +00:00
Kutalmis Bercin	14b328d990	ENH: turbulenceFields: enable custom prefix for output fields	2021-08-04 11:36:53 +00:00
Kutalmis Bercin	2e703543ff	BUG: turbulenceFields: unset duplicate already-registered fields	2021-08-04 11:36:53 +00:00
Kutalmis Bercin	be7d6b3afb	DOC: turbulenceFields: improve header documentation	2021-08-04 11:36:53 +00:00
Kutalmis Bercin	c6c575052f	STYLE: turbulenceFields: apply more recent OpenFOAM-coding practices	2021-08-04 11:36:53 +00:00
Kutalmis Bercin	99291a7ac5	BUG: turbulenceFields: use omega funcs of turbulence models (fixes #2132 )	2021-08-04 11:36:53 +00:00
Andrew Heather	084f89ee94	Merge branch 'issue-2122-dmd-snapshot0' into 'develop' BUG: DMD: write snapshot0 at start times See merge request Development/openfoam!482	2021-08-04 09:23:49 +00:00
Kutalmis Bercin	9c46f86048	ENH: makeSolidReaction: modify macros to allow new models	2021-08-04 09:23:20 +00:00
Bas Nieuwboer	c64c312bae	TUT: hopper: parameterise blockMeshDict content (#2134 )	2021-08-04 09:23:20 +00:00
Kutalmis Bercin	f482f74277	BUG: DMD: write snapshot0 at start times (fixes #2122 )	2021-08-04 09:23:20 +00:00
Andrew Heather	f687a21761	Merge branch 'issue-1615-solver-verbosity' into 'develop' ENH: linear solvers: add variable-specific debug flags See merge request Development/openfoam!478	2021-08-04 09:15:31 +00:00
Kutalmis Bercin	6dedfe078a	ENH: linear solvers: add variable-specific debug flags Introduces a new optional keyword of label type 'log' to linear-solver dictionaries to enable variable-specific debug statements. For example, in fvOptions file: solvers { p { solver GAMG; ... log 2; } U { ... log 0; } } The meanings of values of 'log' are: log 0; <!-- no output log 1; <!-- standard output log 2; <!-- debug output // values higher than 2 are expected to have no effect This keyword does not directly affect the operations of various DebugSwitches and backward compatibility has been ensured in exchange of code cleanness. The related DebugSwitches are: DebugSwitches { SolverPerformance 0; GAMG 0; PCG 0; PBiCG 0; smoothSolver 0; }	2021-08-04 09:15:19 +00:00
Matej Forman	d9b200af28	TUT: rotatingFanInRoom: perturb locationInMesh (fixes #2162 )	2021-08-04 09:15:19 +00:00
Kutalmis Bercin	d6ff22e5e7	ENH: KirchhoffShell: simplification of log output	2021-08-04 09:15:19 +00:00
Andrew Heather	fa03050e03	Merge branch 'feature-implicitCyclics-merge' into 'develop' Implicit treatment of coupled boundary conditions See merge request Development/openfoam!477	2021-08-03 20:08:51 +00:00
Sergio Ferraris	53af23b9fb	Implicit treatment of coupled boundary conditions	2021-08-03 20:08:49 +00:00
mattijs	11e0db96d3	ENH: snappyHexMeshDict: added missing documentation	2021-08-02 14:01:12 +01:00
mattijs	60aeca32a8	ENH: checkMesh: add stabilisation. Fixes #2171	2021-07-28 17:56:48 +01:00
Mark Olesen	9ab8c525cf	CONFIG: set API level to 2107 - liquid film, blockMesh changes, multiWorld communicators	2021-07-28 17:17:41 +02:00
Mark Olesen	9b0051c40d	SUBMODULE: update visualization	2021-07-28 16:54:00 +02:00
Mark Olesen	270f2eecf9	ENH: report API information in header #2126 - previously reported the WM_PROJECT_VERSION in the top banner but this makes it overly sensitive to arbitrary user naming. Change to the current API in the comment banner, include the WM_PROJECT_VERSION as part of the reported "Build: ..." string instead.	2021-07-28 16:53:53 +02:00
Mattijs Janssens	d379f0b61d	Merge branch 'feature-multiWorld' into 'develop' improve handling of multiple world communication See merge request Development/openfoam!480	2021-07-28 14:50:05 +00:00
Mark Olesen	9d26b25c46	TUT: multi-world tutorial setup with circular connectivity	2021-07-28 16:21:14 +02:00
Mark Olesen	2bb91e354b	ENH: improve handling of database mapping for multi-world - cleanup initialization (to avoid stealing tokens)	2021-07-28 16:07:31 +02:00
Mark Olesen	56a5a43910	ENH: simplify mapped patch cloning in extrudeToRegionMesh	2021-07-28 16:07:31 +02:00
Mark Olesen	ec4265343a	ENH: use multi-world connection object for managing mappedPatches	2021-07-28 16:07:30 +02:00
Mark Olesen	b2ec5af9e0	ENH: centralized MeshObject handler for multi-world connections	2021-07-28 16:07:30 +02:00
mattijs	13740de427	BUG: mpi: use group-local communicator allocation. Fixes #2158	2021-07-28 16:07:30 +02:00
Mark Olesen	8f1bad1cb4	BUG: decomposing ConstantField loses information (#2101 ) - previously had a patch-size check, which caused a zero-sized patch to be "sticky" as uniform Zero regardless of what the original input stated (eg, "constant 0.5") This was caused by short-circuiting logic that circumvented reading on zero-size patches.	2021-07-28 16:05:46 +02:00
Mark Olesen	10dd0d34d1	ENH: combine patchDistMethod New methods (code reduction)	2021-07-28 13:14:44 +02:00
Mark Olesen	a530b505ed	TUT: replace m4 versions of angledDuct blockMeshDict with native forms - dictionary parameters to change the bend angle - blockMesh transform parameters to allow different orientations	2021-07-28 13:14:44 +02:00
Mark Olesen	29acee5a50	TUT: use new blockMesh support (faces, transforms)	2021-07-28 13:14:44 +02:00
Mark Olesen	d31f351c15	ENH: blockMesh enhancements - support non-uniform scaling, prescaling and cartesian coordinate transformations. Eg, stretch in one direction and then rotate ``` prescale (1.5 1 1); transform { origin (0 0 0); rotation { type axisAngle; axis (0 0 1); angle 45; } } ``` - support "transformed" versions of blockMesh vertices, topology. With the additional of transformations etc, a simplistic application of a single scale parameter is no longer sufficient. new: blMesh.vertices(true); old: blMesh.vertices() * blMesh.scaleFactor(); new: blMesh.topology(true); old: N/A - add individual edge access for blockDescriptor. Saves copying and duplicate calculations. - handle '(block face)' specification for curved faces, which is ok for external block faces, but likely somewhat questionable if used for internal block faces.	2021-07-28 13:14:44 +02:00
Mark Olesen	4eeff16735	ENH: handle min/span as alternative to min/max for box selection - searchable box, boxToCell, boxToFace, boxToPoint	2021-07-28 11:35:12 +02:00
Mark Olesen	d49566ea34	ENH: improve alignment of tetCell and cellShape - added hexCell for some more encapsulation	2021-07-28 11:35:08 +02:00
Mark Olesen	5b09e59b3f	ENH: support direct specification of coordinate rotation "none" - Can specify ``` transform { origin (1 2 3); rotation none; } ``` instead of the longer form ``` transform { origin (1 2 3); e1 (1 0 0); e3 (0 0 1); } ``` COMPAT: remove old (pre-1812) Euler and STARCD keywords - use "angles", remove old compat name "rotation" as clutter and possible confusion with "coordinate rotation" STYLE: remove coordinate rotation move constructors - an unnecessary holdover from older code.	2021-07-28 10:02:27 +02:00
Mark Olesen	e2be2289a1	ENH: finer granularity for CleanFunctions - separate handling of auxiliary files vs time directories - restore0Dir: avoid removing 0/ if 0.orig/ does not exist	2021-07-28 10:02:27 +02:00
Mark Olesen	eb7d0f0ef0	COMP: add missing includes for runtime-selection	2021-07-28 10:02:27 +02:00
Mark Olesen	847b61d1f9	GIT: fix some merge and edit rubbish	2021-07-28 10:02:27 +02:00
sergio	621d400db2	BUG: Fix to const Uwall helper function	2021-07-27 08:24:33 -07:00
sergio	3e70154ff1	BUG: Fixing movingWall BC mesh flux	2021-07-26 11:03:30 -07:00
Sergio Ferraris	6f20499b57	Merge branch 'feature-liquidFilm' into 'develop' finiteArea: new models and solvers for films See merge request Development/openfoam!475	2021-07-16 16:39:48 +00:00
sergio	1f3c35841f	ENH: Updating tutorials for finite area	2021-07-16 16:38:06 +00:00
sergio	3f670d88ea	ENH: Initialization for area and single layer films for lagrangian	2021-07-16 16:38:06 +00:00
sergio	7d8dba564c	TUT: heatTransfer: new tutorial for thermalShellModel	2021-07-16 16:38:06 +00:00
sergio	5b21050a0b	TUT: incompressible/lagrangian: new tutorials for surfaceFilmModel	2021-07-16 16:38:06 +00:00
sergio	4b05b766d3	ENH: limitVelocity: new faOption	2021-07-16 16:38:06 +00:00
sergio	18ffb1ad07	ENH: KirchhoffShell: replace read() by init()	2021-07-16 16:38:06 +00:00
sergio	d3807a6ee3	ENH: vibrationShellModel and thermalShellModel: remove read() STYLE: vibrationShellModel and thermalShellModel: regroup the access functions	2021-07-16 16:38:06 +00:00
sergio	a9aa6cb0db	ENH: kinematicParcelFoam: new solver and tutorial The original thermoSurfaceFilm sub-models were divided between kinematicSurfaceFilm and thermoSurfaceFilm in order to use the surfaceFilm model in a kinematicCloud. The film interaction models are now in a kinematicSurface class which can be used in a kinematic cloud adding constant thermal properties (p and T) for some sub-models, e.g. drySplashInteraction/wetSplashInteraction. pRef and Tref were added to the kinematicSurfaceFilm as entry to the regionFilm when used with a kinematic cloud. In the finite area surface film model Tref, pRef are stored in filmSubModelBase. TUT: kinematicParcelFoam: new tutorial pitzDailySprinkles	2021-07-16 16:38:06 +00:00
sergio	4476863907	ENH: Adding radiation to thermalShell	2021-07-16 16:38:06 +00:00
sergio	384b9f8100	ENH: reactingParcelFoam: allow to run single species	2021-07-16 16:38:06 +00:00
sergio	f98633e427	ENH: liquidFilm subModels: add force models	2021-07-16 16:38:06 +00:00
sergio	e30cbfc6db	ENH: liquidFilm subModels: add injection models ENH: Several additions to thinFilm models 1) Changing the faSolution entry nonOrthoCorr to nFilmCorr entry. This entry specify the number of iterations on the film thickness Eq which does not correct for laplacian non-orthogonality. 2) Changing laminar turbulent model for the skin friction exerted by the flow to the film. 3) Adding optional relaxation factors for source terms coming from the added mass of the lagragian particles into the film such as rhoSp, pnSp and Usp. Similarly the film pressure (pf) can be relaxed in the inner loops. 4) Adding optional entries to curvatureSepration injection model. minimum (fThreshold) force and minimum curvature (minInvR1) for separation were added to have more control on determining the points of film separation ENH: ThermoSurfaceFilm in lagragian to search for finiteArea film models	2021-07-16 16:38:06 +00:00
sergio	b62306ddc5	ENH: Adding regionFaModels library to solvers	2021-07-16 16:38:06 +00:00
sergio	be09967acf	ENH: regionFaModel: various updates	2021-07-16 16:38:06 +00:00
sergio	efd3f59295	ENH: finiteArea: new models for shallow film modelling - new model framework: liquidFilm - new film turbulence models: filmTurbulenceModel with friction models - new velocity boundary condition: velocityFilmShell to handle/evolve regionFa film - new function object: setTimeStep to control the simulation time-step based on regionFa Courant number - add support for the absorption of Lagrangian particles into films	2021-07-16 16:38:06 +00:00
sergio	b006f5c44e	ENH: movingWallVelocity: add Uwall function	2021-07-16 16:38:06 +00:00
sergio	ce5e2e7e1c	ENH: volSurfaceMapping: new mapToField and mapToSurface funcs	2021-07-16 16:38:06 +00:00
Mark Olesen	db9460f0bc	ENH: copy construct FixedList from fixed subset of input - remove construct from two iterators (#2083)	2021-07-15 16:47:27 +02:00
Mark Olesen	0c15f773e0	STYLE: remove/replace unneeded SortableList	2021-07-15 16:47:27 +02:00
Mark Olesen	aca366d293	TUT: replace remaining templated thermoType TUT: some simplification for multiWorld test STYLE: remove some editor junk	2021-07-15 16:47:22 +02:00
Mark Olesen	09aa3d8b43	COMP: avoid incomplete types in runTime tables, unneeded typedefs etc	2021-07-15 14:31:54 +02:00
Mark Olesen	ccc751e55d	ENH: centralize checks for keyword age into Foam::error STYLE: additional note about include/exclude HashTable.C	2021-07-15 12:53:02 +02:00
sergio	c1c1799051	BUG: Caching alphaPhiUn for VoF for scalarTransport FOcin mixture. #2053	2021-07-12 14:22:08 -07:00
sergio	20ca8e8aac	BUG: Fixing entries for BinghamPlastic model. Fixes #1150	2021-07-08 13:23:39 -07:00
sergio	bda5a87213	BUG: Correcting calculation of phiGas used in coupled BC. Fixes #2149	2021-07-08 12:26:04 -07:00
sergio	2ab81c3b69	BUG: Correcting surface tension calculation from pair order. Fixes #2051	2021-07-08 12:17:05 -07:00
Kutalmis Bercin	70c697fdcb	BUG: momentum: boundary field mismatch on "empty" patches (fixes #2121 )	2021-07-05 15:20:41 +02:00
Kutalmis Bercin	64deb4e4b8	BUG: state error in volRegion::calculateCache() (fixes #2121 ) - need to change requireUpdate_ at the beginning of calculateCache otherwise it is not updated for "ALL" selection (causing unnecessary repeated calculations). For other modes, it would have tripped an error in FULLDEBUG.	2021-07-05 15:17:56 +02:00
Mark Olesen	8cda8e5c3a	ENH: downgrade lookupObject failure from abort() to exit() - lets function objects catch and handle errors (fixes #2144, related to #1779)	2021-07-05 14:55:48 +02:00
Mark Olesen	cc5aa20931	ENH: simplify exit/abort handling and jobinfo (code reduction) - handle failures more robustly - add static shutdown() for similarity with Pstream etc.	2021-07-05 14:55:48 +02:00
Mark Olesen	8a1f667067	COMP: cannot resolve read in UnsortedMeshedSurface constructor - fix linkage for cellModeller compat methods	2021-07-05 14:55:37 +02:00
Mark Olesen	00990698ad	COMP: remove construct SortableList from two iterators (#2083 ) - unused, but missed in `40567b844a`	2021-07-05 14:55:37 +02:00
Mark Olesen	8884d58a52	STYLE: fix/correct some tests. Remove useless tests	2021-07-05 14:55:37 +02:00
Mark Olesen	6a390f6ad1	COMP: remove includes of defunct locations (fixes #2143 )	2021-06-29 13:43:44 +02:00
Mark Olesen	1898395038	ENH: group bools together, relocate override of updateT (basicThermo) - all flags before loaders. Avoids "uninitialized when used here" warning while preserving the field ownership logic. - relocate override of updateT basicThermo into constructor body for clearer logic and initialization STYLE: plain bool instead of Switch for dpdt flag	2021-06-29 13:43:44 +02:00
Mark Olesen	252326df05	STYLE: add fa/fv namespace qualifiers for fa/fv option	2021-06-29 13:43:44 +02:00
Andrew Heather	f75379633c	CONFIG: Updated version to com following release of OpenFOAM v2106	2021-06-28 17:33:22 +01:00
				`@ -0,0 +1 @@`
				`regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();`