diff --git a/applications/solvers/DNS/dnsFoam/createFields.H b/applications/solvers/DNS/dnsFoam/createFields.H
index 64f4ccf4652344f392899588dcaa208bac408e7f..67a17d1677ed167cd7d955907d834d494d7b4f11 100644
--- a/applications/solvers/DNS/dnsFoam/createFields.H
+++ b/applications/solvers/DNS/dnsFoam/createFields.H
@@ -30,6 +30,6 @@ volVectorField U
 
 #include "createPhi.H"
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 #include "readTurbulenceProperties.H"
diff --git a/applications/solvers/basic/potentialFoam/createControls.H b/applications/solvers/basic/potentialFoam/createControls.H
index 1d0c80b47e73468c7145c7235f13775cf57a12ce..81cf83289dc74f7e9579e301f1951b85bcb8988a 100644
--- a/applications/solvers/basic/potentialFoam/createControls.H
+++ b/applications/solvers/basic/potentialFoam/createControls.H
@@ -1,6 +1,6 @@
 const dictionary& potentialFlow
 (
-    mesh.solutionDict().subDict("potentialFlow")
+    mesh.solution().solutionDict("potentialFlow")
 );
 
 const int nNonOrthCorr
diff --git a/applications/solvers/basic/potentialFoam/createFields.H b/applications/solvers/basic/potentialFoam/createFields.H
index 352458d00c04cc89da4391756d1fa9ae741e921d..6808b8f937cb165c5307410ba093fa35673a61c4 100644
--- a/applications/solvers/basic/potentialFoam/createFields.H
+++ b/applications/solvers/basic/potentialFoam/createFields.H
@@ -113,6 +113,6 @@ setRefCell
     PhiRefCell,
     PhiRefValue
 );
-mesh.setFluxRequired(Phi.name());
+mesh.schemes().setFluxRequired(Phi.name());
 
 #include "createMRF.H"
diff --git a/applications/solvers/basic/potentialFoam/overPotentialFoam/createControls.H b/applications/solvers/basic/potentialFoam/overPotentialFoam/createControls.H
index 1d0c80b47e73468c7145c7235f13775cf57a12ce..81cf83289dc74f7e9579e301f1951b85bcb8988a 100644
--- a/applications/solvers/basic/potentialFoam/overPotentialFoam/createControls.H
+++ b/applications/solvers/basic/potentialFoam/overPotentialFoam/createControls.H
@@ -1,6 +1,6 @@
 const dictionary& potentialFlow
 (
-    mesh.solutionDict().subDict("potentialFlow")
+    mesh.solution().solutionDict("potentialFlow")
 );
 
 const int nNonOrthCorr
diff --git a/applications/solvers/combustion/PDRFoam/bEqn.H b/applications/solvers/combustion/PDRFoam/bEqn.H
index e9f1a732799171edcd1a9b1b5e42afd46620a833..446f1046d4d50c106ed619ee9e8a8ae5d0aff17e 100644
--- a/applications/solvers/combustion/PDRFoam/bEqn.H
+++ b/applications/solvers/combustion/PDRFoam/bEqn.H
@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.divScheme("div(phi,ft_b_ha_hau)")
+        mesh.schemes().div("div(phi,ft_b_ha_hau)")
     )
 );
 
diff --git a/applications/solvers/combustion/PDRFoam/createFields.H b/applications/solvers/combustion/PDRFoam/createFields.H
index a0568cd489540558c4e7b67fac39d829420349d6..c4dfc3cf60bd44c56c81be6c9c03463721c05bad 100644
--- a/applications/solvers/combustion/PDRFoam/createFields.H
+++ b/applications/solvers/combustion/PDRFoam/createFields.H
@@ -43,7 +43,7 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::RASModel> turbulence
diff --git a/applications/solvers/combustion/XiFoam/createFields.H b/applications/solvers/combustion/XiFoam/createFields.H
index fa7f36ddf21ab41284b4ffd9f7fbfab35014c524..ba5751f6c54f2f70c34fca78ed0440fc47e7ae80 100644
--- a/applications/solvers/combustion/XiFoam/createFields.H
+++ b/applications/solvers/combustion/XiFoam/createFields.H
@@ -44,7 +44,7 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
diff --git a/applications/solvers/combustion/XiFoam/ftEqn.H b/applications/solvers/combustion/XiFoam/ftEqn.H
index c0f7dcad9c796f13158ee72b901d5a645144d507..b54157988d740cb5c6f5382293d3a89c1192ebee 100644
--- a/applications/solvers/combustion/XiFoam/ftEqn.H
+++ b/applications/solvers/combustion/XiFoam/ftEqn.H
@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.divScheme("div(phi,ft_b_ha_hau)")
+        mesh.schemes().div("div(phi,ft_b_ha_hau)")
     )
 );
 
diff --git a/applications/solvers/combustion/fireFoam/YEEqn.H b/applications/solvers/combustion/fireFoam/YEEqn.H
index e4f6afa17ab1ac609bea4ac0415f1d9eaddf5d53..b8387d42a5c3704390d0cea1d925f7d489d074b3 100644
--- a/applications/solvers/combustion/fireFoam/YEEqn.H
+++ b/applications/solvers/combustion/fireFoam/YEEqn.H
@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.divScheme("div(phi,Yi_h)")
+        mesh.schemes().div("div(phi,Yi_h)")
     )
 );
 {
diff --git a/applications/solvers/combustion/fireFoam/createFields.H b/applications/solvers/combustion/fireFoam/createFields.H
index 7d8fba861972b52bb0dc5fd3e10918f6d5f5a51d..6c59644d015fd5e9213fa2e0f54c4b407930f4fa 100644
--- a/applications/solvers/combustion/fireFoam/createFields.H
+++ b/applications/solvers/combustion/fireFoam/createFields.H
@@ -89,7 +89,7 @@ volScalarField p_rgh
     mesh
 );
 
-mesh.setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
 
 #include "phrghEqn.H"
 
diff --git a/applications/solvers/combustion/reactingFoam/YEqn.H b/applications/solvers/combustion/reactingFoam/YEqn.H
index 5d3ed162072576dd37574a9d754d78ecc68b5255..5d90002717b30cb56f1d56772c579264f1105b83 100644
--- a/applications/solvers/combustion/reactingFoam/YEqn.H
+++ b/applications/solvers/combustion/reactingFoam/YEqn.H
@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.divScheme("div(phi,Yi_h)")
+        mesh.schemes().div("div(phi,Yi_h)")
     )
 );
 
diff --git a/applications/solvers/combustion/reactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/createFields.H
index 4c9f34015db0d34fcf0dac242390743ae5237c44..95cdd06ac5e16270ca858a0e39706a5fb489b921 100644
--- a/applications/solvers/combustion/reactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/createFields.H
@@ -47,7 +47,7 @@ volScalarField& p = thermo.p();
 
 pressureControl pressureControl(p, rho, pimple.dict(), false);
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 Info << "Creating turbulence model.\n" << nl;
 autoPtr<compressible::turbulenceModel> turbulence
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H
index 94b40e98ddf1253e8d2a3b0b670b677e98a85727..f86e1f519c45d6c893f3dcb519d83f1acaf82b53 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H
@@ -48,7 +48,7 @@ volScalarField& p = thermo.p();
 
 #include "compressibleCreatePhi.H"
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 Info << "Creating turbulence model.\n" << nl;
 autoPtr<compressible::turbulenceModel> turbulence
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H
index 2e830f7bcbe60acd5ad5f27c1e4f28de7ac259ae..aafe68fadca1a9612115e1b5be33e64a2520b0c0 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H
@@ -49,7 +49,7 @@ volScalarField& p = thermo.p();
 
 pressureControl pressureControl(p, rho, pimple.dict(), false);
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 
 Info << "Creating turbulence model.\n" << nl;
diff --git a/applications/solvers/compressible/rhoCentralFoam/readFluxScheme.H b/applications/solvers/compressible/rhoCentralFoam/readFluxScheme.H
index e8c28d65c269f444f1b9dfa87ef1a50207f9ff87..4cba5629f7aec305336c2ab7e067f8d8a2d330e2 100644
--- a/applications/solvers/compressible/rhoCentralFoam/readFluxScheme.H
+++ b/applications/solvers/compressible/rhoCentralFoam/readFluxScheme.H
@@ -1,5 +1,5 @@
 word fluxScheme("Kurganov");
-if (mesh.schemesDict().readIfPresent("fluxScheme", fluxScheme))
+if (mesh.schemes().dict().readIfPresent("fluxScheme", fluxScheme))
 {
     if ((fluxScheme == "Tadmor") || (fluxScheme == "Kurganov"))
     {
diff --git a/applications/solvers/compressible/rhoPimpleAdiabaticFoam/createFields.H b/applications/solvers/compressible/rhoPimpleAdiabaticFoam/createFields.H
index 4fd5569a33bc4a43b740bab9532d918dfe3c5f40..f04ac5695fcf047b49642c9452bd2852afe8bc87 100644
--- a/applications/solvers/compressible/rhoPimpleAdiabaticFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleAdiabaticFoam/createFields.H
@@ -79,7 +79,7 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 #include "createMRF.H"
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
index d9fa4879908a6e1e992a04010110eccd671e50ca..531fd411a032592c9cf9c6018d5fac2b413eeb5a 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
@@ -42,7 +42,7 @@ volVectorField U
 
 pressureControl pressureControl(p, rho, pimple.dict(), false);
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
diff --git a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H
index 190763d27936ba6a728d14fca90b8b24dfe40690..d8d7471c032bc7007fcfc9f326defa4fb365818f 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H
@@ -56,7 +56,7 @@ if (mesh.changing())
         pcorrTypes
     );
 
-    mesh.setFluxRequired(pcorr.name());
+    mesh.schemes().setFluxRequired(pcorr.name());
 
 {
     dimensionedScalar rAUf("rAUf", dimTime, 1.0);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createFields.H
index f69bfd980ff58106d8aaf3c1df35959e20fc5659..985c56b0035523f047646a5dbaeead5ea47cd681 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createFields.H
@@ -44,7 +44,7 @@ pressureControl pressureControl(p, rho, pimple.dict(), false);
 const dimensionedScalar rhoMax("rhoMax", dimDensity, GREAT, pimple.dict());
 const dimensionedScalar rhoMin("rhoMin", dimDensity, Zero, pimple.dict());
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 #include "createDpdt.H"
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
index 6822efa11ac4622c2c8bd75d9e7d0b9bebc64341..fd44338a9d801bfbcb98916f96d7780a94231840 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
@@ -40,7 +40,7 @@ volVectorField U
 
 pressureControl pressureControl(p, rho, simple.dict());
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
diff --git a/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createFields.H
index a04acdf9752585d80bcb07a7ccc6db89b021025d..b8765c5f5322023252a040cfb1352b2917bde3e1 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createFields.H
@@ -40,7 +40,7 @@ volVectorField U
 
 pressureControl pressureControl(p, rho, simple.dict());
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
diff --git a/applications/solvers/compressible/sonicFoam/createFields.H b/applications/solvers/compressible/sonicFoam/createFields.H
index 86d85a8c82997bf79f9d8e6aa3854fef787e34c3..a4867f9230d91e1434252759122c532cee522192 100644
--- a/applications/solvers/compressible/sonicFoam/createFields.H
+++ b/applications/solvers/compressible/sonicFoam/createFields.H
@@ -36,7 +36,7 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
diff --git a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/createFields.H b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/createFields.H
index d94b173feec4790996d3303b7929290c2f3e657d..f047e89fff0341ced43697cf90c9f0088c58d1b9 100644
--- a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/createFields.H
+++ b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/createFields.H
@@ -46,4 +46,4 @@ volScalarField rho
 
 #include "compressibleCreatePhi.H"
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
diff --git a/applications/solvers/electromagnetics/mhdFoam/createFields.H b/applications/solvers/electromagnetics/mhdFoam/createFields.H
index 1032504738d463f8c7a2e8357894966dde5d50cd..13a7efe1bc849f52d71205685c42230da3970518 100644
--- a/applications/solvers/electromagnetics/mhdFoam/createFields.H
+++ b/applications/solvers/electromagnetics/mhdFoam/createFields.H
@@ -114,5 +114,5 @@ label pRefCell = 0;
 scalar pRefValue = 0.0;
 setRefCell(p, piso.dict(), pRefCell, pRefValue);
 
-mesh.setFluxRequired(p.name());
-mesh.setFluxRequired(pB.name());
+mesh.schemes().setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(pB.name());
diff --git a/applications/solvers/electromagnetics/mhdFoam/readBPISOControls.H b/applications/solvers/electromagnetics/mhdFoam/readBPISOControls.H
index cb4e3505e95d8ef5e22c150a4c082379c014f75b..d0afaac7f39de5cced694ccb9a3a01a264f98fd6 100644
--- a/applications/solvers/electromagnetics/mhdFoam/readBPISOControls.H
+++ b/applications/solvers/electromagnetics/mhdFoam/readBPISOControls.H
@@ -1,3 +1,3 @@
-    const dictionary& Bpiso = mesh.solutionDict().subDict("BPISO");
+    const dictionary& Bpiso = mesh.solution().solutionDict("BPISO");
 
     const int nBcorr = Bpiso.getOrDefault<int>("nCorrectors", 1);
diff --git a/applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H b/applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H
index 0b4e94b2ca552c8169aa47868e1e7e22e201e5ad..cccdf702258d94f6a5ba533e29ac5847c878fc9b 100644
--- a/applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H
+++ b/applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H
@@ -155,4 +155,4 @@
         dimensionedScalar("one", dimless, 0.01)
     );
 
-    aMesh.setFluxRequired("h");
+    aMesh.schemes().setFluxRequired("h");
diff --git a/applications/solvers/finiteArea/liquidFilmFoam/readSolutionControls.H b/applications/solvers/finiteArea/liquidFilmFoam/readSolutionControls.H
index 3d0c97261c55250fa3124d6cb8fff5b328dfabbf..2aeef3c7dd3687ff0e2a7dd157552cad4aa4e68b 100644
--- a/applications/solvers/finiteArea/liquidFilmFoam/readSolutionControls.H
+++ b/applications/solvers/finiteArea/liquidFilmFoam/readSolutionControls.H
@@ -1 +1 @@
-loopControl iters(runTime, aMesh.solutionDict(), "solution");
+loopControl iters(runTime, aMesh.solution().solutionDict(), "solution");
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
index 908af133dc45275fb9e6a64923072f5d84918c55..87e6e088e6816a206599288ef157aaa5aa6c9d70 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
@@ -119,7 +119,7 @@ if (p_rgh.needReference())
     );
 }
 
-mesh.setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
 
 #include "createMRF.H"
 #include "createIncompressibleRadiationModel.H"
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
index a61f66aea08294af9d04d43a061b4b204469f2ce..fe927dcf52f12384c1c98e61d84096b82bc2ee00 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
@@ -119,7 +119,7 @@ if (p_rgh.needReference())
     );
 }
 
-mesh.setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
 
 #include "createMRF.H"
 #include "createIncompressibleRadiationModel.H"
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
index 813d1c324009bdf8dcc2f4c74d2d0f4975fb8955..508f2c949e4902f1b00b2a28fab2db0c379014ab 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
@@ -74,7 +74,7 @@ volScalarField p_rgh
 // Force p_rgh to be consistent with p
 p_rgh = p - rho*gh;
 
-mesh.setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
 
 label pRefCell = 0;
 scalar pRefValue = 0.0;
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/correctPhi.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/correctPhi.H
index 708a7254d2a4d75041ed825959f26f282bd6bb6b..87902dcfd651d5d34504caad0ae3d1835a7fab88 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/correctPhi.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/correctPhi.H
@@ -56,7 +56,7 @@ if (mesh.changing())
         pcorrTypes
     );
 
-    mesh.setFluxRequired(pcorr.name());
+    mesh.schemes().setFluxRequired(pcorr.name());
 
 {
     dimensionedScalar rAUf("rAUf", dimTime, 1.0);
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/createFields.H
index 050cb8714612b0a53a66230a753f18e5b0aeac39..57d1bb88b63f38f3cdd8c7dc27d7107a430dce96 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/createFields.H
@@ -71,7 +71,7 @@ volScalarField p_rgh
 // Force p_rgh to be consistent with p
 p_rgh = p - rho*gh;
 
-mesh.setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
 
 label pRefCell = 0;
 scalar pRefValue = 0.0;
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
index 33605471f11e778182aa98eea93cfbe7d2ffe477..095abbbcc550d195fe85a7e144762ec3c9f80478 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
@@ -81,7 +81,7 @@ setRefCell
     pRefValue
 );
 
-mesh.setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
 
 dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 dimensionedScalar totalVolume = sum(mesh.V());
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
index ec6a70303ab5a16aa57b1e814d50e20c60d9fe91..29609c2abcd45ea7b2a6788c314a86808297f6ef 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
@@ -174,7 +174,7 @@ forAll(fluidRegions, i)
     // Force p_rgh to be consistent with p
     p_rghFluid[i] = thermoFluid[i].p() - rhoFluid[i]*ghFluid[i];
 
-    fluidRegions[i].setFluxRequired(p_rghFluid[i].name());
+    fluidRegions[i].schemes().setFluxRequired(p_rghFluid[i].name());
 
     radiation.set
     (
@@ -185,7 +185,7 @@ forAll(fluidRegions, i)
     initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
 
     const dictionary& simpleDict =
-        fluidRegions[i].solutionDict().subDict("SIMPLE");
+        fluidRegions[i].solution().solutionDict("SIMPLE");
 
     setRefCell
     (
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H
index 9455864173810658015e99f2519c3b5326a530cd..0864ac2deb60d7310efcc876b7e76f425380056c 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H
@@ -1,4 +1,5 @@
-    const dictionary& simple = fluidRegions[i].solutionDict().subDict("SIMPLE");
+    const dictionary& simple =
+        fluidRegions[i].solution().solutionDict("SIMPLE");
 
     const int nNonOrthCorr =
         simple.getOrDefault<int>("nNonOrthogonalCorrectors", 0);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/readSolidMultiRegionSIMPLEControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/readSolidMultiRegionSIMPLEControls.H
index 1320e4ecc13ab5fcb4611072eb57e7b7a1a1db70..950a80fa776e52a90723723082f5b1e306dee40b 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/readSolidMultiRegionSIMPLEControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/readSolidMultiRegionSIMPLEControls.H
@@ -1,4 +1,5 @@
-const dictionary& simple = mesh.solutionDict().subDict("SIMPLE");
+const dictionary& simple =
+    mesh.solution().solutionDict("SIMPLE");
 
 const int nNonOrthCorr =
     simple.getOrDefault<int>("nNonOrthogonalCorrectors", 0);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/createFluidFields.H
index e7ca574f70d46b9495829703c33f3b6173c784fc..54395b8a6803e27a38d60f43e35ba76ae1fa8ddc 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/createFluidFields.H
@@ -278,14 +278,14 @@ forAll(fluidRegions, i)
         phaseSystemFluid[i].phase1().thermo().p()
      -  phaseSystemFluid[i].phase1().thermo().rho()*ghFluid[i];
 
-    fluidRegions[i].setFluxRequired(p_rghFluid[i].name());
+    fluidRegions[i].schemes().setFluxRequired(p_rghFluid[i].name());
 
     Info<< "    Correcting initialMassFluid\n" << endl;
     initialMassFluid[i] =
         fvc::domainIntegrate(phaseSystemFluid[i].rho()).value();
 
     const dictionary& pimpleDict =
-        fluidRegions[i].solutionDict().subDict("PIMPLE");
+        fluidRegions[i].solution().solutionDict("PIMPLE");
 
     pimpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/readFluidMultiRegionPIMPLEControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/readFluidMultiRegionPIMPLEControls.H
index 9e76ce95743fe5a4c54e1f7398b59cc089b64804..10d0ee8fe4ab536043980dec64c6917acc34c3ca 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/readFluidMultiRegionPIMPLEControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/readFluidMultiRegionPIMPLEControls.H
@@ -1,4 +1,5 @@
-    const dictionary& pimpleDict = mesh.solutionDict().subDict("PIMPLE");
+    const dictionary& pimpleDict =
+        mesh.solution().solutionDict("PIMPLE");
 
     Switch faceMomentum
     (
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H
index 1c17b76eeb7aec7890d53a6fd9869cc325ed59aa..74d7afa259f2ea7151b254e75165a07a3c5a51ee 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H
@@ -9,7 +9,7 @@ if (Y.size())
             mesh,
             fields,
             phi,
-            mesh.divScheme("div(phi,Yi_h)")
+            mesh.schemes().div("div(phi,Yi_h)")
         )
     );
 }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index 2dfd1d314671c104abf1ba7e09cf5e75c7a68936..a012511e52ced53b8fd2dc9e6c5b0194992324e0 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -189,7 +189,7 @@ forAll(fluidRegions, i)
     // Force p_rgh to be consistent with p
     p_rghFluid[i] = thermoFluid[i].p() - rhoFluid[i]*ghFluid[i];
 
