CONFIG: bump 2106 patch level to 211215

BACKPORT: config script bugfixes
- align with fixes for v2112 - align foamConfigurePaths with v2112 modifications CONFIG: scan scotch-64.h as possible fallback (#1904) - For RedHat 8, the /usr/include/scotch.h is actually wrapped with a condition check based on LP64 vs ILP32, so need to check with those headers if the first one failed.
2021-12-15 11:08:43 +01:00 · 2021-12-10 14:00:10 +01:00 · 2021-12-09 18:00:59 +01:00 · 2021-09-08 09:29:27 +02:00
2232 changed files with 18183 additions and 60171 deletions
--- a/29
+++ b/29
@ -8,7 +8,7 @@ wmake -check-dir "$WM_PROJECT_DIR" 2>/dev/null || {
 }
 if [ -f "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments ]
 then  . "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments || \
-    echo "Argument parse error"
+    echo "Argument parse error";
 else
    echo "Error (${0##*/}) : WM_PROJECT_DIR appears to be incorrect"
    echo "    Check your OpenFOAM environment and installation"
@ -23,23 +23,23 @@ else
    echo "mpirun=$(command -v mpirun || true)"
 fi
 echo
-# Report compiler information. First non-blank line from --version output
+# Report compiler information
 compiler="$(wmake -show-path-cxx 2>/dev/null || true)"
 if [ -x "$compiler" ]
 then
    echo "compiler=$compiler"
-    "$compiler" --version 2>/dev/null | sed -e '/^$/d;q'
+    "$compiler" --version 2>/dev/null | sed -ne '1p'
 else
    echo "compiler=unknown"
 fi

 echo
-echo ========================================
+echo "========================================"
 date "+%Y-%m-%d %H:%M:%S %z" 2>/dev/null || echo "date is unknown"
 echo "Starting compile ${WM_PROJECT_DIR##*/} ${0##*/}"
 echo "  $WM_COMPILER ${WM_COMPILER_TYPE:-system} compiler"
 echo "  ${WM_OPTIONS}, with ${WM_MPLIB} ${FOAM_MPI}"
-echo ========================================
+echo "========================================"
 echo

 # Compile tools for wmake
@ -62,21 +62,24 @@ else
    echo "Skip ThirdParty (no directory)"
 fi

-# OpenFOAM libraries
+echo "========================================"
+echo "Compile OpenFOAM libraries"
+echo
 src/Allwmake $targetType $*

-# OpenFOAM applications
+echo "========================================"
+echo "Compile OpenFOAM applications"
+echo
 applications/Allwmake $targetType $*

-
 # Additional components/modules
 if [ "$FOAM_MODULE_PREFIX" = false ] || [ "$FOAM_MODULE_PREFIX" = none ]
 then
-    echo ========================================
+    echo "========================================"
    echo "OpenFOAM modules disabled (prefix=${FOAM_MODULE_PREFIX})"
    echo
-else
-    # Use wmake -all instead of Allwmake to allow for overrides
+elif [ -d "$WM_PROJECT_DIR/modules" ]
+then
    (cd "$WM_PROJECT_DIR/modules" 2>/dev/null && wmake -all)
 fi

@ -90,7 +93,7 @@ _foamCountDirEntries()
 # Some summary information
 echo
 date "+%Y-%m-%d %H:%M:%S %z" 2>/dev/null || echo "date is unknown"
-echo ========================================
+echo "========================================"
 echo "  ${WM_PROJECT_DIR##*/}"
 echo "  $WM_COMPILER ${WM_COMPILER_TYPE:-system} compiler"
 echo "  ${WM_OPTIONS}, with ${WM_MPLIB} ${FOAM_MPI}"
@ -100,6 +103,6 @@ echo "  patch = $(etc/openfoam -show-patch 2>/dev/null)"
 echo "  bin   = $(_foamCountDirEntries "$FOAM_APPBIN") entries"
 echo "  lib   = $(_foamCountDirEntries "$FOAM_LIBBIN") entries"
 echo
-echo ========================================
+echo "========================================"

 #------------------------------------------------------------------------------
--- a/META-INFO/api-info
+++ b/META-INFO/api-info
@ -1,2 +1,2 @@
-api=2108
-patch=0
+api=2106
+patch=211215
--- a/applications/Allwmake
+++ b/applications/Allwmake
@ -8,17 +8,13 @@ wmake -check-dir "$WM_PROJECT_DIR/applications" 2>/dev/null || {
 }
 if [ -f "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments ]
 then  . "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments || \
-    echo "Argument parse error"
+    echo "Argument parse error";
 else
    echo "Error (${0##*/}) : WM_PROJECT_DIR appears to be incorrect"
    echo "    Check your OpenFOAM environment and installation"
    exit 1
 fi

-echo ========================================
-echo Compile OpenFOAM applications
-echo ========================================
-
 #------------------------------------------------------------------------------

 wmake -all $targetType solvers
@ -28,9 +24,3 @@ wmake -all $targetType utilities
 ## wmake -all $targetType tools

 #------------------------------------------------------------------------------
-
-echo ========================================
-echo Done OpenFOAM applications
-echo ========================================
-
-#------------------------------------------------------------------------------
--- a/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModel.H
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModel.H
@ -138,13 +138,6 @@ public:
            const surfaceScalarField& phi
        );

-        //- Return a clone
-        virtual refPtr<regIOobject> deepClone() const
-        {
-            NotImplemented;
-            return nullptr;
-        }
-

    //- Destructor
    virtual ~PDRDragModel();
--- a/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModelNew.C
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModelNew.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -43,9 +43,9 @@ Foam::autoPtr<Foam::PDRDragModel> Foam::PDRDragModel::New

    Info<< "Selecting drag model " << modelType << endl;

-    auto* ctorPtr = dictionaryConstructorTable(modelType);
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);

-    if (!ctorPtr)
+    if (!cstrIter.found())
    {
        FatalIOErrorInLookup
        (
@ -57,9 +57,7 @@ Foam::autoPtr<Foam::PDRDragModel> Foam::PDRDragModel::New
    }

    return autoPtr<PDRDragModel>
-    (
-        ctorPtr(dict, turbulence, rho, U, phi)
-    );
+        (cstrIter()(dict, turbulence, rho, U, phi));
 }


--- a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModelNew.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModelNew.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -42,9 +42,9 @@ Foam::autoPtr<Foam::XiEqModel> Foam::XiEqModel::New

    Info<< "Selecting flame-wrinkling model " << modelType << endl;

-    auto* ctorPtr = dictionaryConstructorTable(modelType);
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);

-    if (!ctorPtr)
+    if (!cstrIter.found())
    {
        FatalIOErrorInLookup
        (
@ -55,7 +55,7 @@ Foam::autoPtr<Foam::XiEqModel> Foam::XiEqModel::New
        ) << exit(FatalIOError);
    }

-    return autoPtr<XiEqModel>(ctorPtr(dict, thermo, turbulence, Su));
+    return autoPtr<XiEqModel>(cstrIter()(dict, thermo, turbulence, Su));
 }


