ENH: improve OpenFOAM I/O for std::vector container

- input: added support (similar to DynamicList)
- output: redirect through UList output, which enables binary etc

- these changes enable support for broadcast and other parallel IO
  for std::vector

ENH: support SubList of std::vector (entire length)

- allows a 'glue' layer for re-casting std::vector to UList etc

.gitmodules

@@ -17,6 +17,3 @@
 [submodule "external-solver"]
 	path = modules/external-solver
 	url = https://develop.openfoam.com/Modules/external-solver.git
-[submodule "turbulence-community"]
-	path = modules/turbulence-community
-	url = https://gitlab.com/openfoam/community/tc-turbulence/turbulence-community.git

CONTRIBUTORS.md

@@ -33,7 +33,6 @@ It is likely incomplete...
 - Alexander Kabat vel Job
 - Thilo Knacke
 - Shannon Leakey
-- Sergey Lesnik
 - Tommaso Lucchini
 - Graham Macpherson
 - Alexey Matveichev
@@ -60,7 +59,6 @@ It is likely incomplete...
 - Vuko Vukcevic
 - Yi Wang
 - Norbert Weber
-- Gregor Weiss
 - Volker Weissmann
 - Henry Weller
 - Niklas Wikstrom

META-INFO/api-info

@@ -1,2 +1,2 @@
-api=2309
+api=2307
 patch=0

applications/solvers/DNS/dnsFoam/createFields.H

@@ -30,6 +30,6 @@ volVectorField U
 
 #include "createPhi.H"
 
-mesh.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 #include "readTurbulenceProperties.H"

applications/solvers/DNS/dnsFoam/dnsFoam.C

@@ -118,7 +118,7 @@ int main(int argc, char *argv[])
                 fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
             );
 
-            pEqn.solve(p.select(piso.finalInnerIter()));
+            pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
             phi = phiHbyA - pEqn.flux();
 

applications/solvers/basic/potentialFoam/createControls.H

@@ -1,6 +1,6 @@
 const dictionary& potentialFlow
 (
-    mesh.solution().solutionDict("potentialFlow")
+    mesh.solutionDict().subDict("potentialFlow")
 );
 
 const int nNonOrthCorr

applications/solvers/basic/potentialFoam/createFields.H

@@ -113,6 +113,6 @@ setRefCell
     PhiRefCell,
     PhiRefValue
 );
-mesh.schemes().setFluxRequired(Phi.name());
+mesh.setFluxRequired(Phi.name());
 
 #include "createMRF.H"

applications/solvers/basic/potentialFoam/overPotentialFoam/createControls.H

@@ -1,6 +1,6 @@
 const dictionary& potentialFlow
 (
-    mesh.solution().solutionDict("potentialFlow")
+    mesh.solutionDict().subDict("potentialFlow")
 );
 
 const int nNonOrthCorr

applications/solvers/combustion/PDRFoam/bEqn.H

@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.schemes().div("div(phi,ft_b_ha_hau)")
+        mesh.divScheme("div(phi,ft_b_ha_hau)")
     )
 );
 

applications/solvers/combustion/PDRFoam/createFields.H

@@ -43,7 +43,7 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-mesh.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::RASModel> turbulence

applications/solvers/combustion/PDRFoam/pEqn.H

@@ -27,7 +27,7 @@ if (pimple.transonic())
             betav*fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -57,7 +57,7 @@ else
             betav*fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {

applications/solvers/combustion/XiFoam/XiDyMFoam/pEqn.H

@@ -30,7 +30,7 @@ if (pimple.transonic())
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -66,7 +66,7 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {

applications/solvers/combustion/XiFoam/XiEngineFoam/pEqn.H

@@ -35,7 +35,7 @@ if (pimple.transonic())
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -71,7 +71,7 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {

applications/solvers/combustion/XiFoam/createFields.H

@@ -44,7 +44,7 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-mesh.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence

applications/solvers/combustion/XiFoam/ftEqn.H

@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.schemes().div("div(phi,ft_b_ha_hau)")
+        mesh.divScheme("div(phi,ft_b_ha_hau)")
     )
 );
 

applications/solvers/combustion/XiFoam/pEqn.H

@@ -29,7 +29,7 @@ if (pimple.transonic())
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -64,7 +64,7 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {

applications/solvers/combustion/chemFoam/YEqn.H

@@ -6,7 +6,7 @@
         solve
         (
             fvm::ddt(rho, Yi) - chemistry.RR(specieI),
-            "Yi"
+            mesh.solver("Yi")
         );
     }
 }