-    fluidRegions[i].setFluxRequired(p_rghFluid[i].name());
+    fluidRegions[i].schemes().setFluxRequired(p_rghFluid[i].name());
 
     Info<< "    Adding to radiationFluid\n" << endl;
     radiation.set
@@ -260,7 +260,8 @@ forAll(fluidRegions, i)
     );
 
     const dictionary& pimpleDict =
-        fluidRegions[i].solutionDict().subDict("PIMPLE");
+        fluidRegions[i].solution().solutionDict("PIMPLE");
+
     pimpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
 
     rhoMaxFluid.set
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H
index d477e4992fef4d442e4ffdf487ef23cac09649dc..a9c699910eb89f99c33f3aa33a3503c9ccfee032 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H
@@ -1,4 +1,5 @@
-    const dictionary& pimple = mesh.solutionDict().subDict("PIMPLE");
+    const dictionary& pimple =
+        mesh.solution().solutionDict("PIMPLE");
 
     const int nCorr =
         pimple.getOrDefault<int>("nCorrectors", 1);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPIMPLEControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPIMPLEControls.H
index 13d35be3cac2eb847c3f8261910f765254369b6c..5469ddfe4976d2bc995993eba509ac0e9a7b9866 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPIMPLEControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPIMPLEControls.H
@@ -1,4 +1,5 @@
-const dictionary& pimple = mesh.solutionDict().subDict("PIMPLE");
+const dictionary& pimple =
+    mesh.solution().solutionDict("PIMPLE");
 
 int nNonOrthCorr =
     pimple.getOrDefault<int>("nNonOrthogonalCorrectors", 0);
diff --git a/applications/solvers/heatTransfer/solidFoam/createFields.H b/applications/solvers/heatTransfer/solidFoam/createFields.H
index 85fc531010f3c96f9913e663d009aaa8bf758061..57006dc5810d66151e541bec78073d90cec1036e 100644
--- a/applications/solvers/heatTransfer/solidFoam/createFields.H
+++ b/applications/solvers/heatTransfer/solidFoam/createFields.H
@@ -92,21 +92,19 @@ else
 
 // Consider mesh flux to correct for mesh deformation
 bool meshFluxCorr(false);
-if (mesh.solutionDict().found("SIMPLE"))
 {
-    meshFluxCorr =
-        mesh.solutionDict().subDict("SIMPLE").getOrDefault<bool>
-        (
-            "meshFluxCorrection", false
-        );
-}
-else if (mesh.solutionDict().found("PIMPLE"))
-{
-    meshFluxCorr =
-        mesh.solutionDict().subDict("PIMPLE").getOrDefault<bool>
-        (
-            "meshFluxCorrection", false
-        );
+    const dictionary& solutionDict = mesh.solution().solutionDict();
+
+    const dictionary* subdict = nullptr;
+
+    if
+    (
+        ((subdict = solutionDict.findDict("SIMPLE")) != nullptr)
+     || ((subdict = solutionDict.findDict("PIMPLE")) != nullptr)
+    )
+    {
+        meshFluxCorr = subdict->getOrDefault("meshFluxCorrection", false);
+    }
 }
 
 #include "createRadiationModel.H"
diff --git a/applications/solvers/heatTransfer/solidFoam/solidFoam.C b/applications/solvers/heatTransfer/solidFoam/solidFoam.C
index 6b75bf7c4dddbd1cecbb04b8362b39608de3a4ab..60df4edcc302e03ad36d448c6e801e110a3d9a78 100644
--- a/applications/solvers/heatTransfer/solidFoam/solidFoam.C
+++ b/applications/solvers/heatTransfer/solidFoam/solidFoam.C
@@ -68,7 +68,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nEvolving thermodynamics\n" << endl;
 
-    if (mesh.solutionDict().found("SIMPLE"))
+    if (mesh.solution().solutionDict().found("SIMPLE"))
     {
         simpleControl simple(mesh);
 
diff --git a/applications/solvers/heatTransfer/thermoFoam/thermoFoam.C b/applications/solvers/heatTransfer/thermoFoam/thermoFoam.C
index 77be2e22601c3a9791816b70ed4f76aab2838698..bf3f2b168543b5a2e85b78298ae81a6cf288a00b 100644
--- a/applications/solvers/heatTransfer/thermoFoam/thermoFoam.C
+++ b/applications/solvers/heatTransfer/thermoFoam/thermoFoam.C
@@ -65,7 +65,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nEvolving thermodynamics\n" << endl;
 
-    if (mesh.solutionDict().found("SIMPLE"))
+    if (mesh.solution().solutionDict().found("SIMPLE"))
     {
         simpleControl simple(mesh);
 
diff --git a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
index be66f4cc539bf1c4a71dd7b8fd17eb068d35e4b1..012917ca90aeb98b5615a39554f4d28f1fffc9a2 100644
--- a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
+++ b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
@@ -106,7 +106,7 @@ int main(int argc, char *argv[])
         //    mesh.relaxationFactor("alpha")
         //   *(lambda*max(Ua & U, zeroSensitivity) - alpha);
         alpha +=
-            mesh.fieldRelaxationFactor("alpha")
+            mesh.solution().fieldRelaxationFactor("alpha")
            *(min(max(alpha + lambda*(Ua & U), zeroAlpha), alphaMax) - alpha);
 
         zeroCells(alpha, inletCells);
diff --git a/applications/solvers/incompressible/adjointShapeOptimizationFoam/createFields.H b/applications/solvers/incompressible/adjointShapeOptimizationFoam/createFields.H
index 4d7b9d881106ae3be279b210dd101c060500ecd4..4f4b6780154c00ca56bfa9f721b9b060e2e1cc71 100644
--- a/applications/solvers/incompressible/adjointShapeOptimizationFoam/createFields.H
+++ b/applications/solvers/incompressible/adjointShapeOptimizationFoam/createFields.H
@@ -32,7 +32,7 @@ volVectorField U
 label pRefCell = 0;
 scalar pRefValue = 0.0;
 setRefCell(p, simple.dict(), pRefCell, pRefValue);
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 
 Info<< "Reading field pa\n" << endl;
@@ -75,7 +75,7 @@ setRefCell
     paRefCell,
     paRefValue
 );
-mesh.setFluxRequired(pa.name());
+mesh.schemes().setFluxRequired(pa.name());
 
 
 singlePhaseTransportModel laminarTransport(U, phi);
diff --git a/applications/solvers/incompressible/icoFoam/createFields.H b/applications/solvers/incompressible/icoFoam/createFields.H
index 02d916bcba9143ad80556b32886ee55bbefc197f..847a1c712479ca021042a0be95efe79aeb60be80 100644
--- a/applications/solvers/incompressible/icoFoam/createFields.H
+++ b/applications/solvers/incompressible/icoFoam/createFields.H
@@ -54,5 +54,11 @@ volVectorField U
 
 label pRefCell = 0;
 scalar pRefValue = 0.0;
-setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-mesh.setFluxRequired(p.name());
+setRefCell
+(
+    p,
+    mesh.solution().solutionDict("PISO"),
+    pRefCell,
+    pRefValue
+);
+mesh.schemes().setFluxRequired(p.name());
diff --git a/applications/solvers/incompressible/nonNewtonianIcoFoam/createFields.H b/applications/solvers/incompressible/nonNewtonianIcoFoam/createFields.H
index d2ec50c8817ebabc59c6d2210729079bb60840e4..f88b03f5cb15f1a15ac08c213e9534004a4f4447 100644
--- a/applications/solvers/incompressible/nonNewtonianIcoFoam/createFields.H
+++ b/applications/solvers/incompressible/nonNewtonianIcoFoam/createFields.H
@@ -36,5 +36,11 @@ singlePhaseTransportModel fluid(U, phi);
 
 label pRefCell = 0;
 scalar pRefValue = 0.0;
-setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-mesh.setFluxRequired(p.name());
+setRefCell
+(
+    p,
+    mesh.solution().solutionDict("PISO"),
+    pRefCell,
+    pRefValue
+);
+mesh.schemes().setFluxRequired(p.name());
diff --git a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/createFields.H b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/createFields.H
index b2ccf32cc8da4c18d55fe6f698b168f75505c6bc..39318782068b0e9321598cd98b3c11a1d626d73f 100644
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/createFields.H
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/createFields.H
@@ -43,7 +43,7 @@ surfaceScalarField phi
 label pRefCell = 0;
 scalar pRefValue = 0.0;
 setRefCell(p, pimple.dict(), pRefCell, pRefValue);
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 Info<< "Creating SRF model\n" << endl;
 autoPtr<SRF::SRFModel> SRF
diff --git a/applications/solvers/incompressible/pimpleFoam/createFields.H b/applications/solvers/incompressible/pimpleFoam/createFields.H
index 4aecb61b06933319b10897af5d9ef470c6c3b7da..b36316cefe1e6bee8b316e0d9ad46961b18358ef 100644
--- a/applications/solvers/incompressible/pimpleFoam/createFields.H
+++ b/applications/solvers/incompressible/pimpleFoam/createFields.H
@@ -34,7 +34,7 @@ volVectorField U
 label pRefCell = 0;
 scalar pRefValue = 0.0;
 setRefCell(p, pimple.dict(), pRefCell, pRefValue);
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 
 singlePhaseTransportModel laminarTransport(U, phi);
diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createFields.H b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createFields.H
index eec6d4e5e6956e25884696aa724f5ade6d282110..55fa20ded34be055eaaffb3f00311753e4d9c271 100644
--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createFields.H
+++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createFields.H
@@ -32,7 +32,7 @@ volVectorField U
 label pRefCell = 0;
 scalar pRefValue = 0.0;
 setRefCell(p, pimple.dict(), pRefCell, pRefValue);
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 
 //- Overset specific
diff --git a/applications/solvers/incompressible/pisoFoam/createFields.H b/applications/solvers/incompressible/pisoFoam/createFields.H
index 6485334fec48f098ca51a8ce8efdb77cc86fc87f..aae3cfecf2bdde50111ba3beb458359428ce68d0 100644
--- a/applications/solvers/incompressible/pisoFoam/createFields.H
+++ b/applications/solvers/incompressible/pisoFoam/createFields.H
@@ -32,7 +32,7 @@ volVectorField U
 label pRefCell = 0;
 scalar pRefValue = 0.0;
 setRefCell(p, piso.dict(), pRefCell, pRefValue);
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 
 singlePhaseTransportModel laminarTransport(U, phi);
diff --git a/applications/solvers/incompressible/shallowWaterFoam/createFields.H b/applications/solvers/incompressible/shallowWaterFoam/createFields.H
index aa121410c817cecc306e4a8c44dfbebbddcb7ff1..7bed9031f9a697ff551d875a51572a56ff6280c0 100644
--- a/applications/solvers/incompressible/shallowWaterFoam/createFields.H
+++ b/applications/solvers/incompressible/shallowWaterFoam/createFields.H
@@ -76,4 +76,4 @@ hTotal.write();
 Info<< "Creating Coriolis Force" << endl;
 const dimensionedVector F("F", ((2.0*Omega) & gHat)*gHat);
 
-mesh.setFluxRequired(h.name());
+mesh.schemes().setFluxRequired(h.name());
diff --git a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/createFields.H b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/createFields.H
index b7b0a5d10b97279c53dd30b9c8db9653310f87e9..bc9657d3b07af5ec12057bce8f790fe1c7085332 100644
--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/createFields.H
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/createFields.H
@@ -43,7 +43,7 @@ surfaceScalarField phi
 label pRefCell = 0;
 scalar pRefValue = 0.0;
 setRefCell(p, simple.dict(), pRefCell, pRefValue);
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 Info<< "Creating SRF model\n" << endl;
 autoPtr<SRF::SRFModel> SRF(SRF::SRFModel::New(Urel));
diff --git a/applications/solvers/incompressible/simpleFoam/createFields.H b/applications/solvers/incompressible/simpleFoam/createFields.H
index 3c4971527b0b544882c2984f6831e80fc5b9fd7e..9f7c3284f57a745eb37ce69c994c8320e7da7f96 100644
--- a/applications/solvers/incompressible/simpleFoam/createFields.H
+++ b/applications/solvers/incompressible/simpleFoam/createFields.H
@@ -32,7 +32,7 @@ volVectorField U
 label pRefCell = 0;
 scalar pRefValue = 0.0;
 setRefCell(p, simple.dict(), pRefCell, pRefValue);
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 
 singlePhaseTransportModel laminarTransport(U, phi);
diff --git a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createFields.H b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createFields.H
index 86789ed65bf8c758c919898f52fa0d94f6dbd0a5..9e67166d00d2578394f269decb53c20eb1b8d7dd 100644
--- a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createFields.H
+++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createFields.H
@@ -32,7 +32,7 @@ volVectorField U
 label pRefCell = 0;
 scalar pRefValue = 0.0;
 setRefCell(p, simple.dict(), pRefCell, pRefValue);
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 
 singlePhaseTransportModel laminarTransport(U, phi);
diff --git a/applications/solvers/lagrangian/DPMFoam/createFields.H b/applications/solvers/lagrangian/DPMFoam/createFields.H
index 256d05eef8de0d1a957f65089cd0302b2c1177e5..a5bcd470af8fd0bbed85dd7c75b44e118ccb418c 100644
--- a/applications/solvers/lagrangian/DPMFoam/createFields.H
+++ b/applications/solvers/lagrangian/DPMFoam/createFields.H
@@ -62,7 +62,7 @@ surfaceScalarField phic
 label pRefCell = 0;
 scalar pRefValue = 0.0;
 setRefCell(p, pimple.dict(), pRefCell, pRefValue);
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 Info<< "Creating turbulence model\n" << endl;
 
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
index 59bb2aae16725529d397a9b630c85513bcc7aa8f..a86ba02f429e32602333fe729ab23b4b8ca825cc 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.divScheme("div(phi,Yi_h)")
+        mesh.schemes().div("div(phi,Yi_h)")
     )
 );
 
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
index 8324b91ae042b87bc9f68ca92adfad7923cf519b..b68d0a2555e997a655cc42d355e323118748c7cb 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
@@ -91,7 +91,7 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
diff --git a/applications/solvers/lagrangian/kinematicParcelFoam/createFields.H b/applications/solvers/lagrangian/kinematicParcelFoam/createFields.H
index 0182144d86e85803af9663c0c508d439831f1119..efe90bca4da9aec2a8b466ed4a8fb2c570a98e7c 100644
--- a/applications/solvers/lagrangian/kinematicParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/createFields.H
@@ -75,7 +75,7 @@ autoPtr<incompressible::turbulenceModel> turbulence
 label pRefCell = 0;
 scalar pRefValue = 0.0;
 setRefCell(p, pimple.dict(), pRefCell, pRefValue);
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 #include "createMRF.H"
 #include "createClouds.H"
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
index c0179eadb860d572bbfa457878efce1061e2595e..87275dd2ac6798136778a93ce347481990eedec6 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.divScheme("div(phi,Yi_h)")
+        mesh.schemes().div("div(phi,Yi_h)")
     )
 );
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
index 6984cc617819b8b05a9d6244b701446e671915ba..be14aedbcad191420c634b2136f16ae1ffe879bd 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
@@ -96,7 +96,7 @@ p_rgh = p - rho*gh;
 
 pressureControl pressureControl(p, rho, pimple.dict(), false);
 
-mesh.setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
 
 Info<< "Creating multi-variate interpolation scheme\n" << endl;
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H
index 925c55df3253f646f5c3e4a29f2c662b9adfd471..f4bebcbcdc5ebac22194dd7d8b9161713f1fd56f 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H
@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.divScheme("div(phi,Yi_h)")
+        mesh.schemes().div("div(phi,Yi_h)")
     )
 );
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createFields.H
index 49f8a1fc42c01e5603067890cb1c86d9651ef5f8..f23a3363870a1782ca7fde855445b59cc52d79cd 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createFields.H
@@ -50,7 +50,7 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 const dimensionedScalar rhoMax("rhoMax", dimDensity, GREAT, simple.dict());
 const dimensionedScalar rhoMin("rhoMin", dimDensity, Zero, simple.dict());
diff --git a/applications/solvers/lagrangian/simpleCoalParcelFoam/YEqn.H b/applications/solvers/lagrangian/simpleCoalParcelFoam/YEqn.H
index 27c1fd46db57e274f10bfcb60eb42032d7a4c501..3c82abe9260f25d24cd4f5cccdf42177c797a408 100644
--- a/applications/solvers/lagrangian/simpleCoalParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/simpleCoalParcelFoam/YEqn.H
@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.divScheme("div(phi,Yi_h)")
+        mesh.schemes().div("div(phi,Yi_h)")
     )
 );
 
diff --git a/applications/solvers/lagrangian/simpleCoalParcelFoam/createFields.H b/applications/solvers/lagrangian/simpleCoalParcelFoam/createFields.H
index 581a52d28537fb4de03eb893f9be84536d831ece..753afc0831959fef255d2d5e8f87a38fda680913 100644
--- a/applications/solvers/lagrangian/simpleCoalParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/simpleCoalParcelFoam/createFields.H
@@ -51,7 +51,7 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 const dimensionedScalar rhoMax("rhoMax", dimDensity, GREAT, simple.dict());
 const dimensionedScalar rhoMin("rhoMin", dimDensity, Zero, simple.dict());
diff --git a/applications/solvers/lagrangian/sprayFoam/YEqn.H b/applications/solvers/lagrangian/sprayFoam/YEqn.H
index f33c64497dd1becff376c826a07b42333010d350..07cdf63cd9d7c6f8a76411ced14f7817a2636821 100644
--- a/applications/solvers/lagrangian/sprayFoam/YEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/YEqn.H
@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.divScheme("div(phi,Yi_h)")
+        mesh.schemes().div("div(phi,Yi_h)")
     )
 );
 
diff --git a/applications/solvers/lagrangian/sprayFoam/createFields.H b/applications/solvers/lagrangian/sprayFoam/createFields.H
index 9eb279fe55b158bd10bc117e29890f3b27ce7fca..d1b3baa82bab626c99ebf6eacbb83ac5fb1cc239 100644
--- a/applications/solvers/lagrangian/sprayFoam/createFields.H
+++ b/applications/solvers/lagrangian/sprayFoam/createFields.H
@@ -50,7 +50,7 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 const dimensionedScalar rhoMax("rhoMax", dimDensity, GREAT, pimple.dict());
 const dimensionedScalar rhoMin("rhoMin", dimDensity, Zero, pimple.dict());