--- a/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/XiGModel/XiGModelNew.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/XiGModel/XiGModelNew.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -42,9 +42,9 @@ Foam::autoPtr<Foam::XiGModel> Foam::XiGModel::New

    Info<< "Selecting flame-wrinkling model " << modelType << endl;

-    auto* ctorPtr = dictionaryConstructorTable(modelType);
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);

-    if (!ctorPtr)
+    if (!cstrIter.found())
    {
        FatalIOErrorInLookup
        (
@ -55,7 +55,7 @@ Foam::autoPtr<Foam::XiGModel> Foam::XiGModel::New
        ) << exit(FatalIOError);
    }

-    return autoPtr<XiGModel>(ctorPtr(dict, thermo, turbulence, Su));
+    return autoPtr<XiGModel>(cstrIter()(dict, thermo, turbulence, Su));
 }


--- a/applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModelNew.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModelNew.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -45,9 +45,9 @@ Foam::autoPtr<Foam::XiModel> Foam::XiModel::New

    Info<< "Selecting flame-wrinkling model " << modelType << endl;

-    auto* ctorPtr = dictionaryConstructorTable(modelType);
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);

-    if (!ctorPtr)
+    if (!cstrIter.found())
    {
        FatalIOErrorInLookup
        (
@ -59,9 +59,7 @@ Foam::autoPtr<Foam::XiModel> Foam::XiModel::New
    }

    return autoPtr<XiModel>
-    (
-        ctorPtr(dict, thermo, turbulence, Su, rho, b, phi)
-    );
+        (cstrIter()(dict, thermo, turbulence, Su, rho, b, phi));
 }


--- a/applications/solvers/combustion/fireFoam/Make/options
+++ b/applications/solvers/combustion/fireFoam/Make/options
@ -1,6 +1,5 @@
 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 \
@ -15,6 +14,7 @@ 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,9 +23,7 @@ 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)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/ODE/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -52,7 +50,4 @@ EXE_LIBS = \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
-    -lODE \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lODE
--- a/applications/solvers/compressible/rhoPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/Make/options
@ -8,8 +8,7 @@ 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)/regionFaModels\lnInclude
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -24,5 +23,4 @@ EXE_LIBS = \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
    -lsampling \
-    -latmosphericModels \
-    -lregionFaModels
+    -latmosphericModels
--- a/applications/solvers/compressible/rhoSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/Make/options
@ -7,7 +7,6 @@ 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 \
@ -19,5 +18,4 @@ EXE_LIBS = \
    -lspecie \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -latmosphericModels \
-    -lregionFaModels
+    -latmosphericModels
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
@ -6,8 +6,7 @@ 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)/regionFaModels\lnInclude
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -20,5 +19,4 @@ EXE_LIBS = \
    -lspecie \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -latmosphericModels \
-    -lregionFaModels
+    -latmosphericModels
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options
@ -6,8 +6,7 @@ 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)/regionFaModels\lnInclude
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -20,5 +19,4 @@ EXE_LIBS = \
    -lradiationModels \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -latmosphericModels \
-    -lregionFaModels
+    -latmosphericModels
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
@ -18,9 +18,7 @@ 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)/regionFaModels\lnInclude
-
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -38,5 +36,4 @@ EXE_LIBS = \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
    -lradiationModels \
-    -lregionModels \
-    -lregionFaModels
+    -lregionModels
--- 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-2019 OpenCFD Ltd.
+    Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -54,7 +54,6 @@ Description
 #include "loopControl.H"
 #include "pressureControl.H"

-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

-    #include "createCoupledRegions.H"
-
    while (runTime.run())
    {
        #include "readTimeControls.H"
@ -112,6 +109,8 @@ 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"
@ -119,41 +118,13 @@ 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)
            {
@ -181,21 +152,6 @@ 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,7 +3,6 @@ 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 \
@ -36,4 +35,3 @@ 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-2019 OpenCFD Ltd.
+    Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -66,7 +66,6 @@ int main(int argc, char *argv[])
    #include "createTime.H"
    #include "createMeshes.H"
    #include "createFields.H"
-    #include "createCoupledRegions.H"
    #include "initContinuityErrs.H"

    while (runTime.loop())
@ -84,32 +83,13 @@ 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");
@ -136,15 +116,6 @@ 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,20 +20,13 @@

    fvOptions.constrain(EEqn);

-    if (coupled)
-    {
-        fvMatrixAssemblyPtr->addFvMatrix(EEqn);
-    }
-    else
-    {
-        EEqn.solve();
+    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);

-    UEqn =
+    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(rho, U)
@ -10,6 +10,7 @@
     ==
        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,8 +21,6 @@ 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
@ -224,11 +222,5 @@ 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,8 +22,6 @@
    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,10 +12,8 @@

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

        fvOptions.constrain(hEqn);

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

-            hEqn.solve();
-            fvOptions.correct(h);
-            thermo.correct();
-
-            Info<< "Min/max T:" << min(thermo.T()).value() << ' '
-                << max(thermo.T()).value() << endl;
-        }
+        fvOptions.correct(h);
    }
+
+    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,7 +7,6 @@ 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
@ -1,39 +0,0 @@
-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,19 +27,13 @@

    fvOptions.constrain(EEqn);

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

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

-        Info<< "Min/max T:" << min(thermo.T()).value() << ' '
-            << max(thermo.T()).value() << endl;
-    }
+    thermo.correct();
+    rad.correct();
+
+    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);

-    UEqn =
+    tmp<fvVectorMatrix> tUEqn
    (
        fvm::ddt(rho, U) + fvm::div(phi, U)
      + MRF.DDt(rho, U)
@ -10,6 +10,7 @@
     ==
        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,8 +16,6 @@ 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);

@ -295,12 +293,6 @@ 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,8 +49,6 @@
    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,20 +18,15 @@ else
    #include "YEqn.H"
    #include "EEqn.H"

-    if (!coupled)
+    // --- PISO loop
+    for (int corr=0; corr<nCorr; corr++)
    {
-        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();
+        #include "pEqn.H"
    }
+
+    turbulence.correct();
+
+    rho = thermo.rho();
 }

 if (finalIter)
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/include/correctThermos.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/include/correctThermos.H
@ -1,20 +0,0 @@
-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
@ -1,69 +0,0 @@
-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;
-    }
-}
-
-if (coupled && !fvMatrixAssemblyPtr)
-{
-    forAll(solidRegions, i)
-    {
-        solidThermo& thermo = thermos[i];
-        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,45 +1,39 @@
+if (finalIter)
 {
-    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)
-    );
+    mesh.data::add("finalIteration", true);
+}

-    hEqn.relax();
-
-    fvOptions.constrain(hEqn);
-
-    if (coupled)
+{
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; ++nonOrth)
    {
-        fvMatrixAssemblyPtr->addFvMatrix(hEqn);
-    }
-    else
-    {
-        Info<< "\nSolving for solid region "<< solidRegions[i].name() << endl;
+        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)
+        );