applications/solvers/combustion/fireFoam/YEEqn.H

@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.schemes().div("div(phi,Yi_h)")
+        mesh.divScheme("div(phi,Yi_h)")
     )
 );
 {
@@ -35,7 +35,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.constrain(YiEqn);
 
-            YiEqn.solve("Yi");
+            YiEqn.solve(mesh.solver("Yi"));
 
             fvOptions.correct(Yi);
 

applications/solvers/combustion/fireFoam/createFields.H

@@ -89,7 +89,7 @@ volScalarField p_rgh
     mesh
 );
 
-mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.setFluxRequired(p_rgh.name());
 
 #include "phrghEqn.H"
 

applications/solvers/combustion/fireFoam/pEqn.H

@@ -36,7 +36,7 @@ while (pimple.correctNonOrthogonal())
       + fvOptions(psi, p_rgh, rho.name())
     );
 
-    p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
     if (pimple.finalNonOrthogonalIter())
     {

applications/solvers/combustion/reactingFoam/YEqn.H

@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.schemes().div("div(phi,Yi_h)")
+        mesh.divScheme("div(phi,Yi_h)")
     )
 );
 
@@ -34,7 +34,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.constrain(YiEqn);
 
-            YiEqn.solve("Yi");
+            YiEqn.solve(mesh.solver("Yi"));
 
             fvOptions.correct(Yi);
 

applications/solvers/combustion/reactingFoam/createFields.H

@@ -47,7 +47,7 @@ volScalarField& p = thermo.p();
 
 pressureControl pressureControl(p, rho, pimple.dict(), false);
 
-mesh.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 Info << "Creating turbulence model.\n" << nl;
 autoPtr<compressible::turbulenceModel> turbulence

applications/solvers/combustion/reactingFoam/pEqn.H

@@ -37,7 +37,7 @@ if (pimple.transonic())
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -72,7 +72,7 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {

applications/solvers/combustion/reactingFoam/pcEqn.H

@@ -49,7 +49,7 @@ if (pimple.transonic())
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -89,7 +89,7 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H

@@ -48,7 +48,7 @@ volScalarField& p = thermo.p();
 
 #include "compressibleCreatePhi.H"
 
-mesh.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 Info << "Creating turbulence model.\n" << nl;
 autoPtr<compressible::turbulenceModel> turbulence

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H

@@ -41,7 +41,7 @@ while (pimple.correctNonOrthogonal())
       - fvm::laplacian(rhorAUf, p_rgh)
     );
 