diff --git a/applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H b/applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H
index 33b6bd75ded220931577fdfca5db0ce24af74604..2cbe3c361885c5716899012085111e341aa4bfd9 100644
--- a/applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H
@@ -10,7 +10,7 @@
             fv::ddtScheme<scalar>::New
             (
                 mesh,
-                mesh.ddtScheme("ddt(alpha)")
+                mesh.schemes().ddt("ddt(alpha)")
             )
         );
         const fv::ddtScheme<scalar>& ddtAlpha = tddtAlpha();
@@ -237,7 +237,7 @@
 
     if
     (
-        word(mesh.ddtScheme("ddt(rho,U)"))
+        word(mesh.schemes().ddt("ddt(rho,U)"))
      == fv::EulerDdtScheme<vector>::typeName
     )
     {
diff --git a/applications/solvers/multiphase/MPPICInterFoam/createFields.H b/applications/solvers/multiphase/MPPICInterFoam/createFields.H
index 1c93103b6e3b4bfc21a7764332b9ebc6c40610df..50d392d3c365540f6e7ef132672a8c59a09b62d3 100644
--- a/applications/solvers/multiphase/MPPICInterFoam/createFields.H
+++ b/applications/solvers/multiphase/MPPICInterFoam/createFields.H
@@ -105,7 +105,7 @@ setRefCell
 (
     p,
     p_rgh,
-    mesh.solutionDict().subDict("PIMPLE"),
+    mesh.solution().solutionDict("PIMPLE"),
     pRefCell,
     pRefValue
 );
@@ -121,8 +121,8 @@ if (p_rgh.needReference())
     p_rgh = p - rho*gh;
 }
 
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha1.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(alpha1.name());
 
 // alphac must be constructed before the cloud
 // so that the drag-models can find it
diff --git a/applications/solvers/multiphase/VoF/alphaEqn.H b/applications/solvers/multiphase/VoF/alphaEqn.H
index c6f318fa2782a7fbb04aa59b8a11b6f1b92b0f28..12a57232da7f4909ea03ef512c769332eb0fbfaf 100644
--- a/applications/solvers/multiphase/VoF/alphaEqn.H
+++ b/applications/solvers/multiphase/VoF/alphaEqn.H
@@ -10,7 +10,7 @@
             fv::ddtScheme<scalar>::New
             (
                 mesh,
-                mesh.ddtScheme("ddt(alpha)")
+                mesh.schemes().ddt("ddt(alpha)")
             )
         );
         const fv::ddtScheme<scalar>& ddtAlpha = tddtAlpha();
@@ -239,9 +239,9 @@
 
     if
     (
-        word(mesh.ddtScheme("ddt(rho,U)"))
+        word(mesh.schemes().ddt("ddt(rho,U)"))
      == fv::EulerDdtScheme<vector>::typeName
-     || word(mesh.ddtScheme("ddt(rho,U)"))
+     || word(mesh.schemes().ddt("ddt(rho,U)"))
      == fv::localEulerDdtScheme<vector>::typeName
     )
     {
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
index 4bde1b55f1cdaf8b059fc58bfe9e8b4507701e1e..7271037aafc7fb1fb5d2cb46d7072a7a98193124 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
@@ -19,7 +19,7 @@ correctUphiBCs(U, phi);
     surfaceScalarField rhof(fvc::interpolate(rho, "div(phi,rho)"));
     dimensionedScalar rAUf("rAUf", dimTime, 1.0);
 
-    mesh.setFluxRequired(pcorr.name());
+    mesh.schemes().setFluxRequired(pcorr.name());
 
     while (pimple.correctNonOrthogonal())
     {
diff --git a/applications/solvers/multiphase/cavitatingFoam/createFields.H b/applications/solvers/multiphase/cavitatingFoam/createFields.H
index 0108b8c56a09800b8186768fd3c466c08069bc98..232d5675c5e372cf02f125628bcc0f6ce837ab89 100644
--- a/applications/solvers/multiphase/cavitatingFoam/createFields.H
+++ b/applications/solvers/multiphase/cavitatingFoam/createFields.H
@@ -43,7 +43,7 @@ volVectorField U
 
 #include "createPhi.H"
 
-mesh.setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(p.name());
 
 // Mass flux (corrected by rhoEqn.H)
 surfaceScalarField rhoPhi
@@ -86,8 +86,8 @@ rho == max
     rhoMin
 );
 
-mesh.setFluxRequired(p.name());
-mesh.setFluxRequired(rho.name());
+mesh.schemes().setFluxRequired(p.name());
+mesh.schemes().setFluxRequired(rho.name());
 
 // Create incompressible turbulence model
 autoPtr<incompressible::turbulenceModel> turbulence
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/alphaControls.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/alphaControls.H
index 2e94f34d752b7465f8eade4425fb76d9d84ab880..3c7ad707be89a57b99af4d173bc30c6f3c31341d 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/alphaControls.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/alphaControls.H
@@ -1,3 +1,3 @@
-const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+const dictionary& alphaControls = mesh.solution().solverDict(alpha1.name());
 
 label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/createFields.H
index 4725a8356da967e410e7c38bff04ab96a062c124..f8db5cefef7cf9f6d379632e1c6ba46ccccb9985 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/createFields.H
@@ -63,8 +63,8 @@ dimensionedScalar pMin
     mixture
 );
 
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha1.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(alpha1.name());
 
 
 #include "readGravitationalAcceleration.H"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index 9ac2034129adc6f0f6bc37a7654aea6150d63b0d..31c7a70c1f1e727f12f898224238fa32e732b195 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@@ -62,8 +62,8 @@ dimensionedScalar pMin
     mixture
 );
 
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha1.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(alpha1.name());
 
 
 #include "readGravitationalAcceleration.H"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/createFields.H
index 57ff76c20a2dd42d14c4ca6a650da6f5343a311c..cc78b9bd06c03e54a000de16ad273e51664153d6 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/createFields.H
@@ -62,8 +62,8 @@ dimensionedScalar pMin
     mixture
 );
 
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha1.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(alpha1.name());
 
 
 #include "readGravitationalAcceleration.H"
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
index 39930ed6386b0871bac395997926ba7f93473565..79e4142029077248977189892082e248a126d87d 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
@@ -45,7 +45,7 @@ volScalarField rho
 
 dimensionedScalar pMin("pMin", dimPressure, mixture);
 
-mesh.setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
 
 
 #include "readGravitationalAcceleration.H"
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
index 77d276973fd32bfd2c40bb27354cc811e8786e89..9f2dce8e69b2df92a8c967bd342ffe8cb8bf6abc 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
@@ -804,7 +804,7 @@ void Foam::multiphaseMixtureThermo::solve()
 {
     const Time& runTime = mesh_.time();
 
-    const dictionary& alphaControls = mesh_.solverDict("alpha");
+    const dictionary& alphaControls = mesh_.solution().solverDict("alpha");
     label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
     scalar cAlpha(alphaControls.get<scalar>("cAlpha"));
 
diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaControls.H b/applications/solvers/multiphase/driftFluxFoam/alphaControls.H
index 133df85af87d184c7876661eb7b40c79d1457a54..2f677b90dae34285a94e5252f0eaa57747b1c3fa 100644
--- a/applications/solvers/multiphase/driftFluxFoam/alphaControls.H
+++ b/applications/solvers/multiphase/driftFluxFoam/alphaControls.H
@@ -1,4 +1,4 @@
-const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+const dictionary& alphaControls = mesh.solution().solverDict(alpha1.name());
 
 label nAlphaCorr(alphaControls.get<label>("nAlphaCorr"));
 
diff --git a/applications/solvers/multiphase/driftFluxFoam/createFields.H b/applications/solvers/multiphase/driftFluxFoam/createFields.H
index d85ba792ea16536b1fc8a8ebd02cdf12bc53d310..80587c681d7ef8ba16ddb3f72eeb194b016194a4 100644
--- a/applications/solvers/multiphase/driftFluxFoam/createFields.H
+++ b/applications/solvers/multiphase/driftFluxFoam/createFields.H
@@ -125,8 +125,8 @@ if (p_rgh.needReference())
     p_rgh = p - rho*gh;
 }
 
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha1.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(alpha1.name());
 
 // MULES Correction
 tmp<surfaceScalarField> talphaPhiCorr0;
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H
index cc4e3c5fb09ce426887262e09ba81bb599ce8c5c..39f31ee936341e50e6801ece3c01dbeeb47415f6 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H
@@ -117,8 +117,8 @@ if (p_rgh.needReference())
     p_rgh = p - rho*gh;
 }
 
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha1.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(alpha1.name());
 
 #include "createMRF.H"
 #include "createFvOptions.H"
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index 9538dccf0da20d8000bf5c40ef49ded7b74313ce..58be476cfb15ca1aae783c6272dc6f57950dbebd 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -123,8 +123,8 @@ if (p_rgh.needReference())
     p_rgh = p - rho*gh;
 }
 
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha1.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(alpha1.name());
 
 #include "createMRF.H"
 #include "createFvOptions.H"
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaControls.H b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaControls.H
index cb83b15a9e01489d56c5acaecfd308677f5128e8..327bdcb72f208f8300941f2da4bff449b8ec9e7b 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaControls.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaControls.H
@@ -1,4 +1,4 @@
-const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+const dictionary& alphaControls = mesh.solution().solverDict(alpha1.name());
 
 label nAlphaCorr(alphaControls.get<label>("nAlphaCorr"));
 
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H b/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
index 7862915150a71f3e6bac50151d950a513a4a844e..12a422cf9db6bb34382fdfabb1acbea95cabb3ae 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
@@ -121,8 +121,8 @@ if (p_rgh.needReference())
     p_rgh = p - rho*gh;
 }
 
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha2.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(alpha2.name());
 
 #include "createMRF.H"
 #include "createFvOptions.H"
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C b/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
index 8ebb7ae852db244fefd4f9e8ee833e468e9002ec..5fe172764af52babc5038d95dffa1e35ee900ac7 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
@@ -157,7 +157,7 @@ Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties
     mixture_(mixture),
     cAlpha_
     (
-        mixture.U().mesh().solverDict
+        mixture.U().mesh().solution().solverDict
         (
             mixture_.alpha1().name()
         ).get<scalar>("cAlpha")
diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/correctPhi.H b/applications/solvers/multiphase/interFoam/overInterDyMFoam/correctPhi.H
index 7df3631b0402c5b0dbce0e90f8309201ad0a2736..125261e6735ce96bc6d4178c119d1810b7acb0e6 100644
--- a/applications/solvers/multiphase/interFoam/overInterDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/correctPhi.H
@@ -60,7 +60,7 @@
         fvc::makeAbsolute(phi, U);
     }
 
-    mesh.setFluxRequired(pcorr.name());
+    mesh.schemes().setFluxRequired(pcorr.name());
 
     dimensionedScalar rAUf("rAUf", dimTime/rho.dimensions(), 1.0);
 
diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/createFields.H b/applications/solvers/multiphase/interFoam/overInterDyMFoam/createFields.H
index 8cfb10551ec861f23b9fb08d3c27c964a723bff6..0c71eeb3ff2c11c9a55e11fde44a0b216f41efbd 100644
--- a/applications/solvers/multiphase/interFoam/overInterDyMFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/createFields.H
@@ -148,7 +148,7 @@ if (p_rgh.needReference())
     p_rgh = p - rho*gh;
 }
 
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha1.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(alpha1.name());
 
 #include "createMRF.H"
diff --git a/applications/solvers/multiphase/interIsoFoam/alphaControls.H b/applications/solvers/multiphase/interIsoFoam/alphaControls.H
index 2e94f34d752b7465f8eade4425fb76d9d84ab880..3c7ad707be89a57b99af4d173bc30c6f3c31341d 100644
--- a/applications/solvers/multiphase/interIsoFoam/alphaControls.H
+++ b/applications/solvers/multiphase/interIsoFoam/alphaControls.H
@@ -1,3 +1,3 @@
-const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+const dictionary& alphaControls = mesh.solution().solverDict(alpha1.name());
 
 label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
diff --git a/applications/solvers/multiphase/interIsoFoam/createFields.H b/applications/solvers/multiphase/interIsoFoam/createFields.H
index 1fdb0243d085b4c91a29265962a566dcd9602b15..2ca4efb2c5b1d96b32b552381a145928384f579c 100644
--- a/applications/solvers/multiphase/interIsoFoam/createFields.H
+++ b/applications/solvers/multiphase/interIsoFoam/createFields.H
@@ -125,8 +125,8 @@ if (p_rgh.needReference())
     p_rgh = p - rho*gh;
 }
 
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha1.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(alpha1.name());
 
 #include "createMRF.H"
 #include "createFvOptions.H"
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaControls.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaControls.H
index e9c832e58493eeb51ce995db56c6add44b75487d..8364c38f9c4619ec676da0d21be2d0493e25fb86 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaControls.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaControls.H
@@ -1,4 +1,4 @@
-const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+const dictionary& alphaControls = mesh.solution().solverDict(alpha1.name());
 
 label nAlphaCorr(alphaControls.get<label>("nAlphaCorr"));
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
index e1b84a64a68b5fb46d1fcbbcf4a28dd04f618f60..61125f4482c3e850da717feb3acdd13bd5164beb 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
@@ -105,8 +105,8 @@ if (p_rgh.needReference())
     p_rgh = p - rho*gh;
 }
 
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha1.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(alpha1.name());
 
 #include "createFvOptions.H"
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/createFields.H b/applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/createFields.H
index d560a9fcbafaf7e19c666f8a70f8222b8fe744aa..c60e356b44d9ad9ad0784e7c7ea91fbe3aa7fb9e 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/createFields.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/createFields.H
@@ -105,8 +105,8 @@ if (p_rgh.needReference())
     p_rgh = p - rho*gh;
 }
 
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha1.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(alpha1.name());
 
 #include "createFvOptions.H"
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H b/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
index 15acdcacb7af30eda99dab760e519e35377a0d8f..048674771a00778820edab8c22181cdb7022d46f 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
@@ -100,7 +100,7 @@ setRefCell
     pRefCell,
     pRefValue
 );
-mesh.setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
 
 
 #include "createMRFZones.H"
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/createFields.H b/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
index 201a4def06cf0393b74f24ffb3dfaee01dfd694d..aafcb82b8dcc4239daa55e92e3384545ae7c8e95 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
@@ -90,7 +90,7 @@ if (p_rgh.needReference())
     );
 }
 
-mesh.setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
 
 #include "createMRF.H"
 #include "createFvOptions.H"
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index 9c6ff5d331497837df56fe988709f44588509190..8a293352a7a95a400f9bef777fd352fdb3fb0816 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@@ -316,7 +316,7 @@ void Foam::multiphaseMixture::solve()
 
     volScalarField& alpha = phases_.first();
 
-    const dictionary& alphaControls = mesh_.solverDict("alpha");
+    const dictionary& alphaControls = mesh_.solution().solverDict("alpha");
     label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
     scalar cAlpha(alphaControls.get<scalar>("cAlpha"));
 
diff --git a/applications/solvers/multiphase/potentialFreeSurfaceFoam/createFields.H b/applications/solvers/multiphase/potentialFreeSurfaceFoam/createFields.H
index a9a7807949867b80a45bd7a6fb19d7f7efadfe7c..2dc9a6d13455f7df8ffa94b7e0f986f3daefe100 100644
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/createFields.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/createFields.H
@@ -74,7 +74,7 @@ p_gh =  p - (g & mesh.C());
 label p_ghRefCell = 0;
 scalar p_ghRefValue = 0.0;
 setRefCell(p_gh, pimple.dict(), p_ghRefCell, p_ghRefValue);
-mesh.setFluxRequired(p_gh.name());
+mesh.schemes().setFluxRequired(p_gh.name());
 
 #include "createMRF.H"
 #include "createFvOptions.H"
diff --git a/applications/solvers/multiphase/reactingMultiphaseEulerFoam/createFields.H b/applications/solvers/multiphase/reactingMultiphaseEulerFoam/createFields.H
index 50f1ed519137e93d0587aa0b452a67cfd68b923f..fc54711a6afa0687fc43bf5d62a033d7ee9892a6 100644
--- a/applications/solvers/multiphase/reactingMultiphaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/reactingMultiphaseEulerFoam/createFields.H
@@ -46,4 +46,4 @@ setRefCell
     pRefCell,
     pRefValue
 );
-mesh.setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createFields.H
index df7c20a798bbc3663575aee20c06e4058fdd2b20..a5bcdb94d1866ba298800baa98e096c2eab1ffb8 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createFields.H
@@ -45,4 +45,4 @@ setRefCell
     pRefCell,
     pRefValue
 );
-mesh.setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaControls.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaControls.H
index 2e94f34d752b7465f8eade4425fb76d9d84ab880..3c7ad707be89a57b99af4d173bc30c6f3c31341d 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaControls.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaControls.H
@@ -1,3 +1,3 @@
-const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+const dictionary& alphaControls = mesh.solution().solverDict(alpha1.name());
 
 label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
index c38f7e2959b7a6b7ab914d8017eaf09f6e7eef36..362062d1204c8e7f5befa4e2babdffb2e6e32530 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
@@ -110,5 +110,5 @@ if (p_rgh.needReference())
     p_rgh = p - rho*gh;
 }
 
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha1.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(alpha1.name());
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
index 1372ae309cd64891103a9f82918d3a833e51ad9a..6189e160139de6b1e6095b136f7b4ed02cc36da4 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
@@ -60,8 +60,8 @@ setRefCell
     pRefCell,
     pRefValue
 );
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha1.name());
+mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.schemes().setFluxRequired(alpha1.name());
 
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
index 2565aa271da5db658ea913919f9359107de1b6b8..6e08627ff91099e3c6895132dce7e5c538a42ed1 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
@@ -74,7 +74,7 @@ int main(int argc, char *argv[])
 
     bool implicitPhasePressure
     (
-        mesh.solverDict(alpha1.name()).getOrDefault
+        mesh.solution().solverDict(alpha1.name()).getOrDefault
         (
             "implicitPhasePressure", false
         )
diff --git a/applications/solvers/stressAnalysis/solidDisplacementFoam/createControl.H b/applications/solvers/stressAnalysis/solidDisplacementFoam/createControl.H
index 40e35f8c0593206d6ee6c308da26c494eede55ce..d1d225d4755fd8605f3450b48e93de9043ff0c29 100644
--- a/applications/solvers/stressAnalysis/solidDisplacementFoam/createControl.H
+++ b/applications/solvers/stressAnalysis/solidDisplacementFoam/createControl.H
@@ -1,3 +1,4 @@
-const dictionary& stressControl = mesh.solutionDict().subDict("stressAnalysis");
+const dictionary& stressControl =
+    mesh.solution().solutionDict("stressAnalysis");
 
 bool compactNormalStress(stressControl.get<bool>("compactNormalStress"));
diff --git a/applications/solvers/stressAnalysis/solidDisplacementFoam/createFields.H b/applications/solvers/stressAnalysis/solidDisplacementFoam/createFields.H
index 8d9ef025c81e69e0dabdf41c96a8a5105e317600..80387f0dd8abb5bb50102ff6c132a70ff51e0fd9 100644
--- a/applications/solvers/stressAnalysis/solidDisplacementFoam/createFields.H
+++ b/applications/solvers/stressAnalysis/solidDisplacementFoam/createFields.H
@@ -76,4 +76,4 @@ else
     divSigmaExp -= fvc::div((2*mu + lambda)*fvc::grad(D), "div(sigmaD)");
 }
 