-        if (finalIter)
-        {
-            mesh.data::add("finalIteration", true);
-        }
+        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");
-        }
    }
+
+    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,8 +7,7 @@ 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)/regionFaModels\lnInclude
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -21,5 +20,4 @@ EXE_LIBS = \
    -ldynamicMesh \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
-    -latmosphericModels \
-    -lregionFaModels
+    -latmosphericModels
--- a/applications/solvers/incompressible/pimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pEqn.H
@ -6,7 +6,10 @@ 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/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options
@ -3,7 +3,6 @@ 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 \
@ -20,9 +19,7 @@ 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)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/dynamicMesh/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -42,7 +39,4 @@ EXE_LIBS = \
    -ldynamicMesh \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
-    -lsampling \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lsampling
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/options
@ -2,7 +2,6 @@ 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,9 +18,7 @@ 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)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -41,7 +38,4 @@ EXE_LIBS = \
    -ldynamicMesh \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
-    -lsampling \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lsampling
--- a/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/DPMIncompressibleTurbulenceModel.C
+++ b/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/DPMIncompressibleTurbulenceModel.C
@ -149,17 +149,6 @@ 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
@ -181,19 +170,6 @@ 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,24 +116,12 @@ 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,7 +2,6 @@ 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,8 +17,6 @@ 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 \
@ -36,7 +33,4 @@ EXE_LIBS = \
    -lDPMTurbulenceModels \
    -lregionModels \
    -lsurfaceFilmModels \
-    -lsampling \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lsampling
--- a/applications/solvers/lagrangian/DPMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/Make/options
@ -1,7 +1,6 @@
 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 \
@ -15,9 +14,7 @@ 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)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -34,7 +31,4 @@ EXE_LIBS = \
    -lDPMTurbulenceModels \
    -lregionModels \
    -lsurfaceFilmModels \
-    -lsampling \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lsampling
--- a/applications/solvers/lagrangian/coalChemistryFoam/Make/options
+++ b/applications/solvers/lagrangian/coalChemistryFoam/Make/options
@ -1,6 +1,5 @@
 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 \
@ -12,14 +11,13 @@ 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 \
@ -47,7 +45,4 @@ EXE_LIBS = \
    -lsurfaceFilmModels \
    -lODE \
    -lcombustionModels \
-    -lsampling \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lsampling
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Make/options
@ -1,6 +1,5 @@
 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 \
@ -14,9 +13,7 @@ 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)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -33,7 +30,4 @@ EXE_LIBS = \
    -lincompressibleTurbulenceModels \
    -lincompressibleTransportModels \
    -lregionModels \
-    -lsurfaceFilmModels \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lsurfaceFilmModels
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options
@ -1,7 +1,6 @@
 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 \
@ -17,9 +16,7 @@ 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)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -39,7 +36,4 @@ EXE_LIBS = \
    -lsurfaceFilmModels \
    -ldynamicMesh \
    -ldynamicFvMesh \
-    -ltopoChangerFvMesh \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -ltopoChangerFvMesh
--- a/applications/solvers/lagrangian/kinematicParcelFoam/Make/files
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/Make/files
@ -1,3 +0,0 @@
-kinematicParcelFoam.C
-
-EXE = $(FOAM_APPBIN)/kinematicParcelFoam
--- a/applications/solvers/lagrangian/kinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/Make/options
@ -1,42 +0,0 @@
-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
@ -1,22 +0,0 @@
-    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
@ -1,16 +0,0 @@
-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
@ -1 +0,0 @@
-regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
--- a/applications/solvers/lagrangian/kinematicParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/createFields.H
@ -1,85 +0,0 @@
-#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
@ -1,155 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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
@ -1,58 +0,0 @@
-
-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,7 +1,6 @@
 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 \
@ -13,14 +12,13 @@ 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 \
@ -51,7 +49,4 @@ EXE_LIBS = \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
    -lODE \
-    -lcombustionModels \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lcombustionModels
--- a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
@ -11,11 +11,7 @@ basicSpecieMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();

 const word inertSpecie(thermo.get<word>("inertSpecie"));
-if
-(
-    !composition.species().found(inertSpecie)
-  && composition.species().size() > 0
-)
+if (!composition.species().found(inertSpecie))
 {
    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,7 +1,6 @@
 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 \
@ -13,14 +12,13 @@ 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 \
@ -50,7 +48,4 @@ EXE_LIBS = \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
    -lODE \
-    -lcombustionModels \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lcombustionModels
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
@ -1,6 +1,5 @@
 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 \
@ -13,6 +12,7 @@ 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,8 +20,6 @@ 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

@ -46,7 +44,4 @@ EXE_LIBS = \
    -lregionModels \
    -lsurfaceFilmModels \
    -lcombustionModels \
-    -lsampling \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lsampling
--- a/applications/solvers/lagrangian/simpleCoalParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/simpleCoalParcelFoam/Make/options
@ -1,6 +1,5 @@
 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 \
@ -13,6 +12,7 @@ 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,9 +21,7 @@ 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)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -46,7 +44,4 @@ EXE_LIBS = \
    -lsurfaceFilmModels \
    -lcombustionModels \
    -lsampling \
-    -lcoalCombustion \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lcoalCombustion
--- a/applications/solvers/lagrangian/sprayFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/Make/options
@ -1,7 +1,6 @@
 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 \
@ -14,6 +13,7 @@ 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,8 +21,6 @@ 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 = \
@ -47,7 +45,4 @@ EXE_LIBS = \
    -lregionModels \
    -lsurfaceFilmModels \
    -lcombustionModels \
-    -lsampling \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lsampling
--- a/applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options
@ -3,7 +3,6 @@ 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 \
@ -16,6 +15,7 @@ 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,9 +24,7 @@ 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)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/combustionModels/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -51,7 +49,4 @@ EXE_LIBS = \
    -lengine \
    -lregionModels \
    -lsurfaceFilmModels \
-    -lcombustionModels \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lcombustionModels
--- a/applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/Make/options
@ -1,7 +1,6 @@
 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 \
@ -14,6 +13,7 @@ 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,9 +21,7 @@ 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)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/combustionModels/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -47,7 +45,4 @@ EXE_LIBS = \
    -lregionModels \
    -lsurfaceFilmModels \
    -lcombustionModels \
-    -lsampling \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lsampling
--- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/Make/options
@ -3,7 +3,6 @@ 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 \
@ -16,6 +15,7 @@ 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,9 +26,7 @@ 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)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -54,7 +52,4 @@ EXE_LIBS = \
    -lcombustionModels \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
-    -ldynamicMesh \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -ldynamicMesh
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/Make/options
@ -1,6 +1,5 @@
 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 \
@ -12,9 +11,7 @@ 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)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -30,7 +27,4 @@ EXE_LIBS = \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
    -lregionModels \
-    -lsurfaceFilmModels \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lsurfaceFilmModels
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
@ -1,7 +1,6 @@
 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 \