-    p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
     if (pimple.finalNonOrthogonalIter())
     {

applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H

@@ -49,7 +49,7 @@ volScalarField& p = thermo.p();
 
 pressureControl pressureControl(p, rho, pimple.dict(), false);
 
-mesh.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 
 Info << "Creating turbulence model.\n" << nl;

applications/solvers/compressible/rhoCentralFoam/readFluxScheme.H

@@ -1,5 +1,5 @@
 word fluxScheme("Kurganov");
-if (mesh.schemes().dict().readIfPresent("fluxScheme", fluxScheme))
+if (mesh.schemesDict().readIfPresent("fluxScheme", fluxScheme))
 {
     if ((fluxScheme == "Tadmor") || (fluxScheme == "Kurganov"))
     {

applications/solvers/compressible/rhoPimpleAdiabaticFoam/createFields.H

@@ -79,7 +79,7 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-mesh.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 #include "createMRF.H"
 

applications/solvers/compressible/rhoPimpleAdiabaticFoam/pEqn.H

@@ -63,7 +63,7 @@
                 fvOptions(psi, p, rho.name())
             );
 
-            pEqn.solve(p.select(pimple.finalInnerIter()));
+            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
             // Rhie & Chow interpolation (part 2)
             if (pimple.finalNonOrthogonalIter())

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -42,7 +42,7 @@ volVectorField U
 
 pressureControl pressureControl(p, rho, pimple.dict(), false);
 
-mesh.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H

@@ -56,7 +56,7 @@ if (mesh.changing())
         pcorrTypes
     );
 
-    mesh.schemes().setFluxRequired(pcorr.name());
+    mesh.setFluxRequired(pcorr.name());
 
 {
     dimensionedScalar rAUf("rAUf", dimTime, 1.0);
@@ -72,7 +72,7 @@ if (mesh.changing())
             divrhoU()
         );
 
-        pcorrEqn.solve(pcorr.select(pimple.finalInnerIter()));
+        pcorrEqn.solve(mesh.solver(pcorr.select(pimple.finalInnerIter())));
         //Bypass virtual layer
         //mesh.fvMesh::solve(pcorrEqn, d);
 

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.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 #include "createDpdt.H"
 

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H

@@ -56,7 +56,7 @@ if (pimple.transonic())
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -78,7 +78,7 @@ else
     {
         fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -54,7 +54,7 @@ if (pimple.transonic())
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -76,7 +76,7 @@ else
     {
         fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {

applications/solvers/compressible/rhoPimpleFoam/pcEqn.H

@@ -65,7 +65,7 @@ if (pimple.transonic())
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -90,7 +90,7 @@ else
     {
         fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
 
-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {

applications/solvers/compressible/rhoSimpleFoam/createFields.H

@@ -40,7 +40,7 @@ volVectorField U
 
 pressureControl pressureControl(p, rho, simple.dict());
 
-mesh.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createFields.H

@@ -40,7 +40,7 @@ volVectorField U
 
 pressureControl pressureControl(p, rho, simple.dict());
 
-mesh.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence

applications/solvers/compressible/sonicFoam/createFields.H

@@ -36,7 +36,7 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-mesh.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence

applications/solvers/compressible/sonicFoam/sonicLiquidFoam/createFields.H

@@ -46,4 +46,4 @@ volScalarField rho
 
 #include "compressibleCreatePhi.H"
 
-mesh.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());

applications/solvers/electromagnetics/mhdFoam/createFields.H

@@ -114,5 +114,5 @@ label pRefCell = 0;
 scalar pRefValue = 0.0;
 setRefCell(p, piso.dict(), pRefCell, pRefValue);
 
-mesh.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
-mesh.schemes().setFluxRequired(pB.name());
+mesh.setFluxRequired(pB.name());

applications/solvers/electromagnetics/mhdFoam/mhdFoam.C

@@ -127,7 +127,7 @@ int main(int argc, char *argv[])
                     );
 
                     pEqn.setReference(pRefCell, pRefValue);
-                    pEqn.solve(p.select(piso.finalInnerIter()));
+                    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
                     if (piso.finalNonOrthogonalIter())
                     {
@@ -167,7 +167,7 @@ int main(int argc, char *argv[])
                     fvm::laplacian(rABf, pB) == fvc::div(phiB)
                 );
 
-                pBEqn.solve(pB.select(bpiso.finalInnerIter()));
+                pBEqn.solve(mesh.solver(pB.select(bpiso.finalInnerIter())));
 
                 if (bpiso.finalNonOrthogonalIter())
                 {

applications/solvers/electromagnetics/mhdFoam/readBPISOControls.H

@@ -1,3 +1,3 @@
-    const dictionary& Bpiso = mesh.solution().solutionDict("BPISO");
+    const dictionary& Bpiso = mesh.solutionDict().subDict("BPISO");
 
     const int nBcorr = Bpiso.getOrDefault<int>("nCorrectors", 1);

applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H

@@ -5,7 +5,7 @@
         (
             "h",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
@@ -20,7 +20,7 @@
         (
             "Us",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
@@ -34,7 +34,7 @@
         (
             "phis",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
@@ -48,7 +48,7 @@
         (
             "phi2s",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
@@ -67,7 +67,7 @@
         (
             "Sm",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
@@ -82,7 +82,7 @@
         (
             "Sd",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
@@ -96,7 +96,7 @@
         (
             "Sg",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
@@ -112,7 +112,7 @@
         (
             "ps",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
@@ -134,7 +134,7 @@
         (
             "manningField",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
@@ -147,7 +147,7 @@
         (
             "frictionFactor",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
@@ -155,4 +155,4 @@
         dimensionedScalar("one", dimless, 0.01)
     );
 
-    aMesh.schemes().setFluxRequired("h");
+    aMesh.setFluxRequired("h");

applications/solvers/finiteArea/liquidFilmFoam/createFvFields.H

@@ -7,7 +7,7 @@ volVectorField U
     (
         "U",
         runTime.timeName(),
-        mesh.thisDb(),
+        mesh,
         IOobject::NO_READ,
         IOobject::AUTO_WRITE
     ),
@@ -22,7 +22,7 @@ volScalarField H
     (
         "H",
         runTime.timeName(),
-        mesh.thisDb(),
+        mesh,
         IOobject::NO_READ,
         IOobject::AUTO_WRITE
     ),

applications/solvers/finiteArea/liquidFilmFoam/readSolutionControls.H

@@ -1 +1 @@
-loopControl iters(runTime, aMesh.solution().solutionDict(), "solution");
+loopControl iters(runTime, aMesh.solutionDict(), "solution");

applications/solvers/finiteArea/sphereSurfactantFoam/createFaFields.H

@@ -5,7 +5,7 @@ areaScalarField Cs
     (
         "Cs",
         runTime.timeName(),
-        aMesh.thisDb(),
+        mesh,
         IOobject::MUST_READ,
         IOobject::AUTO_WRITE
     ),
@@ -26,7 +26,7 @@ areaVectorField Us
     (
         "Us",
         runTime.timeName(),
-        aMesh.thisDb(),
+        mesh,
         IOobject::NO_READ,
         IOobject::NO_WRITE
     ),
@@ -58,7 +58,7 @@ edgeScalarField phis
     (
         "phis",
         runTime.timeName(),
-        aMesh.thisDb(),
+        mesh,
         IOobject::NO_READ,
         IOobject::NO_WRITE
     ),

applications/solvers/finiteArea/sphereSurfactantFoam/createFaMesh.H

@@ -0,0 +1,2 @@
+    // Create Finite Area mesh
+    faMesh aMesh(mesh);

applications/solvers/finiteArea/sphereSurfactantFoam/createVolFields.H

@@ -7,7 +7,7 @@
         (
             "Cvf",
             runTime.timeName(),
-            mesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
@@ -24,7 +24,7 @@
         (
             "U",
             runTime.timeName(),
-            mesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),

applications/solvers/finiteArea/surfactantFoam/createFaFields.H

@@ -5,7 +5,7 @@ areaScalarField Cs
     (
         "Cs",
         runTime.timeName(),
-        aMesh.thisDb(),
+        mesh,
         IOobject::MUST_READ,
         IOobject::AUTO_WRITE
     ),
@@ -37,7 +37,7 @@ areaVectorField Us
     (
         "Us",
         runTime.timeName(),
-        aMesh.thisDb(),
+        mesh,
         IOobject::MUST_READ,
         IOobject::NO_WRITE
     ),
@@ -51,7 +51,7 @@ edgeScalarField phis
     (
         "phis",
         runTime.timeName(),
-        aMesh.thisDb(),
+        mesh,
         IOobject::NO_READ,
         IOobject::NO_WRITE
     ),

applications/solvers/finiteArea/surfactantFoam/createFaMesh.H

@@ -0,0 +1,2 @@
+    // Create Finite Area mesh
+    faMesh aMesh(mesh);

applications/solvers/finiteArea/surfactantFoam/createVolFields.H

@@ -7,7 +7,7 @@
         (
             "Cvf",
             runTime.timeName(),
-            mesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
@@ -24,7 +24,7 @@
         (
             "U",
             runTime.timeName(),
-            mesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H

@@ -119,7 +119,7 @@ if (p_rgh.needReference())
     );
 }
 
-mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.setFluxRequired(p_rgh.name());
 
 #include "createMRF.H"
 #include "createIncompressibleRadiationModel.H"

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H

@@ -27,7 +27,7 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H

@@ -119,7 +119,7 @@ if (p_rgh.needReference())
     );
 }
 
-mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.setFluxRequired(p_rgh.name());
 
 #include "createMRF.H"
 #include "createIncompressibleRadiationModel.H"

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.schemes().setFluxRequired(p_rgh.name());
+mesh.setFluxRequired(p_rgh.name());
 
 label pRefCell = 0;
 scalar pRefValue = 0.0;

applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/correctPhi.H

@@ -56,7 +56,7 @@ if (mesh.changing())
         pcorrTypes
     );
 
-    mesh.schemes().setFluxRequired(pcorr.name());
+    mesh.setFluxRequired(pcorr.name());
 
 {
     dimensionedScalar rAUf("rAUf", dimTime, 1.0);
@@ -72,7 +72,7 @@ if (mesh.changing())
             divrhoU()
         );
 
-        //pcorrEqn.solve(pcorr.select(pimple.finalInnerIter()));
+        //pcorrEqn.solve(mesh.solver(pcorr.select(pimple.finalInnerIter())));
         //Bypass virtual layer
         const dictionary& d = mesh.solver
             (

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.schemes().setFluxRequired(p_rgh.name());
+mesh.setFluxRequired(p_rgh.name());
 
 label pRefCell = 0;
 scalar pRefValue = 0.0;

applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/pEqn.H

@@ -50,7 +50,7 @@ while (pimple.correctNonOrthogonal())
       - fvm::laplacian(rhorAUf, p_rgh)
     );
 
-    p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
     if (pimple.finalNonOrthogonalIter())
     {

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -52,7 +52,7 @@ while (pimple.correctNonOrthogonal())
         compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
     );
 
-    p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
     if (pimple.finalNonOrthogonalIter())
     {

applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H

@@ -81,7 +81,7 @@ setRefCell
     pRefValue
 );
 
-mesh.schemes().setFluxRequired(p_rgh.name());
+mesh.setFluxRequired(p_rgh.name());
 
 dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 dimensionedScalar totalVolume = sum(mesh.V());

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].schemes().setFluxRequired(p_rghFluid[i].name());
+    fluidRegions[i].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].solution().solutionDict("SIMPLE");
+        fluidRegions[i].solutionDict().subDict("SIMPLE");
 
     setRefCell
     (

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

@@ -1,5 +1,4 @@
-    const dictionary& simple =
+    const dictionary& simple = fluidRegions[i].solutionDict().subDict("SIMPLE");
-        fluidRegions[i].solution().solutionDict("SIMPLE");
 
     const int nNonOrthCorr =
         simple.getOrDefault<int>("nNonOrthogonalCorrectors", 0);

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

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

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].schemes().setFluxRequired(p_rghFluid[i].name());
+    fluidRegions[i].setFluxRequired(p_rghFluid[i].name());
 
     Info<< "    Correcting initialMassFluid\n" << endl;
     initialMassFluid[i] =
         fvc::domainIntegrate(phaseSystemFluid[i].rho()).value();
 
     const dictionary& pimpleDict =
-        fluidRegions[i].solution().solutionDict("PIMPLE");
+        fluidRegions[i].solutionDict().subDict("PIMPLE");
 
     pimpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
 

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

@@ -1,5 +1,4 @@
-    const dictionary& pimpleDict =
+    const dictionary& pimpleDict = mesh.solutionDict().subDict("PIMPLE");
-        mesh.solution().solutionDict("PIMPLE");
 
     Switch faceMomentum
     (

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/solveFluid.H

@@ -1,6 +1,6 @@
 if (finalIter)
 {
-    mesh.data().setFinalIteration(true);
+    mesh.data::add("finalIteration", true);
 }
 
 if (frozenFlow)
@@ -35,5 +35,5 @@ else
 
 if (finalIter)
 {
-    mesh.data().setFinalIteration(false);
+    mesh.data::remove("finalIteration");
 }

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

@@ -1,6 +1,6 @@
 if (finalIter)
 {
-    mesh.data().setFinalIteration(true);
+    mesh.data::add("finalIteration", true);
 }
 
 {
@@ -22,7 +22,7 @@ if (finalIter)
 
         fvOptions.constrain(hEqn);
 
-        hEqn.solve(h.select(finalIter));
+        hEqn.solve(mesh.solver(h.select(finalIter)));
 
         fvOptions.correct(h);
     }
@@ -35,5 +35,5 @@ if (finalIter)
 
 if (finalIter)
 {
-    mesh.data().setFinalIteration(false);
+    mesh.data::remove("finalIteration");
 }

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

@@ -33,7 +33,7 @@
     }
     else
     {
-        EEqn.solve(he.select(finalIter));
+        EEqn.solve(mesh.solver(he.select(finalIter)));
         fvOptions.correct(he);
 
         thermo.correct();

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

@@ -28,7 +28,7 @@
                   - fvc::snGrad(p_rgh)
                 )*mesh.magSf()
             ),
-            U.select(finalIter)
+            mesh.solver(U.select(finalIter))
         );
 
         fvOptions.correct(U);

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

@@ -9,7 +9,7 @@ if (Y.size())
             mesh,
             fields,
             phi,
-            mesh.schemes().div("div(phi,Yi_h)")
+            mesh.divScheme("div(phi,Yi_h)")
         )
     );
 }
@@ -44,7 +44,7 @@ if (Y.size())
 
             fvOptions.constrain(YiEqn);
 
-            YiEqn.solve("Yi");
+            YiEqn.solve(mesh.solver("Yi"));
 
             fvOptions.correct(Yi);
 

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].schemes().setFluxRequired(p_rghFluid[i].name());
+    fluidRegions[i].setFluxRequired(p_rghFluid[i].name());
 
     Info<< "    Adding to radiationFluid\n" << endl;
     radiation.set
@@ -260,8 +260,7 @@ forAll(fluidRegions, i)
     );
 
     const dictionary& pimpleDict =
-        fluidRegions[i].solution().solutionDict("PIMPLE");
+        fluidRegions[i].solutionDict().subDict("PIMPLE");
-
     pimpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
 
     rhoMaxFluid.set

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

@@ -52,12 +52,15 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
         p_rghEqn.solve
         (
-            p_rgh.select
+            mesh.solver
             (
+                p_rgh.select
                 (
-                    oCorr == nOuterCorr-1
+                    (
-                 && corr == nCorr-1
+                        oCorr == nOuterCorr-1
-                 && nonOrth == nNonOrthCorr
+                     && corr == nCorr-1
+                     && nonOrth == nNonOrthCorr
+                    )
                 )
             )
         );

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

@@ -1,5 +1,4 @@
-    const dictionary& pimple =
+    const dictionary& pimple = mesh.solutionDict().subDict("PIMPLE");
-        mesh.solution().solutionDict("PIMPLE");
 
     const int nCorr =
         pimple.getOrDefault<int>("nCorrectors", 1);

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

@@ -1,6 +1,6 @@
 if (finalIter)
 {
-    mesh.data().setFinalIteration(true);
+    mesh.data::add("finalIteration", true);
 }
 
 if (frozenFlow)
@@ -36,5 +36,5 @@ else
 
 if (finalIter)
 {
-    mesh.data().setFinalIteration(false);
+    mesh.data::remove("finalIteration");
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/include/createCoupledRegions.H

@@ -1,8 +1,6 @@
-bool coupled = false;
+fvSolution solutionDict(runTime);
-{
+
-    fvSolution solutionDict(runTime);
+bool coupled(solutionDict.getOrDefault("coupledEnergyField", false));
-    solutionDict.readIfPresent("coupledEnergyField", coupled);
-}
 
 autoPtr<fvMatrix<scalar>> fvMatrixAssemblyPtr;
 

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

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

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

@@ -25,10 +25,10 @@
 
         if (finalIter)
         {
-            mesh.data().setFinalIteration(true);
+            mesh.data::add("finalIteration", true);
         }
 
-        hEqn.solve(h.select(finalIter));
+        hEqn.solve(mesh.solver(h.select(finalIter)));
 
         fvOptions.correct(h);
 
@@ -39,7 +39,7 @@
 
         if (finalIter)
         {
-            mesh.data().setFinalIteration(false);
+            mesh.data::remove("finalIteration");
         }
     }
 }

applications/solvers/heatTransfer/solidFoam/createFields.H

@@ -92,19 +92,21 @@ else
 
 // Consider mesh flux to correct for mesh deformation
 bool meshFluxCorr(false);
+if (mesh.solutionDict().found("SIMPLE"))
 {
-    const dictionary& solutionDict = mesh.solution().solutionDict();
+    meshFluxCorr =
-
+        mesh.solutionDict().subDict("SIMPLE").getOrDefault<bool>
-    const dictionary* subdict = nullptr;
+        (
-
+            "meshFluxCorrection", false
-    if
+        );
-    (
+}
-        ((subdict = solutionDict.findDict("SIMPLE")) != nullptr)
+else if (mesh.solutionDict().found("PIMPLE"))
-     || ((subdict = solutionDict.findDict("PIMPLE")) != nullptr)
+{
-    )
+    meshFluxCorr =
-    {
+        mesh.solutionDict().subDict("PIMPLE").getOrDefault<bool>
-        meshFluxCorr = subdict->getOrDefault("meshFluxCorrection", false);
+        (
-    }
+            "meshFluxCorrection", false
+        );
 }
 
 #include "createRadiationModel.H"

applications/solvers/heatTransfer/solidFoam/solidFoam.C

@@ -68,7 +68,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nEvolving thermodynamics\n" << endl;
 
-    if (mesh.solution().solutionDict().found("SIMPLE"))
+    if (mesh.solutionDict().found("SIMPLE"))
     {
         simpleControl simple(mesh);
 

applications/solvers/heatTransfer/thermoFoam/thermoFoam.C

@@ -65,7 +65,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nEvolving thermodynamics\n" << endl;
 
-    if (mesh.solution().solutionDict().found("SIMPLE"))
+    if (mesh.solutionDict().found("SIMPLE"))
     {
         simpleControl simple(mesh);
 

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.solution().fieldRelaxationFactor("alpha")
+            mesh.fieldRelaxationFactor("alpha")
            *(min(max(alpha + lambda*(Ua & U), zeroAlpha), alphaMax) - alpha);
 
         zeroCells(alpha, inletCells);

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.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 
 Info<< "Reading field pa\n" << endl;
@@ -75,7 +75,7 @@ setRefCell
     paRefCell,
     paRefValue
 );
-mesh.schemes().setFluxRequired(pa.name());
+mesh.setFluxRequired(pa.name());
 
 
 singlePhaseTransportModel laminarTransport(U, phi);

applications/solvers/incompressible/icoFoam/createFields.H

@@ -54,11 +54,5 @@ volVectorField U
 
 label pRefCell = 0;
 scalar pRefValue = 0.0;
-setRefCell
+setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-(
+mesh.setFluxRequired(p.name());
-    p,
-    mesh.solution().solutionDict("PISO"),
-    pRefCell,
-    pRefValue
-);
-mesh.schemes().setFluxRequired(p.name());

applications/solvers/incompressible/icoFoam/icoFoam.C

@@ -141,7 +141,7 @@ int main(int argc, char *argv[])
 
                 pEqn.setReference(pRefCell, pRefValue);
 
-                pEqn.solve(p.select(piso.finalInnerIter()));
+                pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
                 if (piso.finalNonOrthogonalIter())
                 {

applications/solvers/incompressible/nonNewtonianIcoFoam/createFields.H

@@ -36,11 +36,5 @@ singlePhaseTransportModel fluid(U, phi);
 
 label pRefCell = 0;
 scalar pRefValue = 0.0;
-setRefCell
+setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-(
+mesh.setFluxRequired(p.name());
-    p,
-    mesh.solution().solutionDict("PISO"),
-    pRefCell,
-    pRefValue
-);
-mesh.schemes().setFluxRequired(p.name());

applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C

@@ -114,7 +114,7 @@ int main(int argc, char *argv[])
 
                 pEqn.setReference(pRefCell, pRefValue);
 
-                pEqn.solve(p.select(piso.finalInnerIter()));
+                pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
                 if (piso.finalNonOrthogonalIter())
                 {

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.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 Info<< "Creating SRF model\n" << endl;
 autoPtr<SRF::SRFModel> SRF

applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H

@@ -40,7 +40,7 @@ while (pimple.correctNonOrthogonal())
 
     pEqn.setReference(pRefCell, pRefValue);
 
-    pEqn.solve(p.select(pimple.finalInnerIter()));
+    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
     if (pimple.finalNonOrthogonalIter())
     {

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.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 
 singlePhaseTransportModel laminarTransport(U, phi);

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.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 
 //- Overset specific

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H

@@ -32,7 +32,7 @@ while (pimple.correctNonOrthogonal())
 
     pEqn.setReference(pRefCell, pRefValue);
 
-    pEqn.solve(p.select(pimple.finalInnerIter()));
+    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
     if (pimple.finalNonOrthogonalIter())
     {

applications/solvers/incompressible/pimpleFoam/pEqn.H

@@ -45,7 +45,7 @@ while (pimple.correctNonOrthogonal())
 
     pEqn.setReference(pRefCell, pRefValue);
 
-    pEqn.solve(p.select(pimple.finalInnerIter()));
+    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
     if (pimple.finalNonOrthogonalIter())
     {

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.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 
 singlePhaseTransportModel laminarTransport(U, phi);

applications/solvers/incompressible/pisoFoam/pEqn.H

@@ -26,7 +26,7 @@ while (piso.correctNonOrthogonal())
 
     pEqn.setReference(pRefCell, pRefValue);
 
-    pEqn.solve(p.select(piso.finalInnerIter()));
+    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
     if (piso.finalNonOrthogonalIter())
     {

applications/solvers/incompressible/pisoFoam/pisoFoam.C

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -106,9 +105,6 @@ int main(int argc, char *argv[])
 
         #include "CourantNo.H"
 
-        // Update settings from the control dictionary
-        piso.read();
-
         // Pressure-velocity PISO corrector
         {
             #include "UEqn.H"

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.schemes().setFluxRequired(h.name());
+mesh.setFluxRequired(h.name());

applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C

@@ -135,7 +135,7 @@ int main(int argc, char *argv[])
                       - fvm::laplacian(ghrAUf, h)
                     );
 
-                    hEqn.solve(h.select(pimple.finalInnerIter()));
+                    hEqn.solve(mesh.solver(h.select(pimple.finalInnerIter())));
 
                     if (pimple.finalNonOrthogonalIter())
                     {

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.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 Info<< "Creating SRF model\n" << endl;
 autoPtr<SRF::SRFModel> SRF(SRF::SRFModel::New(Urel));

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.schemes().setFluxRequired(p.name());
+mesh.setFluxRequired(p.name());
 
 
 singlePhaseTransportModel laminarTransport(U, phi);