-mesh.setFluxRequired(D.name());
+mesh.schemes().setFluxRequired(D.name());
diff --git a/applications/solvers/stressAnalysis/solidEquilibriumDisplacementFoam/createControls.H b/applications/solvers/stressAnalysis/solidEquilibriumDisplacementFoam/createControls.H
index 95d74370cc8b3da65fa219a56e6a4bc883040c5a..b1bc8560e1cb929f844314351f8fe0328626ba03 100644
--- a/applications/solvers/stressAnalysis/solidEquilibriumDisplacementFoam/createControls.H
+++ b/applications/solvers/stressAnalysis/solidEquilibriumDisplacementFoam/createControls.H
@@ -1 +1,2 @@
-const dictionary& stressControl = mesh.solutionDict().subDict("stressAnalysis");
+const dictionary& stressControl =
+    mesh.solution().solutionDict("stressAnalysis");
diff --git a/applications/test/parallelOverset/createFields.H b/applications/test/parallelOverset/createFields.H
index fc66a21b162a5710cc040d6118d2d3c78ba7010b..f4ec3455b723eee6e7f6bb7da29044da19314703 100644
--- a/applications/test/parallelOverset/createFields.H
+++ b/applications/test/parallelOverset/createFields.H
@@ -20,7 +20,7 @@
 
         const_cast<dictionary&>
         (
-            mesh.schemesDict()
+            mesh.schemes().dict()
         ).add
         (
             "oversetInterpolationRequired",
diff --git a/applications/test/reconstructedDistanceFunction/Test-reconstructedDistanceFunction.C b/applications/test/reconstructedDistanceFunction/Test-reconstructedDistanceFunction.C
index 55d421d0ce1b7d323066761335b38d5a7c3f5f7d..c6c9b3903647e15b629cc2dfccbd5861e27d2c89 100644
--- a/applications/test/reconstructedDistanceFunction/Test-reconstructedDistanceFunction.C
+++ b/applications/test/reconstructedDistanceFunction/Test-reconstructedDistanceFunction.C
@@ -76,7 +76,7 @@ int main(int argc, char *argv[])
     Info<< "create field phi\n" << endl;
     surfaceScalarField phi(fvc::interpolate(U) & mesh.Sf());
 
-    dictionary dict = mesh.solverDict(alpha1.name());
+    dictionary dict = mesh.solution().solverDict(alpha1.name());
 
     autoPtr<reconstructionSchemes> surf =
         reconstructionSchemes::New(alpha1,phi,U,dict);
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
index 048999ae4d24c036474bd4a1898497832cbc28a1..b1739cd538a5f08e813647c6eb23907770c8dd33 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
@@ -1160,7 +1160,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::relax()
 
     scalar relaxCoeff = 1;
 
-    if (this->mesh().relaxField(name, relaxCoeff))
+    if (this->mesh().solution().relaxField(name, relaxCoeff))
     {
         relax(relaxCoeff);
     }
diff --git a/src/dynamicFaMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C b/src/dynamicFaMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C
index ef431b1d9e40eadaef70b1ef70b0456a1b324d48..07fa8ac82af7eff7aef587b3f4e07ea32475c7ba 100644
--- a/src/dynamicFaMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C
+++ b/src/dynamicFaMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C
@@ -2112,10 +2112,7 @@ bool Foam::interfaceTrackingFvMesh::update()
 
         word ddtScheme
         (
-            mesh().ddtScheme
-            (
-                "ddt(" + U().name() + ')'
-            )
+            mesh().schemes().ddt("ddt(" + U().name() + ')')
         );
 
         if
diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrix.C b/src/finiteArea/faMatrices/faMatrix/faMatrix.C
index 55aa863bf71724c4fcecbb41c0c21b080b183645..975b89970bba7c965e91bdd5d9a7d20db02f4462 100644
--- a/src/finiteArea/faMatrices/faMatrix/faMatrix.C
+++ b/src/finiteArea/faMatrices/faMatrix/faMatrix.C
@@ -597,7 +597,7 @@ void Foam::faMatrix<Type>::relax()
 {
     scalar relaxCoeff = 0;
 
-    if (psi_.mesh().relaxEquation(psi_.name(), relaxCoeff))
+    if (psi_.mesh().solution().relaxEquation(psi_.name(), relaxCoeff))
     {
         relax(relaxCoeff);
     }
@@ -676,7 +676,7 @@ template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::faePatchField, Foam::edgeMesh>>
 Foam::faMatrix<Type>::flux() const
 {
-    if (!psi_.mesh().fluxRequired(psi_.name()))
+    if (!psi_.mesh().schemes().fluxRequired(psi_.name()))
     {
         FatalErrorInFunction
             << "flux requested but " << psi_.name()
@@ -749,14 +749,14 @@ const Foam::dictionary& Foam::faMatrix<Type>::solverDict
     const word& name
 ) const
 {
-    return psi_.mesh().solverDict(name);
+    return psi_.mesh().solution().solverDict(name);
 }
 
 
 template<class Type>
 const Foam::dictionary& Foam::faMatrix<Type>::solverDict() const
 {
-    return psi_.mesh().solverDict
+    return psi_.mesh().solution().solverDict
     (
         psi_.name()
     );
diff --git a/src/finiteArea/finiteArea/fac/facD2dt2.C b/src/finiteArea/finiteArea/fac/facD2dt2.C
index eb680c3b3f8a53f5decefc8f7b2885cecf403305..231e1ec415864f52fd825411afbc804800358786 100644
--- a/src/finiteArea/finiteArea/fac/facD2dt2.C
+++ b/src/finiteArea/finiteArea/fac/facD2dt2.C
@@ -51,7 +51,7 @@ tmp<GeometricField<Type, faPatchField, areaMesh>> d2dt2
     return fa::faD2dt2Scheme<Type>::New
     (
         mesh,
-        mesh.d2dt2Scheme("d2dt2(" + dt.name() + ')')
+        mesh.schemes().d2dt2("d2dt2(" + dt.name() + ')')
     ).ref().facD2dt2(dt);
 }
 
@@ -65,7 +65,7 @@ tmp<GeometricField<Type, faPatchField, areaMesh>> d2dt2
     return fa::faD2dt2Scheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().d2dt2Scheme("d2dt2(" + vf.name() + ')')
+        vf.mesh().schemes().d2dt2("d2dt2(" + vf.name() + ')')
     ).ref().facD2dt2(vf);
 }
 
@@ -80,7 +80,7 @@ tmp<GeometricField<Type, faPatchField, areaMesh>> d2dt2
     return fa::faD2dt2Scheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().d2dt2Scheme
+        vf.mesh().schemes().d2dt2
         (
             "d2dt2(" + rho.name() + ',' + vf.name() + ')'
         )
@@ -98,7 +98,7 @@ tmp<GeometricField<Type, faPatchField, areaMesh>> d2dt2
     return fa::faD2dt2Scheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().d2dt2Scheme
+        vf.mesh().schemes().d2dt2
         (
             "d2dt2(" + rho.name() + ',' + vf.name() + ')'
         )
diff --git a/src/finiteArea/finiteArea/fac/facDdt.C b/src/finiteArea/finiteArea/fac/facDdt.C
index fde27bdfeb532b1859b3b053a708665d61f1b1cd..fd5df0f64fdd040ab0ba4630a67daec941c97aba 100644
--- a/src/finiteArea/finiteArea/fac/facDdt.C
+++ b/src/finiteArea/finiteArea/fac/facDdt.C
@@ -52,7 +52,7 @@ ddt
     return fa::faDdtScheme<Type>::New
     (
         mesh,
-        mesh.ddtScheme("ddt(" + dt.name() + ')')
+        mesh.schemes().ddt("ddt(" + dt.name() + ')')
     ).ref().facDdt(dt);
 }
 
@@ -67,7 +67,7 @@ ddt
     return fa::faDdtScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme("ddt(" + vf.name() + ')')
+        vf.mesh().schemes().ddt("ddt(" + vf.name() + ')')
     ).ref().facDdt(vf);
 }
 
@@ -83,7 +83,7 @@ ddt
     return fa::faDdtScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme
+        vf.mesh().schemes().ddt
         (
             "ddt(" + rho.name() + ',' + vf.name() + ')'
         )
@@ -102,7 +102,7 @@ ddt
     return fa::faDdtScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme
+        vf.mesh().schemes().ddt
         (
             "ddt(" + rho.name() + ',' + vf.name() + ')'
         )
diff --git a/src/finiteArea/finiteArea/fac/facDiv.C b/src/finiteArea/finiteArea/fac/facDiv.C
index 12b95407ba361652d0f695d5861fcb5312d707a1..53a30008af9e0e291f65c9792170078d9c566b9a 100644
--- a/src/finiteArea/finiteArea/fac/facDiv.C
+++ b/src/finiteArea/finiteArea/fac/facDiv.C
@@ -107,7 +107,8 @@ div
     (
         fa::divScheme<Type>::New
         (
-            vf.mesh(), vf.mesh().divScheme(name)
+            vf.mesh(),
+            vf.mesh().schemes().div(name)
         ).ref().facDiv(vf)
     );
     GeometricField
@@ -207,7 +208,7 @@ div
         (
             vf.mesh(),
             flux,
-            vf.mesh().divScheme(name)
+            vf.mesh().schemes().div(name)
         ).ref().facDiv(flux, vf)
     );
     GeometricField<Type, faPatchField, areaMesh>& Div = tDiv.ref();
diff --git a/src/finiteArea/finiteArea/fac/facGrad.C b/src/finiteArea/finiteArea/fac/facGrad.C
index 262b20e1299e8de50444aca3e503455727cfd154..72cab10d473df35415c12a19cfe37c1a5e45145f 100644
--- a/src/finiteArea/finiteArea/fac/facGrad.C
+++ b/src/finiteArea/finiteArea/fac/facGrad.C
@@ -108,7 +108,7 @@ grad
         fa::gradScheme<Type>::New
         (
             vf.mesh(),
-            vf.mesh().gradScheme(name)
+            vf.mesh().schemes().grad(name)
         )
         .cref()  // const ref to tmp contents
         .grad(vf).ptr()  // steal ptr or deep copy of cached gradient
diff --git a/src/finiteArea/finiteArea/fac/facLaplacian.C b/src/finiteArea/finiteArea/fac/facLaplacian.C
index 97a0e9f3e2bedc649542e5762c8f48be0596605c..53d38640d7c1b4179b5ec78623cf0250731cd06f 100644
--- a/src/finiteArea/finiteArea/fac/facLaplacian.C
+++ b/src/finiteArea/finiteArea/fac/facLaplacian.C
@@ -52,7 +52,7 @@ laplacian
     return fa::laplacianScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().laplacianScheme(name)
+        vf.mesh().schemes().laplacian(name)
     ).ref().facLaplacian(vf);
 }
 
@@ -172,7 +172,7 @@ laplacian
     return fa::laplacianScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().laplacianScheme(name)
+        vf.mesh().schemes().laplacian(name)
     ).ref().facLaplacian(gamma, vf);
 }
 
@@ -306,7 +306,7 @@ laplacian
     return fa::laplacianScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().laplacianScheme(name)
+        vf.mesh().schemes().laplacian(name)
     ).ref().facLaplacian(gamma, vf);
 }
 
diff --git a/src/finiteArea/finiteArea/fac/facLnGrad.C b/src/finiteArea/finiteArea/fac/facLnGrad.C
index 951d6bb0e274f9c09e5a005ee43bb700fe2f0b5d..f1fb9e4add5793aca34c09acd27128b91145ad6b 100644
--- a/src/finiteArea/finiteArea/fac/facLnGrad.C
+++ b/src/finiteArea/finiteArea/fac/facLnGrad.C
@@ -52,7 +52,7 @@ lnGrad
     return fa::lnGradScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().lnGradScheme(name)
+        vf.mesh().schemes().lnGrad(name)
     ).ref().lnGrad(vf);
 }
 
diff --git a/src/finiteArea/finiteArea/fac/facNDiv.C b/src/finiteArea/finiteArea/fac/facNDiv.C
index e0956bac0bc3494ccd176e268b2d49fddcad458e..955629209187ec34e710f7fd25244ce44832e76c 100644
--- a/src/finiteArea/finiteArea/fac/facNDiv.C
+++ b/src/finiteArea/finiteArea/fac/facNDiv.C
@@ -97,9 +97,10 @@ ndiv
     tmp<GeometricField<Type, faPatchField, areaMesh>> tDiv
     (
         fa::divScheme<Type>::New
-            (
-                vf.mesh(), vf.mesh().divScheme(name)
-            ).ref().facDiv(vf)
+        (
+            vf.mesh(),
+            vf.mesh().schemes().div(name)
+        ).ref().facDiv(vf)
     );
     GeometricField<Type, faPatchField, areaMesh>& Div = tDiv.ref();
 
@@ -192,7 +193,7 @@ ndiv
         (
             vf.mesh(),
             flux,
-            vf.mesh().divScheme(name)
+            vf.mesh().schemes().div(name)
         ).ref().facDiv(flux, vf)
     );
 
diff --git a/src/finiteArea/finiteArea/fac/facNGrad.C b/src/finiteArea/finiteArea/fac/facNGrad.C
index fe7f572e6592fbef292696834926275fef342063..6cf70146aadde306c401300afc158ed71f02e8fc 100644
--- a/src/finiteArea/finiteArea/fac/facNGrad.C
+++ b/src/finiteArea/finiteArea/fac/facNGrad.C
@@ -118,7 +118,7 @@ ngrad
         fa::gradScheme<Type>::New
         (
             vf.mesh(),
-            vf.mesh().gradScheme(name)
+            vf.mesh().schemes().grad(name)
         ).ref().grad(vf);
 
     GeometricField<GradType, faPatchField, areaMesh>& gGrad = tgGrad.ref();
diff --git a/src/finiteArea/finiteArea/fam/famD2dt2.C b/src/finiteArea/finiteArea/fam/famD2dt2.C
index 3743c5953aaf58958db8db7a5d141487cd78f581..9a2191033e67a0b0bf8e0e880c88cd4fecfda11d 100644
--- a/src/finiteArea/finiteArea/fam/famD2dt2.C
+++ b/src/finiteArea/finiteArea/fam/famD2dt2.C
@@ -51,7 +51,7 @@ tmp<faMatrix<Type>> d2dt2
     return fa::faD2dt2Scheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().d2dt2Scheme("d2dt2(" + vf.name() + ')')
+        vf.mesh().schemes().d2dt2("d2dt2(" + vf.name() + ')')
     ).ref().famD2dt2(vf);
 }
 
@@ -66,7 +66,7 @@ tmp<faMatrix<Type>> d2dt2
     return fa::faD2dt2Scheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().d2dt2Scheme
+        vf.mesh().schemes().d2dt2
         (
             "d2dt2(" + rho.name() + ',' + vf.name() + ')'
         )
@@ -84,7 +84,7 @@ tmp<faMatrix<Type>> d2dt2
     return fa::faD2dt2Scheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme
+        vf.mesh().schemes().d2dt2
         (
             "d2dt2(" + rho.name() + ',' + vf.name() + ')'
         )
diff --git a/src/finiteArea/finiteArea/fam/famDdt.C b/src/finiteArea/finiteArea/fam/famDdt.C
index 94ec07e6f2d1a0e16653ff8f68091f4e7896e057..13dde85af3042eb08f2ddb9cdf2175b690c4ec80 100644
--- a/src/finiteArea/finiteArea/fam/famDdt.C
+++ b/src/finiteArea/finiteArea/fam/famDdt.C
@@ -52,7 +52,7 @@ ddt
     return fa::faDdtScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme("ddt(" + vf.name() + ')')
+        vf.mesh().schemes().ddt("ddt(" + vf.name() + ')')
     ).ref().famDdt(vf);
 }
 
@@ -68,7 +68,7 @@ ddt
     return fa::faDdtScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme
+        vf.mesh().schemes().ddt
         (
             "ddt(" + rho.name() + ',' + vf.name() + ')'
         )
@@ -87,7 +87,7 @@ ddt
     return fa::faDdtScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme
+        vf.mesh().schemes().ddt
         (
             "ddt(" + rho.name() + ',' + vf.name() + ')'
         )
diff --git a/src/finiteArea/finiteArea/fam/famDiv.C b/src/finiteArea/finiteArea/fam/famDiv.C
index db4f4a449ffcc4aa4daa2f37c9724c11e7a49eaf..bd4e8b07a5df8b73c932ccadbb53d9660ee9163d 100644
--- a/src/finiteArea/finiteArea/fam/famDiv.C
+++ b/src/finiteArea/finiteArea/fam/famDiv.C
@@ -59,7 +59,7 @@ div
         (
             vf.mesh(),
             flux,
-            vf.mesh().divScheme(name)
+            vf.mesh().schemes().div(name)
         ).ref().famDiv(flux, vf)
     );
     faMatrix<Type>& M = tM.ref();
@@ -70,7 +70,7 @@ div
         (
             vf.mesh(),
             flux,
-            vf.mesh().divScheme(name)
+            vf.mesh().schemes().div(name)
         ).ref().facDiv(flux, vf)
     );
 
diff --git a/src/finiteArea/finiteArea/fam/famLaplacian.C b/src/finiteArea/finiteArea/fam/famLaplacian.C
index cd2b4687baf2283bd79fc92dbfffc7ca2ca006f5..2eae375945e9395d35901121d4a8a4dc85aed710 100644
--- a/src/finiteArea/finiteArea/fam/famLaplacian.C
+++ b/src/finiteArea/finiteArea/fam/famLaplacian.C
@@ -172,7 +172,7 @@ laplacian
     return fa::laplacianScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().laplacianScheme(name)
+        vf.mesh().schemes().laplacian(name)
     ).ref().famLaplacian(gamma, vf);
 }
 
@@ -218,7 +218,7 @@ laplacian
     return fa::laplacianScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().laplacianScheme(name)
+        vf.mesh().schemes().laplacian(name)
     ).ref().famLaplacian(gamma, vf);
 }
 
diff --git a/src/finiteArea/finiteArea/fam/famNDiv.C b/src/finiteArea/finiteArea/fam/famNDiv.C
index d0e359cfd4ab3318b0e9ceb6568d51877e39809a..b61baed963ca36e715f985b938b4d8f3d4c002cf 100644
--- a/src/finiteArea/finiteArea/fam/famNDiv.C
+++ b/src/finiteArea/finiteArea/fam/famNDiv.C
@@ -55,7 +55,7 @@ ndiv
     (
         vf.mesh(),
         flux,
-        vf.mesh().divScheme(name)
+        vf.mesh().schemes().div(name)
     ).ref().famDiv(flux, vf);//TODO calculate normal
 }
 
diff --git a/src/finiteArea/finiteArea/fam/vectorFamDiv.C b/src/finiteArea/finiteArea/fam/vectorFamDiv.C
index 337c47f9481d91d3d410f5d1570840ddaab53ce8..ee0e2ec038380b4237c25ef965cb4e71030c1084 100644
--- a/src/finiteArea/finiteArea/fam/vectorFamDiv.C
+++ b/src/finiteArea/finiteArea/fam/vectorFamDiv.C
@@ -55,7 +55,7 @@ div
     (
         vf.mesh(),
         flux,
-        vf.mesh().divScheme(name)
+        vf.mesh().schemes().div(name)
     ).ref().famDiv(flux, vf);
 }
 