@ -15,9 +14,7 @@ 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)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -36,7 +33,4 @@ EXE_LIBS = \
    -lsurfaceFilmModels \
    -ldynamicMesh \
    -ldynamicFvMesh \
-    -ltopoChangerFvMesh \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -ltopoChangerFvMesh
--- a/applications/solvers/multiphase/MPPICInterFoam/Make/options
+++ b/applications/solvers/multiphase/MPPICInterFoam/Make/options
@ -2,7 +2,6 @@ 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 \
@ -21,9 +20,7 @@ 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)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude

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

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

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

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

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

            MULES::explicitSolve
            (
--- a/applications/solvers/multiphase/VoF/createAlphaFluxes.H
+++ b/applications/solvers/multiphase/VoF/createAlphaFluxes.H
@ -24,18 +24,3 @@ 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,6 +14,7 @@ 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-2021 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -122,7 +122,7 @@ Foam::fv::VoFSolidificationMeltingSource::VoFSolidificationMeltingSource
    const fvMesh& mesh
 )
 :
-    fv::cellSetOption(sourceName, modelType, dict, mesh),
+    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 fv::cellSetOption
+    public 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 (fv::cellSetOption::read(dict))
+    if (cellSetOption::read(dict))
    {
        alphaSolidT_ = Function1<scalar>::New("alphaSolidT", coeffs_);
        coeffs_.readEntry("L", L_);
--- a/applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/BinghamPlastic/BinghamPlastic.C
+++ b/applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/BinghamPlastic/BinghamPlastic.C
@ -6,7 +6,6 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2014-2015 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -131,9 +130,9 @@ bool Foam::mixtureViscosityModels::BinghamPlastic::read
 {
    plastic::read(viscosityProperties);

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

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

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

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

--- a/applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/mixtureViscosityModel/mixtureViscosityModelNew.C
+++ b/applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/mixtureViscosityModel/mixtureViscosityModelNew.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2014-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -44,9 +44,9 @@ Foam::autoPtr<Foam::mixtureViscosityModel> Foam::mixtureViscosityModel::New

    Info<< "Selecting incompressible transport model " << modelType << endl;

-    auto* ctorPtr = dictionaryConstructorTable(modelType);
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);

-    if (!ctorPtr)
+    if (!cstrIter.found())
    {
        FatalIOErrorInLookup
        (
@ -58,9 +58,7 @@ Foam::autoPtr<Foam::mixtureViscosityModel> Foam::mixtureViscosityModel::New
    }

    return autoPtr<mixtureViscosityModel>
-    (
-        ctorPtr(name, dict, U, phi)
-    );
+        (cstrIter()(name, dict, U, phi));
 }


--- a/applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/relativeVelocityModel/relativeVelocityModel.C
+++ b/applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/relativeVelocityModel/relativeVelocityModel.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2014-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -111,9 +111,9 @@ Foam::autoPtr<Foam::relativeVelocityModel> Foam::relativeVelocityModel::New

    Info<< "Selecting relative velocity model " << modelType << endl;

-    auto* ctorPtr = dictionaryConstructorTable(modelType);
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);

-    if (!ctorPtr)
+    if (!cstrIter.found())
    {
        FatalIOErrorInLookup
        (
@ -127,7 +127,7 @@ Foam::autoPtr<Foam::relativeVelocityModel> Foam::relativeVelocityModel::New
    return
        autoPtr<relativeVelocityModel>
        (
-            ctorPtr
+            cstrIter()
            (
                dict.optionalSubDict(modelType + "Coeffs"),
                mixture
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/reflectionModel/reflectionModel/reflectionModelNew.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/reflectionModel/reflectionModel/reflectionModelNew.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -41,9 +41,9 @@ Foam::radiation::reflectionModel::New

    Info<< "Selecting reflectionModel " << modelType << endl;

-    auto* ctorPtr = dictionaryConstructorTable(modelType);
+    const auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);

-    if (!ctorPtr)
+    if (!cstrIter.found())
    {
        FatalIOErrorInLookup
        (
@ -54,7 +54,7 @@ Foam::radiation::reflectionModel::New
        ) << exit(FatalIOError);
    }

-    return autoPtr<reflectionModel>(ctorPtr(dict, mesh));
+    return autoPtr<reflectionModel>(cstrIter()(dict, mesh));
 }


--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H
@ -120,6 +120,21 @@ 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,6 +1,5 @@
 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 \
@ -10,7 +9,6 @@ EXE_INC = \

 LIB_LIBS = \
    -lfiniteVolume \
-    -lmeshTools \
    -lgeometricVoF \
    -ltwoPhaseMixture \
    -linterfaceProperties \
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.C
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
    Copyright (C) 2020 Henning Scheufler
 -------------------------------------------------------------------------------
 License
@ -269,14 +269,14 @@ correct()

    // interface heat resistance
    mDotc_ = interfaceArea_*R_*max(TSat - T, T0)/L;
-    mDote_ = interfaceArea_*R_*max(T - TSat, T0)/L;
+    mDote_ = interfaceArea_*R_*(T - TSat)/L;

-    // Limiting max condensation
    forAll(mDotc_, celli)
    {
        scalar rhobyDt = mixture_.rho1().value()/mesh_.time().deltaTValue();
        scalar maxEvap = mixture_.alpha1()[celli]*rhobyDt; // positive
        scalar maxCond = -mixture_.alpha2()[celli]*rhobyDt; // negative
+        mDote_[celli] = min(max(mDote_[celli], maxCond), maxEvap);
        mDotc_[celli] = min(max(mDotc_[celli], maxCond), maxEvap);
    }

@ -326,6 +326,14 @@ correct()
 void Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
 updateInterface()
 {
+    const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
+    const twoPhaseMixtureEThermo& thermo =
+    refCast<const twoPhaseMixtureEThermo>
+    (
+        mesh_.lookupObject<basicThermo>(basicThermo::dictName)
+    );
+
+    const dimensionedScalar& TSat = thermo.TSat();

    // interface heat resistance
    // Interpolating alpha1 cell centre values to mesh points (vertices)
@ -346,6 +354,30 @@ updateInterface()
                mag(cutCell.faceArea())/mesh_.V()[celli];
        }
    }
+
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
+    forAll(pbm, patchi)
+    {
+        if (isA<wallPolyPatch>(pbm[patchi]))
+        {
+            const polyPatch& pp = pbm[patchi];
+            forAll(pp.faceCells(),i)
+            {
+                const label pCelli = pp.faceCells()[i];
+
+                if
+                (
+                    (TSat.value() - T[pCelli]) > 0
+                  && mixture_.alpha1()[pCelli] < 0.9
+                )
+                {
+                    interfaceArea_[pCelli] =
+                        mag(pp.faceAreas()[i])/mesh_.V()[pCelli];
+                }
+            }
+        }
+    }
 }


--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixtureNew.C
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixtureNew.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -54,9 +54,9 @@ Foam::temperaturePhaseChangeTwoPhaseMixture::New

    Info<< "Selecting phaseChange model " << modelType << endl;