diff --git a/src/finiteArea/finiteArea/gradSchemes/faGradScheme/faGradScheme.C b/src/finiteArea/finiteArea/gradSchemes/faGradScheme/faGradScheme.C
index 597355f7147a9c85d94e382339238fbefcf1b15c..8b640be8068b935087c7f8d1dcb04ef38753c40e 100644
--- a/src/finiteArea/finiteArea/gradSchemes/faGradScheme/faGradScheme.C
+++ b/src/finiteArea/finiteArea/gradSchemes/faGradScheme/faGradScheme.C
@@ -109,7 +109,7 @@ gradScheme<Type>::grad
     GradFieldType* pgGrad =
         mesh().thisDb().template getObjectPtr<GradFieldType>(name);
 
-    if (!this->mesh().cache(name)) // || this->mesh().changing()
+    if (!this->mesh().solution().cache(name)) // || this->mesh().changing()
     {
         // Delete any old occurrences to avoid double registration
         if (pgGrad && pgGrad->ownedByRegistry())
diff --git a/src/finiteArea/finiteArea/laplacianSchemes/gaussFaLaplacianScheme/gaussFaLaplacianScheme.C b/src/finiteArea/finiteArea/laplacianSchemes/gaussFaLaplacianScheme/gaussFaLaplacianScheme.C
index 46e7d4bb4f3502bb177523dc985d30905cc2b5d1..898914c4fde76180b4d615a6e894b9bff4edf441 100644
--- a/src/finiteArea/finiteArea/laplacianSchemes/gaussFaLaplacianScheme/gaussFaLaplacianScheme.C
+++ b/src/finiteArea/finiteArea/laplacianSchemes/gaussFaLaplacianScheme/gaussFaLaplacianScheme.C
@@ -80,7 +80,7 @@ gaussLaplacianScheme<Type>::famLaplacian
 
     if (this->tlnGradScheme_().corrected())
     {
-        if (this->mesh().fluxRequired(vf.name()))
+        if (this->mesh().schemes().fluxRequired(vf.name()))
         {
             fam.faceFluxCorrectionPtr() = new
             GeometricField<Type, faePatchField, edgeMesh>
diff --git a/src/finiteArea/finiteArea/lnGradSchemes/correctedLnGrad/correctedLnGrad.C b/src/finiteArea/finiteArea/lnGradSchemes/correctedLnGrad/correctedLnGrad.C
index 808fd89648f635540ed85a66f4abf2ff8dac62f1..7afd4581acad234f84a90934a9c3c5f248937dd1 100644
--- a/src/finiteArea/finiteArea/lnGradSchemes/correctedLnGrad/correctedLnGrad.C
+++ b/src/finiteArea/finiteArea/lnGradSchemes/correctedLnGrad/correctedLnGrad.C
@@ -90,7 +90,7 @@ correctedLnGrad<Type>::correction
                 gradScheme<typename pTraits<Type>::cmptType>::New
                 (
                     mesh,
-                    mesh.gradScheme(ssf.name())
+                    mesh.schemes().grad(ssf.name())
                 )()
                .grad(vf.component(cmpt))
             )
diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolate.C b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolate.C
index 54dcd271b1460ce156813086109b87c999756221..6808171353aeaf112e3c866cea9165c1310fdda3 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolate.C
+++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolate.C
@@ -56,7 +56,7 @@ Foam::tmp<Foam::edgeInterpolationScheme<Type>> Foam::fac::scheme
     (
         faceFlux.mesh(),
         faceFlux,
-        faceFlux.mesh().interpolationScheme(name)
+        faceFlux.mesh().schemes().interpolation(name)
     );
 }
 
@@ -86,7 +86,7 @@ Foam::tmp<Foam::edgeInterpolationScheme<Type>> Foam::fac::scheme
     return edgeInterpolationScheme<Type>::New
     (
         mesh,
-        mesh.interpolationScheme(name)
+        mesh.schemes().interpolation(name)
     );
 }
 
diff --git a/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C b/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C
index 06dcfb73077c30d985ee177facc5a52c839a1187..3f35996652ee6de637844c8246b3e2a143904824 100644
--- a/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C
+++ b/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C
@@ -92,7 +92,7 @@ void Foam::CorrectPhi
         fvc::makeAbsolute(phi, U);
     }
 
-    mesh.setFluxRequired(pcorr.name());
+    mesh.schemes().setFluxRequired(pcorr.name());
 
     while (pimple.correctNonOrthogonal())
     {
@@ -165,7 +165,7 @@ void Foam::CorrectPhi
         pcorrTypes
     );
 
-    mesh.setFluxRequired(pcorr.name());
+    mesh.schemes().setFluxRequired(pcorr.name());
 
     while (pimple.correctNonOrthogonal())
     {
diff --git a/src/finiteVolume/cfdTools/general/SRF/derivedFvPatchFields/SRFFreestreamVelocityFvPatchVectorField/SRFFreestreamVelocityFvPatchVectorField.C b/src/finiteVolume/cfdTools/general/SRF/derivedFvPatchFields/SRFFreestreamVelocityFvPatchVectorField/SRFFreestreamVelocityFvPatchVectorField.C
index b2f876dffe942b4f9806518a3801f94883f51e79..ea30b0a113bb407b8b5e494c3a193ccf8ba17bae 100644
--- a/src/finiteVolume/cfdTools/general/SRF/derivedFvPatchFields/SRFFreestreamVelocityFvPatchVectorField/SRFFreestreamVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/cfdTools/general/SRF/derivedFvPatchFields/SRFFreestreamVelocityFvPatchVectorField/SRFFreestreamVelocityFvPatchVectorField.C
@@ -123,7 +123,7 @@ void Foam::SRFFreestreamVelocityFvPatchVectorField::updateCoeffs()
     word ddtScheme
     (
         this->internalField().mesh()
-       .ddtScheme(this->internalField().name())
+       .schemes().ddt(this->internalField().name())
     );
 
     if (ddtScheme == fv::steadyStateDdtScheme<scalar>::typeName)
diff --git a/src/finiteVolume/cfdTools/general/include/alphaControls.H b/src/finiteVolume/cfdTools/general/include/alphaControls.H
index 8b2179c1f7fe836e17d2c928c71ff0d405497c39..514ad80fe1ef9cb759fe29db66ae6b037515b88f 100644
--- a/src/finiteVolume/cfdTools/general/include/alphaControls.H
+++ b/src/finiteVolume/cfdTools/general/include/alphaControls.H
@@ -1,4 +1,4 @@
-const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+const dictionary& alphaControls = mesh.solution().solverDict(alpha1.name());
 
 label nAlphaCorr(alphaControls.get<label>("nAlphaCorr"));
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
index a1630c943428dc32e049df98bb458317cec15ed6..a6c36c37d93345846be6687b3d74a15e28c1e3ca 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
@@ -299,7 +299,7 @@ Foam::solutionControl::solutionControl(fvMesh& mesh, const word& algorithmName)
 
 const Foam::dictionary Foam::solutionControl::dict() const
 {
-    return mesh_.solutionDict().subOrEmptyDict(algorithmName_);
+    return mesh_.solution().dict().subOrEmptyDict(algorithmName_);
 }
 
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlTemplates.C b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlTemplates.C
index eac74386367af7c3a66601b14b6c55ff29eb4835..78baa622ab8968e73734ac8886e09223a46a655f 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlTemplates.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlTemplates.C
@@ -41,7 +41,11 @@ void Foam::solutionControl::storePrevIter() const
     {
         const word& fldName = fld.name();
 
-        if (!fldName.contains("PrevIter") && mesh_.relaxField(fldName))
+        if
+        (
+            !fldName.contains("PrevIter")
+         && mesh_.solution().relaxField(fldName)
+        )
         {
             DebugInfo
                 << algorithmName_ << ": storing previous iter for "
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C
index 0d6dac5c1564422a277a4c5aa3fd942835817a0d..c0a529d3e13536c1184bbfeb47d87b8c602bc097 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C
@@ -295,7 +295,13 @@ void Foam::cyclicFvPatchField<Type>::manipulateMatrix
             }
         }
 
-        if (matrix.psi(mat).mesh().fluxRequired(this->internalField().name()))
+        if
+        (
+            matrix.psi(mat).mesh().schemes().fluxRequired
+            (
+                this->internalField().name()
+            )
+        )
         {
             matrix.internalCoeffs().set
             (
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C
index 99a12cbc27d810603f79fba23b0e7930d10bac89..fa4c236934ce938d004c49323e4e80c6c7e59315 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C
@@ -844,7 +844,13 @@ void Foam::cyclicACMIFvPatchField<Type>::manipulateMatrix
         // Set internalCoeffs and boundaryCoeffs in the assembly matrix
         // on cyclicACMI patches to be used in the individual matrix by
         // matrix.flux()
-        if (matrix.psi(mat).mesh().fluxRequired(this->internalField().name()))
+        if
+        (
+            matrix.psi(mat).mesh().schemes().fluxRequired
+            (
+                this->internalField().name()
+            )
+        )
         {
             matrix.internalCoeffs().set
             (
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.C
index 9f3d146b1c9e2198d8747f2ceec3692f2069a49e..f4dded80117677c7a2d281ed96214453f72ac928 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.C
@@ -758,7 +758,13 @@ void Foam::cyclicAMIFvPatchField<Type>::manipulateMatrix
         // Set internalCoeffs and boundaryCoeffs in the assembly matrix
         // on clyclicAMI patches to be used in the individual matrix by
         // matrix.flux()
-        if (matrix.psi(mat).mesh().fluxRequired(this->internalField().name()))
+        if
+        (
+            matrix.psi(mat).mesh().schemes().fluxRequired
+            (
+                this->internalField().name()
+            )
+        )
         {
             matrix.internalCoeffs().set
             (
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/advective/advectiveFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/advective/advectiveFvPatchField.C
index 250b46a0eace0db06f83dc41d58662dfe18cece6..1a56bdfe86f9a1cfcb635327583fe959ee770551 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/advective/advectiveFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/advective/advectiveFvPatchField.C
@@ -179,7 +179,7 @@ void Foam::advectiveFvPatchField<Type>::updateCoeffs()
 
     word ddtScheme
     (
-        mesh.ddtScheme(this->internalField().name())
+        mesh.schemes().ddt(this->internalField().name())
     );
     scalar deltaT = this->db().time().deltaTValue();
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/pressurePIDControlInletVelocity/pressurePIDControlInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/pressurePIDControlInletVelocity/pressurePIDControlInletVelocityFvPatchVectorField.C
index 84c265f9908fd50d927625ef42c523b3c61a65b5..d5af9523cba95ae48f8323cf46e71b3294b6c3c0 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/pressurePIDControlInletVelocity/pressurePIDControlInletVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/pressurePIDControlInletVelocity/pressurePIDControlInletVelocityFvPatchVectorField.C
@@ -299,7 +299,7 @@ void Foam::pressurePIDControlInletVelocityFvPatchVectorField::updateCoeffs()
     // Target and measured flow rates
     scalar QTarget, QMeasured;
     const scalar a = (1/sqr(Ab) - 1/sqr(Aa))/(2*rho);
-    if (!mesh.steady() && db().foundObject<volScalarField>(pName_))
+    if (!mesh.schemes().steady() && db().foundObject<volScalarField>(pName_))
     {
         const scalar b = LbyA/deltaT;
         const scalar c = - LbyA/deltaT*oldQ_ /* - deltaP */;
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
index 0eb44f87fa69ddbc73fe328411f22f8f375b89f0..81406d2b52dac9e5bc55ef7015fcb7c161a062c6 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
@@ -151,7 +151,7 @@ void Foam::waveSurfacePressureFvPatchScalarField::updateCoeffs()
     vectorField& zetap = zeta.boundaryFieldRef()[patchi];
 
     // Lookup d/dt scheme from database for zeta
-    const word ddtSchemeName(zeta.mesh().ddtScheme(zeta.name()));
+    const word ddtSchemeName(zeta.mesh().schemes().ddt(zeta.name()));
     ddtSchemeType ddtScheme(ddtSchemeTypeNames_[ddtSchemeName]);
 
     // Retrieve the flux field from the database
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdt.C b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdt.C
index 77542fdddc6d43950cf8fa8ce03e733b174db537..392075572b86aab18ff9e89e108f26c870e97c56 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdt.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdt.C
@@ -39,7 +39,7 @@ Foam::word Foam::fv::localEulerDdt::rSubDeltaTName("rSubDeltaT");
 bool Foam::fv::localEulerDdt::enabled(const fvMesh& mesh)
 {
     return
-        word(mesh.ddtScheme("default"))
+        word(mesh.schemes().ddt("default"))
      == fv::localEulerDdtScheme<scalar>::typeName;
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcCorrectAlpha.C b/src/finiteVolume/finiteVolume/fvc/fvcCorrectAlpha.C
index 2f0feefc64bf66ed6338163301e2843a933433c9..cfef3413c926592e221c47a8bc2408452c2c8a12 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcCorrectAlpha.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcCorrectAlpha.C
@@ -61,7 +61,7 @@ tmp<GeometricField<scalar, fvsPatchField, surfaceMesh>> alphaCorr
     const word fieldName = U.select(finalIter);
 
     scalar alpha = 1;
-    mesh.relaxEquation(fieldName, alpha);
+    mesh.solution().relaxEquation(fieldName, alpha);
 
     return
         (1 - alpha)
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcDdt.C b/src/finiteVolume/finiteVolume/fvc/fvcDdt.C
index 2ae8cc422a17f426173a717b68b81716b52a597f..40eec4f43b29dea27e0d4938ddda3520837affbe 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcDdt.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcDdt.C
@@ -52,7 +52,7 @@ ddt
     return fv::ddtScheme<Type>::New
     (
         mesh,
-        mesh.ddtScheme("ddt(" + dt.name() + ')')
+        mesh.schemes().ddt("ddt(" + dt.name() + ')')
     ).ref().fvcDdt(dt);
 }
 
@@ -67,7 +67,7 @@ ddt
     return fv::ddtScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme("ddt(" + vf.name() + ')')
+        vf.mesh().schemes().ddt("ddt(" + vf.name() + ')')
     ).ref().fvcDdt(vf);
 }
 
@@ -83,7 +83,7 @@ ddt
     return fv::ddtScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme("ddt(" + rho.name() + ',' + vf.name() + ')')
+        vf.mesh().schemes().ddt("ddt(" + rho.name() + ',' + vf.name() + ')')
     ).ref().fvcDdt(rho, vf);
 }
 
@@ -99,7 +99,7 @@ ddt
     return fv::ddtScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme("ddt(" + rho.name() + ',' + vf.name() + ')')
+        vf.mesh().schemes().ddt("ddt(" + rho.name() + ',' + vf.name() + ')')
     ).ref().fvcDdt(rho, vf);
 }
 
@@ -116,7 +116,7 @@ ddt
     return fv::ddtScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme
+        vf.mesh().schemes().ddt
         (
             "ddt("
           + alpha.name() + ','
@@ -137,7 +137,7 @@ ddt
     return fv::ddtScheme<Type>::New
     (
         sf.mesh(),
-        sf.mesh().ddtScheme("ddt(" + sf.name() + ')')
+        sf.mesh().schemes().ddt("ddt(" + sf.name() + ')')
     ).ref().fvcDdt(sf);
 }
 
@@ -177,7 +177,7 @@ ddtCorr
     return fv::ddtScheme<Type>::New
     (
         U.mesh(),
-        U.mesh().ddtScheme("ddt(" + U.name() + ')')
+        U.mesh().schemes().ddt("ddt(" + U.name() + ')')
     ).ref().fvcDdtUfCorr(U, Uf);
 }
 
@@ -198,7 +198,7 @@ ddtCorr
     return fv::ddtScheme<Type>::New
     (
         U.mesh(),
-        U.mesh().ddtScheme("ddt(" + U.name() + ')')
+        U.mesh().schemes().ddt("ddt(" + U.name() + ')')
     ).ref().fvcDdtPhiCorr(U, phi);
 }
 
@@ -240,7 +240,7 @@ ddtCorr
     return fv::ddtScheme<Type>::New
     (
         U.mesh(),
-        U.mesh().ddtScheme("ddt(" + U.name() + ')')
+        U.mesh().schemes().ddt("ddt(" + U.name() + ')')
     ).ref().fvcDdtUfCorr(rho, U, Uf);
 }
 
@@ -262,7 +262,7 @@ ddtCorr
     return fv::ddtScheme<Type>::New
     (
         U.mesh(),
-        U.mesh().ddtScheme("ddt(" + rho.name() + ',' + U.name() + ')')
+        U.mesh().schemes().ddt("ddt(" + rho.name() + ',' + U.name() + ')')
     ).ref().fvcDdtPhiCorr(rho, U, phi);
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcDiv.C b/src/finiteVolume/finiteVolume/fvc/fvcDiv.C
index 1424f643f5c0ab0c79c2d14dc18b347e8a199980..09aefc911fec0342031bb5488cd51b6d9c3cc407 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcDiv.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcDiv.C
@@ -90,7 +90,7 @@ div
 {
     return fv::divScheme<Type>::New
     (
-        vf.mesh(), vf.mesh().divScheme(name)
+        vf.mesh(), vf.mesh().schemes().div(name)
     ).ref().fvcDiv(vf);
 }
 
@@ -168,7 +168,7 @@ div
     (
         vf.mesh(),
         flux,
-        vf.mesh().divScheme(name)
+        vf.mesh().schemes().div(name)
     ).ref().fvcDiv(flux, vf);
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcFluxTemplates.C b/src/finiteVolume/finiteVolume/fvc/fvcFluxTemplates.C
index 16aeb5b405411a6381da2cdd94614f3a41766f79..ea8ce76476d1803713ac68d309fa80a41d56ea77 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcFluxTemplates.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcFluxTemplates.C
@@ -54,7 +54,7 @@ flux
     (
         vf.mesh(),
         phi,
-        vf.mesh().divScheme(name)
+        vf.mesh().schemes().div(name)
     )().flux(phi, vf);
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcGrad.C b/src/finiteVolume/finiteVolume/fvc/fvcGrad.C
index 287888ad76020699b34793750b74dcc1017bc06b..ab5f927a75a4c635f6333dd3996c225eb2b180c9 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcGrad.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcGrad.C
@@ -99,7 +99,7 @@ grad
     return fv::gradScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().gradScheme(name)
+        vf.mesh().schemes().grad(name)
     )().grad(vf, name);
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcLaplacian.C b/src/finiteVolume/finiteVolume/fvc/fvcLaplacian.C
index b99bb59d60c137f0f4767271672c50228307fb4a..76754cbb53e801d80084b76d0d3c01f139ca68e2 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcLaplacian.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcLaplacian.C
@@ -52,7 +52,7 @@ laplacian
     return fv::laplacianScheme<Type, scalar>::New
     (
         vf.mesh(),
-        vf.mesh().laplacianScheme(name)
+        vf.mesh().schemes().laplacian(name)
     ).ref().fvcLaplacian(vf);
 }
 
@@ -203,7 +203,7 @@ laplacian
     return fv::laplacianScheme<Type, GType>::New
     (
         vf.mesh(),
-        vf.mesh().laplacianScheme(name)
+        vf.mesh().schemes().laplacian(name)
     ).ref().fvcLaplacian(gamma, vf);
 }
 
@@ -345,7 +345,7 @@ laplacian
     return fv::laplacianScheme<Type, GType>::New
     (
         vf.mesh(),
-        vf.mesh().laplacianScheme(name)
+        vf.mesh().schemes().laplacian(name)
     ).ref().fvcLaplacian(gamma, vf);
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C b/src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C
index 42becab798ab12db3637d704470656a79281f794..7001d74c9f7eb591e50545792165fd31ab4a34bb 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C
@@ -40,7 +40,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::fvc::meshPhi
     return fv::ddtScheme<vector>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme("ddt(" + vf.name() + ')')