-    auto* ctorPtr = componentsConstructorTable(modelType);
+    auto cstrIter = componentsConstructorTablePtr_->cfind(modelType);

-    if (!ctorPtr)
+    if (!cstrIter.found())
    {
        FatalIOErrorInLookup
        (
@ -70,7 +70,7 @@ Foam::temperaturePhaseChangeTwoPhaseMixture::New
    return
        autoPtr<temperaturePhaseChangeTwoPhaseMixture>
        (
-            ctorPtr(thermo, mesh)
+            cstrIter()(thermo, mesh)
        );
 }

--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@ -126,5 +126,20 @@ 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/createFields.H
+++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/createFields.H
@ -147,4 +147,19 @@ 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/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/options
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/options
@ -1,7 +1,6 @@
 EXE_INC = \
    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
    -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude

@ -9,5 +8,4 @@ LIB_LIBS = \
    -ltwoPhaseMixture \
    -ltwoPhaseProperties \
    -lincompressibleTransportModels \
-    -lfiniteVolume \
-    -lmeshTools
+    -lfiniteVolume
--- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
@ -26,7 +26,6 @@ License
 \*---------------------------------------------------------------------------*/

 #include "phaseChangeTwoPhaseMixture.H"
-#include "surfaceMesh.H"

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

--- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixtureNew.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixtureNew.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -28,7 +28,6 @@ License

 #include "phaseChangeTwoPhaseMixture.H"
 #include "incompressibleTwoPhaseMixture.H"
-#include "surfaceMesh.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -56,9 +55,9 @@ Foam::phaseChangeTwoPhaseMixture::New

    Info<< "Selecting phaseChange model " << modelType << endl;

-    auto* ctorPtr = componentsConstructorTable(modelType);
+    auto cstrIter = componentsConstructorTablePtr_->cfind(modelType);

-    if (!ctorPtr)
+    if (!cstrIter.found())
    {
        FatalIOErrorInLookup
        (
@ -69,7 +68,7 @@ Foam::phaseChangeTwoPhaseMixture::New
        ) << exit(FatalIOError);
    }

-    return autoPtr<phaseChangeTwoPhaseMixture>(ctorPtr(U, phi));
+    return autoPtr<phaseChangeTwoPhaseMixture>(cstrIter()(U, phi));
 }


--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@ -6,7 +6,6 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -27,6 +26,7 @@ License
 \*---------------------------------------------------------------------------*/

 #include "multiphaseMixture.H"
+#include "alphaContactAngleFvPatchScalarField.H"
 #include "Time.H"
 #include "subCycle.H"
 #include "MULES.H"
@ -36,7 +36,6 @@ License
 #include "fvcDiv.H"
 #include "fvcFlux.H"
 #include "unitConversion.H"
-#include "alphaContactAngleFvPatchScalarField.H"

 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

@ -417,123 +416,102 @@ void Foam::multiphaseMixture::correctContactAngle
    surfaceVectorField::Boundary& nHatb
 ) const
 {
-    const volScalarField::Boundary& gb1f = alpha1.boundaryField();
-    const volScalarField::Boundary& gb2f = alpha2.boundaryField();
+    const volScalarField::Boundary& gbf
+        = alpha1.boundaryField();

    const fvBoundaryMesh& boundary = mesh_.boundary();

    forAll(boundary, patchi)
    {
-        if (isA<alphaContactAngleFvPatchScalarField>(gb1f[patchi]))
+        if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
        {
            const alphaContactAngleFvPatchScalarField& acap =
-                refCast<const alphaContactAngleFvPatchScalarField>(gb1f[patchi]);
+                refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);

-            correctBoundaryContactAngle(acap, patchi, alpha1, alpha2, nHatb);
-        }
-        else if (isA<alphaContactAngleFvPatchScalarField>(gb2f[patchi]))
-        {
-            const alphaContactAngleFvPatchScalarField& acap =
-                refCast<const alphaContactAngleFvPatchScalarField>(gb2f[patchi]);
+            vectorField& nHatPatch = nHatb[patchi];

-            correctBoundaryContactAngle(acap, patchi, alpha2, alpha1, nHatb);
+            vectorField AfHatPatch
+            (
+                mesh_.Sf().boundaryField()[patchi]
+               /mesh_.magSf().boundaryField()[patchi]
+            );
+
+            const auto tp =
+                acap.thetaProps().cfind(interfacePair(alpha1, alpha2));
+
+            if (!tp.found())
+            {
+                FatalErrorInFunction
+                    << "Cannot find interface " << interfacePair(alpha1, alpha2)
+                    << "\n    in table of theta properties for patch "
+                    << acap.patch().name()
+                    << exit(FatalError);
+            }
+
+            const bool matched = (tp.key().first() == alpha1.name());
+
+            const scalar theta0 = degToRad(tp().theta0(matched));
+            scalarField theta(boundary[patchi].size(), theta0);
+
+            const scalar uTheta = tp().uTheta();
+
+            // Calculate the dynamic contact angle if required
+            if (uTheta > SMALL)
+            {
+                const scalar thetaA = degToRad(tp().thetaA(matched));
+                const scalar thetaR = degToRad(tp().thetaR(matched));
+
+                // Calculated the component of the velocity parallel to the wall
+                vectorField Uwall
+                (
+                    U_.boundaryField()[patchi].patchInternalField()
+                  - U_.boundaryField()[patchi]
+                );
+                Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
+
+                // Find the direction of the interface parallel to the wall
+                vectorField nWall
+                (
+                    nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
+                );
+
+                // Normalise nWall
+                nWall /= (mag(nWall) + SMALL);
+
+                // Calculate Uwall resolved normal to the interface parallel to
+                // the interface
+                scalarField uwall(nWall & Uwall);
+
+                theta += (thetaA - thetaR)*tanh(uwall/uTheta);
+            }
+
+
+            // Reset nHatPatch to correspond to the contact angle
+
+            scalarField a12(nHatPatch & AfHatPatch);
+
+            scalarField b1(cos(theta));
+
+            scalarField b2(nHatPatch.size());
+
+            forAll(b2, facei)
+            {
+                b2[facei] = cos(acos(a12[facei]) - theta[facei]);
+            }
+
+            scalarField det(1.0 - a12*a12);
+
+            scalarField a((b1 - a12*b2)/det);
+            scalarField b((b2 - a12*b1)/det);
+
+            nHatPatch = a*AfHatPatch + b*nHatPatch;
+
+            nHatPatch /= (mag(nHatPatch) + deltaN_.value());
        }
    }
 }