+        vf.mesh().schemes().ddt("ddt(" + vf.name() + ')')
     ).ref().meshPhi(vf);
 }
 
@@ -54,7 +54,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::fvc::meshPhi
     return fv::ddtScheme<vector>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme("ddt(" + rho.name() + ',' + vf.name() + ')')
+        vf.mesh().schemes().ddt("ddt(" + rho.name() + ',' + vf.name() + ')')
     ).ref().meshPhi(vf);
 }
 
@@ -68,7 +68,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::fvc::meshPhi
     return fv::ddtScheme<vector>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme("ddt(" + rho.name() + ',' + vf.name() + ')')
+        vf.mesh().schemes().ddt("ddt(" + rho.name() + ',' + vf.name() + ')')
     ).ref().meshPhi(vf);
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSnGrad.C b/src/finiteVolume/finiteVolume/fvc/fvcSnGrad.C
index e373cdcb2e72903cc51430c945c347b8c3599283..4dc84a50433052a1c5772d0da5fd265d93ac582f 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcSnGrad.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSnGrad.C
@@ -52,7 +52,7 @@ snGrad
     return fv::snGradScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().snGradScheme(name)
+        vf.mesh().schemes().snGrad(name)
     )().snGrad(vf);
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.C b/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.C
index a0d64a120c9afffa6b0962eb5dd220c663e01b7d..666a5960484cf8f0e587a51c81ea3b7ed7ed723c 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.C
@@ -52,7 +52,7 @@ d2dt2
     return fv::d2dt2Scheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().d2dt2Scheme("d2dt2(" + vf.name() + ')')
+        vf.mesh().schemes().d2dt2("d2dt2(" + vf.name() + ')')
     ).ref().fvmD2dt2(vf);
 }
 
@@ -68,7 +68,7 @@ d2dt2
     return fv::d2dt2Scheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().d2dt2Scheme("d2dt2(" + rho.name() + ',' + vf.name() + ')')
+        vf.mesh().schemes().d2dt2("d2dt2(" + rho.name() + ',' + vf.name() + ')')
     ).ref().fvmD2dt2(rho, vf);
 }
 
@@ -84,7 +84,7 @@ d2dt2
     return fv::d2dt2Scheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().d2dt2Scheme("d2dt2(" + rho.name() + ',' + vf.name() + ')')
+        vf.mesh().schemes().d2dt2("d2dt2(" + rho.name() + ',' + vf.name() + ')')
     ).ref().fvmD2dt2(rho, vf);
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDdt.C b/src/finiteVolume/finiteVolume/fvm/fvmDdt.C
index 2c76dca9d031d6443382dbe46d8aaf68f6743628..eac67e6a3cef152dd41f753cadd4290a28ff5425 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmDdt.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmDdt.C
@@ -52,7 +52,7 @@ ddt
     return fv::ddtScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme("ddt(" + vf.name() + ')')
+        vf.mesh().schemes().ddt("ddt(" + vf.name() + ')')
     ).ref().fvmDdt(vf);
 }
 
@@ -80,7 +80,7 @@ ddt
     return fv::ddtScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme("ddt(" + rho.name() + ',' + vf.name() + ')')
+        vf.mesh().schemes().ddt("ddt(" + rho.name() + ',' + vf.name() + ')')
     ).ref().fvmDdt(rho, vf);
 }
 
@@ -96,7 +96,7 @@ ddt
     return fv::ddtScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme("ddt(" + rho.name() + ',' + vf.name() + ')')
+        vf.mesh().schemes().ddt("ddt(" + rho.name() + ',' + vf.name() + ')')
     ).ref().fvmDdt(rho, vf);
 }
 
@@ -113,7 +113,7 @@ ddt
     return fv::ddtScheme<Type>::New
     (
         vf.mesh(),
-        vf.mesh().ddtScheme
+        vf.mesh().schemes().ddt
         (
             "ddt("
           + alpha.name() + ','
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDiv.C b/src/finiteVolume/finiteVolume/fvm/fvmDiv.C
index f1e6b234b7117512e29266f2336f7bf2658bb1cd..88d1a5c1939d3c86c96ae3da3bd40e252d914555 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmDiv.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmDiv.C
@@ -55,7 +55,7 @@ div
     (
         vf.mesh(),
         flux,
-        vf.mesh().divScheme(name)
+        vf.mesh().schemes().div(name)
     )().fvmDiv(flux, vf);
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C b/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C
index d1729c417d5ac6b917fa753264d23ce80b13f044..459c06705664e117fe5581246a6f0d7a9c48b45a 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C
@@ -217,7 +217,7 @@ laplacian
     return fv::laplacianScheme<Type, GType>::New
     (
         vf.mesh(),
-        vf.mesh().laplacianScheme(name)
+        vf.mesh().schemes().laplacian(name)
     ).ref().fvmLaplacian(gamma, vf);
 }
 
@@ -282,7 +282,7 @@ laplacian
     return fv::laplacianScheme<Type, GType>::New
     (
         vf.mesh(),
-        vf.mesh().laplacianScheme(name)
+        vf.mesh().schemes().laplacian(name)
     ).ref().fvmLaplacian(gamma, vf);
 }
 
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C
index 37bd371ff88b9c1b6f5f1570f741cbbb3565ccbb..1de3bfb3181f16935a8de20304809b2c13235d1a 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C
@@ -96,7 +96,7 @@ Foam::fv::gradScheme<Type>::grad
     GradFieldType* pgGrad =
         mesh().objectRegistry::template getObjectPtr<GradFieldType>(name);
 
-    if (!this->mesh().cache(name) || this->mesh().changing())
+    if (!this->mesh().solution().cache(name) || this->mesh().changing())
     {
         // Delete any old occurrences to avoid double registration
         if (pgGrad && pgGrad->ownedByRegistry())
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/iterativeGaussGrad/iterativeGaussGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/iterativeGaussGrad/iterativeGaussGrad.C
index b46c064352303bf7327b4428369ffce77ee62fc6..0c59fef7a30ccc86f96a2bc64303d0b82f8022e9 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/iterativeGaussGrad/iterativeGaussGrad.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/iterativeGaussGrad/iterativeGaussGrad.C
@@ -62,7 +62,7 @@ Foam::fv::iterativeGaussGrad<Type>::calcGrad
 
     scalar relax = 1;
     const bool useRelax =
-        vsf.mesh().relaxField("grad(" + vsf.name() + ")", relax);
+        vsf.mesh().solution().relaxField("grad(" + vsf.name() + ")", relax);
 
     for (label i = 0; i < nIter_; ++i)
     {
diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C b/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C
index 20f20e80e1b4df8b11b79b96804819de65f402b7..d995a569aaa5a3fbd02f7ebec14d554e9b2384c6 100644
--- a/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C
+++ b/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C
@@ -192,7 +192,7 @@ gaussLaplacianScheme<Type, GType>::fvmLaplacian
 
     fvm.source() -= mesh.V()*fvc::div(tfaceFluxCorrection())().primitiveField();
 
-    if (mesh.fluxRequired(vf.name()))
+    if (mesh.schemes().fluxRequired(vf.name()))
     {
         fvm.faceFluxCorrectionPtr() = tfaceFluxCorrection.ptr();
     }
diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianSchemes.C b/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianSchemes.C
index 1201bad1fe2f66ae24587653e2da1b8d49d94a6b..5b957b482931f52d60a4d974fbd2034a0f0ce9fe 100644
--- a/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianSchemes.C
+++ b/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianSchemes.C
@@ -59,7 +59,7 @@ Foam::fv::gaussLaplacianScheme<Foam::Type, Foam::scalar>::fvmLaplacian         \
                                                                                \
     if (this->tsnGradScheme_().corrected())                                    \
     {                                                                          \
-        if (mesh.fluxRequired(vf.name()))                                      \
+        if (mesh.schemes().fluxRequired(vf.name()))                            \
         {                                                                      \
             fvm.faceFluxCorrectionPtr() = new                                  \
             GeometricField<Type, fvsPatchField, surfaceMesh>                   \
diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.C b/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.C
index d0732bf12414e737e983563f026fb960ac864b67..d38c8f873516ce2fc35476b9d6a7f5ca013a5128 100644
--- a/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.C
+++ b/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.C
@@ -195,7 +195,10 @@ relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvmLaplacian
     tmp<SType> trelaxedCorrection(new SType(tfaceFluxCorrection()));
 
     const word oldName(corrName + "_0");
-    const scalar relax(vf.mesh().equationRelaxationFactor(oldName));
+    const scalar relax
+    (
+        vf.mesh().solution().equationRelaxationFactor(oldName)
+    );
 
     const objectRegistry& obr = vf.db();
     if (obr.foundObject<SType>(oldName))
@@ -220,7 +223,7 @@ relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvmLaplacian
             trelaxedCorrection()
         )().primitiveField();
 
-    if (mesh.fluxRequired(vf.name()))
+    if (mesh.schemes().fluxRequired(vf.name()))
     {
         fvm.faceFluxCorrectionPtr() = trelaxedCorrection.ptr();
     }
diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianSchemes.C b/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianSchemes.C
index 7aa0027ec98c39be9dd6a7e249fda1636e5644ec..a119a75fb1f10d0983037dc3f245078e1949ca20 100644
--- a/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianSchemes.C
+++ b/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianSchemes.C
@@ -69,7 +69,10 @@ fvmLaplacian                                                                   \
         tmp<SType> trelaxedCorrection(new SType(tfaceFluxCorrection()));       \
                                                                                \
         const word oldName(corrName + "_0");                                   \
-        const scalar relax(vf.mesh().equationRelaxationFactor(corrName));      \
+        const scalar relax                                                     \
+        (                                                                      \
+            vf.mesh().solution().equationRelaxationFactor(corrName)            \
+        );                                                                     \
         const objectRegistry& obr = vf.db();                                   \
         if (obr.foundObject<SType>(oldName))                                   \
         {                                                                      \
@@ -96,7 +99,7 @@ fvmLaplacian                                                                   \
                                                                                \
         fvm.source() -= tcorr();                                               \
                                                                                \
-        if (mesh.fluxRequired(vf.name()))                                      \
+        if (mesh.schemes().fluxRequired(vf.name()))                            \
         {                                                                      \
             fvm.faceFluxCorrectionPtr() = trelaxedCorrection.ptr();            \
         }                                                                      \
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C b/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C
index ab5d7f406334c07ff0be26eed131193a5a1e1400..8e6f6897d5996fd3a536c1e13fedd8e5278b3150 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C
@@ -52,7 +52,7 @@ Foam::fv::correctedSnGrad<Type>::fullGradCorrection
             gradScheme<Type>::New
             (
                 mesh,
-                mesh.gradScheme("grad(" + vf.name() + ')')
+                mesh.schemes().grad("grad(" + vf.name() + ')')
             )().grad(vf, "grad(" + vf.name() + ')')
         );
     tssf.ref().rename("snGradCorr(" + vf.name() + ')');
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/relaxedSnGrad/relaxedSnGrad.C b/src/finiteVolume/finiteVolume/snGradSchemes/relaxedSnGrad/relaxedSnGrad.C
index 0b018126770a811d634eab17d209b5e9491627e5..1f77c180298c9e79204b6cd998ea82a6214ddb23 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/relaxedSnGrad/relaxedSnGrad.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/relaxedSnGrad/relaxedSnGrad.C
@@ -48,7 +48,7 @@ Foam::fv::relaxedSnGrad<Type>::correction
     const word fieldName(vf.name());
     const word oldFieldName(fieldName + "_0");
     const scalar relax =
-        vf.mesh().fieldRelaxationFactor("snGrad("+fieldName+")");
+        vf.mesh().solution().fieldRelaxationFactor("snGrad("+fieldName+")");
 
     // Return explicit correction field if
     // previous-time step correction is unavailable
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/skewCorrectedSnGrad/skewCorrectedSnGrad.C b/src/finiteVolume/finiteVolume/snGradSchemes/skewCorrectedSnGrad/skewCorrectedSnGrad.C
index 6d1c26c5143b0714c2399b25829fbf06d928ef9e..06275cc85223aef1a180094e39a6e0124fdebfa9 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/skewCorrectedSnGrad/skewCorrectedSnGrad.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/skewCorrectedSnGrad/skewCorrectedSnGrad.C
@@ -155,7 +155,7 @@ Foam::fv::skewCorrectedSnGrad<Type>::fullGradCorrection
             gradScheme<typename pTraits<Type>::cmptType>::New
             (
                 mesh,
-                mesh.gradScheme("grad(" + vf.name() + ')')
+                mesh.schemes().grad("grad(" + vf.name() + ')')
             )()
            .grad(vf.component(cmpt))
         );
@@ -208,7 +208,7 @@ Foam::fv::skewCorrectedSnGrad<Type>::fullGradCorrection
     //         gradScheme<Type>::New
     //         (
     //             mesh,
-    //             mesh.gradScheme("grad(" + vf.name() + ')')
+    //             mesh.schemes().grad("grad(" + vf.name() + ')')
     //         )().grad(vf, "grad(" + vf.name() + ')')
     //     );
     //
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 6b15f83645ad712ad3820fe17e074193b5a0d217..ca8bac1517643d8b429dea09c4a0a0ca541fddc1 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -1262,7 +1262,7 @@ void Foam::fvMatrix<Type>::relax()
 
     scalar relaxCoeff = 0;
 
-    if (psi_.mesh().relaxEquation(name, relaxCoeff))
+    if (psi_.mesh().solution().relaxEquation(name, relaxCoeff))
     {
         relax(relaxCoeff);
     }
@@ -1428,7 +1428,7 @@ Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>>
 Foam::fvMatrix<Type>::
 flux() const
 {
-    if (!psi_.mesh().fluxRequired(psi_.name()))
+    if (!psi_.mesh().schemes().fluxRequired(psi_.name()))
     {
         FatalErrorInFunction
             << "flux requested but " << psi_.name()
@@ -1534,14 +1534,14 @@ const Foam::dictionary& Foam::fvMatrix<Type>::solverDict
     const word& name
 ) const
 {
-    return psi_.mesh().solverDict(name);
+    return psi_.mesh().solution().solverDict(name);
 }
 
 
 template<class Type>
 const Foam::dictionary& Foam::fvMatrix<Type>::solverDict() const
 {
-    return psi_.mesh().solverDict
+    return psi_.mesh().solution().solverDict
     (
         psi_.select(psi_.mesh().data().isFinalIteration())
     );
diff --git a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
index a8d87158f8f5ebe0bbf994dfa2b3008bc06f4221..b5305f4aaca54276f028968b77713fdaf21ac60f 100644
--- a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
@@ -190,7 +190,7 @@ Foam::solverPerformance Foam::fvMatrix<Foam::scalar>::solveSegregated
     {
         createOrUpdateLduPrimitiveAssembly();
 
-        if (psi_.mesh().fluxRequired(psi_.name()))
+        if (psi_.mesh().schemes().fluxRequired(psi_.name()))
         {
             // Save lower/upper for flux calculation
             if (asymmetric())
@@ -292,7 +292,7 @@ Foam::solverPerformance Foam::fvMatrix<Foam::scalar>::solveSegregated
 
     if (useImplicit_)
     {
-        if (psi_.mesh().fluxRequired(psi_.name()))
+        if (psi_.mesh().schemes().fluxRequired(psi_.name()))
         {
             // Restore lower/upper
             if (asymmetric())
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
index 070e9b557ba497ca9fde91060b7043ecbde03c76..a51de37e544f67fdec601b152320d0f9cd437228 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
@@ -219,7 +219,7 @@ void Foam::MULES::limiterCorr
 
     const fvMesh& mesh = psi.mesh();
 
-    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
+    const dictionary& MULEScontrols = mesh.solution().solverDict(psi.name());
 
     const label nLimiterIter
     (
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index 30ad56a9d38e2b31094aff4c2f7e24050dd5982c..16123894b154fafa5e5918f6f4a7482fddf9b398 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -210,7 +210,7 @@ void Foam::MULES::limiter
 
     const fvMesh& mesh = psi.mesh();
 
-    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
+    const dictionary& MULEScontrols = mesh.solution().solverDict(psi.name());
 
     const label nLimiterIter
     (
diff --git a/src/finiteVolume/fvMesh/fvMeshGeometry.C b/src/finiteVolume/fvMesh/fvMeshGeometry.C
index e7afa733b174673484c3c668121b95b92214dcad..bb8d4255765e4eb5379db13695d16ffe419a69e4 100644
--- a/src/finiteVolume/fvMesh/fvMeshGeometry.C
+++ b/src/finiteVolume/fvMesh/fvMeshGeometry.C
@@ -261,7 +261,7 @@ const Foam::volScalarField::Internal& Foam::fvMesh::V00() const
 
 Foam::tmp<Foam::volScalarField::Internal> Foam::fvMesh::Vsc() const
 {
-    if (!steady() && moving() && time().subCycling())
+    if (!schemes().steady() && moving() && time().subCycling())
     {
         const TimeState& ts = time();
         const TimeState& ts0 = time().prevTimeState();
@@ -283,7 +283,7 @@ Foam::tmp<Foam::volScalarField::Internal> Foam::fvMesh::Vsc() const
 
 Foam::tmp<Foam::volScalarField::Internal> Foam::fvMesh::Vsc0() const
 {
-    if (!steady() && moving() && time().subCycling())
+    if (!schemes().steady() && moving() && time().subCycling())
     {
         const TimeState& ts = time();
         const TimeState& ts0 = time().prevTimeState();
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C
index aed48475eb88e3d4b0def1c0659c56a8e0ef5ba9..1fd75194970cc3caa09d10d05e2cd4f550806134 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C
@@ -149,7 +149,7 @@ Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter
 
     const word limiterFieldName(type() + "Limiter(" + phi.name() + ')');
 
-    if (this->mesh().cache("limiter"))
+    if (this->mesh().solution().cache("limiter"))
     {
         auto* fldptr = mesh.getObjectPtr<surfaceScalarField>(limiterFieldName);
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.C
index 7848088c17306fb0a51cb7ad07a322ee2f2d5a4d..377996f4a46b02cc5e50b593240d3cc0f8466c4a 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.C
@@ -72,7 +72,7 @@ Foam::linearUpwind<Type>::correction
         fv::gradScheme<scalar>::New
         (
             mesh,
-            mesh.gradScheme(gradSchemeName_)
+            mesh.schemes().grad(gradSchemeName_)
         )
     );
 
@@ -183,7 +183,7 @@ Foam::linearUpwind<Foam::vector>::correction
         fv::gradScheme<vector>::New
         (
             mesh,
-            mesh.gradScheme(gradSchemeName_)
+            mesh.schemes().grad(gradSchemeName_)
         )
     );
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwindV.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwindV.H
index 793a4ab115930329173c51ee01509afbf0042a63..14111b3f5e754217bc3e7769e494d82cfbfd6d4e 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwindV.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwindV.H
@@ -109,7 +109,7 @@ public:
                 fv::gradScheme<Type>::New
                 (
                     mesh,
-                    mesh.gradScheme(gradSchemeName_)
+                    mesh.schemes().grad(gradSchemeName_)
                 )
             )
         {}
@@ -129,7 +129,7 @@ public:
                 fv::gradScheme<Type>::New
                 (
                     mesh,
-                    mesh.gradScheme(gradSchemeName_)
+                    mesh.schemes().grad(gradSchemeName_)
                 )
             )
         {}
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C
index edfdf8a9921eeb3b1e72908a475d0e13f1f4381f..72d24defcb58067897a092c4953c1b10bc735cc2 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C
@@ -57,7 +57,7 @@ Foam::tmp<Foam::surfaceInterpolationScheme<Type>> Foam::fvc::scheme
     (
         faceFlux.mesh(),
         faceFlux,
-        faceFlux.mesh().interpolationScheme(name)
+        faceFlux.mesh().schemes().interpolation(name)
     );
 }
 
@@ -87,7 +87,7 @@ Foam::tmp<Foam::surfaceInterpolationScheme<Type>> Foam::fvc::scheme
     return surfaceInterpolationScheme<Type>::New
     (
         mesh,
-        mesh.interpolationScheme(name)
+        mesh.schemes().interpolation(name)
     );
 }
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
index acd90e7ff5a0e2fbb88cea4f3e9a299a5ce219ec..509b9934690435197f4c3177bab0b72a39d1d9f2 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
@@ -87,7 +87,7 @@ const Foam::fvGeometryScheme& Foam::surfaceInterpolation::geometry() const
         geometryPtr_ = fvGeometryScheme::New
         (
             mesh_,
-            mesh_.schemesDict().subOrEmptyDict("geometry"),
+            mesh_.schemes().dict().subOrEmptyDict("geometry"),
             basicFvGeometryScheme::typeName
         );
     }
diff --git a/src/functionObjects/field/age/age.C b/src/functionObjects/field/age/age.C
index 72cf677d2ff68696d19291f4de140b92c0de578b..64264095d9b1e95d0a62499ca55827ff91c8b10f 100644
--- a/src/functionObjects/field/age/age.C
+++ b/src/functionObjects/field/age/age.C
@@ -173,7 +173,7 @@ bool Foam::functionObjects::age::execute()
 
     // Set under-relaxation coeff
     scalar relaxCoeff = 0;
-    mesh_.relaxEquation(schemesField_, relaxCoeff);
+    mesh_.solution().relaxEquation(schemesField_, relaxCoeff);
 
     Foam::fv::options& fvOptions(Foam::fv::options::New(mesh_));
 
diff --git a/src/functionObjects/field/blendingFactor/blendingFactorTemplates.C b/src/functionObjects/field/blendingFactor/blendingFactorTemplates.C
index e238f4677ec54fda535a97add708f0b254b3c5fd..41afeb2e2597f7e76ea6b985d42f571162a6b735 100644
--- a/src/functionObjects/field/blendingFactor/blendingFactorTemplates.C
+++ b/src/functionObjects/field/blendingFactor/blendingFactorTemplates.C
@@ -102,7 +102,7 @@ bool Foam::functionObjects::blendingFactor::calcScheme()
     const FieldType& field = lookupObject<FieldType>(fieldName_);
 
     const word divScheme("div(" + phiName_ + ',' + fieldName_ + ')');
-    ITstream& its = mesh_.divScheme(divScheme);
+    ITstream& its = mesh_.schemes().div(divScheme);
 
     const surfaceScalarField& phi = lookupObject<surfaceScalarField>(phiName_);
 
diff --git a/src/functionObjects/field/ddt2/ddt2.C b/src/functionObjects/field/ddt2/ddt2.C
index 0cdc3252977c2af98cfe6cd6f7e16a06cb10be38..fd58ac3841c13ea3ab3259cfcfc43ad1486f6f48 100644
--- a/src/functionObjects/field/ddt2/ddt2.C
+++ b/src/functionObjects/field/ddt2/ddt2.C
@@ -129,7 +129,8 @@ bool Foam::functionObjects::ddt2::read(const dictionary& dict)
 {
     fvMeshFunctionObject::read(dict);
 
-    if (word(mesh_.ddtScheme("default")) == "steadyState")
+    word ddtScheme(mesh_.schemes().ddt("default"));
+    if (ddtScheme == "steady" || ddtScheme == "steadyState")
     {
         WarningInFunction
             << typeName << " function object not appropriate for steady-state"
diff --git a/src/functionObjects/field/grad/gradTemplates.C b/src/functionObjects/field/grad/gradTemplates.C
index 46229086353d12e2cf62017f1c491507bdce5a77..9505a55dce677234a1e53253e3b866fcfcc3d9c2 100644
--- a/src/functionObjects/field/grad/gradTemplates.C
+++ b/src/functionObjects/field/grad/gradTemplates.C
@@ -41,7 +41,7 @@ bool Foam::functionObjects::grad::calcGrad()
         (
             resultName_,
             fvc::grad(lookupObject<VolFieldType>(fieldName_)),
-            mesh_.changing() && mesh_.cache(resultName_)
+            mesh_.changing() && mesh_.solution().cache(resultName_)
         );
     }
     else if (foundObject<SurfaceFieldType>(fieldName_, false))
@@ -50,7 +50,7 @@ bool Foam::functionObjects::grad::calcGrad()
         (
             resultName_,
             fvc::grad(lookupObject<SurfaceFieldType>(fieldName_)),
-            mesh_.changing() && mesh_.cache(resultName_)
+            mesh_.changing() && mesh_.solution().cache(resultName_)
         );
     }
 
diff --git a/src/functionObjects/solvers/energyTransport/energyTransport.C b/src/functionObjects/solvers/energyTransport/energyTransport.C
index ea5b022f117853cb267d42d18c7ca1c25b41a7cf..95ee13625958b3457c51d329b5fad26d77f2e8c7 100644
--- a/src/functionObjects/solvers/energyTransport/energyTransport.C
+++ b/src/functionObjects/solvers/energyTransport/energyTransport.C
@@ -372,7 +372,7 @@ bool Foam::functionObjects::energyTransport::execute()
 
     // Set under-relaxation coeff
     scalar relaxCoeff = 0;
-    mesh_.relaxEquation(schemesField_, relaxCoeff);
+    mesh_.solution().relaxEquation(schemesField_, relaxCoeff);
 
     if (phi.dimensions() == dimMass/dimTime)
     {
diff --git a/src/functionObjects/solvers/scalarTransport/scalarTransport.C b/src/functionObjects/solvers/scalarTransport/scalarTransport.C
index ec36978276013c0150671fadddad2263deb44362..1d5e88cf93c70533f3eecdfe99bf08b4411e9051 100644
--- a/src/functionObjects/solvers/scalarTransport/scalarTransport.C
+++ b/src/functionObjects/solvers/scalarTransport/scalarTransport.C
@@ -78,7 +78,7 @@ Foam::volScalarField& Foam::functionObjects::scalarTransport::transportedField()
 
         if (phaseName_ != "none")
         {
-            mesh_.setFluxRequired(fieldName_);
+            mesh_.schemes().setFluxRequired(fieldName_);
         }
     }
 
@@ -264,7 +264,7 @@ bool Foam::functionObjects::scalarTransport::execute()
 
     // Set under-relaxation coeff
     scalar relaxCoeff = 0;
-    mesh_.relaxEquation(schemesField_, relaxCoeff);
+    mesh_.solution().relaxEquation(schemesField_, relaxCoeff);
 
     // Two phase scalar transport
     if (phaseName_ != "none")
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
index e3298fe8a57d2e8a59d7d928862f4466094c46d9..9225893add80b1b8a6cfd955c17b76311df9ea12 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
@@ -142,7 +142,7 @@ void Foam::cloudSolution::read()
             transient_ = false;
         }
 
-        if (mesh_.steady())
+        if (mesh_.schemes().steady())
         {
             IOWarningInFunction(dict_)
                 << "Transient tracking is not supported for steady-state"
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
index c811306bc4cab56107b6282ff97a2d96b2b2ef0d..65ab932cd27b0e14d78ee66181a68b8c2fcda23e 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
@@ -128,7 +128,7 @@ void Foam::PackingModels::Implicit<CloudType>::cacheFields(const bool store)
                 cloudName + ":uSqrAverage"
             );
 
-        mesh.setFluxRequired(alpha_.name());
+        mesh.schemes().setFluxRequired(alpha_.name());
 
         // Property fields
         // ~~~~~~~~~~~~~~~
diff --git a/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointBoundaryCondition/adjointBoundaryCondition.C b/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointBoundaryCondition/adjointBoundaryCondition.C
index 9523408d2d2855d854e0efbc8ad660de76875787..0805287fe98147bf9f71a34c7fcf6d8f00f5f73c 100644
--- a/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointBoundaryCondition/adjointBoundaryCondition.C
+++ b/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointBoundaryCondition/adjointBoundaryCondition.C
@@ -69,7 +69,7 @@ adjointBoundaryCondition<Type>::computePatchGrad(word name)
     // Interpolation scheme is now read through interpolation schemes.
     /*
     word gradSchemeName("grad(" + name + ')');
-    ITstream& is = mesh.gradScheme(gradSchemeName);
+    ITstream& is = mesh.schemes().grad(gradSchemeName);
     word schemeData(is);
     */
 
@@ -89,7 +89,8 @@ adjointBoundaryCondition<Type>::computePatchGrad(word name)
         surfFieldPtr =
             surfaceInterpolationScheme<Type2>::New
             (
-                mesh, mesh.interpolationScheme("interpolate(" + name + ")")
+                mesh,
+                mesh.schemes().interpolation("interpolate(" + name + ")")
             )().interpolate(field).ptr();
         surfFieldPtr->rename(surfFieldName);
         regIOobject::store(surfFieldPtr);
@@ -108,7 +109,8 @@ adjointBoundaryCondition<Type>::computePatchGrad(word name)
             surfFieldPtr =
                 surfaceInterpolationScheme<Type2>::New
                 (
-                    mesh, mesh.interpolationScheme("interpolate(" + name + ")")
+                    mesh,
+                    mesh.schemes().interpolation("interpolate(" + name + ")")
                 )().interpolate(field).ptr();
             surfFieldPtr->rename(surfFieldName);
             regIOobject::store(surfFieldPtr);
diff --git a/src/optimisation/adjointOptimisation/adjoint/displacementMethod/displacementMethod/displacementMethod.C b/src/optimisation/adjointOptimisation/adjoint/displacementMethod/displacementMethod/displacementMethod.C
index d0da1a72118b629d2fb6567d1d66e78bf8a5caaa..bdb9cf95d0bf72f187266bd39a084c32c378504c 100644
--- a/src/optimisation/adjointOptimisation/adjoint/displacementMethod/displacementMethod/displacementMethod.C
+++ b/src/optimisation/adjointOptimisation/adjoint/displacementMethod/displacementMethod/displacementMethod.C
@@ -128,7 +128,7 @@ void Foam::displacementMethod::update()
     mesh_.movePoints(tnewPoints());
     scalar timeAft = mesh_.time().elapsedCpuTime();
     Info<< "Mesh movement took " << timeAft - timeBef << " seconds" << endl;
-    if (!mesh_.steady())
+    if (!mesh_.schemes().steady())
     {
         mesh_.moving(false);
     }
diff --git a/src/optimisation/adjointOptimisation/adjoint/finiteVolume/interpolation/surfaceInterpolation/schemes/limitedSchemes/linearUpwindNormal/linearUpwindNormal.H b/src/optimisation/adjointOptimisation/adjoint/finiteVolume/interpolation/surfaceInterpolation/schemes/limitedSchemes/linearUpwindNormal/linearUpwindNormal.H
index 93f48dd598ef16f4ea67f44c558c580e8048f032..1b11e364d9c9b413083dac14db68d9f5e51258d6 100644
--- a/src/optimisation/adjointOptimisation/adjoint/finiteVolume/interpolation/surfaceInterpolation/schemes/limitedSchemes/linearUpwindNormal/linearUpwindNormal.H
+++ b/src/optimisation/adjointOptimisation/adjoint/finiteVolume/interpolation/surfaceInterpolation/schemes/limitedSchemes/linearUpwindNormal/linearUpwindNormal.H
@@ -113,7 +113,7 @@ public:
                 fv::gradScheme<Type>::New
                 (
                     mesh,
-                    mesh.gradScheme(gradSchemeName_)
+                    mesh.schemes().grad(gradSchemeName_)
                 )
             )
         {}
@@ -133,7 +133,7 @@ public:
                 fv::gradScheme<Type>::New
                 (
                     mesh,
-                    mesh.gradScheme(gradSchemeName_)
+                    mesh.schemes().grad(gradSchemeName_)
                 )
             )
         {}
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/adjointEikonalSolver/adjointEikonalSolverIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/adjointEikonalSolver/adjointEikonalSolverIncompressible.C
index bdffd2682eb2bacec758f5189c1e5a8d97b2b818..6c5d87671bc2732bc9879a8410191b421986dcde 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/adjointEikonalSolver/adjointEikonalSolverIncompressible.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/adjointEikonalSolver/adjointEikonalSolverIncompressible.C
@@ -68,7 +68,7 @@ void adjointEikonalSolver::read()
     tolerance_ = dict_.getOrDefault<scalar>("tolerance", 1e-6);
     epsilon_ = dict_.getOrDefault<scalar>("epsilon", 0.1);
     const scalar defaultEps =
-        mesh_.schemesDict().subDict("wallDist").
+        mesh_.schemes().dict().subDict("wallDist").
             subOrEmptyDict("advectionDiffusionCoeffs").
                 getOrDefault<scalar>("epsilon", 0.1);
     epsilon_ = dict_.getOrDefault<scalar>("epsilon", defaultEps);
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/adjointMeshMovementSolver/adjointMeshMovementSolverIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/adjointMeshMovementSolver/adjointMeshMovementSolverIncompressible.C
index 663db50998daa3e5bbeac91d1872e577541a49a0..4614aa99eda1fd7f77f23de2f8c7bde3d59cabf8 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/adjointMeshMovementSolver/adjointMeshMovementSolverIncompressible.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/adjointMeshMovementSolver/adjointMeshMovementSolverIncompressible.C
@@ -155,7 +155,11 @@ void adjointMeshMovementSolver::solve()
 
         //scalar residual = max(maEqn.solve().initialResidual());
         scalar residual =
-            mag(Foam::solve(maEqn, mesh_.solverDict("ma")).initialResidual());
+            mag
+            (
+                Foam::solve(maEqn, mesh_.solution().solverDict("ma"))
+                .initialResidual()
+            );
 
         Info<< "Max ma " << gMax(mag(ma_)()) << endl;
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVars.C b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVars.C
index 0fe83537d5ab48643484bb68763264b38923fe38..55f88deea40aba336d637fb37a2841d90d520736 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVars.C
+++ b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVars.C
@@ -55,7 +55,7 @@ void incompressibleVars::setFields()
         useSolverNameForFields_
     );
 
-    mesh_.setFluxRequired(pPtr_->name());
+    mesh_.schemes().setFluxRequired(pPtr_->name());
 
     // if required, correct boundary conditions of mean flow fields here in
     // order to have the correct bcs for e.g. turbulence models that follow.
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVars.C b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVars.C
index 7a23aba9254ffa9f3bbc863822bea61e6b39398c..b9e27d2efca7a4e5903a680dfaeda34be27f6673 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVars.C
+++ b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVars.C
@@ -54,7 +54,7 @@ void incompressibleAdjointMeanFlowVars::setFields()
         useSolverNameForFields_
     );
 
-    mesh_.setFluxRequired(paPtr_->name());
+    mesh_.schemes().setFluxRequired(paPtr_->name());
 }
 
 void incompressibleAdjointMeanFlowVars::setMeanFields()
diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSST.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSST.C
index 862cea770d408e2fd214543925878cbd4daef2c3..7c080c31fc346b69d3bde601b1213469d29ccb95 100644
--- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSST.C
+++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSST.C
@@ -1187,7 +1187,7 @@ tmp<surfaceInterpolationScheme<scalar>> adjointkOmegaSST::convectionScheme
     const surfaceScalarField& phi = primalVars_.phi();
     const surfaceScalarField& phiInst = primalVars_.phiInst();
     word divEntry("div(" + phiInst.name() + ',' + varName +')');
-    ITstream& divScheme = mesh_.divScheme(divEntry);
+    ITstream& divScheme = mesh_.schemes().div(divEntry);
     // Skip the first entry which might be 'bounded' or 'Gauss'.
     // If it is 'bounded', skip the second entry as well
     word discarded(divScheme);
diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSSTTemplates.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSSTTemplates.C
index 467cd8b01f64da3f6aca0035bc403888d0e5a3c3..343ec646cac86f0e2744f03836c2d210b12262ff 100644
--- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSSTTemplates.C
+++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSSTTemplates.C
@@ -48,11 +48,11 @@ tmp<surfaceInterpolationScheme<Type>> adjointkOmegaSST::interpolationScheme
     return
         tmp<surfaceInterpolationScheme<Type>>
         (
-            mesh_.interpolationSchemes().found(schemeName) ?
+            mesh_.schemes().interpolationSchemes().found(schemeName) ?
             surfaceInterpolationScheme<Type>::New
             (
                 mesh_,
-                mesh_.interpolationScheme(schemeName)
+                mesh_.schemes().interpolation(schemeName)
             ) :
             tmp<surfaceInterpolationScheme<Type>>
                 (new reverseLinear<Type>(mesh_))
diff --git a/src/overset/cellCellStencil/cellCellStencil/cellCellStencilObject.H b/src/overset/cellCellStencil/cellCellStencil/cellCellStencilObject.H
index 4ecd06006f32abbd229512b4c3ab1004ca98eb55..de394707106810df7958517e36700748280cb9d3 100644
--- a/src/overset/cellCellStencil/cellCellStencil/cellCellStencilObject.H
+++ b/src/overset/cellCellStencil/cellCellStencil/cellCellStencilObject.H
@@ -90,10 +90,7 @@ public:
                 cellCellStencil::New
                 (
                     mesh,
-                    mesh.schemesDict().subDict
-                    (
-                        "oversetInterpolation"
-                    ),
+                    mesh.schemes().dict().subDict("oversetInterpolation"),
                     update
                 )
             )
diff --git a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C
index b0e999c0f88855ab14d359dfafa3d54b7533974a..4b19d46976a810a4dc892c6922019129aadfffca 100644
--- a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C
+++ b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C
@@ -2183,14 +2183,14 @@ bool Foam::cellCellStencils::inverseDistance::update()
 
 
     label useLayer = dict_.getOrDefault("useLayer", -1);
-    const dictionary& fvSchemes = mesh_.schemesDict();
-    if (fvSchemes.found("oversetInterpolation"))
     {
-        const dictionary& intDict = fvSchemes.subDict
-        (
-            "oversetInterpolation"
-        );
-        useLayer = intDict.getOrDefault("useLayer", -1);
+        const dictionary* interpDict =
+            mesh_.schemes().dict().findDict("oversetInterpolation");
+
+        if (interpDict)
+        {
+            useLayer = interpDict->getOrDefault("useLayer", -1);
+        }
     }
 
     walkFront
diff --git a/src/overset/oversetFvMesh/oversetFvMeshBase.C b/src/overset/oversetFvMesh/oversetFvMeshBase.C
index b5a95780282a6e8520486190bc52da75866cf164..0ceb4376addd46bf69ed8a0b3650ad27e6b1a332 100644
--- a/src/overset/oversetFvMesh/oversetFvMeshBase.C
+++ b/src/overset/oversetFvMesh/oversetFvMeshBase.C
@@ -606,7 +606,7 @@ bool Foam::oversetFvMeshBase::interpolateFields()
     // Use whatever the solver has set up as suppression list
     const dictionary* dictPtr
     (
-        mesh_.schemesDict().findDict("oversetInterpolationSuppressed")
+        mesh_.schemes().dict().findDict("oversetInterpolationSuppressed")
     );
     if (dictPtr)
     {
diff --git a/src/overset/oversetPolyPatch/oversetFvPatchField.C b/src/overset/oversetPolyPatch/oversetFvPatchField.C
index 3beba4a72f3c05bf72e76fceebd70ee225bcd590..6ae659bbd9288bd04a781d544ec2ba2341363efd 100644
--- a/src/overset/oversetPolyPatch/oversetFvPatchField.C
+++ b/src/overset/oversetPolyPatch/oversetFvPatchField.C
@@ -602,7 +602,7 @@ void Foam::oversetFvPatchField<Type>::initEvaluate
     {
         // Trigger interpolation
         const fvMesh& mesh = this->internalField().mesh();
-        const dictionary& fvSchemes = mesh.schemesDict();
+        const dictionary& fvSchemes = mesh.schemes().dict();
         const word& fldName = this->internalField().name();
 
         if (&mesh.lduAddr() != &mesh.fvMesh::lduAddr())
diff --git a/src/phaseSystemModels/multiphaseEuler/multiphaseSystem/multiphaseSystem/multiphaseSystem.C b/src/phaseSystemModels/multiphaseEuler/multiphaseSystem/multiphaseSystem/multiphaseSystem.C
index c8142628ac59f39d0909458d0413686dc17c79b5..825e12751437b160de08dbd6184d0fedf94e3527 100644
--- a/src/phaseSystemModels/multiphaseEuler/multiphaseSystem/multiphaseSystem/multiphaseSystem.C
+++ b/src/phaseSystemModels/multiphaseEuler/multiphaseSystem/multiphaseSystem/multiphaseSystem.C
@@ -839,7 +839,7 @@ void Foam::multiphaseSystem::solve()
 
     const Time& runTime = mesh_.time();
 
-    const dictionary& alphaControls = mesh_.solverDict("alpha");
+    const dictionary& alphaControls = mesh_.solution().solverDict("alpha");
     label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
 
     if (nAlphaSubCycles > 1)
diff --git a/src/phaseSystemModels/multiphaseInter/phasesSystem/derivedFvPatchFields/timeVaryingMassSorption/timeVaryingMassSorptionFvPatchScalarField.C b/src/phaseSystemModels/multiphaseInter/phasesSystem/derivedFvPatchFields/timeVaryingMassSorption/timeVaryingMassSorptionFvPatchScalarField.C
index 771d7e0b77c7947d3a55a0a6bf76a91e1bd0716a..f7790b21411af87405ef9614c08a58ea0e117210 100644
--- a/src/phaseSystemModels/multiphaseInter/phasesSystem/derivedFvPatchFields/timeVaryingMassSorption/timeVaryingMassSorptionFvPatchScalarField.C
+++ b/src/phaseSystemModels/multiphaseInter/phasesSystem/derivedFvPatchFields/timeVaryingMassSorption/timeVaryingMassSorptionFvPatchScalarField.C
@@ -189,7 +189,7 @@ void Foam::timeVaryingMassSorptionFvPatchScalarField::updateCoeffs()
     const volScalarField& fld0 = fld.oldTime();
 
     // Lookup d/dt scheme from database
-    const word ddtSchemeName(fld.mesh().ddtScheme(fld.name()));
+    const word ddtSchemeName(fld.mesh().schemes().ddt(fld.name()));
     const ddtSchemeType& ddtScheme = ddtSchemeTypeNames_[ddtSchemeName];
 
     const scalarField cp(*this);
diff --git a/src/phaseSystemModels/multiphaseInter/phasesSystem/multiphaseSystem/multiphaseSystem.C b/src/phaseSystemModels/multiphaseInter/phasesSystem/multiphaseSystem/multiphaseSystem.C
index 4cf91a30e38faa9f6492a36d4716718eadfc99ee..0a6b93d2fb6015e2a27c93ee464684ae621844d3 100644
--- a/src/phaseSystemModels/multiphaseInter/phasesSystem/multiphaseSystem/multiphaseSystem.C
+++ b/src/phaseSystemModels/multiphaseInter/phasesSystem/multiphaseSystem/multiphaseSystem.C
@@ -77,7 +77,7 @@ Foam::multiphaseInter::multiphaseSystem::multiphaseSystem
         phases_.set(phasei++, &pm);
     }
 
-    mesh.solverDict("alpha").readEntry("cAlphas", cAlphas_);
+    mesh.solution().solverDict("alpha").readEntry("cAlphas", cAlphas_);
 
     // Initiate Su and Sp
     forAllConstIters(phaseModels_, iter)
@@ -129,7 +129,7 @@ void Foam::multiphaseInter::multiphaseSystem::calculateSuSp()
 
 void Foam::multiphaseInter::multiphaseSystem::solve()
 {
-    const dictionary& alphaControls = mesh_.solverDict("alpha");
+    const dictionary& alphaControls = mesh_.solution().solverDict("alpha");
     label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
 
     volScalarField& alpha = phases_.first();
@@ -172,7 +172,7 @@ void Foam::multiphaseInter::multiphaseSystem::solve()
 void Foam::multiphaseInter::multiphaseSystem::solveAlphas()
 {
 
-    const dictionary& alphaControls = mesh_.solverDict("alpha");
+    const dictionary& alphaControls = mesh_.solution().solverDict("alpha");
     alphaControls.readEntry("cAlphas", cAlphas_);
     label nAlphaCorr(alphaControls.get<label>("nAlphaCorr"));
 
diff --git a/src/phaseSystemModels/multiphaseInter/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C b/src/phaseSystemModels/multiphaseInter/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
index e6ab27f26804d5013cac6f39d6281afa4072679e..39745d48d5a2aa3c03cb682fe5e87e983a7d5c21 100644
--- a/src/phaseSystemModels/multiphaseInter/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
+++ b/src/phaseSystemModels/multiphaseInter/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
@@ -213,7 +213,7 @@ void Foam::MultiComponentPhaseModel<BasePhaseModel, phaseThermo>::solveYi
 
     const fvMesh& mesh = alpha1.mesh();
 
-    const dictionary& MULEScontrols = mesh.solverDict(alpha1.name());
+    const dictionary& MULEScontrols = mesh.solution().solverDict(alpha1.name());
 
     scalar cAlpha(MULEScontrols.get<scalar>("cYi"));
 
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
index 7623bd2d7029e49ae7ee475a15baa2f7a4fce043..60a4ea89d665f07c71131f662df713d6487c654d 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
@@ -575,12 +575,8 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiFs
         this->addField(phase, "phiF", pPrimeByAf*snGradAlpha1, phiFs);
 
         const bool implicitPhasePressure =
-            this->mesh_.solverDict(phase.volScalarField::name()).
-            template getOrDefault<Switch>
-            (
-                "implicitPhasePressure",
-                false
-            );
+            this->mesh_.solution().solverDict(phase.volScalarField::name())
+                .getOrDefault("implicitPhasePressure", false);
 
         if (implicitPhasePressure)
         {
@@ -738,12 +734,8 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiFfs
         this->addField(phase, "phiFf", pPrimeByAf*snGradAlpha1, phiFfs);
 
         const bool implicitPhasePressure =
-            this->mesh_.solverDict(phase.volScalarField::name()).
-            template getOrDefault<Switch>
-            (
-                "implicitPhasePressure",
-                false
-            );
+            this->mesh_.solution().solverDict(phase.volScalarField::name())
+                .getOrDefault("implicitPhasePressure", false);
 
         if (implicitPhasePressure)
         {
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
index e3a2be4b9ab3e2f4e0cfe420978722f177ddf985..91d7d0b3f24f82324461bef3edb54ca1e8c37af4 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -436,7 +436,10 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
             << ", max = " << gMax(Tf.primitiveField())
             << endl;
 
-        scalar iDmdtRelax(this->mesh().fieldRelaxationFactor("iDmdt"));
+        scalar iDmdtRelax
+        (
+            this->mesh().solution().fieldRelaxationFactor("iDmdt")
+        );
         iDmdt = (1 - iDmdtRelax)*iDmdt + iDmdtRelax*iDmdtNew;
 
         if (phaseChange_)
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/diameterModels/velocityGroup/velocityGroup.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/diameterModels/velocityGroup/velocityGroup.C
index 15e50d8b8f1b1883820450aa6d96021593206efe..b92147daf98a5a0aa218cdcb53802781093c513d 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/diameterModels/velocityGroup/velocityGroup.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/diameterModels/velocityGroup/velocityGroup.C
@@ -130,7 +130,7 @@ Foam::diameterModels::velocityGroup::mvconvection() const
             phase_.mesh(),
             fields_,
             phase_.alphaRhoPhi(),
-            phase_.mesh().divScheme
+            phase_.mesh().schemes().div
             (
                 "div(" + phase_.alphaRhoPhi()().name() + ",f)"
             )
@@ -222,17 +222,11 @@ Foam::diameterModels::velocityGroup::velocityGroup
 {
     if
     (
-        phase_.mesh().solverDict(popBalName_).getOrDefault<Switch>
-        (
-            "renormalizeAtRestart",
-            false
-        )
+        phase_.mesh().solution().solverDict(popBalName_)
+        .getOrDefault("renormalizeAtRestart", false)
      ||
-        phase_.mesh().solverDict(popBalName_).getOrDefault<Switch>
-        (
-            "renormalize",
-            false
-        )
+        phase_.mesh().solution().solverDict(popBalName_)
+        .getOrDefault("renormalize", false)
     )
     {
         renormalize();
@@ -307,11 +301,8 @@ void Foam::diameterModels::velocityGroup::postSolve()
 
     if
     (
-        phase_.mesh().solverDict(popBalName_).getOrDefault<Switch>
-        (
-            "renormalize",
-            false
-        )
+        phase_.mesh().solution().solverDict(popBalName_)
+        .getOrDefault("renormalize", false)
     )
     {
         renormalize();
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/multiphaseSystem/multiphaseSystem.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/multiphaseSystem/multiphaseSystem.C
index e232ec1e39738a4db475bcf35cbeabac45ebfccf..62238f58f6b3c93bfe7ca0ea28de2a448ee61c70 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/multiphaseSystem/multiphaseSystem.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/multiphaseSystem/multiphaseSystem.C
@@ -545,7 +545,7 @@ Foam::multiphaseSystem::multiphaseSystem
     forAll(phases(), phasei)
     {
         volScalarField& alphai = phases()[phasei];
-        mesh_.setFluxRequired(alphai.name());
+        mesh_.schemes().setFluxRequired(alphai.name());
     }
 }
 
@@ -627,7 +627,7 @@ void Foam::multiphaseSystem::solve()
 {
     const Time& runTime = mesh_.time();
 
-    const dictionary& alphaControls = mesh_.solverDict("alpha");
+    const dictionary& alphaControls = mesh_.solution().solverDict("alpha");
     label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
 
     bool LTS = fv::localEulerDdt::enabled(mesh_);
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
index 9482f0244cbf1bb86173005d4cf0e0f8c1c4ee6e..5afd8448946100073a71cf203de038ae26abb5fd 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
@@ -59,7 +59,7 @@ Foam::MultiComponentPhaseModel<BasePhaseModel>::MultiComponentPhaseModel
     (
         "residualAlpha",
         dimless,
-        fluid.mesh().solverDict("Yi")
+        fluid.mesh().solution().solverDict("Yi")
     ),
     inertIndex_(-1)
 {
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/populationBalanceModel/populationBalanceModel.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
index db4ccbfac4d3b9053de7f2e24421b72fe1c5ad2b..4e3fd97445e885f03877d3d9a70327aaace29bc2 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
@@ -1233,7 +1233,7 @@ Foam::diameterModels::populationBalanceModel::continuousTurbulence() const
 
 void Foam::diameterModels::populationBalanceModel::solve()
 {
-    const dictionary& solutionControls = mesh_.solverDict(name_);
+    const dictionary& solutionControls = mesh_.solution().solverDict(name_);
     const bool solveOnFinalIterOnly =
         solutionControls.getOrDefault("solveOnFinalIterOnly", false);
 
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H b/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H
index 95f85381de5231a8d4a554f9aa5e3e48e0eb5430..4ba2bc94aa928ce145b85ca99eb1d3f74a795896 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H
@@ -30,7 +30,7 @@ License
 
 inline Foam::label Foam::diameterModels::populationBalanceModel::nCorr() const
 {
-    return mesh_.solverDict(name_).get<label>("nCorr");
+    return mesh_.solution().solverDict(name_).get<label>("nCorr");
 }
 
 
@@ -38,7 +38,7 @@ inline Foam::label
 Foam::diameterModels::populationBalanceModel::sourceUpdateInterval() const
 {
     return
-        mesh_.solverDict(name_)
+        mesh_.solution().solverDict(name_)
        .getOrDefault<label>("sourceUpdateInterval", 1);
 }
 
diff --git a/src/phaseSystemModels/reactingEuler/twoPhaseSystem/twoPhaseSystem.C b/src/phaseSystemModels/reactingEuler/twoPhaseSystem/twoPhaseSystem.C
index 51f8dd9898ccc067a4d5e06f05dd6444ad727e12..87926186f47f32943f594873fa80f656fef58990 100644
--- a/src/phaseSystemModels/reactingEuler/twoPhaseSystem/twoPhaseSystem.C
+++ b/src/phaseSystemModels/reactingEuler/twoPhaseSystem/twoPhaseSystem.C
@@ -67,7 +67,7 @@ Foam::twoPhaseSystem::twoPhaseSystem
     phase2_.volScalarField::operator=(scalar(1) - phase1_);
 
     volScalarField& alpha1 = phase1_;
-    mesh.setFluxRequired(alpha1.name());
+    mesh.schemes().setFluxRequired(alpha1.name());
 }
 
 
@@ -126,7 +126,8 @@ void Foam::twoPhaseSystem::solve()
     volScalarField& alpha1 = phase1_;
     volScalarField& alpha2 = phase2_;
 
-    const dictionary& alphaControls = mesh_.solverDict(alpha1.name());
+    const dictionary& alphaControls =
+        mesh_.solution().solverDict(alpha1.name());
 
     label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
     label nAlphaCorr(alphaControls.get<label>("nAlphaCorr"));
diff --git a/src/phaseSystemModels/twoPhaseEuler/twoPhaseSystem/twoPhaseSystem/twoPhaseSystem.C b/src/phaseSystemModels/twoPhaseEuler/twoPhaseSystem/twoPhaseSystem/twoPhaseSystem.C
index fdb52fea29c308129913b8aca7ac4cddf14fb118..6168f2bc37d9e8f60c2b743d430ec761c9f047e6 100644
--- a/src/phaseSystemModels/twoPhaseEuler/twoPhaseSystem/twoPhaseSystem/twoPhaseSystem.C
+++ b/src/phaseSystemModels/twoPhaseEuler/twoPhaseSystem/twoPhaseSystem/twoPhaseSystem.C
@@ -359,7 +359,7 @@ void Foam::twoPhaseSystem::solve()
     const surfaceScalarField& phi1 = phase1_.phi();
     const surfaceScalarField& phi2 = phase2_.phi();
 
-    const dictionary& alphaControls = mesh_.solverDict
+    const dictionary& alphaControls = mesh_.solution().solverDict
     (
         alpha1.name()
     );
diff --git a/src/regionFaModels/regionFaModel/regionFaModelI.H b/src/regionFaModels/regionFaModel/regionFaModelI.H
index 7e8a407e7f9581ea5f39a02ef810c5a01af0f022..8fef5562d21144fd9e1c8596f97ceafa66877bd9 100644
--- a/src/regionFaModels/regionFaModel/regionFaModelI.H
+++ b/src/regionFaModels/regionFaModel/regionFaModelI.H
@@ -85,7 +85,7 @@ Foam::regionModels::regionFaModel::outputProperties()
 inline const Foam::dictionary&
 Foam::regionModels::regionFaModel::solution() const
 {
-    return regionMesh().solutionDict();
+    return regionMesh().solution().dict();
 }
 
 
diff --git a/src/regionModels/regionModel/regionModel/regionModelI.H b/src/regionModels/regionModel/regionModel/regionModelI.H
index 91590d87b60c748c776ad57d2f989b552a4a0663..e18dea7265a86799f5c17d85488aa04949cab136 100644
--- a/src/regionModels/regionModel/regionModel/regionModelI.H
+++ b/src/regionModels/regionModel/regionModel/regionModelI.H
@@ -61,7 +61,7 @@ inline Foam::fvMesh& Foam::regionModels::regionModel::regionMesh()
 inline const Foam::dictionary&
 Foam::regionModels::regionModel::solution() const
 {
-    return regionMesh().solutionDict();
+    return regionMesh().solution().dict();
 }
 
 
diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
index 58a6beeb1b5bf2e992aa77a53d99e25581e0e119..84ba8ffeb4483a5040f2da97ac9b2aca20b1533c 100644
--- a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
@@ -384,7 +384,7 @@ void kinematicSingleLayer::solveThickness
         fvc::interpolate(delta_)*deltarUAf*rhof*fvc::interpolate(pp)
     );
 
-    regionMesh().setFluxRequired(delta_.name());
+    regionMesh().schemes().setFluxRequired(delta_.name());
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr_; nonOrth++)
     {
diff --git a/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C b/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C
index 9cf0a65f6ba2bfaa81e393224db32facad4fe0de..e92c6ba892761e40be0485b400e6aff85e4b5bf0 100644
--- a/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C
+++ b/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C
@@ -172,7 +172,9 @@ void Foam::solarCalculator::initialise()
         }
         case mSunDirTracking:
         {
-            if (word(mesh_.ddtScheme("default")) == "steadyState")
+            const word ddtScheme(mesh_.schemes().ddt("default"));
+
+            if (ddtScheme == "steady" || ddtScheme == "steadyState")
             {
                 FatalErrorInFunction
                     << " Sun direction model can not be sunDirtracking if the "
diff --git a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C
index 1b54ff2ca1a3eda72675618ba98733a1abd2c2db..5645779cd38d4c6aa11efcaa5694891b95a46865 100644
--- a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C
+++ b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C
@@ -60,7 +60,7 @@ Foam::isoAdvection::isoAdvection
 )
 :
     mesh_(alpha1.mesh()),
-    dict_(mesh_.solverDict(alpha1.name())),
+    dict_(mesh_.solution().solverDict(alpha1.name())),
     alpha1_(alpha1),
     alpha1In_(alpha1.primitiveFieldRef()),
     phi_(phi),
diff --git a/src/transportModels/interfaceProperties/interfaceProperties.C b/src/transportModels/interfaceProperties/interfaceProperties.C
index 6d13da17ef73cfe4918cddb53707960aaec42f50..cc41d7952ff09f6ce80d59dae6c380823599b9dc 100644
--- a/src/transportModels/interfaceProperties/interfaceProperties.C
+++ b/src/transportModels/interfaceProperties/interfaceProperties.C
@@ -178,12 +178,12 @@ Foam::interfaceProperties::interfaceProperties
     transportPropertiesDict_(dict),
     nAlphaSmoothCurvature_
     (
-        alpha1.mesh().solverDict(alpha1.name()).
+        alpha1.mesh().solution().solverDict(alpha1.name()).
             getOrDefault<int>("nAlphaSmoothCurvature", 0)
     ),
     cAlpha_
     (
-        alpha1.mesh().solverDict(alpha1.name()).get<scalar>("cAlpha")
+        alpha1.mesh().solution().solverDict(alpha1.name()).get<scalar>("cAlpha")
     ),
 
     sigmaPtr_(surfaceTensionModel::New(dict, alpha1.mesh())),
@@ -256,7 +256,8 @@ void Foam::interfaceProperties::correct()
 
 bool Foam::interfaceProperties::read()
 {
-    alpha1_.mesh().solverDict(alpha1_.name()).readEntry("cAlpha", cAlpha_);
+    alpha1_.mesh().solution()
+        .solverDict(alpha1_.name()).readEntry("cAlpha", cAlpha_);
     sigmaPtr_->readDict(transportPropertiesDict_);
 
     return true;