-void Foam::multiphaseMixture::correctBoundaryContactAngle
-(
-    const alphaContactAngleFvPatchScalarField& acap,
-    label patchi,
-    const phase& alpha1,
-    const phase& alpha2,
-    surfaceVectorField::Boundary& nHatb
-) const
-{
-    const fvBoundaryMesh& boundary = mesh_.boundary();
-
-    vectorField& nHatPatch = nHatb[patchi];
-
-    vectorField AfHatPatch
-    (
-        mesh_.Sf().boundaryField()[patchi]
-       /mesh_.magSf().boundaryField()[patchi]
-    );
-
-    const auto tp = acap.thetaProps().cfind(interfacePair(alpha1, alpha2));
-
-    if (!tp.found())
-    {
-        FatalErrorInFunction
-            << "Cannot find interface " << interfacePair(alpha1, alpha2)
-            << "\n    in table of theta properties for patch "
-            << acap.patch().name()
-            << exit(FatalError);
-    }
-
-    const bool matched = (tp.key().first() == alpha1.name());
-
-    const scalar theta0 = degToRad(tp().theta0(matched));
-    scalarField theta(boundary[patchi].size(), theta0);
-
-    const scalar uTheta = tp().uTheta();
-
-    // Calculate the dynamic contact angle if required
-    if (uTheta > SMALL)
-    {
-        const scalar thetaA = degToRad(tp().thetaA(matched));
-        const scalar thetaR = degToRad(tp().thetaR(matched));
-
-        // Calculated the component of the velocity parallel to the wall
-        vectorField Uwall
-        (
-            U_.boundaryField()[patchi].patchInternalField()
-          - U_.boundaryField()[patchi]
-        );
-        Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
-
-        // Find the direction of the interface parallel to the wall
-        vectorField nWall
-        (
-            nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
-        );
-
-        // Normalise nWall
-        nWall /= (mag(nWall) + SMALL);
-
-        // Calculate Uwall resolved normal to the interface parallel to
-        // the interface
-        scalarField uwall(nWall & Uwall);
-
-        theta += (thetaA - thetaR)*tanh(uwall/uTheta);
-    }
-
-
-    // Reset nHatPatch to correspond to the contact angle
-
-    scalarField a12(nHatPatch & AfHatPatch);
-
-    scalarField b1(cos(theta));
-
-    scalarField b2(nHatPatch.size());
-
-    forAll(b2, facei)
-    {
-        b2[facei] = cos(acos(a12[facei]) - theta[facei]);
-    }
-
-    scalarField det(1.0 - a12*a12);
-
-    scalarField a((b1 - a12*b2)/det);
-    scalarField b((b2 - a12*b1)/det);
-
-    nHatPatch = a*AfHatPatch + b*nHatPatch;
-
-    nHatPatch /= (mag(nHatPatch) + deltaN_.value());
-}
-
-
 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
 (
    const phase& alpha1,
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H
@ -57,9 +57,6 @@ SourceFiles
 namespace Foam
 {

-// Class forward declarations
-class alphaContactAngleFvPatchScalarField;
-
 /*---------------------------------------------------------------------------*\
                      Class multiphaseMixture Declaration
 \*---------------------------------------------------------------------------*/
@ -173,15 +170,6 @@ private:
            surfaceVectorField::Boundary& nHatb
        ) const;

-        void correctBoundaryContactAngle
-        (
-            const alphaContactAngleFvPatchScalarField& acap,
-            label patchi,
-            const phase& alpha1,
-            const phase& alpha2,
-            surfaceVectorField::Boundary& nHatb
-        ) const;
-
        tmp<volScalarField> K(const phase& alpha1, const phase& alpha2) const;


--- a/applications/test/DynamicList2/Make/files
+++ b/applications/test/DynamicList2/Make/files
@ -1,3 +0,0 @@
-Test-DynamicList2.C
-
-EXE = $(FOAM_USER_APPBIN)/Test-DynamicList2
--- a/applications/test/DynamicList2/Make/options
+++ b/applications/test/DynamicList2/Make/options
@ -1,2 +0,0 @@
-/* EXE_INC = */
-/* EXE_LIBS = */
--- a/applications/test/DynamicList2/Test-DynamicList2.C
+++ b/applications/test/DynamicList2/Test-DynamicList2.C
@ -1,155 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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/>.
-
-Description
-    Test allocation patterns when reading into an existing list.
-
-\*---------------------------------------------------------------------------*/
-
-#include "DynamicList.H"
-#include "DynamicField.H"
-#include "IOstreams.H"
-#include "ITstream.H"
-#include "OTstream.H"
-#include "StringStream.H"
-#include "FlatOutput.H"
-#include "ListOps.H"
-#include "labelRange.H"
-#include "labelIndList.H"
-
-using namespace Foam;
-
-template<class T, int SizeMin>
-void printInfo
-(
-    const word& tag,
-    const DynamicList<T, SizeMin>& list,
-    const bool showSize = true
-)
-{
-    Info<< '<' << tag;
-    if (showSize)
-    {
-        Info<< " size=\"" << list.size()
-            << "\" capacity=\"" << list.capacity() << "\"";
-        if (list.cdata())
-        {
-            Info<< " ptr=\"" << name(list.cdata()) << "\"";
-        }
-        else
-        {
-            Info<< " ptr=\"nullptr\"";
-        }
-    }
-    Info<< '>' << nl << flatOutput(list) << nl
-        << "</" << tag << ">\n" << endl;
-}
-
-
-template<class T, int SizeMin>
-void printInfo
-(
-    const word& tag,
-    const DynamicField<T, SizeMin>& list,
-    const bool showSize = true
-)
-{
-    Info<< '<' << tag;
-    if (showSize)
-    {
-        Info<< " size=\"" << list.size()
-            << "\" capacity=\"" << list.capacity() << "\"";
-        if (list.cdata())
-        {
-            Info<< " ptr=\"" << name(list.cdata()) << "\"";
-        }
-        else
-        {
-            Info<< " ptr=\"nullptr\"";
-        }
-    }
-    Info<< '>' << nl << flatOutput(list) << nl
-        << "</" << tag << ">\n" << endl;
-}
-
-
-template<class T, int SizeMin>
-void readList
-(
-    DynamicList<T, SizeMin>& output,
-    const UList<T>& input
-)
-{
-    OTstream os;
-    os << input;
-    ITstream is("input", os.tokens());
-
-    is >> output;
-}
-
-template<class T, int SizeMin>
-void readList
-(
-    DynamicField<T, SizeMin>& output,
-    const UList<T>& input
-)
-{
-    OTstream os;
-    os << input;
-    ITstream is("input", os.tokens());
-
-    is >> output;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-// Main program:
-
-int main(int argc, char *argv[])
-{
-    //
-    {
-        DynamicList<label, 64> list1;
-
-        list1.resize(4);
-        ListOps::identity(list1);
-
-        list1.resize(3);
-        printInfo("", list1);
-
-        // list1.clear();
-        // printInfo("", list1);
-
-        list1.setCapacity(3);
-        printInfo("", list1);
-    }
-
-    Info<< "\nEnd\n";
-
-    return 0;
-}
-
-
-// ************************************************************************* //
--- a/applications/test/FixedList/Test-FixedList.C
+++ b/applications/test/FixedList/Test-FixedList.C
@ -280,24 +280,16 @@ int main(int argc, char *argv[])


    List<label> list3{0, 1, 2, 3};
-    FixedList<label, 4> list4(list3);
+    FixedList<label, 4> list4(list3.begin(), list3.end());
    Info<< "list3: " << list3 << nl
        << "list4: " << list4 << nl;

-    list4 = {1, 20, 3, 40};
+    list4 = {1, 2, 3, 5};
    Info<< "list4: " << list4 << nl;

    FixedList<label, 5> list5{0, 1, 2, 3, 4};
    Info<< "list5: " << list5 << nl;

-    {
-        const FixedList<label, 2> indices({3, 1});
-        FixedList<label, 2> list4b(list4, indices);
-
-        Info<< "subset " << list4 << " with " << indices << " -> "
-            << list4b << nl;
-    }
-
    List<FixedList<label, 2>> list6{{0, 1}, {2, 3}};
    Info<< "list6: " << list6 << nl;

--- a/applications/test/HashTable1/Test-HashTable1.C
+++ b/applications/test/HashTable1/Test-HashTable1.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -28,12 +28,11 @@ License

 #include "HashTable.H"
 #include "List.H"
+#include "SortableList.H"
 #include "DynamicList.H"
 #include "FlatOutput.H"
 #include "IOstreams.H"
 #include "StringStream.H"
-#include "ListOps.H"
-#include "flipOp.H"

 using namespace Foam;

@ -210,16 +209,12 @@ int main()
            << "/" << table1.keys().size()
            << nl;

-        List<word> sortKeys
+        SortableList<word> sortKeys
+        // DynamicList<word> sortKeys
        (
-            ListOps::create<word>
-            (
-                table1.keys().begin(),
-                table1.keys().end(),
-                noOp{}
-            )
+            table1.keys().begin(),
+            table1.keys().end()
        );
-        sort(sortKeys);
        Info<<"sortKeys: " << flatOutput(sortKeys) << nl;
    }

--- a/applications/test/ITstream/Test-ITstream.C
+++ b/applications/test/ITstream/Test-ITstream.C
@ -33,8 +33,6 @@ Description
 #include "IOstreams.H"
 #include "argList.H"
 #include "ITstream.H"
-#include "ListOps.H"
-#include "flipOp.H"

 using namespace Foam;

@ -111,7 +109,7 @@ void doTest
        << nl
        << "====" << nl << endl;

-    ITstream its(input);
+    ITstream its(name, input);
    Info<< "got " << its.size() << " tokens - index at "
        << its.tokenIndex() << endl;

@ -166,15 +164,7 @@ int main(int argc, char *argv[])

    string stringInput("( string     ;    input \"string\" to tokenize )");

-    List<char> listInput
-    (
-        ListOps::create<char>
-        (
-            stringInput.cbegin(),
-            stringInput.cend(),
-            Foam::noOp{}
-        )
-    );
+    List<char> listInput(stringInput.cbegin(), stringInput.cend());

    doTest("empty", "", true, true);

--- a/applications/test/List/Test-List.C
+++ b/applications/test/List/Test-List.C
@ -48,9 +48,7 @@ See also
 #include "scalarList.H"
 #include "HashOps.H"
 #include "ListOps.H"
-#include "IndirectList.H"
 #include "SubList.H"
-#include "SliceList.H"
 #include "ListPolicy.H"

 #include <list>
@ -346,7 +344,7 @@ int main(int argc, char *argv[])
        labelList longLabelList = identity(25);
        reverse(longLabelList);

-        FixedList<label, 6> fixedLabelList({0,1,2,3,4,5});
+        FixedList<label, 6> fixedLabelList{0,1,2,3,4,5};
        const labelList constLabelList = identity(25);

        Info<< "full-list: " << flatOutput(longLabelList) << nl;
@ -355,36 +353,6 @@ int main(int argc, char *argv[])
        Info<<"sub range:" << range1 << "=";
        Info<< SubList<label>(longLabelList, range1) << nl;

-        {
-            // A valid range
-            const labelRange subset(4, 5);
-
-            // Assign some values
-            longLabelList.slice(subset) = identity(subset.size());
-
-            Info<<"assigned identity in range:" << subset
-                << "=> " << flatOutput(longLabelList) << nl;
-
-            labelList someList(identity(24));
-
-            longLabelList.slice(subset) =
-                SliceList<label>(someList, sliceRange(8, subset.size(), 2));
-
-            Info<<"assigned sliced/stride in range:" << subset
-                << "=> " << flatOutput(longLabelList) << nl;
-
-            // Does not work - need a reference, not a temporary
-            // Foam::reverse(longLabelList[subset]);
-
-            {
-                auto sub(longLabelList.slice(subset));
-                Foam::reverse(sub);
-            }
-
-            Info<<"reversed range:" << subset
-                << "=> " << flatOutput(longLabelList) << nl;
-        }
-
        labelRange range2(7, 8);
        Info<<"sub range:" << range2 << "=";
        Info<< SubList<label>(longLabelList, range2) << nl;
@ -400,29 +368,29 @@ int main(int argc, char *argv[])
        // > constLabelList[labelRange(23,5)] = 5;

        // Check correct overlaps
-        longLabelList.slice(labelRange(-10, 12)) = 200;
-        longLabelList.slice({18,3}) = 100;
-        longLabelList.slice({23,3}) = 400;
+        longLabelList[labelRange(-10, 12)] = 200;
+        longLabelList[{18,3}] = 100;
+        longLabelList[{23,3}] = 400;
        // and complete misses
-        longLabelList.slice({500,50}) = 100;
+        longLabelList[{500,50}] = 100;

        // -ve size suppressed by internal 'validateRange' = no-op
-        longLabelList.slice({5,-5}) = 42;
-        longLabelList.slice({21,100}) = 42;
+        longLabelList[{5,-5}] = 42;
+        longLabelList[{21,100}] = 42;

        //Good: does not compile
-        longLabelList.slice(labelRange(20,50)) = constLabelList;
+        //> longLabelList[labelRange(20,50)] = constLabelList;

        //Good: does not compile
        // longLabelList[labelRange(20,50)] = fixedLabelList;

-        Info<< "updated: " << constLabelList.slice(labelRange(23,5)) << nl;
+        Info<< "updated: " << constLabelList[labelRange(23,5)] << nl;
        Info<< "updated: " << flatOutput(longLabelList) << nl;

-        //Nope: sort(longLabelList.slice(labelRange(18,5)));
+        //Nope: sort(longLabelList[labelRange(18,5)]);
        {
            // Instead
-            auto sub = longLabelList.slice(labelRange(8));
+            UList<label> sub = longLabelList[labelRange(8)];
            sort(sub);
        }
        Info<< "sub-sorted: " << flatOutput(longLabelList) << nl;
--- a/applications/test/PDRblockMesh/Test-PDRblockMesh.C
+++ b/applications/test/PDRblockMesh/Test-PDRblockMesh.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -180,7 +180,7 @@ int main(int argc, char *argv[])

    // Fatal with FULLDEBUG
    {
-        const bool oldThrowingError = FatalError.throwing(true);
+        const bool throwingError = FatalError.throwExceptions();

        const label nx = mesh.sizes().x();
        const label ny = mesh.sizes().y();
@ -230,7 +230,7 @@ int main(int argc, char *argv[])
                << err.message().c_str() << nl;
        }

-        FatalError.throwing(oldThrowingError);
+        FatalError.throwExceptions(throwingError);
    }


--- a/applications/test/PackedList1/Test-PackedList1.C
+++ b/applications/test/PackedList1/Test-PackedList1.C
@ -6,7 +6,6 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -172,7 +171,7 @@ int main(int argc, char *argv[])
        Info<< "\ntest operator[] non-const with out-of-range index\n";

        // Expect failure
-        const bool oldThrowingError = FatalError.throwing(true);
+        const bool throwingError = FatalError.throwExceptions();

        try
        {
@ -191,7 +190,7 @@ int main(int argc, char *argv[])
            Info<< "Failed (expected) " << err << nl << endl;
        }

-        FatalError.throwing(oldThrowingError);
+        FatalError.throwExceptions(throwingError);
        report(list1);
    }

@ -199,7 +198,7 @@ int main(int argc, char *argv[])
        Info<< "\ntest operator[] with out-of-range index\n";

        // Expect failure
-        const bool oldThrowingError = FatalError.throwing(true);
+        const bool throwingError = FatalError.throwExceptions();

        try
        {
@ -213,7 +212,7 @@ int main(int argc, char *argv[])
            Info<< "Failed (expected) " << err << nl << endl;
        }

-        FatalError.throwing(oldThrowingError);
+        FatalError.throwExceptions(throwingError);
        report(list1);
    }

@ -269,7 +268,7 @@ int main(int argc, char *argv[])

    // Expect failure
    {
-        const bool oldThrowingError = FatalError.throwing(true);
+        const bool throwingError = FatalError.throwExceptions();

        Info<< "\ntest operator[] assignment with auto-vivify\n";

@ -287,7 +286,7 @@ int main(int argc, char *argv[])
            list1.set(36, list1.max_value);
        }

-        FatalError.throwing(oldThrowingError);
+        FatalError.throwExceptions(throwingError);
        report(list1);
    }

--- a/applications/test/PatchFunction1/Make/options
+++ b/applications/test/PatchFunction1/Make/options
@ -2,6 +2,7 @@ EXE_INC = \
    -DFULLDEBUG -g -O0 \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude

--- a/applications/test/PrecisionAdaptor/Make/files
+++ b/applications/test/PrecisionAdaptor/Make/files
@ -1,3 +0,0 @@
-Test-PrecisionAdaptor.C
-
-EXE = $(FOAM_USER_APPBIN)/Test-PrecisionAdaptor
--- a/applications/test/PrecisionAdaptor/Make/options
+++ b/applications/test/PrecisionAdaptor/Make/options
@ -1,2 +0,0 @@
-/* EXE_INC = */
-/* EXE_LIBS = */
--- a/applications/test/PrecisionAdaptor/Test-PrecisionAdaptor.C
+++ b/applications/test/PrecisionAdaptor/Test-PrecisionAdaptor.C
@ -1,111 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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
-    Test-PrecisionAdaptor
-
-Description
-
-\*---------------------------------------------------------------------------*/
-
-#include "argList.H"
-#include "primitiveFields.H"
-#include "PrecisionAdaptor.H"
-
-using namespace Foam;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    Field<double> content1(8);
-    Field<double> content2(8);
-
-    forAll(content1, i)
-    {
-        content1[i] = 10 * i;
-        content2[i] = 10 * i;
-    }
-
-    Foam::reverse(content2);
-
-    ConstPrecisionAdaptor<float, double, Field> cadaptor1;
-
-    cadaptor1.set(content1);
-    cadaptor1.commit();  // This is a no-op
-    Info<< "wrapped: " << cadaptor1() << nl;
-
-    cadaptor1.set(content2);
-    Info<< "wrapped: " << cadaptor1() << nl;
-
-    Info<< nl;
-
-    PrecisionAdaptor<float, double, Field> adaptor2;
-
-    adaptor2.set(content1);
-    adaptor2.ref() *= 2;
-    adaptor2.commit();  // Propagate changes back to input now
-
-    Info<< "modified wrapped: " << adaptor2() << nl;
-
-    adaptor2.set(content2);
-    adaptor2.ref() *= 2;
-    adaptor2.commit();  // Propagate changes back to input now
-
-    Info<< "modified wrapped: " << adaptor2() << nl;
-    Info<< "source: " << content1 << nl;
-    Info<< "source: " << content2 << nl;
-
-
-    content2 *= 2;
-    Info<< nl
-        << "set with " << content2 << nl;
-    Info<< "wrapped was " << adaptor2() << nl;
-    adaptor2.set(content2);
-    Info<< "wrapped now " << adaptor2() << nl;
-    Info<< "source: " << content2 << nl;
-
-    // Can even do this
-    Foam::reverse(adaptor2.ref());
-
-    adaptor2.ref() *= 2;
-    adaptor2.set(content1);  // implicit commit
-    Info<< "updated: " << content2 << nl;
-
-    Info<< nl
-        << "input: " << content1 << nl;
-
-    adaptor2.ref() *= 2;
-    adaptor2.clear();  // discard
-    adaptor2.commit();
-
-    Info<< "unchanged: " << content1 << nl;
-
-    Info<< nl << "Done" << nl << endl;
-    return 0;
-}
-
-
-// ************************************************************************* //
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
Mark Olesen	127f4268af	CONFIG: bump 2106 patch level to 211215	2021-12-15 11:08:43 +01:00
Mark Olesen	0bd1b0feac	BACKPORT: config script bugfixes - align with fixes for v2112 - align foamConfigurePaths with v2112 modifications CONFIG: scan scotch-64.h as possible fallback (#1904) - For RedHat 8, the /usr/include/scotch.h is actually wrapped with a condition check based on LP64 vs ILP32, so need to check with those headers if the first one failed.	2021-12-10 14:00:10 +01:00
Mark Olesen	42178fdd5d	BACKPORT: simpler hashing of string-types includes: ENH: simplify hashing overloads of string-types - this revises the changes made in `95cd8ee75c` to replace the SFINAE-type of handling of string hashes with direct definitions. This places a bit more burden on the developer if creating hashable classes derived from std::string or variants of Foam::string, but improves reliability when linking. STYLE: drop template key defaulting from HashSet - this was never used and `HashSet<>` is much less transparent than writing `HashSet<word>` or `wordHashSet` STYLE: use Hash<word> instead of string::hasher for runTimeSelectionTables	2021-12-09 18:00:59 +01:00
Mark Olesen	f815a12bba	BUG: error with empty distributed roots specification (fixes #2196 ) - previously used the size of distributed roots to transmit if the case was running in distributed mode, but this behaves rather poorly with bad input. Specifically, the following questionable setup: distributed true; roots ( /none/ ); Now transmit the ParRunControl distributed() value instead, and also emit a gentle warning for the user: WARNING: running distributed but did not specify roots!	2021-09-08 09:29:27 +02:00
				`@ -1 +0,0 @@`
				`regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();`