Merge branch 'nextRelease' of ssh://git.code.sf.net/p/foam-extend/foam-extend-4.0 into nextRelease

.gitignore

@@ -45,6 +45,7 @@ lnInclude
 linux*Gcc*/
 linux*Icc*/
 mingw*Gcc*/
+linuxARM*/
 darwin*Gcc*/
 darwin*Intel*/
 linuxming*/

CMakeLists.txt

@@ -1,7 +1,7 @@
 # /*-------------------------------------------------------------------------*\
 #   =========                 |
 #   \\      /  F ield         | foam-extend: Open Source CFD
-#    \\    /   O peration     | Version:     4.0
+#    \\    /   O peration     | Version:     4.1
 #     \\  /    A nd           | Web:         http://www.foam-extend.org
 #      \\/     M anipulation  | For copyright notice see file Copyright
 # -----------------------------------------------------------------------------
@@ -23,7 +23,7 @@
 #
 # Description
 #     CMakeLists.txt file for implementing a test harness for the compilation
-#     and test of foam-extend-3.2 using Kitware CTest/CMake/CDash
+#     and test of foam-extend-4.1 using Kitware CTest/CMake/CDash
 #
 #     The results will be submitted to the CDash server identified by the file
 #     CTestConfig.cmake
@@ -36,7 +36,7 @@
 
 cmake_minimum_required (VERSION 2.8)
 
-PROJECT(foam-extend-4.0)
+PROJECT(foam-extend-4.1)
 
 #-----------------------------------------------------------------------------
 # Utility functions
@@ -52,13 +52,16 @@ function(GetHostName var)
         execute_process(
             COMMAND         hostname
             OUTPUT_VARIABLE thisHostname
+            OUTPUT_STRIP_TRAILING_WHITESPACE
         )
     else()
         execute_process(
             COMMAND         hostname -f
             OUTPUT_VARIABLE thisHostname
+            OUTPUT_STRIP_TRAILING_WHITESPACE
         )
     endif()
+
     set(${var} ${thisHostname} PARENT_SCOPE)
 endfunction()
 

ReleaseNotes

@@ -1,6 +1,7 @@
 #                            -*- mode: org; -*-
 #
 #+TITLE:             *Release notes for foam-extend-4.0*
+#+TITLE:                  *Version 4.0 - Guimaraes*
 #+AUTHOR:                foam-extend administrators:
 #+AUTHOR:                        Hrvoje Jasak
 #+AUTHOR:                       Håkan Nilsson

applications/solvers/basic/potentialDyMFoam/potentialDyMFoam.C

@@ -72,6 +72,11 @@ int main(int argc, char *argv[])
             U.internalField() = vector::zero;
         }
 
+        p.correctBoundaryConditions();
+        U.correctBoundaryConditions();
+
+        phi == (fvc::interpolate(U) & mesh.Sf());
+
 #       include "volContinuity.H"
 #       include "meshCourantNo.H"
 

applications/solvers/compressible/sonicDyMFoam/UEqn.H

@@ -7,7 +7,10 @@
 
     UEqn.relax
     (
-        mesh.solutionDict().equationRelaxationFactor(U.select(pimple.finalIter()))
+        mesh.solutionDict().equationRelaxationFactor
+        (
+            U.select(pimple.finalIter())
+        )
     );
 
     solve(UEqn == -fvc::grad(p));

applications/solvers/compressible/sonicFoam/pEqn.H

@@ -27,8 +27,6 @@
             fvm::ddt(psis, p)
           + fvm::div(phid, p)
             // Convective flux relaxation terms
-          + fvm::SuSp(-divPhid, p)
-          + divPhid*p
           + fvc::div(phid2)
           - fvm::laplacian(rho*rUA, p)
         );

applications/solvers/compressible/steadyCompressibleFoam/UEqn.H

@@ -3,8 +3,7 @@
 
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::div(phi, U)
-      + fvm::div(phi, U)
       + turbulence->divDevRhoReff()
     );
 

applications/solvers/compressible/steadyCompressibleFoam/hEqn.H

@@ -1,13 +1,14 @@
 {
     // Solve the enthalpy equation in total enthalpy formulation (see K)
 
+    T.storePrevIter();
+    
     K = 0.5*(magSqr(U));
 
     fvScalarMatrix hEqn
     (
-        fvm::ddt(rho, h)
+        fvm::div(phi, h)
-      + fvc::ddt(rho, K)
+      + fvm::SuSp(-fvc::div(phi), h)
-      + fvm::div(phi, h)
       + fvc::div(phi, K)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
@@ -16,7 +17,6 @@
     );
 
     hEqn.relax();
-
     hEqn.solve();
 
     // Bounding of enthalpy taken out

applications/solvers/compressible/steadyCompressibleFoam/pEqn.H

@@ -21,21 +21,11 @@
 
         p.storePrevIter();
 
-        volScalarField divPhid
-        (
-            "divPhid",
-            fvc::div(phid)
-        );
-
         while (pimple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
             (
-                fvm::ddt(psis, p)
+                fvm::div(phid, p)
-              + fvm::div(phid, p)
-                // Convective flux relaxation terms
-              + fvm::SuSp(-divPhid, p)
-              + divPhid*p
               + fvc::div(phid2)
               - fvm::laplacian(rho*rUA, p)
             );
@@ -49,7 +39,8 @@
             }
         }
 
-#       include "compressibleContinuityErrs.H"
+        // Use incompressible continuity error check: div(rho U) = 0
+#       include "continuityErrs.H"
 
         // Relax the pressure
         p.relax();

applications/solvers/compressible/steadyCompressibleMRFFoam/Make/options

@@ -2,7 +2,7 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
-    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/solvers/compressible/steadyCompressibleMRFFoam/UEqn.H

@@ -1,10 +1,6 @@
-    // Solve the momentum equation
-    U.storePrevIter();
-
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::div(phi, U)
-      + fvm::div(phi, U)
       + turbulence->divDevRhoReff()
     );
 

applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H

@@ -24,6 +24,7 @@
     );
 
     hEqn.relax();
+    hEqn.solve();
 
     // Bounding of enthalpy taken out
     thermo.correct();

applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H

@@ -5,13 +5,30 @@
     Urel == U;
     mrfZones.relativeVelocity(Urel);
 
+    // Bound the relative velocity to preserve the i to h conversion bound
+    // HJ, 22/Mar/2017
+    volScalarField magUrel = mag(Urel);
+
+    if (max(magUrel) > UMax)
+    {
+        volScalarField UrelLimiter = pos(magUrel - UMax)*UMax/(magUrel + smallU)
+            + neg(magUrel - UMax);
+        UrelLimiter.max(scalar(0));
+        UrelLimiter.min(scalar(1));
+
+        Urel *= UrelLimiter;
+        Urel.correctBoundaryConditions();
+    }
+
     // Create rotational velocity (= omega x r)
     Urot == U - Urel;
 
+    T.storePrevIter();
+
     fvScalarMatrix iEqn
     (
-        fvm::ddt(rho, i)
+        fvm::div(phi, i)
-      + fvm::div(phi, i)
+      + fvm::SuSp(-fvc::div(phi), i)
       - fvm::laplacian(turbulence->alphaEff(), i)
      ==
         // Viscous heating: note sign (devRhoReff has a minus in it)
@@ -19,7 +36,6 @@
     );
 
     iEqn.relax();
-
     iEqn.solve();
 
     // Calculate enthalpy out of rothalpy

applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H

@@ -16,7 +16,7 @@
         mrfZones.correctBoundaryVelocity(U);
 
         // Calculate phi for boundary conditions
-        phi = rhof*fvc::interpolate(U) & mesh.Sf();
+        phi = rhof*(fvc::interpolate(U) & mesh.Sf());
 
         surfaceScalarField phid2 = rhoReff/rhof*phi;
 
@@ -29,21 +29,11 @@
 
         p.storePrevIter();
 
-        volScalarField divPhid
-        (
-            "divPhid",
-            fvc::div(phid)
-        );
-
         while (pimple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
             (
-                fvm::ddt(psis, p)
+                fvm::div(phid, p)
-              + fvm::div(phid, p)
-                // Convective flux relaxation terms
-              + fvm::SuSp(-divPhid, p)
-              + divPhid*p
               + fvc::div(phid2)
               - fvm::laplacian(rho*rUA, p)
             );
@@ -57,7 +47,8 @@
             }
         }
 
-#       include "compressibleContinuityErrs.H"
+        // Use incompressible continuity error check: div(rho U) = 0
+#       include "continuityErrs.H"
 
         // Relax the pressure
         p.relax();

applications/solvers/compressible/steadyCompressibleSRFFoam/UEqn.H

@@ -3,8 +3,7 @@
 
     fvVectorMatrix UrelEqn
     (
-        fvm::ddt(rho, Urel)
+        fvm::div(phi, Urel)
-      + fvm::div(phi, Urel)
       + turbulence->divDevRhoReff()
       + rho*SRF->Su()
     );

applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H

@@ -1,10 +1,11 @@
 {
     // Solve the rothalpy equation
+    T.storePrevIter();
 
     fvScalarMatrix iEqn
     (
-        fvm::ddt(rho, i)
+        fvm::div(phi, i)
-      + fvm::div(phi, i)
+      + fvm::SuSp(-fvc::div(phi), i)
       - fvm::laplacian(turbulence->alphaEff(), i)
      ==
         // Viscous heating: note sign (devRhoReff has a minus in it)
@@ -12,7 +13,6 @@
     );
 
     iEqn.relax();
-
     iEqn.solve();
 
     // Calculate enthalpy out of rothalpy

applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H

@@ -20,21 +20,11 @@
 
         p.storePrevIter();
 
-        volScalarField divPhid
-        (
-            "divPhid",
-            fvc::div(phid)
-        );
-
         while (pimple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
             (
-                fvm::ddt(psis, p)
+                fvm::div(phid, p)
-              + fvm::div(phid, p)
-                // Convective flux relaxation terms
-              + fvm::SuSp(-divPhid, p)
-              + divPhid*p
               + fvc::div(phid2)
               - fvm::laplacian(rho*rUrelA, p)
             );
@@ -48,7 +38,8 @@
             }
         }
 
-#       include "compressibleContinuityErrs.H"
+        // Use incompressible continuity error check: div(rho U) = 0
+#       include "continuityErrs.H"
 
         // Relax the pressure
         p.relax();

applications/solvers/compressible/steadyUniversalFoam/UEqn.H

@@ -3,8 +3,7 @@
 
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::div(phi, U)
-      + fvm::div(phi, U)
       + turbulence->divDevRhoReff()
     );
 

applications/solvers/compressible/steadyUniversalFoam/hEqn.H

@@ -10,8 +10,8 @@
 
     fvScalarMatrix hEqn
     (
-        fvm::ddt(rho, h)
+        fvm::div(phi, h)
-      + fvm::div(phi, h)
+      + fvm::SuSp(-fvc::div(phi), h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
         fvc::div(faceU, p, "div(U,p)")
@@ -22,7 +22,6 @@
     );
 
     hEqn.relax();
-
     hEqn.solve();
 
     // Bounding of enthalpy taken out

applications/solvers/compressible/steadyUniversalFoam/pEqn.H

@@ -24,8 +24,7 @@
         {
             fvScalarMatrix pEqn
             (
-                fvm::ddt(psis, p)
+                fvm::div(phid, p)
-              + fvm::div(phid, p)
               + fvc::div(phid2)
               - fvm::laplacian(rho*rUA, p)
             );
@@ -39,8 +38,8 @@
             }
         }
 
-        // Use custom continuity error check
+        // Use incompressible continuity error check: div(rho U) = 0
-#       include "universalContinuityErrs.H"
+#       include "continuityErrs.H"
 
         // Relax the pressure
         p.relax();

applications/solvers/compressible/steadyUniversalFoam/rhoFromP.H

@@ -12,4 +12,5 @@
     rho = Foam::min(rho, rhoMax);
     rho = Foam::max(rho, rhoMin);
     rho.relax();
+    rho.correctBoundaryConditions();
 }

applications/solvers/compressible/steadyUniversalMRFFoam/UEqn.H

@@ -3,8 +3,7 @@
 
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::div(phi, U)
-      + fvm::div(phi, U)
       + turbulence->divDevRhoReff()
     );
 

applications/solvers/compressible/steadyUniversalMRFFoam/iEqn.H

@@ -3,14 +3,29 @@
     Urel == U;
     mrfZones.relativeVelocity(Urel);
 
+    // Bound the relative velocity to preserve the i to h conversion bound
+    // HJ, 22/Mar/2017
+    volScalarField magUrel = mag(Urel);
+
+    if (max(magUrel) > UMax)
+    {
+        volScalarField UrelLimiter = pos(magUrel - UMax)*UMax/(magUrel + smallU)
+            + neg(magUrel - UMax);
+        UrelLimiter.max(scalar(0));
+        UrelLimiter.min(scalar(1));
+
+        Urel *= UrelLimiter;
+        Urel.correctBoundaryConditions();
+    }
+
     // Create rotational velocity (= omega x r)
     Urot == U - Urel;
 
     fvScalarMatrix iEqn
     (
-          fvm::ddt(rho, i)
+        fvm::div(phi, i)
-        + fvm::div(phi, i)
+      + fvm::SuSp(-fvc::div(phi), i)
-        - fvm::laplacian(turbulence->alphaEff(), i)
+      - fvm::laplacian(turbulence->alphaEff(), i)
      ==
       // Viscous heating: note sign (devRhoReff has a minus in it)
        - (turbulence->devRhoReff() && fvc::grad(U))
@@ -23,8 +38,8 @@
     // From rothalpy, calculate enthalpy after solution of rothalpy equation
     h = i + 0.5*(magSqr(Urot) - magSqr(Urel));
     h.correctBoundaryConditions();
+
     // Update thermo for new h
     thermo.correct();
     psis = thermo.psi()/thermo.Cp()*thermo.Cv();
 }
-

applications/solvers/compressible/steadyUniversalMRFFoam/pEqn.H

@@ -33,8 +33,7 @@
         {
             fvScalarMatrix pEqn
             (
-                fvm::ddt(psis, p)
+                fvm::div(phid, p)
-              + fvm::div(phid, p)
               + fvc::div(phid2)
               - fvm::laplacian(rho*rUA, p)
             );
@@ -48,8 +47,8 @@
             }
         }
 
-        // Use custom continuity error check
+        // Use incompressible continuity error check: div(rho U) = 0
-#       include "universalContinuityErrs.H"
+#       include "continuityErrs.H"
 
         // Relax the pressure
         p.relax();

applications/solvers/compressible/steadyUniversalMRFFoam/rhoFromP.H

@@ -12,4 +12,5 @@
     rho = Foam::min(rho, rhoMax);
     rho = Foam::max(rho, rhoMin);
     rho.relax();
+    rho.correctBoundaryConditions();
 }

applications/solvers/compressible/steadyUniversalMRFFoam/universalContinuityErrs.H

@@ -1,49 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Global
-    continuityErrs
-
-Description
-    Calculates and prints the continuity errors.
-
-\*---------------------------------------------------------------------------*/
-
-{
-    volScalarField contErr = fvc::ddt(rho) + fvc::div(phi);
-
-    sumLocalContErr = runTime.deltaT().value()*
-        mag(contErr)().weightedAverage(mesh.V()).value();
-
-    globalContErr = runTime.deltaT().value()*
-        contErr.weightedAverage(mesh.V()).value();
-
-    cumulativeContErr += globalContErr;
-
-    Info<< "time step continuity errors : sum local = " << sumLocalContErr
-        << ", global = " << globalContErr
-        << ", cumulative = " << cumulativeContErr
-        << endl;
-}
-
-// ************************************************************************* //

applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C

@@ -94,7 +94,7 @@ int main(int argc, char *argv[])
             
             solve
             (
-                tpEqn
+                tpEqn()
              ==
               - fvc::div(U)
             );
@@ -112,7 +112,9 @@ int main(int argc, char *argv[])
             U.correctBoundaryConditions();
             p.correctBoundaryConditions();
 
-            phi = (fvc::interpolate(U) & mesh.Sf()) + tpEqn().flux() + tpresSource;
+            phi = (fvc::interpolate(U) & mesh.Sf())
+                + tpEqn().flux()
+                + tpresSource;
         }
         
         // Make flux relative in rotating zones

applications/solvers/coupled/MRFPorousFoam/UEqn.H

@@ -1,10 +1,4 @@
 {
-    volScalarField divPhi
-    (
-        "divPhi",
-        fvc::div(phi)
-    );
-
     // Update boundary velocity for consistency with the flux
     mrfZones.correctBoundaryVelocity(U);
 
@@ -15,13 +9,10 @@
       + turbulence->divDevReff()
     );
 
-    // Add MRF sources
+    // Add MRF and porous sources implicitly.  HJ, 18/Nov/2017
-    mrfZones.addCoriolis(UEqn);
-
-    // Add porous sources
     tmp<volTensorField> tTU;
 
-    if (addPorosity)
+    if (addMRF || addPorosity)
     {
         tTU = tmp<volTensorField>
         (
@@ -41,6 +32,10 @@
         );
         volTensorField& TU = tTU();
 
+        // Add implicit MRF source as a Hodge dual of the rotational velocity
+        TU += *mrfZones.omega();
+        
+        // Add implicit resistance
         pZones.addResistance(UEqn, TU);
 
         trTU = inv(TU + tensor(I)*UEqn.A());
@@ -58,7 +53,7 @@
     // Insert momentum equation
     UpEqn.insertEquation(0, UEqn);
 
-    if (addPorosity)
+    if (addMRF || addPorosity)
     {
         // Manually over-ride the 3x3 block to handle the off-diagonal
         // part of the Ap coefficient
@@ -68,11 +63,43 @@
 
         const scalarField& V = mesh.V().field();
 
-        // Note: insertion should only happen in porous cell zones
+        // Note: insertion should only happen in MRF and porous cell zones
         // HJ, 14/Mar/2016
 
+        // Warning.  Possible problem with a zone which is both MRF and porous
+        // The solution is to do the loop below everywhere
+        // Currently only inserting zones for efficiency.  HJ, 18/Nov/2017
         register label cellI;
 
+        forAll (mrfZones, mrfZoneI)
+        {
+            const labelList& curZoneCells = mrfZones[mrfZoneI].zone();
+
+            // Loop over all cells in the zone
+            forAll (curZoneCells, zcI)
+            {
+                cellI = curZoneCells[zcI];
+
+                const scalar& cellV = V[cellI];
+
+                const tensor& cellTU = TUIn[cellI];
+
+                CoeffField<vector4>::squareType& cellDD = DD[cellI];
+
+                cellDD(0, 0) += cellV*cellTU.xx();
+                cellDD(0, 1) += cellV*cellTU.xy();
+                cellDD(0, 2) += cellV*cellTU.xz();
+
+                cellDD(1, 0) += cellV*cellTU.yx();
+                cellDD(1, 1) += cellV*cellTU.yy();
+                cellDD(2, 2) += cellV*cellTU.yz();
+
+                cellDD(2, 0) += cellV*cellTU.zx();
+                cellDD(2, 1) += cellV*cellTU.zy();
+                cellDD(2, 2) += cellV*cellTU.zz();
+            }
+        }
+
         forAll (pZones, pZoneI)
         {
             const labelList& curZoneCells = pZones[pZoneI].zone();

applications/solvers/coupled/MRFPorousFoam/createMRF.H

@@ -1,2 +1,9 @@
     MRFZones mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U);
+
+    bool addMRF = false;
+
+    if (!mrfZones.empty())
+    {
+        addMRF = true;
+    }

applications/solvers/coupled/MRFPorousFoam/pEqn.H

@@ -2,7 +2,7 @@
 tmp<fvScalarMatrix> tpEqn;
 tmp<surfaceScalarField> tpresSource;
 
-if (addPorosity)
+if (addMRF || addPorosity)
 {
     // Collect pressure source with tensorial 1/Ap
     const volTensorField& rTU = trTU();
@@ -10,7 +10,7 @@ if (addPorosity)
     tpresSource =
     (
         (mesh.Sf() & fvc::interpolate(rTU))
-        & fvc::interpolate(fvc::grad(p))
+      & fvc::interpolate(fvc::grad(p))
     );
     const surfaceScalarField& presSource = tpresSource();
 

applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H

@@ -1,5 +1,6 @@
     label pRefCell = 0;
     scalar pRefValue = 0;
+
     setRefCell
     (
         p,

applications/solvers/coupled/pUCoupledFoam/UEqn.H

@@ -1,10 +1,4 @@
 {
-    volScalarField divPhi
-    (
-        "divPhi",
-        fvc::div(phi)
-    );
-
     // Momentum equation
     fvVectorMatrix UEqn
     (

applications/solvers/engine/sonicTurbDyMEngineFoam/UEqn.H

@@ -8,7 +8,10 @@
 
     UEqn.relax
     (
-        mesh.solutionDict().equationRelaxationFactor(U.select(pimple.finalIter()))
+        mesh.solutionDict().equationRelaxationFactor
+        (
+            U.select(pimple.finalIter())
+        )
     );
 
     solve(UEqn == -fvc::grad(p));

applications/solvers/engine/sonicTurbDyMEngineFoam/hEqn.H

@@ -12,4 +12,8 @@
     hEqn.solve();
 
     thermo.correct();
+
+    // Recalculate density
+    rho = thermo.rho();
+    rho.correctBoundaryConditions();
 }

applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H

@@ -1,6 +1,4 @@
 {
-    rho = thermo.rho();
-
     rUA = 1.0/UEqn.A();
     U = rUA*UEqn.H();
 
@@ -56,12 +54,10 @@
         p.relax();
     }
 
-#   include "rhoEqn.H"
 #   include "compressibleContinuityErrs.H"
 
     U -= rUA*fvc::grad(p);
     U.correctBoundaryConditions();
 
     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
-    dpdt = fvc::ddt(p);
 }

applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C

@@ -79,14 +79,15 @@ int main(int argc, char *argv[])
 
         Info<< "Crank angle = " << runTime.theta() << " CA-deg" << endl;
 
-        // Make flux absolute
+        // Make the fluxes absolute (using the ddt(rho, U) scheme)
-        phi += meshFlux;
+        phi += fvc::interpolate(rho)*fvc::meshPhi(rho, U);
 
         bool meshChanged = mesh.update();
 
 #       include "volContinuity.H"
 
-        mesh.setBoundaryVelocity(U);
+        // Make the fluxes relative (using the ddt(rho, U) scheme)
+        phi -= fvc::interpolate(rho)*fvc::meshPhi(rho, U);
 
         if (meshChanged)
         {
@@ -95,14 +96,7 @@ int main(int argc, char *argv[])
             rho.correctBoundaryConditions();
         }
 
-        meshFlux = fvc::interpolate(rho)*fvc::meshPhi(rho, U);
+        if (meshChanged)
-
-        phi = fvc::interpolate(rho)
-            *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
-
-        DpDt = dpdt + fvc::div(phi/fvc::interpolate(rho), p)
-            - fvc::div(phi/fvc::interpolate(rho) + fvc::meshPhi(U))*p;
-
         {
 #           include "compressibleCourantNo.H"
         }
@@ -111,22 +105,20 @@ int main(int argc, char *argv[])
         while (pimple.loop())
         {
 #           include "rhoEqn.H"
+#           include "hEqn.H"
 #           include "UEqn.H"
 
             // --- PISO loop
             while (pimple.correct())
             {
 #               include "pEqn.H"
-#               include "hEqn.H"
             }
+
+            turbulence->correct();
         }
 
-        turbulence->correct();
-
 #       include "logSummary.H"
 
-        rho = thermo.rho();
-
         runTime.write();
 
 #       include "infoDataOutput.H"

applications/solvers/immersedBoundary/SOLVER_COOKBOOK

@@ -1,49 +0,0 @@
-Make/options: add
-
-    -I$(LIB_SRC)/triSurface/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I../immersedBoundary/lnInclude
-
-
-    -ltriSurface \
-    -lmeshTools \
-    -L$(FOAM_USER_LIBBIN) -limmersedBoundary \
-
-
-createFields.H: add
-
-#   include "createIbMasks.H"
-
-
-solver:
-
-1) add immersed boundary headers
-
-#include "immersedBoundaryFvPatch.H"
-#include "immersedBoundaryAdjustPhi.H"
-
-2) on calculation of face fluxes (or their components), mask with faceIbMask
-
-3) before to adjustPhi, add:
-
-        // Adjust immersed boundary fluxes
-        immersedBoundaryAdjustPhi(phi, U);
-
-4) on explicit updates, add correctBoundaryConditions();
-
-eg.
-
-        U = fvc::reconstruct(phi);
-        U.correctBoundaryConditions();
-
-5) On reports of continuity error, add masking:
-
-        Info<< "IB-masked continuity error = "
-            << mag(cellIbMask*fvc::div(phi))().weightedAverage(mesh.V()).value()
-            << endl;
-
-  or use immersedBoundaryContinuityErrs.H
-
-5) Chenge Courant number check to be IB-sensitive, using
-
-  immersedBoundaryCourantNo.H

applications/solvers/immersedBoundary/icoDyMIbFoam/Make/files

@@ -1,3 +0,0 @@
-icoDyMIbFoam.C
-
-EXE = $(FOAM_APPBIN)/icoDyMIbFoam

applications/solvers/immersedBoundary/icoDyMIbFoam/createFields.H

@@ -1,56 +0,0 @@
-    Info<< "Reading transportProperties\n" << endl;
-
-    IOdictionary transportProperties
-    (
-        IOobject
-        (
-            "transportProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    );
-
-    dimensionedScalar nu
-    (
-        transportProperties.lookup("nu")
-    );
-
-
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-#   include "createPhi.H"
-
-
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, pimple.dict(), pRefCell, pRefValue);
-    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/immersedBoundary/icoDyMIbFoam/icoDyMIbFoam.C

@@ -1,152 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    icoDyMOversetFoam
-
-Description
-    Transient solver for incompressible, laminar flow of Newtonian fluids
-    with dynamic mesh and immersed boundary mesh support.
-
-Author
-    Hrvoje Jasak, Wikki Ltd.  All rights reserved
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "dynamicFvMesh.H"
-#include "immersedBoundaryFvPatch.H"
-#include "immersedBoundaryAdjustPhi.H"
-#include "pimpleControl.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createDynamicFvMesh.H"
-
-    pimpleControl pimple(mesh);
-
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.loop())
-    {
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        // Make the fluxes absolute
-        fvc::makeAbsolute(phi, U);
-
-        bool meshChanged = mesh.update();
-        reduce(meshChanged, orOp<bool>());
-        Info<< "Mesh update" << meshChanged << endl;
-#       include "createIbMasks.H"
-
-        // Make the fluxes relative to the mesh motion
-        fvc::makeRelative(phi, U);
-
-#       include "CourantNo.H"
-
-        // Pressure-velocity corrector
-        while (pimple.loop())
-        {
-            fvVectorMatrix UEqn
-            (
-                fvm::ddt(U)
-              + fvm::div(phi, U)
-              - fvm::laplacian(nu, U)
-            );
-
-            if (pimple.momentumPredictor())
-            {
-                solve(UEqn == -fvc::grad(p));
-            }
-
-            // --- PISO loop
-            while (pimple.correct())
-            {
-                volScalarField rUA = 1.0/UEqn.A();
-
-                U = rUA*UEqn.H();
-                // Immersed boundary update
-                U.correctBoundaryConditions();
-
-                phi = faceIbMask*(fvc::interpolate(U) & mesh.Sf());
-
-                // Adjust immersed boundary fluxes
-                immersedBoundaryAdjustPhi(phi, U);
-                adjustPhi(phi, U, p);
-
-                // Non-orthogonal pressure corrector loop
-                while (pimple.correctNonOrthogonal())
-                {
-                    fvScalarMatrix pEqn
-                    (
-                        fvm::laplacian(rUA, p) == fvc::div(phi)
-                    );
-
-                    pEqn.setReference(pRefCell, pRefValue);
-                    pEqn.solve
-                    (
-                        mesh.solutionDict().solver
-                        (
-                            p.select(pimple.finalInnerIter())
-                        )
-                    );
-
-                    if (pimple.finalNonOrthogonalIter())
-                    {
-                        phi -= pEqn.flux();
-                    }
-                }
-
-#               include "immersedBoundaryContinuityErrs.H"
-
-                // Make the fluxes relative to the mesh motion
-                fvc::makeRelative(phi, U);
-
-                U -= rUA*fvc::grad(p);
-                U.correctBoundaryConditions();
-            }
-        }
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return(0);
-}
-
-
-// ************************************************************************* //

applications/solvers/immersedBoundary/icoIbFoam/Make/files

@@ -1,3 +0,0 @@
-icoIbFoam.C
-
-EXE = $(FOAM_APPBIN)/icoIbFoam

applications/solvers/immersedBoundary/icoIbFoam/Make/options

@@ -1,13 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude
-
-EXE_LIBS = \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lsurfMesh \
-    -lsampling \
-    -ldynamicMesh \
-    -limmersedBoundary \
-    -llduSolvers

applications/solvers/immersedBoundary/icoIbFoam/icoIbFoam.C

@@ -1,139 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    icoIbFoam
-
-Description
-    Transient solver for incompressible, laminar flow of Newtonian fluids
-    with immersed boundary support.
-
-Author
-    Hrvoje Jasak, Wikki Ltd.  All rights reserved
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "immersedBoundaryFvPatch.H"
-#include "immersedBoundaryAdjustPhi.H"
-#include "pimpleControl.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-#   include "setRootCase.H"
-
-#   include "createTime.H"
-#   include "createMesh.H"
-
-    pimpleControl pimple(mesh);
-
-#   include "createIbMasks.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.loop())
-    {
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-#       include "CourantNo.H"
-
-        // Pressure-velocity corrector
-        while (pimple.loop())
-        {
-            fvVectorMatrix UEqn
-            (
-                fvm::ddt(U)
-              + fvm::div(phi, U)
-              - fvm::laplacian(nu, U)
-            );
-
-            if (pimple.momentumPredictor())
-            {
-                solve(UEqn == -fvc::grad(p));
-            }
-
-            // --- PISO loop
-            while (pimple.correct())
-            {
-                volScalarField rUA = 1.0/UEqn.A();
-
-                U = rUA*UEqn.H();
-                // Immersed boundary update
-                U.correctBoundaryConditions();
-
-                phi = faceIbMask*(fvc::interpolate(U) & mesh.Sf());
-
-                // Adjust immersed boundary fluxes
-                immersedBoundaryAdjustPhi(phi, U);
-                adjustPhi(phi, U, p);
-
-                // Non-orthogonal pressure corrector loop
-                while (pimple.correctNonOrthogonal())
-                {
-                    fvScalarMatrix pEqn
-                    (
-                        fvm::laplacian(rUA, p) == fvc::div(phi)
-                    );
-
-                    pEqn.setReference(pRefCell, pRefValue);
-                    pEqn.solve
-                    (
-                        mesh.solutionDict().solver
-                        (
-                            p.select(pimple.finalInnerIter())
-                        )
-                    );
-
-                    if (pimple.finalNonOrthogonalIter())
-                    {
-                        phi -= pEqn.flux();
-                    }
-                }
-
-#               include "immersedBoundaryContinuityErrs.H"
-
-                U -= rUA*fvc::grad(p);
-                U.correctBoundaryConditions();
-            }
-        }
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/immersedBoundary/interIbFoam/Make/files

@@ -1,3 +0,0 @@
-interIbFoam.C
-
-EXE = $(FOAM_APPBIN)/interIbFoam

applications/solvers/immersedBoundary/interIbFoam/UEqn.H

@@ -1,27 +0,0 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U)
-      + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muEff, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    );
-
-    UEqn.relax();
-
-    if (pimple.momentumPredictor())
-    {
-        solve
-        (
-            UEqn
-         ==
-          - gh*fvc::grad(rho)
-          - fvc::grad(pd)
-        );
-    }

applications/solvers/immersedBoundary/interIbFoam/alphaEqn.H

@@ -1,56 +0,0 @@
-{
-    word alphaScheme("div(phi,alpha)");
-    word alpharScheme("div(phirb,alpha)");
-
-    // New formulation of phir: Co condition
-    surfaceScalarField phir
-    (
-        "phir",
-        faceIbMask*interface.cAlpha()*interface.nHatf()*
-        min
-        (
-            scalar(0.5)/    // This is ref compression Co number for cAlpha = 1
-            (
-                mesh.surfaceInterpolation::deltaCoeffs()*
-                mesh.time().deltaT()
-            ),
-            mag(phi/mesh.magSf())
-          + dimensionedScalar("small", dimVelocity, 1e-3)
-        )
-    );
-
-    fvScalarMatrix alpha1Eqn
-    (
-        fvm::ddt(alpha1)
-      + fvm::div(phi, alpha1, alphaScheme)
-      + fvm::div
-        (
-            -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
-            alpha1,
-            alpharScheme
-        )
-    );
-
-    alpha1Eqn.relax();
-    alpha1Eqn.solve();
-
-    Info<< "Liquid phase volume fraction = "
-        << alpha1.weightedAverage(mesh.V()).value()
-        << "  Min(alpha1) = " << min(cellIbMask*alpha1).value()
-        << "  Max(alpha1) = " << max(cellIbMask*alpha1).value()
-        << endl;
-
-    rhoPhi = faceIbMask*(alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2);
-
-    // Limit alpha to handle IB cells
-//     alpha1.max(0);
-//     alpha1.min(1);
-
-    // Update of interface and density
-    interface.correct();
-    twoPhaseProperties.correct();
-
-    // Calculate density using limited alpha1
-    rho = twoPhaseProperties.rho();
-    rho.correctBoundaryConditions();
-}

applications/solvers/immersedBoundary/interIbFoam/correctPhi.H

@@ -1,58 +0,0 @@
-{
-#   include "continuityErrs.H"
-
-    wordList pcorrTypes
-    (
-        pd.boundaryField().size(),
-        zeroGradientFvPatchScalarField::typeName
-    );
-
-    for (label i=0; i<pd.boundaryField().size(); i++)
-    {
-        if (pd.boundaryField()[i].fixesValue())
-        {
-            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
-        }
-    }
-
-    volScalarField pcorr
-    (
-        IOobject
-        (
-            "pcorr",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
-        pcorrTypes
-    );
-
-    dimensionedScalar rUAf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
-
-    // Adjust immersed boundary fluxes
-    immersedBoundaryAdjustPhi(phi, U);
-    adjustPhi(phi, U, pcorr);
-
-    mesh.schemesDict().setFluxRequired(pcorr.name());
-
-    while(pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pcorrEqn
-        (
-            fvm::laplacian(rUAf, pcorr) == fvc::div(phi)
-        );
-
-        pcorrEqn.setReference(pdRefCell, pdRefValue);
-        pcorrEqn.solve();
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi -= pcorrEqn.flux();
-        }
-    }
-
-#   include "immersedBoundaryContinuityErrs.H"
-}

applications/solvers/immersedBoundary/interIbFoam/createFields.H

@@ -1,118 +0,0 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
-    (
-        IOobject
-        (
-            "pd",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading field alpha1\n" << endl;
-    volScalarField alpha1
-    (
-        IOobject
-        (
-            "alpha1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-#   include "createPhi.H"
-
-
-    Info<< "Reading transportProperties\n" << endl;
-    twoPhaseMixture twoPhaseProperties(U, phi, "alpha1");
-
-    const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
-    const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
-
-
-    // Need to store rho for ddt(rho, U)
-    volScalarField rho
-    (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT
-        ),
-        alpha1*rho1 + (scalar(1) - alpha1)*rho2,
-        alpha1.boundaryField().types()
-    );
-    rho.oldTime();
-
-
-    // Mass flux
-    // Initialisation does not matter because rhoPhi is reset after the
-    // alpha1 solution before it is used in the U equation.
-    surfaceScalarField rhoPhi
-    (
-        IOobject
-        (
-            "rho*phi",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        rho1*phi
-    );
-
-
-    Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("gh", g & mesh.Cf());
-
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pd + rho*gh
-    );
-
-
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
-    setRefCell(pd, pimple.dict(), pdRefCell, pdRefValue);
-
-    mesh.schemesDict().setFluxRequired(pd.name());
-    mesh.schemesDict().setFluxRequired(alpha1.name());
-
-    // Construct interface from alpha1 distribution
-    interfaceProperties interface(alpha1, U, twoPhaseProperties);
-
-     // Construct incompressible turbulence model
-    autoPtr<incompressible::turbulenceModel> turbulence
-    (
-        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
-    );

applications/solvers/immersedBoundary/interIbFoam/limitU.H

@@ -1,24 +0,0 @@
-{
-    scalar limitMagU = readScalar(pimple.dict().lookup("limitMagU"));
-
-    volScalarField magU(mag(U));
-
-    scalar maxMagU = max(magU).value();
-
-    Info<< "mag(U): max: " << maxMagU
-        << " avg: " << magU.weightedAverage(mesh.V()).value();
-
-    if (maxMagU > limitMagU)
-    {
-        U.internalField() *=
-            neg(magU.internalField() - limitMagU)
-          + pos(magU.internalField() - limitMagU)*
-            limitMagU/(magU.internalField() + SMALL);
-        U.correctBoundaryConditions();
-        Info << " ...limiting" << endl;
-    }
-    else
-    {
-        Info << endl;
-    }
-}

applications/solvers/immersedBoundary/interIbFoam/pEqn.H

@@ -1,53 +0,0 @@
-{
-    volScalarField rUA = 1.0/UEqn.A();
-    surfaceScalarField rUAf = fvc::interpolate(rUA);
-
-    U = rUA*UEqn.H();
-    // Immersed boundary update
-    U.correctBoundaryConditions();
-#   include "limitU.H"
-
-    surfaceScalarField phiU
-    (
-        "phiU",
-        faceIbMask*(fvc::interpolate(U) & mesh.Sf())
-    );
-
-    phi = phiU
-      + faceIbMask*
-        (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-          - ghf*fvc::snGrad(rho)
-        )*rUAf*mesh.magSf();
-
-    // Adjust immersed boundary fluxes
-    immersedBoundaryAdjustPhi(phi, U);
-    adjustPhi(phi, U, pd);
-
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pdEqn
-        (
-            fvm::laplacian(rUAf, pd, "laplacian(rAU,pd)") == fvc::div(phi)
-        );
-
-        pdEqn.setReference(pdRefCell, pdRefValue);
-
-        pdEqn.solve
-        (
-            mesh.solutionDict().solver(pd.select(pimple.finalInnerIter()))
-        );
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi -= pdEqn.flux();
-        }
-    }
-
-    // Explicitly relax pressure
-    pd.relax();
-
-    U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
-    U.correctBoundaryConditions();
-}

applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/files

@@ -0,0 +1,3 @@
+pimpleDyMIbFoam.C
+
+EXE = $(FOAM_APPBIN)/pimpleDyMIbFoam

applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/options

@@ -1,22 +1,23 @@
 EXE_INC = \
-    -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel
 
 EXE_LIBS = \
-    -linterfaceProperties \
+    -lfiniteVolume \
+    -limmersedBoundary \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -ldynamicMesh \
+    -lmeshTools \
     -lincompressibleTransportModels \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
     -lincompressibleLESModels \
-    -lfiniteVolume \
+    -llduSolvers \
-    -lmeshTools \
+    -L$(MESQUITE_LIB_DIR) -lmesquite
-    -lsurfMesh \
-    -lsampling \
-    -ldynamicMesh \
-    -limmersedBoundary \
-    -llduSolvers

applications/solvers/immersedBoundary/pimpleDyMIbFoam/UEqn.H

@@ -0,0 +1,26 @@
+    // Convection-diffusion matrix
+    fvVectorMatrix HUEqn
+    (
+        fvm::div(phi, U)
+      + turbulence->divDevReff()
+    );
+
+    // Time derivative matrix
+    fvVectorMatrix ddtUEqn(fvm::ddt(U));
+
+    // Get under-relaxation factor
+    scalar UUrf =
+        mesh.solutionDict().equationRelaxationFactor(U.select(pimple.finalIter()));
+
+    if (pimple.momentumPredictor())
+    {
+        // Solve momentum predictor
+        solve
+        (
+            ddtUEqn
+          + relax(HUEqn, UUrf)
+         ==
+          - fvc::grad(p),
+            mesh.solutionDict().solver((U.select(pimple.finalIter())))
+        );
+    }

applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctPhi.H

@@ -0,0 +1,31 @@
+{
+    volScalarField pcorr("pcorr", p);
+
+    // Initialise flux with interpolated velocity
+    phi = fvc::interpolate(U) & mesh.Sf();
+
+    adjustPhi(phi, U, pcorr);
+
+    mesh.schemesDict().setFluxRequired(pcorr.name());
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pcorrEqn
+        (
+            fvm::laplacian(1/aU, pcorr) == fvc::div(phi)
+        );
+
+        pcorrEqn.setReference(pRefCell, pRefValue);
+        pcorrEqn.solve();
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi -= pcorrEqn.flux();
+        }
+
+        // Fluxes are corrected to absolute velocity and further corrected
+        // later.  HJ, 6/Feb/2009
+    }
+#   include "continuityErrs.H"
+}
+

applications/solvers/immersedBoundary/pimpleDyMIbFoam/createControls.H

@@ -0,0 +1,11 @@
+#include "createTimeControls.H"
+
+bool correctPhi
+(
+    pimple.dict().lookupOrDefault("correctPhi", false)
+);
+
+bool checkMeshCourantNo
+(
+    pimple.dict().lookupOrDefault("checkMeshCourantNo", false)
+);

applications/solvers/immersedBoundary/pimpleDyMIbFoam/createFields.H

@@ -1,23 +1,3 @@
-    Info<< "Reading transportProperties\n" << endl;
-
-    IOdictionary transportProperties
-    (
-        IOobject
-        (
-            "transportProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    );
-
-    dimensionedScalar nu
-    (
-        transportProperties.lookup("nu")
-    );
-
-
     Info<< "Reading field p\n" << endl;
     volScalarField p
     (
@@ -54,3 +34,26 @@
     scalar pRefValue = 0.0;
     setRefCell(p, pimple.dict(), pRefCell, pRefValue);
     mesh.schemesDict().setFluxRequired(p.name());
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );
+
+    Info<< "Reading field aU if present\n" << endl;
+    volScalarField aU
+    (
+        IOobject
+        (
+            "aU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        1/runTime.deltaT(),
+        zeroGradientFvPatchScalarField::typeName
+    );

applications/solvers/immersedBoundary/pimpleDyMIbFoam/pEqn.H

@@ -0,0 +1,58 @@
+{
+    p.boundaryField().updateCoeffs();
+
+    // Prepare clean Ap without time derivative contribution and
+    // without contribution from under-relaxation
+    // HJ, 26/Oct/2015
+    aU = HUEqn.A();
+
+    // Store velocity under-relaxation point before using U for the flux
+    // precursor
+    U.storePrevIter();
+
+    U = HUEqn.H()/aU;
+    phi = (fvc::interpolate(U) & mesh.Sf());
+
+    // Global flux balance
+    adjustPhi(phi, U, p);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::laplacian(1/aU, p) == fvc::div(phi)
+        );
+
+        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.solve
+        (
+            mesh.solutionDict().solver(p.select(pimple.finalInnerIter()))
+        );
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi -= pEqn.flux();
+        }
+    }
+
+    // Explicitly relax pressure for momentum corrector except for last corrector
+    if (!pimple.finalIter())
+    {
+        p.relax();
+    }
+
+    // Make the fluxes relative to the mesh motion
+    fvc::makeRelative(phi, U);
+
+#   include "movingMeshContinuityErrs.H"
+
+    U = UUrf*
+        (
+            1.0/(aU + ddtUEqn.A())*
+            (
+                U*aU - fvc::grad(p) + ddtUEqn.H()
+            )
+        )
+      + (1 - UUrf)*U.prevIter();
+    U.correctBoundaryConditions();
+}

applications/solvers/immersedBoundary/pimpleDyMIbFoam/pimpleDyMIbFoam.C

@@ -22,19 +22,17 @@ License
     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    interIbFoam
+    pimpleDyMFoam.C
 
 Description
-    Solver for 2 incompressible, isothermal immiscible fluids using a VOF
+    Transient solver for incompressible, flow of Newtonian fluids
-    (volume of fluid) phase-fraction based interface capturing approach,
+    with dynamic mesh using the PIMPLE (merged PISO-SIMPLE) algorithm.
-    with immersed boundary support
-
-    The momentum and other fluid properties are of the "mixture" and a single
-    momentum equation is solved.
 
     Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
 
-    For a two-fluid approach see twoPhaseEulerFoam.
+    Consistent formulation without time-step and relaxation dependence by Jasak
+
+    Support for immersed boundary
 
 Author
     Hrvoje Jasak, Wikki Ltd.  All rights reserved
@@ -42,32 +40,30 @@ Author
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "interfaceProperties.H"
+#include "singlePhaseTransportModel.H"
-#include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
+#include "dynamicFvMesh.H"
 #include "pimpleControl.H"
 
+#include "immersedBoundaryPolyPatch.H"
 #include "immersedBoundaryFvPatch.H"
-#include "immersedBoundaryAdjustPhi.H"
+#include "emptyFvPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
 #   include "setRootCase.H"
+
 #   include "createTime.H"
-#   include "createMesh.H"
+#   include "createDynamicFvMesh.H"
 
     pimpleControl pimple(mesh);
 
-#   include "readGravitationalAcceleration.H"
 #   include "initContinuityErrs.H"
-#   include "createFields.H"
 #   include "createIbMasks.H"
-#   include "createTimeControls.H"
+#   include "createFields.H"
-#   include "correctPhi.H"
+#   include "createControls.H"
-#   include "CourantNo.H"
-#   include "setInitialDeltaT.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -75,19 +71,46 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readTimeControls.H"
+#       include "readControls.H"
-#       include "immersedBoundaryCourantNo.H"
+#       include "CourantNo.H"
 #       include "setDeltaT.H"
 
         runTime++;
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        // Pressure-velocity corrector
+        bool meshChanged = mesh.update();
+
+        U.correctBoundaryConditions();
+        p.correctBoundaryConditions();
+        aU.correctBoundaryConditions();
+        phi.correctBoundaryConditions();
+        turbulence->correct();
+
+#       include "updateIbMasks.H"
+#       include "volContinuity.H"
+
+        if (runTime.outputTime())
+        {
+            volScalarField divMeshPhi("divMeshPhi", mag(fvc::surfaceIntegrate(mesh.phi())));
+            divMeshPhi.write();
+        }
+
+        if (checkMeshCourantNo)
+        {
+#           include "meshCourantNo.H"
+        }
+
+        // Fluxes will be corrected to absolute velocity
+        // HJ, 6/Feb/2009
+#       include "correctPhi.H"
+
+        // Make the fluxes relative to the mesh motion
+        fvc::makeRelative(phi, U);
+
+        // --- PIMPLE loop
         while (pimple.loop())
         {
-            twoPhaseProperties.correct();
-
 #           include "UEqn.H"
 
             // --- PISO loop
@@ -96,15 +119,6 @@ int main(int argc, char *argv[])
 #               include "pEqn.H"
             }
 
-#           include "immersedBoundaryContinuityErrs.H"
-
-#           include "limitU.H"
-
-            p = pd + rho*gh;
-            p.correctBoundaryConditions();
-
-#           include "alphaEqn.H"
-
             turbulence->correct();
         }
 

applications/solvers/immersedBoundary/pimpleDyMIbFoam/readControls.H

@@ -0,0 +1,5 @@
+#include "readTimeControls.H"
+
+correctPhi = pimple.dict().lookupOrDefault("correctPhi", false);
+
+checkMeshCourantNo = pimple.dict().lookupOrDefault("checkMeshCourantNo", false);

applications/solvers/immersedBoundary/porousSimpleIbFoam/Make/files

@@ -1,3 +0,0 @@
-porousSimpleIbFoam.C
-
-EXE = $(FOAM_APPBIN)/porousSimpleIbFoam

applications/solvers/immersedBoundary/porousSimpleIbFoam/Make/options

@@ -1,21 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
-    -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude
-
-EXE_LIBS = \
-    -lincompressibleRASModels \
-    -lincompressibleTransportModels \
-    -lincompressibleTurbulenceModel \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lsurfMesh \
-    -lsampling \
-    -ldynamicMesh \
-    -limmersedBoundary \
-    -limmersedBoundaryTurbulence \
-    -llduSolvers

applications/solvers/immersedBoundary/porousSimpleIbFoam/UEqn.H

@@ -1,12 +0,0 @@
-    // Solve the Momentum equation
-    tmp<fvVectorMatrix> UEqn
-    (
-        fvm::div(phi, U)
-      + turbulence->divDevReff()
-    );
-
-    UEqn().relax();
-
-    pZones.addResistance(UEqn());
-
-    solve(UEqn() == -fvc::grad(p));

applications/solvers/immersedBoundary/porousSimpleIbFoam/createFields.H

@@ -1,44 +0,0 @@
-    Info << "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info << "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-#   include "createPhi.H"
-
-
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, simple.dict(), pRefCell, pRefValue);
-    mesh.schemesDict().setFluxRequired(p.name());
-
-    singlePhaseTransportModel laminarTransport(U, phi);
-
-    autoPtr<incompressible::RASModel> turbulence
-    (
-        incompressible::RASModel::New(U, phi, laminarTransport)
-    );
-
-    porousZones pZones(mesh);

applications/solvers/immersedBoundary/porousSimpleIbFoam/pEqn.H

@@ -1,41 +0,0 @@
-{
-    volScalarField AU = UEqn().A();
-
-    U = UEqn().H()/AU;
-    // Immersed boundary update
-    U.correctBoundaryConditions();
-
-    UEqn.clear();
-    phi = faceIbMask*(fvc::interpolate(U) & mesh.Sf());
-
-    // Adjust immersed boundary fluxes
-    immersedBoundaryAdjustPhi(phi, U);
-    adjustPhi(phi, U, p);
-
-    // Non-orthogonal pressure corrector loop
-    while (simple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::laplacian(1.0/AU, p) == fvc::div(phi)
-        );
-
-        pEqn.setReference(pRefCell, pRefValue);
-
-        pEqn.solve();
-
-        if (simple.finalNonOrthogonalIter())
-        {
-            phi -= pEqn.flux();
-        }
-    }
-
-#   include "immersedBoundaryContinuityErrs.H"
-
-    // Explicitly relax pressure for momentum corrector
-    p.relax();
-
-    // Momentum corrector
-    U -= fvc::grad(p)/AU;
-    U.correctBoundaryConditions();
-}

applications/solvers/immersedBoundary/porousSimpleIbFoam/porousSimpleIbFoam.C

@@ -1,88 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    simpleIbFoam
-
-Description
-    Steady-state solver for incompressible, turbulent flow
-    with porous zones and immersed boundary support.
-
-Author
-    Hrvoje Jasak, Wikki Ltd.  All rights reserved
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "singlePhaseTransportModel.H"
-#include "RASModel.H"
-#include "porousZones.H"
-#include "simpleControl.H"
-
-#include "immersedBoundaryFvPatch.H"
-#include "immersedBoundaryAdjustPhi.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-
-    simpleControl simple(mesh);
-
-#   include "createIbMasks.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (simple.loop())
-    {
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        // Pressure-velocity SIMPLE corrector
-        {
-#           include "UEqn.H"
-#           include "pEqn.H"
-        }
-
-        turbulence->correct();
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/immersedBoundary/potentialIbFoam/Make/files

@@ -1,3 +0,0 @@
-potentialIbFoam.C
-
-EXE = $(FOAM_APPBIN)/potentialIbFoam

applications/solvers/immersedBoundary/potentialIbFoam/Make/options

@@ -1,16 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude
-
-EXE_LIBS = \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lsurfMesh \
-    -lsampling \
-    -ldynamicMesh \
-    -ldynamicFvMesh \
-    -limmersedBoundary \
-    -llduSolvers

applications/solvers/immersedBoundary/potentialIbFoam/createFields.H

@@ -1,50 +0,0 @@
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    // Do not reset p
-
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    // Do not reset U
-
-    surfaceScalarField phi
-    (
-        IOobject
-        (
-            "phi",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        fvc::interpolate(U) & mesh.Sf()
-    );
-
-
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, simple.dict(), pRefCell, pRefValue);
-    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/immersedBoundary/potentialIbFoam/potentialIbFoam.C

@@ -1,231 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    potentialIbFoam
-
-Description
-    Potential flow solver with immersed boundary support.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "immersedBoundaryFvPatch.H"
-#include "immersedBoundaryAdjustPhi.H"
-#include "simpleControl.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    argList::validOptions.insert("writep", "");
-
-#   include "setRootCase.H"
-
-#   include "createTime.H"
-#   include "createMesh.H"
-
-    simpleControl simple(mesh);
-
-#   include "createIbMasks.H"
-#   include "createFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< nl << "Calculating potential flow" << endl;
-
-    // Do correctors over the complete set
-    while (simple.correctNonOrthogonal())
-    {
-        phi = faceIbMask*(linearInterpolate(U) & mesh.Sf());
-
-        // Adjust immersed boundary fluxes
-        immersedBoundaryAdjustPhi(phi, U);
-
-        // Adjust fluxes
-        adjustPhi(phi, U, p);
-
-        p.storePrevIter();
-
-        fvScalarMatrix pEqn
-        (
-            fvm::laplacian
-            (
-                dimensionedScalar
-                (
-                    "1",
-                    dimTime/p.dimensions()*dimensionSet(0, 2, -2, 0, 0),
-                    1
-                ),
-                p
-            )
-         ==
-            fvc::div(phi)
-        );
-
-        pEqn.setReference(pRefCell, pRefValue);
-        pEqn.solve();
-
-        // Correct the flux
-        phi -= pEqn.flux();
-
-        if (!simple.finalNonOrthogonalIter())
-        {
-            p.relax();
-        }
-
-        Info<< "p min " << gMin(p.internalField())
-            << " max " << gMax(p.internalField())
-            << " masked min "
-            << gMin(cellIbMask.internalField()*p.internalField())
-            << " max "
-            << gMax(cellIbMask.internalField()*p.internalField())
-            << endl;
-
-        Info<< "continuity error = "
-            << mag
-               (
-                    fvc::div(faceIbMask*phi)
-               )().weightedAverage(mesh.V()).value()
-            << endl;
-
-        Info<< "Contour continuity error = "
-            << mag(sum(phi.boundaryField()))
-            << endl;
-
-        U = fvc::reconstruct(phi);
-        U.correctBoundaryConditions();
-
-        Info<< "Interpolated U error = "
-            << (
-                   sqrt
-                   (
-                       sum
-                       (
-                           sqr
-                           (
-                               faceIbMask*
-                               (
-                                   fvc::interpolate(U) & mesh.Sf()
-                               )
-                             - phi
-                           )
-                       )
-                   )/sum(mesh.magSf())
-               ).value()
-            << endl;
-    }
-
-    // Calculate velocity magnitude
-    {
-        volScalarField magU = cellIbMask*mag(U);
-
-        Info << "IB-masked mag(U): max: " << gMax(magU.internalField())
-            << " min: " << gMin(magU.internalField()) << endl;
-    }
-
-    // Force the write
-    U.write();
-    phi.write();
-
-    cellIbMask.write();
-    cellIbMaskExt.write();
-
-    if (args.optionFound("writep"))
-    {
-        // Find reference patch
-        label refPatch = -1;
-        scalar maxMagU = 0;
-
-        // Go through all velocity patches and find the one that fixes
-        // velocity to the largest value
-
-        forAll (U.boundaryField(), patchI)
-        {
-            const fvPatchVectorField& Upatch = U.boundaryField()[patchI];
-
-            if (Upatch.fixesValue())
-            {
-                // Calculate mean velocity
-                scalar u = sum(mag(Upatch));
-                label patchSize = Upatch.size();
-
-                reduce(u, sumOp<scalar>());
-                reduce(patchSize, sumOp<label>());
-
-                if (patchSize > 0)
-                {
-                    scalar curMag = u/patchSize;
-
-                    if (curMag > maxMagU)
-                    {
-                        refPatch = patchI;
-
-                        maxMagU = curMag;
-                    }
-                }
-            }
-        }
-
-        if (refPatch > -1)
-        {
-            // Calculate reference pressure
-            const fvPatchVectorField& Upatch = U.boundaryField()[refPatch];
-            const fvPatchScalarField& pPatch = p.boundaryField()[refPatch];
-
-            scalar patchE = sum(mag(pPatch + 0.5*magSqr(Upatch)));
-            label patchSize = Upatch.size();
-
-            reduce(patchE, sumOp<scalar>());
-            reduce(patchSize, sumOp<label>());
-
-            scalar e = patchE/patchSize;
-
-
-            Info<< "Using reference patch " << refPatch
-                << " with mag(U) = " << maxMagU
-                << " p + 0.5*U^2 = " << e << endl;
-
-            p.internalField() = e - 0.5*magSqr(U.internalField());
-            p.correctBoundaryConditions();
-        }
-        else
-        {
-            Info<< "No reference patch found.  Writing potential function"
-                << endl;
-        }
-
-        p.write();
-    }
-
-    Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-        << nl << endl;
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/immersedBoundary/simpleIbFoam/Make/files

@@ -1,3 +0,0 @@
-simpleIbFoam.C
-
-EXE = $(FOAM_APPBIN)/simpleIbFoam

applications/solvers/immersedBoundary/simpleIbFoam/Make/options

@@ -1,21 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
-    -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude
-
-EXE_LIBS = \
-    -lincompressibleRASModels \
-    -lincompressibleTransportModels \
-    -lincompressibleTurbulenceModel \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lsurfMesh \
-    -lsampling \
-    -ldynamicMesh \
-    -limmersedBoundary \
-    -limmersedBoundaryTurbulence \
-    -llduSolvers

applications/solvers/immersedBoundary/simpleIbFoam/UEqn.H

@@ -1,10 +0,0 @@
-    // Solve the Momentum equation
-    tmp<fvVectorMatrix> UEqn
-    (
-        fvm::div(phi, U)
-      + turbulence->divDevReff()
-    );
-
-    UEqn().relax();
-
-    solve(UEqn() == -fvc::grad(p));

applications/solvers/immersedBoundary/simpleIbFoam/createFields.H

@@ -1,42 +0,0 @@
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-#   include "createPhi.H"
-
-
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, simple.dict(), pRefCell, pRefValue);
-    mesh.schemesDict().setFluxRequired(p.name());
-
-    singlePhaseTransportModel laminarTransport(U, phi);
-
-    autoPtr<incompressible::RASModel> turbulence
-    (
-        incompressible::RASModel::New(U, phi, laminarTransport)
-    );

applications/solvers/immersedBoundary/simpleIbFoam/pEqn.H

@@ -1,41 +0,0 @@
-{
-    volScalarField AU = UEqn().A();
-
-    U = UEqn().H()/AU;
-    // Immersed boundary update
-    U.correctBoundaryConditions();
-
-    UEqn.clear();
-    phi = faceIbMask*(fvc::interpolate(U) & mesh.Sf());
-
-    // Adjust immersed boundary fluxes
-    immersedBoundaryAdjustPhi(phi, U);
-    adjustPhi(phi, U, p);
-
-    // Non-orthogonal pressure corrector loop
-    while (simple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::laplacian(1.0/AU, p) == fvc::div(phi)
-        );
-
-        pEqn.setReference(pRefCell, pRefValue);
-
-        pEqn.solve();
-
-        if (simple.finalNonOrthogonalIter())
-        {
-            phi -= pEqn.flux();
-        }
-    }
-
-#   include "immersedBoundaryContinuityErrs.H"
-
-    // Explicitly relax pressure for momentum corrector
-    p.relax();
-
-    // Momentum corrector
-    U -= fvc::grad(p)/AU;
-    U.correctBoundaryConditions();
-}

applications/solvers/immersedBoundary/simpleIbFoam/simpleIbFoam.C

@@ -1,87 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    simpleIbFoam
-
-Description
-    Steady-state solver for incompressible, turbulent flow
-    with immersed boundary support.
-
-Author
-    Hrvoje Jasak, Wikki Ltd.  All rights reserved
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "singlePhaseTransportModel.H"
-#include "RASModel.H"
-#include "simpleControl.H"
-
-#include "immersedBoundaryFvPatch.H"
-#include "immersedBoundaryAdjustPhi.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-
-    simpleControl simple(mesh);
-
-#   include "createIbMasks.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (simple.loop())
-    {
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        // Pressure-velocity SIMPLE corrector
-        {
-#           include "UEqn.H"
-#           include "pEqn.H"
-        }
-
-        turbulence->correct();
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/incompressible/MRFSimpleFoam/createFields.H

@@ -12,7 +12,7 @@
         mesh
     );
 
-    Info << "Reading field U\n" << endl;
+    Info<< "Reading field U\n" << endl;
     volVectorField U
     (
         IOobject

applications/solvers/incompressible/pimpleDyMFoam/correctPhi.H

@@ -1,34 +1,6 @@
 {
-#   include "continuityErrs.H"
+    volScalarField pcorr("pcorr", p);
-
+    pcorr *= 0;
-    wordList pcorrTypes
-    (
-        p.boundaryField().size(),
-        zeroGradientFvPatchScalarField::typeName
-    );
-
-    for (label i=0; i<p.boundaryField().size(); i++)
-    {
-        if (p.boundaryField()[i].fixesValue())
-        {
-            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
-        }
-    }
-
-    volScalarField pcorr
-    (
-        IOobject
-        (
-            "pcorr",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("pcorr", p.dimensions(), 0.0),
-        pcorrTypes
-    );
 
     // Initialise flux with interpolated velocity
     phi = fvc::interpolate(U) & mesh.Sf();

applications/solvers/multiphase/interPhaseChangeFoam/createFields.H

@@ -104,7 +104,6 @@
     label pdRefCell = 0;
     scalar pdRefValue = 0.0;
     setRefCell(pd, pimple.dict(), pdRefCell, pdRefValue);
-    mesh.schemesDict().setFluxRequired(pd.name());
 
     scalar pRefValue = 0.0;
 
@@ -128,3 +127,6 @@
     (
         incompressible::turbulenceModel::New(U, phi, twoPhaseProperties())
     );
+
+    mesh.schemesDict().setFluxRequired(pd.name());
+    mesh.schemesDict().setFluxRequired(alpha1.name());

applications/solvers/solidMechanics/elasticThermalSolidFoam/Make/options

@@ -5,7 +5,8 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/topoChangerFvMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/solidModels/lnInclude
+    -I$(LIB_SRC)/solidModels/lnInclude \
+    -I$(LIB_SRC)/conjugateHeatTransfer/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
@@ -13,4 +14,5 @@ EXE_LIBS = \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -ltopoChangerFvMesh \
-    -lsolidModels
+    -lsolidModels \
+    -lconjugateHeatTransfer

applications/solvers/surfaceTracking/freeSurface/freeSurface.C

@@ -187,13 +187,13 @@ freeSurface::freeSurface
     // Set point normal correction patches
     boolList& correction = aMesh().correctPatchPointNormals();
 
-    forAll(pointNormalsCorrectionPatches_, patchI)
+    forAll (pointNormalsCorrectionPatches_, patchI)
     {
         word patchName = pointNormalsCorrectionPatches_[patchI];
 
         label patchID = aMesh().boundary().findPatchID(patchName);
 
-        if(patchID == -1)
+        if (patchID == -1)
         {
             FatalErrorIn
             (
@@ -212,7 +212,7 @@ freeSurface::freeSurface
     // Detect the free surface patch
     forAll (mesh().boundary(), patchI)
     {
-        if(mesh().boundary()[patchI].name() == "freeSurface")
+        if (mesh().boundary()[patchI].name() == "freeSurface")
         {
             aPatchID_ = patchI;
 
@@ -221,7 +221,7 @@ freeSurface::freeSurface
         }
     }
 
-    if(aPatchID() == -1)
+    if (aPatchID() == -1)
     {
         FatalErrorIn("freeSurface::freeSurface(...)")
             << "Free surface patch not defined.  Please make sure that "
@@ -235,7 +235,7 @@ freeSurface::freeSurface
     {
         forAll (mesh().boundary(), patchI)
         {
-            if(mesh().boundary()[patchI].name() == "freeSurfaceShadow")
+            if (mesh().boundary()[patchI].name() == "freeSurfaceShadow")
             {
                 bPatchID_ = patchI;
 
@@ -244,7 +244,7 @@ freeSurface::freeSurface
             }
         }
 
-        if(bPatchID() == -1)
+        if (bPatchID() == -1)
         {
             FatalErrorIn("freeSurface::freeSurface(...)")
                 << "Free surface shadow patch not defined. "
@@ -257,7 +257,7 @@ freeSurface::freeSurface
 
     // Mark free surface boundary points
     // which belonge to processor patches
-    forAll(aMesh().boundary(), patchI)
+    forAll (aMesh().boundary(), patchI)
     {
         if
         (
@@ -268,7 +268,7 @@ freeSurface::freeSurface
             const labelList& patchPoints =
                 aMesh().boundary()[patchI].pointLabels();
 
-            forAll(patchPoints, pointI)
+            forAll (patchPoints, pointI)
             {
                 motionPointsMask()[patchPoints[pointI]] = -1;
             }
@@ -277,7 +277,7 @@ freeSurface::freeSurface
 
 
     // Mark fixed free surface boundary points
-    forAll(fixedFreeSurfacePatches_, patchI)
+    forAll (fixedFreeSurfacePatches_, patchI)
     {
         label fixedPatchID =
             aMesh().boundary().findPatchID
@@ -285,7 +285,7 @@ freeSurface::freeSurface
                 fixedFreeSurfacePatches_[patchI]
             );
 
-        if(fixedPatchID == -1)
+        if (fixedPatchID == -1)
         {
             FatalErrorIn("freeSurface::freeSurface(...)")
                 << "Wrong faPatch name in the fixedFreeSurfacePatches list"
@@ -296,7 +296,7 @@ freeSurface::freeSurface
         const labelList& patchPoints =
             aMesh().boundary()[fixedPatchID].pointLabels();
 
-        forAll(patchPoints, pointI)
+        forAll (patchPoints, pointI)
         {
             motionPointsMask()[patchPoints[pointI]] = 0;
         }
@@ -305,7 +305,7 @@ freeSurface::freeSurface
 
     // Mark free-surface boundary point
     // at the axis of 2-D axisymmetic cases
-    forAll(aMesh().boundary(), patchI)
+    forAll (aMesh().boundary(), patchI)
     {
         if
         (
@@ -316,7 +316,7 @@ freeSurface::freeSurface
             const wedgeFaPatch& wedgePatch =
                 refCast<const wedgeFaPatch>(aMesh().boundary()[patchI]);
 
-            if(wedgePatch.axisPoint() > -1)
+            if (wedgePatch.axisPoint() > -1)
             {
                 motionPointsMask()[wedgePatch.axisPoint()] = 0;
 
@@ -357,7 +357,7 @@ freeSurface::~freeSurface()
 
 void freeSurface::updateDisplacementDirections()
 {
-    if(normalMotionDir())
+    if (normalMotionDir())
     {
         // Update point displacement correction
         pointsDisplacementDir() = aMesh().pointAreaNormals();
@@ -365,9 +365,9 @@ void freeSurface::updateDisplacementDirections()
         // Correcte point displacement direction
         // at the "centerline" symmetryPlane which represents the axis
         // of an axisymmetric case
-        forAll(aMesh().boundary(), patchI)
+        forAll (aMesh().boundary(), patchI)
         {
-            if(aMesh().boundary()[patchI].type() == wedgeFaPatch::typeName)
+            if (aMesh().boundary()[patchI].type() == wedgeFaPatch::typeName)
             {
                 const wedgeFaPatch& wedgePatch =
                     refCast<const wedgeFaPatch>(aMesh().boundary()[patchI]);
@@ -377,12 +377,12 @@ void freeSurface::updateDisplacementDirections()
                 label centerLinePatchID =
                     aMesh().boundary().findPatchID("centerline");
 
-                if(centerLinePatchID != -1)
+                if (centerLinePatchID != -1)
                 {
                     const labelList& pointLabels =
                         aMesh().boundary()[centerLinePatchID].pointLabels();
 
-                    forAll(pointLabels, pointI)
+                    forAll (pointLabels, pointI)
                     {
                         vector dir =
                             pointsDisplacementDir()[pointLabels[pointI]];
@@ -449,8 +449,8 @@ bool freeSurface::correctPoints()
 {
     for
     (
-        int freeSurfCorr=0;
+        int freeSurfCorr = 0;
-        freeSurfCorr<nFreeSurfCorr_;
+        freeSurfCorr < nFreeSurfCorr_;
         freeSurfCorr++
     )
     {
@@ -463,7 +463,7 @@ bool freeSurface::correctPoints()
 
 bool freeSurface::movePoints(const scalarField& interfacePhi)
 {
-    pointField newMeshPoints = mesh().points();
+    pointField newMeshPoints = mesh().allPoints();
 
     scalarField sweptVolCorr =
         interfacePhi
@@ -535,7 +535,7 @@ bool freeSurface::movePoints(const scalarField& interfacePhi)
         newMeshPoints[meshPointsA[pointI]] += displacement[pointI];
     }
 
-    if(twoFluids_)
+    if (twoFluids_)
     {
         const labelList& meshPointsB =
             mesh().boundaryMesh()[bPatchID_].meshPoints();
@@ -554,12 +554,12 @@ bool freeSurface::movePoints(const scalarField& interfacePhi)
 
     // Update total displacement field
 
-    if(totalDisplacementPtr_ && (curTimeIndex_ < DB().timeIndex()))
+    if (totalDisplacementPtr_ && (curTimeIndex_ < DB().timeIndex()))
     {
         FatalErrorIn("freeSurface::movePoints()")
             << "Total displacement of free surface points "
-                << "from previous time step is not absorbed by the mesh."
+            << "from previous time step is not absorbed by the mesh."
-                << abort(FatalError);
+            << abort(FatalError);
     }
     else if (curTimeIndex_ < DB().timeIndex())
     {
@@ -585,7 +585,7 @@ bool freeSurface::movePoints(const scalarField& interfacePhi)
 
     // Move correctedFvPatchField fvSubMeshes
 
-    forAll(U().boundaryField(), patchI)
+    forAll (U().boundaryField(), patchI)
     {
         if
         (
@@ -615,7 +615,7 @@ bool freeSurface::movePoints(const scalarField& interfacePhi)
         }
     }
 
-    forAll(p().boundaryField(), patchI)
+    forAll (p().boundaryField(), patchI)
     {
         if
         (
@@ -651,7 +651,7 @@ bool freeSurface::movePoints(const scalarField& interfacePhi)
 
 bool freeSurface::moveMeshPointsForOldFreeSurfDisplacement()
 {
-    if(totalDisplacementPtr_)
+    if (totalDisplacementPtr_)
     {
         pointField newPoints = mesh().points();
 
@@ -708,7 +708,7 @@ bool freeSurface::moveMeshPointsForOldFreeSurfDisplacement()
                     totalDisplacement()/DB().deltaT().value()
                 );
 
-            if(twoFluids_)
+            if (twoFluids_)
             {
                 const labelList& meshPointsB =
                     mesh().boundaryMesh()[bPatchID()].meshPoints();
@@ -764,7 +764,7 @@ bool freeSurface::moveMeshPointsForOldFreeSurfDisplacement()
             motionUaPatch ==
                 totalDisplacement()/DB().deltaT().value();
 
-            if(twoFluids_)
+            if (twoFluids_)
             {
                 const labelList& meshPointsB =
                     mesh().boundaryMesh()[bPatchID()].meshPoints();
@@ -808,7 +808,7 @@ bool freeSurface::moveMeshPointsForOldFreeSurfDisplacement()
 
         // Move correctedFvPatchField fvSubMeshes
 
-        forAll(U().boundaryField(), patchI)
+        forAll (U().boundaryField(), patchI)
         {
             if
             (
@@ -838,7 +838,7 @@ bool freeSurface::moveMeshPointsForOldFreeSurfDisplacement()
             }
         }
 
-        forAll(p().boundaryField(), patchI)
+        forAll (p().boundaryField(), patchI)
         {
             if
             (
@@ -994,7 +994,7 @@ bool freeSurface::moveMeshPoints()
 
         // Move correctedFvPatchField fvSubMeshes
 
-        forAll(U().boundaryField(), patchI)
+        forAll (U().boundaryField(), patchI)
         {
             if
             (
@@ -1024,7 +1024,7 @@ bool freeSurface::moveMeshPoints()
             }
         }
 
-        forAll(p().boundaryField(), patchI)
+        forAll (p().boundaryField(), patchI)
         {
             if
             (
@@ -1068,7 +1068,7 @@ void freeSurface::updateBoundaryConditions()
 
 void freeSurface::updateVelocity()
 {
-    if(twoFluids())
+    if (twoFluids())
     {
         vectorField nA = mesh().boundary()[aPatchID()].nf();
 
@@ -1141,7 +1141,7 @@ void freeSurface::updateVelocity()
 
         vectorField tangentialSurfaceTensionForce(nA.size(), vector::zero);
 
-        if(!cleanInterface())
+        if (!cleanInterface())
         {
             tangentialSurfaceTensionForce =
                 surfaceTensionGrad()().internalField();
@@ -1262,7 +1262,7 @@ void freeSurface::updateVelocity()
 
         vectorField tangentialSurfaceTensionForce(nA.size(), vector::zero);
 
-        if(!cleanInterface())
+        if (!cleanInterface())
         {
             tangentialSurfaceTensionForce =
                 surfaceTensionGrad()().internalField();
@@ -1357,7 +1357,7 @@ void freeSurface::updatePressure()
 
     vectorField nA = mesh().boundary()[aPatchID()].nf();
 
-    if(twoFluids())
+    if (twoFluids())
     {
         scalarField pA =
             interpolatorBA().faceInterpolate
@@ -1371,7 +1371,7 @@ void freeSurface::updatePressure()
             << ", max = " << gMax(K)
             << ", average = " << gAverage(K) << endl << flush;
 
-        if(cleanInterface())
+        if (cleanInterface())
         {
 //             pA -= cleanInterfaceSurfTension().value()*(K - gAverage(K));
             pA -= cleanInterfaceSurfTension().value()*K;
@@ -1422,7 +1422,7 @@ void freeSurface::updatePressure()
             << ", max = " << gMax(K) << ", average = " << gAverage(K)
             << endl;
 
-        if(cleanInterface())
+        if (cleanInterface())
         {
 //             pA -= cleanInterfaceSurfTension().value()*(K - gAverage(K));
             pA -= cleanInterfaceSurfTension().value()*K;
@@ -1474,7 +1474,7 @@ void freeSurface::updateSurfaceFlux()
 
 void freeSurface::updateSurfactantConcentration()
 {
-    if(!cleanInterface())
+    if (!cleanInterface())
     {
         Info << "Correct surfactant concentration" << endl << flush;
 
@@ -1493,7 +1493,7 @@ void freeSurface::updateSurfactantConcentration()
         );
 
 
-        if(surfactant().soluble())
+        if (surfactant().soluble())
         {
             const scalarField& C =
                 mesh().boundary()[aPatchID()]
@@ -1541,7 +1541,7 @@ void freeSurface::updateSurfactantConcentration()
            *log(1.0 - surfactantConcentration()
            /surfactant().surfactSaturatedConc());
 
-        if(neg(min(surfaceTension().internalField())))
+        if (neg(min(surfaceTension().internalField())))
         {
             FatalErrorIn
             (
@@ -1555,7 +1555,7 @@ void freeSurface::updateSurfactantConcentration()
 
 void freeSurface::correctUsBoundaryConditions()
 {
-    forAll(Us().boundaryField(), patchI)
+    forAll (Us().boundaryField(), patchI)
     {
         if
         (
@@ -1570,7 +1570,7 @@ void freeSurface::correctUsBoundaryConditions()
             label ngbPolyPatchID =
                 aMesh().boundary()[patchI].ngbPolyPatchIndex();
 
-            if(ngbPolyPatchID != -1)
+            if (ngbPolyPatchID != -1)
             {
                 if
                 (
@@ -1646,7 +1646,7 @@ vector freeSurface::totalSurfaceTensionForce() const
 
     vectorField surfTensionForces(n.size(), vector::zero);
 
-    if(cleanInterface())
+    if (cleanInterface())
     {
         surfTensionForces =
             S*cleanInterfaceSurfTension().value()
@@ -1678,7 +1678,7 @@ void freeSurface::initializeControlPointsPosition()
 
     scalarField sweptVol(faces.size(), 0.0);
 
-    forAll(faces, faceI)
+    forAll (faces, faceI)
     {
         sweptVol[faceI] = -faces[faceI].sweptVol(points, newPoints);
     }
@@ -1690,7 +1690,7 @@ void freeSurface::initializeControlPointsPosition()
         faceArea[faceI] = faces[faceI].normal(newPoints);
     }
 
-    forAll(deltaH, faceI)
+    forAll (deltaH, faceI)
     {
         deltaH[faceI] = sweptVol[faceI]/
             (faceArea[faceI] & facesDisplacementDir()[faceI]);
@@ -1704,9 +1704,9 @@ scalar freeSurface::maxCourantNumber()
 {
     scalar CoNum = 0;
 
-    if(cleanInterface())
+    if (cleanInterface())
     {
-        const scalarField& dE =aMesh().lPN();
+        const scalarField& dE = aMesh().lPN();
 
         CoNum = gMax
         (
@@ -1781,7 +1781,7 @@ void freeSurface::writeVTKControlPoints()
         << "DATASET POLYDATA" << nl
         << "POINTS " << controlPoints().size() << " float" << nl;
 
-    forAll(controlPoints(), pointI)
+    forAll (controlPoints(), pointI)
     {
         mps << controlPoints()[pointI].x() << ' '
             << controlPoints()[pointI].y() << ' '
@@ -1792,7 +1792,7 @@ void freeSurface::writeVTKControlPoints()
     mps << "VERTICES " << controlPoints().size() << ' '
         << controlPoints().size()*2 << nl;
 
-    forAll(controlPoints(), pointI)
+    forAll (controlPoints(), pointI)
     {
         mps << 1 << ' ' << pointI << nl;
     }

applications/solvers/surfaceTracking/freeSurface/freeSurface.H

@@ -135,6 +135,7 @@ class freeSurface
         //- Interface smoothing at the begining of time step
         Switch smoothing_;
 
+
     // Demand-driven data
 
         //- Patch to patch interpolation object which deals with
@@ -229,6 +230,7 @@ class freeSurface
             const vector& axis
         ) const;
 
+    
 public:
 
     // Declare name of the class and it's debug switch
@@ -247,9 +249,8 @@ public:
     );
 
 
-    // Destructor
+    //- Destructor
-
+    virtual ~freeSurface();
-    ~freeSurface();
 
 
     // Member Functions

applications/solvers/surfaceTracking/interTrackFoam/interTrackFoam.C

@@ -41,24 +41,30 @@ Description
 int main(int argc, char *argv[])
 {
 #   include "setRootCase.H"
-
 #   include "createTime.H"
-
 #   include "createDynamicFvMesh.H"
 
     pimpleControl pimple(mesh);
 
 #   include "createFields.H"
-
+#   include "createTimeControls.H"
 #   include "initContinuityErrs.H"
+#   include "CourantNo.H"
+#   include "setInitialDeltaT.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info << "\nStarting time loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.run())
     {
-        Info << "Time = " << runTime.value() << endl << endl;
+#       include "readTimeControls.H"
+#       include "CourantNo.H"
+        //#       include "setSurfaceStabilityDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
 
         interface.moveMeshPointsForOldFreeSurfDisplacement();
 
@@ -116,6 +122,8 @@ int main(int argc, char *argv[])
                     {
                         phi -= pEqn.flux();
                     }
+                    
+                    p.relax();
                 }
 
 #               include "continuityErrs.H"

applications/solvers/surfaceTracking/interTrackFoam/setSurfaceStabilityDeltaT.H

@@ -22,28 +22,34 @@ License
     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
 
 Global
-    continuityErrs
+    setSurfaceStabilityDeltaT
 
 Description
-    Calculates and prints the continuity errors.
+    Reset the timestep to maintain a constant maximum courant Number.
+    Reduction of time-step is immediate, but increase is damped to avoid
+    unstable oscillations.
 
 \*---------------------------------------------------------------------------*/
 
+if (adjustTimeStep)
 {
-    volScalarField contErr = fvc::ddt(rho) + fvc::div(phi);
+    scalar maxDeltaTFact = maxCo/(CoNum + SMALL);
+    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
-    sumLocalContErr = runTime.deltaT().value()*
+    runTime.setDeltaT
-        mag(contErr)().weightedAverage(mesh.V()).value();
+    (
+        Foam::min
+        (
+            interface.surfStabCrit().value(),
+            Foam::min
+            (
+                deltaTFact*runTime.deltaT().value(),
+                maxDeltaT
+            )
+        )
+    );
 
-    globalContErr = runTime.deltaT().value()*
+    Info<< "deltaT = " <<  runTime.deltaT().value() << endl;
-        contErr.weightedAverage(mesh.V()).value();
-
-    cumulativeContErr += globalContErr;
-
-    Info<< "time step continuity errors : sum local = " << sumLocalContErr
-        << ", global = " << globalContErr
-        << ", cumulative = " << cumulativeContErr
-        << endl;
 }
 
 // ************************************************************************* //

applications/utilities/immersedBoundary/refineImmersedBoundaryMesh/Make/options

@@ -1,13 +1,13 @@
 EXE_INC = \
-    -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude \
+    -I$(LIB_SRC)/immersedBoundary/immersedBoundaryDynamicMesh/lnInclude \
 
 EXE_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
-    -lsurfMesh \
-    -lsampling \
     -ldynamicMesh \
-    -limmersedBoundary
+    -limmersedBoundary \
+    -limmersedBoundaryDynamicMesh

applications/utilities/immersedBoundary/refineImmersedBoundaryMesh/refineImmersedBoundaryMesh.C

@@ -38,7 +38,6 @@ int main(int argc, char *argv[])
 {
     argList::validOptions.insert("ibCells", "");
     argList::validOptions.insert("ibCellCells", "");
-    argList::validOptions.insert("ibCellCellFaces", "");
 
 #   include "setRootCase.H"
 
@@ -63,13 +62,6 @@ int main(int argc, char *argv[])
             refineImmersedBoundaryMesh::IB_CELL_CELLS
         );
     }
-    else if (args.optionFound("ibCellCellFaces"))
-    {
-        rc = rib.refinementCells
-        (
-            refineImmersedBoundaryMesh::IB_CELL_CELL_FACES
-        );
-    }
     else
     {
         FatalErrorIn(args.executable())

applications/utilities/immersedBoundary/writeIbMasks/Make/options

@@ -1,9 +1,11 @@
 EXE_INC = \
+    -I$(LIB_SRC)/postProcessing/postCalc \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude
 
 EXE_LIBS = \
+    $(FOAM_LIBBIN)/postCalc.o \
     -lfiniteVolume \
     -lmeshTools \
     -lsurfMesh \

applications/utilities/immersedBoundary/writeIbMasks/writeIbMasks.C

@@ -29,32 +29,128 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "fvCFD.H"
+#include "calc.H"
+#include "fvc.H"
+#include "fvMatrices.H"
 #include "immersedBoundaryFvPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-int main(int argc, char *argv[])
+void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
 {
-#   include "setRootCase.H"
+    Info<< nl << "Calculating gamma" << endl;
-#   include "createTime.H"
+    volScalarField gamma
-#   include "createMesh.H"
+    (
-#   include "createIbMasks.H"
+        IOobject
+        (
+            "gamma",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("one", dimless, 1)
+    );
+    gamma.internalField() = mesh.V()/mesh.cellVolumes();
+    gamma.write();
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // Report minimal live cell volume
+    scalar minLiveGamma = GREAT;
+    label minLiveCell = -1;
+    const scalarField& gammaIn = gamma.internalField();
+    
+    forAll (mesh.boundary(), patchI)
+    {
+        if (isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
+        {
+            const immersedBoundaryFvPatch& ibPatch =
+                refCast<const immersedBoundaryFvPatch>
+                (
+                    mesh.boundary()[patchI]
+                );
 
-    Info<< "Time = " << runTime.timeName() << nl << endl;
+            const labelList& ibCells = ibPatch.ibPolyPatch().ibCells();
 
-    cellIbMask.write();
+            forAll (ibCells, dcI)
-    cellIbMaskExt.write();
+            {
+                if (gammaIn[ibCells[dcI]] < minLiveGamma)
+                {
+                    minLiveGamma = gammaIn[ibCells[dcI]];
+                    minLiveCell = ibCells[dcI];
+                }
+            }
+        }
+    }
 
-    Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+    Info<< "Min live cell " << minLiveCell
-        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+        << " gamma = " << minLiveGamma
-        << nl << endl;
+        << endl;
 
-    Info<< "End\n" << endl;
+    Info<< nl << "Calculating sGamma" << endl;
+    surfaceScalarField sGamma
+    (
+        IOobject
+        (
+            "sGamma",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("one", dimless, 0)
+    );
 
-    return 0;
+    const surfaceScalarField& magSf = mesh.magSf();
+    const scalarField magFaceAreas = mag(mesh.faceAreas());
+
+    sGamma.internalField() =
+        magSf.internalField()/
+        scalarField::subField(magFaceAreas, mesh.nInternalFaces());
+
+    forAll (mesh.boundary(), patchI)
+    {
+        if (!isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
+        {
+            sGamma.boundaryField()[patchI] =
+                magSf.boundaryField()[patchI]/
+                mesh.boundary()[patchI].patchSlice(magFaceAreas);
+
+            gamma.boundaryField()[patchI] =
+                sGamma.boundaryField()[patchI];
+        }
+    }
+
+    sGamma.write();
+
+    // Check consistency of face area vectors
+
+    Info<< nl << "Calculating divSf" << endl;
+    volVectorField divSf
+    (
+        "divSf",
+        fvc::div(mesh.Sf())
+    );
+    divSf.write();
+
+    // Check divergence of face area vectors
+    scalarField magDivSf = mag(divSf)().internalField();
+
+    Info<< "Face areas divergence (min, max, average): "
+        << "(" << min(magDivSf) << " " << max(magDivSf)
+        << " " << average(magDivSf) << ")"
+        << endl;
+
+    if (max(magDivSf) > 1e-9)
+    {
+        WarningIn("writeIbMasks")
+            << "Possible problem with immersed boundary face area vectors: "
+            << max(magDivSf)
+            << endl;
+    }
+
+    Info<< endl;
 }
 
 

applications/utilities/mesh/manipulation/createPatch/createPatch.C

@@ -344,7 +344,7 @@ void syncPoints
             (
                 isA<processorPolyPatch>(pp)
              && pp.nPoints() > 0
-             && refCast<const processorPolyPatch>(pp).owner()
+             && refCast<const processorPolyPatch>(pp).master()
             )
             {
                 const processorPolyPatch& procPatch =
@@ -381,7 +381,7 @@ void syncPoints
             (
                 isA<processorPolyPatch>(pp)
              && pp.nPoints() > 0
-             && !refCast<const processorPolyPatch>(pp).owner()
+             && !refCast<const processorPolyPatch>(pp).master()
             )
             {
                 const processorPolyPatch& procPatch =

applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C

@@ -82,7 +82,11 @@ int main(int argc, char *argv[])
                 )
             );
 
-            magMeshCo.write();
+            if (runTime.outputTime())
+            {
+                Info<< "Writing mesh motion Co number" << endl;
+                magMeshCo.write();
+            }
         }
 
         runTime.write();

applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C

@@ -338,30 +338,39 @@ autoPtr<mapPolyMesh> reorderMesh
     (
         new mapPolyMesh
         (
-            mesh,                       //const polyMesh& mesh,
+            mesh,                       // const polyMesh& mesh,
             mesh.nPoints(),             // nOldPoints,
             mesh.nFaces(),              // nOldFaces,
             mesh.nCells(),              // nOldCells,
+
             identity(mesh.nPoints()),   // pointMap,
             List<objectMap>(0),         // pointsFromPoints,
+
             faceOrder,                  // faceMap,
             List<objectMap>(0),         // facesFromPoints,
             List<objectMap>(0),         // facesFromEdges,
             List<objectMap>(0),         // facesFromFaces,
+
             cellOrder,                  // cellMap,
             List<objectMap>(0),         // cellsFromPoints,
             List<objectMap>(0),         // cellsFromEdges,
             List<objectMap>(0),         // cellsFromFaces,
             List<objectMap>(0),         // cellsFromCells,
+
             identity(mesh.nPoints()),   // reversePointMap,
             reverseFaceOrder,           // reverseFaceMap,
             reverseCellOrder,           // reverseCellMap,
+
             labelHashSet(0),            // flipFaceFlux,
+
             patchPointMap,              // patchPointMap,
             labelListList(0),           // pointZoneMap,
             labelListList(0),           // faceZonePointMap,
             labelListList(0),           // faceZoneFaceMap,
             labelListList(0),           // cellZoneMap,
+
+            boolList(mesh.boundaryMesh().size(), false), // resetPatchFlag
+
             pointField(0),              // preMotionPoints,
             patchStarts,                // oldPatchStarts,
             oldPatchNMeshPoints         // oldPatchNMeshPoints

applications/utilities/parallelProcessing/decomposePar/Make/files

@@ -1,12 +1,3 @@
-decomposeMesh.C
 decomposePar.C
-domainDecomposition.C
-distributeCells.C
-faMeshDecomposition.C
-fvFieldDecomposer.C
-faFieldDecomposer.C
-pointFieldDecomposer.C
-tetPointFieldDecomposer.C
-lagrangianFieldDecomposer.C
 
 EXE = $(FOAM_APPBIN)/decomposePar

applications/utilities/parallelProcessing/decomposePar/Make/options

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/decompositionMethods/decompositionMethods/lnInclude \
+    -I$(LIB_SRC)/decompositionMethods/decomposeReconstruct/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/finiteArea/lnInclude \
@@ -8,6 +9,7 @@ EXE_INC = \
 
 EXE_LIBS = \
     -ldecompositionMethods \
+    -ldecomposeReconstruct \
     -lmeshTools \
     -lfiniteVolume \
     -lfiniteArea \

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -235,26 +235,45 @@ int main(int argc, char *argv[])
     }
 
     Info<< "Create mesh for region " << regionName << endl;
-    domainDecomposition mesh
+    fvMesh mesh
     (
         IOobject
         (
             regionName,
             runTime.timeName(),
-            runTime
+            runTime,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
         )
     );
 
+    domainDecomposition meshDecomp
+    (
+        mesh,
+        IOdictionary
+        (
+            IOobject
+            (
+                "decomposeParDict",
+                runTime.system(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+        )
+    );
+
+
     // Decompose the mesh
     if (!decomposeFieldsOnly)
     {
-        mesh.decomposeMesh(filterPatches);
+        meshDecomp.decomposeMesh(filterPatches);
 
-        mesh.writeDecomposition();
+        meshDecomp.writeDecomposition();
 
         if (writeCellDist)
         {
-            const labelList& procIds = mesh.cellToProc();
+            const labelList& procIds = meshDecomp.cellToProc();
 
             // Write the decomposition as labelList for use with 'manual'
             // decomposition method.
@@ -345,7 +364,7 @@ int main(int argc, char *argv[])
 
     // Construct the point fields
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-    pointMesh pMesh(mesh);
+    const pointMesh& pMesh = pointMesh::New(mesh);
 
     PtrList<pointScalarField> pointScalarFields;
     readFields(pMesh, objects, pointScalarFields);
@@ -595,7 +614,7 @@ int main(int argc, char *argv[])
     Info<< endl;
 
     // Split the fields over processors
-    for (label procI = 0; procI < mesh.nProcs(); procI++)
+    for (label procI = 0; procI < meshDecomp.nProcs(); procI++)
     {
         Info<< "Processor " << procI << ": field transfer" << endl;
 
@@ -632,6 +651,33 @@ int main(int argc, char *argv[])
                 processorDb
             )
         );
+        procMesh.syncUpdateMesh();
+
+        labelIOList pointProcAddressing
+        (
+            IOobject
+            (
+                "pointProcAddressing",
+                procMesh.facesInstance(),
+                procMesh.meshSubDir,
+                procMesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+        );
+
+        labelIOList faceProcAddressing
+        (
+            IOobject
+            (
+                "faceProcAddressing",
+                procMesh.facesInstance(),
+                procMesh.meshSubDir,
+                procMesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+        );
 
         labelIOList cellProcAddressing
         (
@@ -674,19 +720,6 @@ int main(int argc, char *argv[])
          || surfaceTensorFields.size()
         )
         {
-            labelIOList faceProcAddressing
-            (
-                IOobject
-                (
-                    "faceProcAddressing",
-                    procMesh.facesInstance(),
-                    procMesh.meshSubDir,
-                    procMesh,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                )
-            );
-
             fvFieldDecomposer fieldDecomposer
             (
                 mesh,
@@ -720,20 +753,7 @@ int main(int argc, char *argv[])
          || pointTensorFields.size()
         )
         {
-            labelIOList pointProcAddressing
+            const pointMesh& procPMesh = pointMesh::New(procMesh);
-            (
-                IOobject
-                (
-                    "pointProcAddressing",
-                    procMesh.facesInstance(),
-                    procMesh.meshSubDir,
-                    procMesh,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                )
-            );
-
-            pointMesh procPMesh(procMesh, true);
 
             pointFieldDecomposer fieldDecomposer
             (
@@ -757,34 +777,6 @@ int main(int argc, char *argv[])
             const tetPolyMesh& tetMesh = *tetMeshPtr;
             tetPolyMesh procTetMesh(procMesh);
 
-            // Read the point addressing information
-            labelIOList pointProcAddressing
-            (
-                IOobject
-                (
-                    "pointProcAddressing",
-                    procMesh.facesInstance(),
-                    procMesh.meshSubDir,
-                    procMesh,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                )
-            );
-
-            // Read the point addressing information
-            labelIOList faceProcAddressing
-            (
-                IOobject
-                (
-                    "faceProcAddressing",
-                    procMesh.facesInstance(),
-                    procMesh.meshSubDir,
-                    procMesh,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                )
-            );
-
             tetPointFieldDecomposer fieldDecomposer
             (
                 tetMesh,
@@ -872,7 +864,7 @@ int main(int argc, char *argv[])
         {
             const fileName timePath = processorDb.timePath();
 
-            if (copyUniform || mesh.distributed())
+            if (copyUniform || meshDecomp.distributed())
             {
                 cp
                 (
@@ -962,7 +954,7 @@ int main(int argc, char *argv[])
         Info << endl;
 
         // Split the fields over processors
-        for (label procI = 0; procI < mesh.nProcs(); procI++)
+        for (label procI = 0; procI < meshDecomp.nProcs(); procI++)
         {
             Info<< "Processor " << procI
                 << ": finite area field transfer" << endl;

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -1,694 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "domainDecomposition.H"
-#include "foamTime.H"
-#include "dictionary.H"
-#include "labelIOList.H"
-#include "processorPolyPatch.H"
-#include "fvMesh.H"
-#include "OSspecific.H"
-#include "Map.H"
-#include "globalMeshData.H"
-#include "DynamicList.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-domainDecomposition::domainDecomposition(const IOobject& io)
-:
-    fvMesh(io),
-    decompositionDict_
-    (
-        IOobject
-        (
-            "decomposeParDict",
-            time().system(),
-            *this,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    ),
-    nProcs_(readInt(decompositionDict_.lookup("numberOfSubdomains"))),
-    distributed_(false),
-    cellToProc_(nCells()),
-    procPointAddressing_(nProcs_),
-    procFaceAddressing_(nProcs_),
-    nInternalProcFaces_(nProcs_),
-    nLiveProcFaces_(nProcs_),
-    procCellAddressing_(nProcs_),
-    procBoundaryAddressing_(nProcs_),
-    procPatchSize_(nProcs_),
-    procPatchStartIndex_(nProcs_),
-    procNeighbourProcessors_(nProcs_),
-    procProcessorPatchSize_(nProcs_),
-    procProcessorPatchStartIndex_(nProcs_),
-    globallySharedPoints_(0),
-    cyclicParallel_(false)
-{
-    if (decompositionDict_.found("distributed"))
-    {
-        distributed_ = Switch(decompositionDict_.lookup("distributed"));
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-domainDecomposition::~domainDecomposition()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool domainDecomposition::writeDecomposition()
-{
-    Info<< "\nConstructing processor meshes" << endl;
-
-    // Make a lookup map for globally shared points
-    Map<label> sharedPointLookup(2*globallySharedPoints_.size());
-
-    forAll (globallySharedPoints_, pointi)
-    {
-        sharedPointLookup.insert(globallySharedPoints_[pointi], pointi);
-    }
-
-
-    // Mark point/faces/cells that are in zones.  Bad coding - removed
-    // HJ, 31/Mar/2009
-
-    label totProcFaces = 0;
-    label maxProcPatches = 0;
-    label maxProcFaces = 0;
-
-
-    // Note: get cellLevel and pointLevel. Avoid checking whether they exist or
-    // not by hand. If they don't exist, simpy assume that the level is 0
-    const labelIOList globalCellLevel
-    (
-        IOobject
-        (
-            "cellLevel",
-            this->facesInstance(),
-            polyMesh::meshSubDir,
-            *this,
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
-        labelList(nCells(), 0)
-    );
-
-    const labelIOList globalPointLevel
-    (
-        IOobject
-        (
-            "pointLevel",
-            this->facesInstance(),
-            polyMesh::meshSubDir,
-            *this,
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
-        labelList(nPoints(), 0)
-    );
-
-
-    // Write out the meshes
-    for (label procI = 0; procI < nProcs_; procI++)
-    {
-        // Create processor points
-        const labelList& curPointLabels = procPointAddressing_[procI];
-
-        // Access list of all points in the mesh.  HJ, 27/Mar/2009
-        const pointField& meshPoints = allPoints();
-
-        labelList pointLookup(meshPoints.size(), -1);
-
-        pointField procPoints(curPointLabels.size());
-
-        forAll (curPointLabels, pointi)
-        {
-            procPoints[pointi] = meshPoints[curPointLabels[pointi]];
-
-            pointLookup[curPointLabels[pointi]] = pointi;
-        }
-
-        // Create processor faces
-        const labelList& curFaceLabels = procFaceAddressing_[procI];
-
-        // Access list of all faces in the mesh.  HJ, 27/Mar/2009
-        const faceList& meshFaces = allFaces();
-
-        labelList faceLookup(meshFaces.size(), -1);
-
-        faceList procFaces(curFaceLabels.size());
-
-        forAll (curFaceLabels, facei)
-        {
-            // Mark the original face as used
-            // Remember to decrement the index by one (turning index)
-            // HJ, 5/Dec/2001
-            label curF = mag(curFaceLabels[facei]) - 1;
-
-            faceLookup[curF] = facei;
-
-            // get the original face
-            labelList origFaceLabels;
-
-            if (curFaceLabels[facei] >= 0)
-            {
-                // face not turned
-                origFaceLabels = meshFaces[curF];
-            }
-            else
-            {
-                origFaceLabels = meshFaces[curF].reverseFace();
-            }
-
-            // translate face labels into local point list
-            face& procFaceLabels = procFaces[facei];
-
-            procFaceLabels.setSize(origFaceLabels.size());
-
-            forAll (origFaceLabels, pointi)
-            {
-                procFaceLabels[pointi] = pointLookup[origFaceLabels[pointi]];
-            }
-        }
-
-        // Create cell lookup
-        labelList cellLookup(nCells(), -1);
-        const labelList& curCellLabels = procCellAddressing_[procI];
-
-        forAll (curCellLabels, cellI)
-        {
-            cellLookup[curCellLabels[cellI]] = cellI;
-        }
-
-        // Get complete owner-neighour addressing in the mesh
-        const labelList& own = faceOwner();
-        const labelList& nei = faceNeighbour();
-
-        // Calculate owner and neighbour list
-        // Owner list is sized to number of live faces
-        // Neighbour list is sized to number of internal faces
-
-        labelList procOwner(nLiveProcFaces_[procI]);
-
-        // Note: loop over owner, not all faces: sizes are different
-        forAll (procOwner, faceI)
-        {
-            // Remember to decrement the index by one (turning index)
-            // HJ, 28/Mar/2009
-            label curF = mag(curFaceLabels[faceI]) - 1;
-
-            if (curFaceLabels[faceI] >= 0)
-            {
-                procOwner[faceI] = cellLookup[own[curF]];
-            }
-            else
-            {
-                procOwner[faceI] = cellLookup[nei[curF]];
-            }
-        }
-
-        labelList procNeighbour(nInternalProcFaces_[procI]);
-
-        // Note: loop over neighbour, not all faces: sizes are different
-        forAll (procNeighbour, faceI)
-        {
-            // Remember to decrement the index by one (turning index)
-            // HJ, 28/Mar/2009
-            label curF = mag(curFaceLabels[faceI]) - 1;
-
-            if (curFaceLabels[faceI] >= 0)
-            {
-                procNeighbour[faceI] = cellLookup[nei[curF]];
-            }
-            else
-            {
-                procNeighbour[faceI] = cellLookup[own[curF]];
-            }
-        }
-
-        // Create processor cells.  No longer needed: using owner and neighbour
-        // HJ, 28/Mar/2009
-//         const cellList& meshCells = cells();
-
-//         cellList procCells(curCellLabels.size());
-
-//         forAll (curCellLabels, cellI)
-//         {
-//             const labelList& origCellLabels = meshCells[curCellLabels[cellI]];
-
-//             cell& curCell = procCells[cellI];
-
-//             curCell.setSize(origCellLabels.size());
-
-//             forAll (origCellLabels, cellFaceI)
-//             {
-//                 curCell[cellFaceI] = faceLookup[origCellLabels[cellFaceI]];
-//             }
-//         }
-
-        // Create processor mesh without a boundary
-
-        fileName processorCasePath
-        (
-            time().caseName()/fileName(word("processor") + Foam::name(procI))
-        );
-
-        // make the processor directory
-        mkDir(time().rootPath()/processorCasePath);
-
-        // create a database
-        Time processorDb
-        (
-            Time::controlDictName,
-            time().rootPath(),
-            processorCasePath,
-            "system",
-            "constant",
-            true
-        );
-
-        // Create the mesh
-        polyMesh procMesh
-        (
-            IOobject
-            (
-                this->polyMesh::name(),  // region name of undecomposed mesh
-                pointsInstance(),
-                processorDb
-            ),
-            xferMove(procPoints),
-            xferMove(procFaces),
-            xferMove(procOwner),
-            xferMove(procNeighbour),
-            false          // Do not sync par
-//   xferMove(procCells)   // Old-fashioned mesh creation using cells.
-                           // Deprecated: using face owner/neighbour
-                           // HJ, 30/Mar/2009
-        );
-
-        // Create processor boundary patches
-        const labelList& curBoundaryAddressing = procBoundaryAddressing_[procI];
-
-        const labelList& curPatchSizes = procPatchSize_[procI];
-
-        const labelList& curPatchStarts = procPatchStartIndex_[procI];
-
-        const labelList& curNeighbourProcessors =
-            procNeighbourProcessors_[procI];
-
-        const labelList& curProcessorPatchSizes =
-            procProcessorPatchSize_[procI];
-
-        const labelList& curProcessorPatchStarts =
-            procProcessorPatchStartIndex_[procI];
-
-        const polyPatchList& meshPatches = boundaryMesh();
-
-        List<polyPatch*> procPatches
-        (
-            curPatchSizes.size()
-          + curProcessorPatchSizes.size(),
-            reinterpret_cast<polyPatch*>(0)
-        );
-
-        label nPatches = 0;
-
-        forAll (curPatchSizes, patchi)
-        {
-            procPatches[nPatches] =
-                meshPatches[curBoundaryAddressing[patchi]].clone
-                (
-                    procMesh.boundaryMesh(),
-                    nPatches,
-                    curPatchSizes[patchi],
-                    curPatchStarts[patchi]
-                ).ptr();
-
-            nPatches++;
-        }
-
-        forAll (curProcessorPatchSizes, procPatchI)
-        {
-            procPatches[nPatches] =
-                new processorPolyPatch
-                (
-                    word("procBoundary") + Foam::name(procI)
-                  + word("to")
-                  + Foam::name(curNeighbourProcessors[procPatchI]),
-                    curProcessorPatchSizes[procPatchI],
-                    curProcessorPatchStarts[procPatchI],
-                    nPatches,
-                    procMesh.boundaryMesh(),
-                    procI,
-                    curNeighbourProcessors[procPatchI]
-            );
-
-            nPatches++;
-        }
-
-        // Add boundary patches
-        procMesh.addPatches(procPatches);
-
-        // Create and add zones
-
-        // Note:
-        // This coding was all wrong, as each point/face/cell may only
-        //  belong to a single zone.
-        // Additionally, ordering of points/faces/cells in the processor mesh
-        // needs to match the ordering in global mesh zones.  Full rewrite.
-        // HJ, 30/Mar/2009
-
-        // Create zones if needed
-        if
-        (
-            pointZones().size() > 0
-         || faceZones().size() > 0
-         || cellZones().size() > 0
-        )
-        {
-            // Point zones
-            List<pointZone*> procPz(pointZones().size());
-
-            {
-                const pointZoneMesh& pz = pointZones();
-
-                // Go through all the zoned points and find out if they
-                // belong to a processor.  If so, add it to the zone as
-                // necessary
-                forAll (pz, zoneI)
-                {
-                    const labelList& zonePoints = pz[zoneI];
-
-                    labelList procZonePoints(zonePoints.size());
-                    label nZonePoints = 0;
-
-                    forAll (zonePoints, pointI)
-                    {
-                        const label localIndex =
-                            pointLookup[zonePoints[pointI]];
-
-                        if (localIndex >= 0)
-                        {
-                            // Point live on processor: add to zone
-                            procZonePoints[nZonePoints] = localIndex;
-                            nZonePoints++;
-                        }
-                    }
-
-                    // Add the zone
-                    procZonePoints.setSize(nZonePoints);
-
-                    procPz[zoneI] = new pointZone
-                    (
-                        pz[zoneI].name(),
-                        procZonePoints,
-                        zoneI,
-                        procMesh.pointZones()
-                    );
-                }
-            }
-
-
-            // Face zones
-            List<faceZone*> procFz(faceZones().size());
-
-            {
-                const faceZoneMesh& fz = faceZones();
-
-                forAll (fz, zoneI)
-                {
-                    const labelList& zoneFaces = fz[zoneI];
-                    const boolList& flipMap = fz[zoneI].flipMap();
-
-                    // Go through all the zoned faces and find out if they
-                    // belong to a processor.  If so, add it to the zone as
-                    // necessary
-
-                    labelList procZoneFaces(zoneFaces.size());
-                    boolList procZoneFaceFlips(zoneFaces.size());
-                    label nZoneFaces = 0;
-
-                    forAll (zoneFaces, faceI)
-                    {
-                        const label localIndex = faceLookup[zoneFaces[faceI]];
-
-                        if (localIndex >= 0)
-                        {
-                            // Face is present on the processor
-
-                            // Add the face to the zone
-                            procZoneFaces[nZoneFaces] = localIndex;
-
-                            // Grab the flip
-                            bool flip = flipMap[faceI];
-
-                            if (curFaceLabels[localIndex] < 0)
-                            {
-                                flip = !flip;
-                            }
-
-                            procZoneFaceFlips[nZoneFaces] = flip;
-                            nZoneFaces++;
-                        }
-                    }
-
-                    // Add the zone
-                    procZoneFaces.setSize(nZoneFaces);
-                    procZoneFaceFlips.setSize(nZoneFaces);
-
-                    procFz[zoneI] = new faceZone
-                    (
-                        fz[zoneI].name(),
-                        procZoneFaces,
-                        procZoneFaceFlips,
-                        zoneI,
-                        procMesh.faceZones()
-                    );
-                }
-            }
-
-            // Cell zones
-            List<cellZone*> procCz(cellZones().size());
-
-            {
-                const cellZoneMesh& cz = cellZones();
-
-                // Go through all the zoned cells and find out if they
-                // belong to a processor.  If so, add it to the zone as
-                // necessary
-
-                forAll (cz, zoneI)
-                {
-                    const labelList& zoneCells = cz[zoneI];
-
-                    labelList procZoneCells(zoneCells.size());
-                    label nZoneCells = 0;
-
-                    forAll (zoneCells, cellI)
-                    {
-                        const label localIndex = cellLookup[zoneCells[cellI]];
-
-                        if (localIndex >= 0)
-                        {
-                            procZoneCells[nZoneCells] = localIndex;
-                            nZoneCells++;
-                        }
-                    }
-
-                    // Add the zone
-                    procZoneCells.setSize(nZoneCells);
-
-                    procCz[zoneI] = new cellZone
-                    (
-                        cz[zoneI].name(),
-                        procZoneCells,
-                        zoneI,
-                        procMesh.cellZones()
-                    );
-                }
-            }
-
-            // Add zones
-            procMesh.addZones(procPz, procFz, procCz);
-        }
-
-        // Set the precision of the points data to 10
-        IOstream::defaultPrecision(10);
-
-        procMesh.write();
-
-        Info<< endl
-            << "Processor " << procI << nl
-            << "    Number of cells = " << procMesh.nCells()
-            << endl;
-
-        label nBoundaryFaces = 0;
-        label nProcPatches = 0;
-        label nProcFaces = 0;
-
-        forAll (procMesh.boundaryMesh(), patchi)
-        {
-            if
-            (
-                procMesh.boundaryMesh()[patchi].type()
-             == processorPolyPatch::typeName
-            )
-            {
-                const processorPolyPatch& ppp =
-                refCast<const processorPolyPatch>
-                (
-                    procMesh.boundaryMesh()[patchi]
-                );
-
-                Info<< "    Number of faces shared with processor "
-                    << ppp.neighbProcNo() << " = " << ppp.size() << endl;
-
-                nProcPatches++;
-                nProcFaces += ppp.size();
-            }
-            else
-            {
-                nBoundaryFaces += procMesh.boundaryMesh()[patchi].size();
-            }
-        }
-
-        Info<< "    Number of processor patches = " << nProcPatches << nl
-            << "    Number of processor faces = " << nProcFaces << nl
-            << "    Number of boundary faces = " << nBoundaryFaces << endl;
-
-        totProcFaces += nProcFaces;
-        maxProcPatches = max(maxProcPatches, nProcPatches);
-        maxProcFaces = max(maxProcFaces, nProcFaces);
-
-        // create and write the addressing information
-        labelIOList pointProcAddressing
-        (
-            IOobject
-            (
-                "pointProcAddressing",
-                procMesh.facesInstance(),
-                procMesh.meshSubDir,
-                procMesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            procPointAddressing_[procI]
-        );
-        pointProcAddressing.write();
-
-        labelIOList faceProcAddressing
-        (
-            IOobject
-            (
-                "faceProcAddressing",
-                procMesh.facesInstance(),
-                procMesh.meshSubDir,
-                procMesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            procFaceAddressing_[procI]
-        );
-        faceProcAddressing.write();
-
-        labelIOList cellProcAddressing
-        (
-            IOobject
-            (
-                "cellProcAddressing",
-                procMesh.facesInstance(),
-                procMesh.meshSubDir,
-                procMesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            procCellAddressing_[procI]
-        );
-        cellProcAddressing.write();
-
-        labelIOList boundaryProcAddressing
-        (
-            IOobject
-            (
-                "boundaryProcAddressing",
-                procMesh.facesInstance(),
-                procMesh.meshSubDir,
-                procMesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            procBoundaryAddressing_[procI]
-        );
-        boundaryProcAddressing.write();
-
-        // Create and write cellLevel and pointLevel information
-        const unallocLabelList& cellMap = cellProcAddressing;
-        labelIOList procCellLevel
-        (
-            IOobject
-            (
-                "cellLevel",
-                procMesh.facesInstance(),
-                procMesh.meshSubDir,
-                procMesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            labelList(globalCellLevel, cellMap)
-        );
-        procCellLevel.write();
-
-        const unallocLabelList& pointMap = pointProcAddressing;
-        labelIOList procPointLevel
-        (
-            IOobject
-            (
-                "pointLevel",
-                procMesh.facesInstance(),
-                procMesh.meshSubDir,
-                procMesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            labelList(globalPointLevel, pointMap)
-        );
-        procPointLevel.write();
-    }
-
-    Info<< nl
-        << "Number of processor faces = " << totProcFaces/2 << nl
-        << "Max number of processor patches = " << maxProcPatches << nl
-        << "Max number of faces between processors = " << maxProcFaces
-        << endl;
-
-    return true;
-}
-
-
-// ************************************************************************* //

applications/utilities/parallelProcessing/decomposePar/readFields.C

@@ -48,7 +48,7 @@ void Foam::readFields
     // Construct the vol scalar fields
     fields.setSize(fieldObjects.size());
 
-    label fieldi=0;
+    label fieldI = 0;
     for
     (
         IOobjectList::iterator iter = fieldObjects.begin();
@@ -58,7 +58,7 @@ void Foam::readFields
     {
         fields.set
         (
-            fieldi++,
+            fieldI++,
             new GeoField
             (
                 *iter(),

applications/utilities/parallelProcessing/reconstructPar/Make/files

@@ -1,10 +1,4 @@
 processorMeshes.C
-processorFaMeshes.C
-fvFieldReconstructor.C
-faFieldReconstructor.C
-pointFieldReconstructor.C
-tetPointFieldReconstructor.C
-reconstructLagrangianPositions.C
 reconstructPar.C
 
 EXE = $(FOAM_APPBIN)/reconstructPar

applications/utilities/parallelProcessing/reconstructPar/Make/options

@@ -1,12 +1,17 @@
 EXE_INC = \
+    -I$(LIB_SRC)/decompositionMethods/decompositionMethods/lnInclude \
+    -I$(LIB_SRC)/decompositionMethods/decomposeReconstruct/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/tetFiniteElement/lnInclude
 
 EXE_LIBS = \
+    -ldecompositionMethods \
+    -ldecomposeReconstruct \
+    -lmeshTools \
     -lfiniteVolume \
     -lfiniteArea \
     -llagrangian \
-    -lmeshTools \
     -ltetFiniteElement

applications/utilities/parallelProcessing/reconstructPar/faFieldReconstructorReconstructFields.C

@@ -1,640 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "faFieldReconstructor.H"
-#include "foamTime.H"
-#include "PtrList.H"
-#include "faPatchFields.H"
-#include "emptyFaPatch.H"
-#include "emptyFaPatchField.H"
-#include "emptyFaePatchField.H"
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::tmp<Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh> >
-Foam::faFieldReconstructor::reconstructFaAreaField
-(
-    const IOobject& fieldIoObject
-)
-{
-    // Read the field for all the processors
-    PtrList<GeometricField<Type, faPatchField, areaMesh> > procFields
-    (
-        procMeshes_.size()
-    );
-
-    forAll (procMeshes_, procI)
-    {
-        procFields.set
-        (
-            procI,
-            new GeometricField<Type, faPatchField, areaMesh>
-            (
-                IOobject
-                (
-                    fieldIoObject.name(),
-                    procMeshes_[procI].time().timeName(),
-                    procMeshes_[procI](),
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                ),
-                procMeshes_[procI]
-            )
-        );
-    }
-
-    // Create the internalField
-    Field<Type> internalField(mesh_.nFaces());
-
-    // Create the patch fields
-    PtrList<faPatchField<Type> > patchFields(mesh_.boundary().size());
-
-
-    // Create global mesh patchs starts
-
-    labelList gStarts(mesh_.boundary().size(), -1);
-
-    if (mesh_.boundary().size() > 0)
-    {
-        gStarts[0] = mesh_.nInternalEdges();
-    }
-
-    for(label i=1; i<mesh_.boundary().size(); i++)
-    {
-        gStarts[i] = gStarts[i-1] + mesh_.boundary()[i-1].labelList::size();
-    }
-
-    forAll (procMeshes_, procI)
-    {
-        const GeometricField<Type, faPatchField, areaMesh>& procField =
-            procFields[procI];
-
-        // Set the face values in the reconstructed field
-        internalField.rmap
-        (
-            procField.internalField(),
-            faceProcAddressing_[procI]
-        );
-
-
-
-        // Set the boundary patch values in the reconstructed field
-
-        labelList starts(procMeshes_[procI].boundary().size(), -1);
-
-        if(procMeshes_[procI].boundary().size() > 0)
-        {
-            starts[0] = procMeshes_[procI].nInternalEdges();
-        }
-
-        for(label i=1; i<procMeshes_[procI].boundary().size(); i++)
-        {
-            starts[i] =
-                starts[i-1]
-              + procMeshes_[procI].boundary()[i-1].labelList::size();
-        }
-
-        forAll(boundaryProcAddressing_[procI], patchI)
-        {
-            // Get patch index of the original patch
-            const label curBPatch = boundaryProcAddressing_[procI][patchI];
-
-            // Get addressing slice for this patch
-
-//             const labelList::subList cp =
-//                 procMeshes_[procI].boundary()[patchI].patchSlice
-//                 (
-//                     edgeProcAddressing_[procI]
-//                 );
-
-            const labelList::subList cp =
-                labelList::subList
-                (
-                    edgeProcAddressing_[procI],
-                    procMeshes_[procI].boundary()[patchI].size(),
-                    starts[patchI]
-                );
-
-            // check if the boundary patch is not a processor patch
-            if (curBPatch >= 0)
-            {
-                // Regular patch. Fast looping
-
-                if (!patchFields(curBPatch))
-                {
-                    patchFields.set
-                    (
-                        curBPatch,
-                        faPatchField<Type>::New
-                        (
-                            procField.boundaryField()[patchI],
-                            mesh_.boundary()[curBPatch],
-                            DimensionedField<Type, areaMesh>::null(),
-                            faPatchFieldReconstructor
-                            (
-                                mesh_.boundary()[curBPatch].size(),
-                                procField.boundaryField()[patchI].size()
-                            )
-                        )
-                    );
-                }
-
-                const label curPatchStart = gStarts[curBPatch];
-//                     mesh_.boundary()[curBPatch].start();
-
-                labelList reverseAddressing(cp.size());
-
-                forAll(cp, edgeI)
-                {
-                    // Subtract one to take into account offsets for
-                    // face direction.
-//                     reverseAddressing[edgeI] = cp[edgeI] - 1 - curPatchStart;
-                    reverseAddressing[edgeI] = cp[edgeI] - curPatchStart;
-                }
-
-                patchFields[curBPatch].rmap
-                (
-                    procField.boundaryField()[patchI],
-                    reverseAddressing
-                );
-            }
-            else
-            {
-                const Field<Type>& curProcPatch =
-                    procField.boundaryField()[patchI];
-
-                // In processor patches, there's a mix of internal faces (some
-                // of them turned) and possible cyclics. Slow loop
-                forAll(cp, edgeI)
-                {
-                    // Subtract one to take into account offsets for
-                    // face direction.
-//                     label curE = cp[edgeI] - 1;
-                    label curE = cp[edgeI];
-
-                    // Is the face on the boundary?
-                    if (curE >= mesh_.nInternalEdges())
-                    {
-//                         label curBPatch = mesh_.boundary().whichPatch(curE);
-                        label curBPatch = -1;
-
-                        forAll (mesh_.boundary(), pI)
-                        {
-                            if
-                            (
-                                curE >= gStarts[pI]
-                             && curE <
-                                (
-                                    gStarts[pI]
-                                  + mesh_.boundary()[pI].labelList::size()
-                                )
-                            )
-                            {
-                                curBPatch = pI;
-                            }
-                        }
-
-                        if (!patchFields(curBPatch))
-                        {
-                            patchFields.set
-                            (
-                                curBPatch,
-                                faPatchField<Type>::New
-                                (
-                                    mesh_.boundary()[curBPatch].type(),
-                                    mesh_.boundary()[curBPatch],
-                                    DimensionedField<Type, areaMesh>::null()
-                                )
-                            );
-                        }
-
-                        // add the edge
-//                         label curPatchEdge =
-//                             mesh_.boundary()
-//                                 [curBPatch].whichEdge(curE);
-
-                        label curPatchEdge = curE - gStarts[curBPatch];
-
-                        patchFields[curBPatch][curPatchEdge] =
-                            curProcPatch[edgeI];
-                    }
-                }
-            }
-        }
-    }
-
-    forAll(mesh_.boundary(), patchI)
-    {
-        // add empty patches
-        if
-        (
-            typeid(mesh_.boundary()[patchI]) == typeid(emptyFaPatch)
-         && !patchFields(patchI)
-        )
-        {
-            patchFields.set
-            (
-                patchI,
-                faPatchField<Type>::New
-                (
-                    emptyFaPatchField<Type>::typeName,
-                    mesh_.boundary()[patchI],
-                    DimensionedField<Type, areaMesh>::null()
-                )
-            );
-        }
-    }
-
-
-    // Now construct and write the field
-    // setting the internalField and patchFields
-    return tmp<GeometricField<Type, faPatchField, areaMesh> >
-    (
-        new GeometricField<Type, faPatchField, areaMesh>
-        (
-            IOobject
-            (
-                fieldIoObject.name(),
-                mesh_.time().timeName(),
-                mesh_(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            procFields[0].dimensions(),
-            internalField,
-            patchFields
-        )
-    );
-}
-
-
-template<class Type>
-Foam::tmp<Foam::GeometricField<Type, Foam::faePatchField, Foam::edgeMesh> >
-Foam::faFieldReconstructor::reconstructFaEdgeField
-(
-    const IOobject& fieldIoObject
-)
-{
-    // Read the field for all the processors
-    PtrList<GeometricField<Type, faePatchField, edgeMesh> > procFields
-    (
-        procMeshes_.size()
-    );
-
-    forAll (procMeshes_, procI)
-    {
-        procFields.set
-        (
-            procI,
-            new GeometricField<Type, faePatchField, edgeMesh>
-            (
-                IOobject
-                (
-                    fieldIoObject.name(),
-                    procMeshes_[procI].time().timeName(),
-                    procMeshes_[procI](),
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                ),
-                procMeshes_[procI]
-            )
-        );
-    }
-
-
-    // Create the internalField
-    Field<Type> internalField(mesh_.nInternalEdges());
-
-    // Create the patch fields
-    PtrList<faePatchField<Type> > patchFields(mesh_.boundary().size());
-
-
-    labelList gStarts(mesh_.boundary().size(), -1);
-
-    if(mesh_.boundary().size() > 0)
-    {
-        gStarts[0] = mesh_.nInternalEdges();
-    }
-
-    for(label i=1; i<mesh_.boundary().size(); i++)
-    {
-        gStarts[i] = gStarts[i-1] + mesh_.boundary()[i-1].labelList::size();
-    }
-
-
-    forAll (procMeshes_, procI)
-    {
-        const GeometricField<Type, faePatchField, edgeMesh>& procField =
-            procFields[procI];
-
-        // Set the face values in the reconstructed field
-
-        // It is necessary to create a copy of the addressing array to
-        // take care of the face direction offset trick.
-        //
-        {
-            labelList curAddr(edgeProcAddressing_[procI]);
-
-//             forAll (curAddr, addrI)
-//             {
-//                 curAddr[addrI] -= 1;
-//             }
-
-            internalField.rmap
-            (
-                procField.internalField(),
-                curAddr
-            );
-        }
-
-        // Set the boundary patch values in the reconstructed field
-
-        labelList starts(procMeshes_[procI].boundary().size(), -1);
-
-        if(procMeshes_[procI].boundary().size() > 0)
-        {
-            starts[0] = procMeshes_[procI].nInternalEdges();
-        }
-
-        for(label i=1; i<procMeshes_[procI].boundary().size(); i++)
-        {
-            starts[i] =
-                starts[i-1]
-              + procMeshes_[procI].boundary()[i-1].labelList::size();
-        }
-
-        forAll(boundaryProcAddressing_[procI], patchI)
-        {
-            // Get patch index of the original patch
-            const label curBPatch = boundaryProcAddressing_[procI][patchI];
-
-            // Get addressing slice for this patch
-
-//             const labelList::subList cp =
-//                 procMeshes_[procI].boundary()[patchI].patchSlice
-//                 (
-//                     faceProcAddressing_[procI]
-//                 );
-
-            const labelList::subList cp =
-                labelList::subList
-                (
-                    edgeProcAddressing_[procI],
-                    procMeshes_[procI].boundary()[patchI].size(),
-                    starts[patchI]
-                );
-
-            // check if the boundary patch is not a processor patch
-            if (curBPatch >= 0)
-            {
-                // Regular patch. Fast looping
-
-                if (!patchFields(curBPatch))
-                {
-                    patchFields.set
-                    (
-                        curBPatch,
-                        faePatchField<Type>::New
-                        (
-                            procField.boundaryField()[patchI],
-                            mesh_.boundary()[curBPatch],
-                            DimensionedField<Type, edgeMesh>::null(),
-                            faPatchFieldReconstructor
-                            (
-                                mesh_.boundary()[curBPatch].size(),
-                                procField.boundaryField()[patchI].size()
-                            )
-                        )
-                    );
-                }
-
-                const label curPatchStart = gStarts[curBPatch];
-//                     mesh_.boundary()[curBPatch].start();
-
-                labelList reverseAddressing(cp.size());
-
-                forAll(cp, edgeI)
-                {
-                    // Subtract one to take into account offsets for
-                    // face direction.
-//                     reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart;
-                    reverseAddressing[edgeI] = cp[edgeI] - curPatchStart;
-                }
-
-                patchFields[curBPatch].rmap
-                (
-                    procField.boundaryField()[patchI],
-                    reverseAddressing
-                );
-            }
-            else
-            {
-                const Field<Type>& curProcPatch =
-                    procField.boundaryField()[patchI];
-
-                // In processor patches, there's a mix of internal faces (some
-                // of them turned) and possible cyclics. Slow loop
-                forAll(cp, edgeI)
-                {
-//                     label curF = cp[edgeI] - 1;
-                    label curE = cp[edgeI];
-
-                    // Is the face turned the right side round
-                    if (curE >= 0)
-                    {
-                        // Is the face on the boundary?
-                        if (curE >= mesh_.nInternalEdges())
-                        {
-//                             label curBPatch =
-//                                 mesh_.boundary().whichPatch(curF);
-
-                            label curBPatch = -1;
-
-                            forAll (mesh_.boundary(), pI)
-                            {
-                                if
-                                (
-                                    curE >= gStarts[pI]
-                                 && curE <
-                                    (
-                                        gStarts[pI]
-                                      + mesh_.boundary()[pI].labelList::size()
-                                    )
-                                )
-                                {
-                                    curBPatch = pI;
-                                }
-                            }
-
-                            if (!patchFields(curBPatch))
-                            {
-                                patchFields.set
-                                (
-                                    curBPatch,
-                                    faePatchField<Type>::New
-                                    (
-                                        mesh_.boundary()[curBPatch].type(),
-                                        mesh_.boundary()[curBPatch],
-                                        DimensionedField<Type, edgeMesh>
-                                           ::null()
-                                    )
-                                );
-                            }
-
-                            // add the face
-//                             label curPatchFace =
-//                                 mesh_.boundary()
-//                                 [curBPatch].whichEdge(curF);
-
-                            label curPatchEdge = curE - gStarts[curBPatch];
-
-                            patchFields[curBPatch][curPatchEdge] =
-                                curProcPatch[edgeI];
-                        }
-                        else
-                        {
-                            // Internal face
-                            internalField[curE] = curProcPatch[edgeI];
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-    forAll(mesh_.boundary(), patchI)
-    {
-        // add empty patches
-        if
-        (
-            typeid(mesh_.boundary()[patchI]) == typeid(emptyFaPatch)
-         && !patchFields(patchI)
-        )
-        {
-            patchFields.set
-            (
-                patchI,
-                faePatchField<Type>::New
-                (
-                    emptyFaePatchField<Type>::typeName,
-                    mesh_.boundary()[patchI],
-                    DimensionedField<Type, edgeMesh>::null()
-                )
-            );
-        }
-    }
-
-
-    // Now construct and write the field
-    // setting the internalField and patchFields
-    return tmp<GeometricField<Type, faePatchField, edgeMesh> >
-    (
-        new GeometricField<Type, faePatchField, edgeMesh>
-        (
-            IOobject
-            (
-                fieldIoObject.name(),
-                mesh_.time().timeName(),
-                mesh_(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            procFields[0].dimensions(),
-            internalField,
-            patchFields
-        )
-    );
-}
-
-
-// Reconstruct and write all area fields
-template<class Type>
-void Foam::faFieldReconstructor::reconstructFaAreaFields
-(
-    const IOobjectList& objects
-)
-{
-    const word& fieldClassName =
-        GeometricField<Type, faPatchField, areaMesh>::typeName;
-
-    IOobjectList fields = objects.lookupClass(fieldClassName);
-
-    if (fields.size())
-    {
-        Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
-
-        for
-        (
-            IOobjectList::const_iterator fieldIter = fields.begin();
-            fieldIter != fields.end();
-            ++fieldIter
-        )
-        {
-            Info << "        " << fieldIter()->name() << endl;
-
-            reconstructFaAreaField<Type>(*fieldIter())().write();
-        }
-
-        Info<< endl;
-    }
-}
-
-// Reconstruct and write all edge fields
-template<class Type>
-void Foam::faFieldReconstructor::reconstructFaEdgeFields
-(
-    const IOobjectList& objects
-)
-{
-    const word& fieldClassName =
-        GeometricField<Type, faePatchField, edgeMesh>::typeName;
-
-    IOobjectList fields = objects.lookupClass(fieldClassName);
-
-    if (fields.size())
-    {
-        Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
-
-        for
-        (
-            IOobjectList::const_iterator fieldIter = fields.begin();
-            fieldIter != fields.end();
-            ++fieldIter
-        )
-        {
-            Info<< "        " << fieldIter()->name() << endl;
-
-            reconstructFaEdgeField<Type>(*fieldIter())().write();
-        }
-
-        Info<< endl;
-    }
-}
-
-
-// ************************************************************************* //

applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H

@@ -1,196 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::fvFieldReconstructor
-
-Description
-    FV volume and surface field reconstructor.
-
-SourceFiles
-    fvFieldReconstructor.C
-    fvFieldReconstructorReconstructFields.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef fvFieldReconstructor_H
-#define fvFieldReconstructor_H
-
-#include "PtrList.H"
-#include "fvMesh.H"
-#include "IOobjectList.H"
-#include "fvPatchFieldMapper.H"
-#include "labelIOList.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-
-/*---------------------------------------------------------------------------*\
-                 Class fvFieldReconstructor Declaration
-\*---------------------------------------------------------------------------*/
-
-class fvFieldReconstructor
-{
-    // Private data
-
-        //- Reconstructed mesh reference
-        fvMesh& mesh_;
-
-        //- List of processor meshes
-        const PtrList<fvMesh>& procMeshes_;
-
-        //- List of processor face addressing lists
-        const PtrList<labelIOList>& faceProcAddressing_;
-
-        //- List of processor cell addressing lists
-        const PtrList<labelIOList>& cellProcAddressing_;
-
-        //- List of processor boundary addressing lists
-        const PtrList<labelIOList>& boundaryProcAddressing_;
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        fvFieldReconstructor(const fvFieldReconstructor&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const fvFieldReconstructor&);
-
-
-public:
-
-        class fvPatchFieldReconstructor
-        :
-            public fvPatchFieldMapper
-        {
-            label size_;
-            label sizeBeforeMapping_;
-
-        public:
-
-            // Constructors
-
-                //- Construct given size
-                fvPatchFieldReconstructor
-                (
-                    const label size,
-                    const label sizeBeforeMapping
-                )
-                :
-                    size_(size),
-                    sizeBeforeMapping_(sizeBeforeMapping)
-                {}
-
-
-            // Member functions
-
-                virtual label size() const
-                {
-                    return size_;
-                }
-
-                virtual label sizeBeforeMapping() const
-                {
-                    return sizeBeforeMapping_;
-                }
-
-                virtual bool direct() const
-                {
-                    return true;
-                }
-
-                virtual const unallocLabelList& directAddressing() const
-                {
-                    return unallocLabelList::null();
-                }
-        };
-
-
-    // Constructors
-
-        //- Construct from components
-        fvFieldReconstructor
-        (
-            fvMesh& mesh,
-            const PtrList<fvMesh>& procMeshes,
-            const PtrList<labelIOList>& faceProcAddressing,
-            const PtrList<labelIOList>& cellProcAddressing,
-            const PtrList<labelIOList>& boundaryProcAddressing
-        );
-
-
-    // Member Functions
-
-        //- Reconstruct volume field
-        template<class Type>
-        tmp<GeometricField<Type, fvPatchField, volMesh> >
-        reconstructFvVolumeField
-        (
-            const IOobject& fieldIoObject
-        );
-
-        //- Reconstruct surface field
-        template<class Type>
-        tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
-        reconstructFvSurfaceField
-        (
-            const IOobject& fieldIoObject
-        );
-
-        //- Reconstruct and write all/selected volume fields
-        template<class Type>
-        void reconstructFvVolumeFields
-        (
-            const IOobjectList& objects,
-            const HashSet<word>& selectedFields
-        );
-
-        //- Reconstruct and write all/selected volume fields
-        template<class Type>
-        void reconstructFvSurfaceFields
-        (
-            const IOobjectList& objects,
-            const HashSet<word>& selectedFields
-        );
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "fvFieldReconstructorReconstructFields.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C

@@ -1,526 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvFieldReconstructor.H"
-#include "foamTime.H"
-#include "PtrList.H"
-#include "fvPatchFields.H"
-#include "emptyFvPatch.H"
-#include "emptyFvPatchField.H"
-#include "emptyFvsPatchField.H"
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
-Foam::fvFieldReconstructor::reconstructFvVolumeField
-(
-    const IOobject& fieldIoObject
-)
-{
-    // Read the field for all the processors
-    PtrList<GeometricField<Type, fvPatchField, volMesh> > procFields
-    (
-        procMeshes_.size()
-    );
-
-    forAll (procMeshes_, procI)
-    {
-        procFields.set
-        (
-            procI,
-            new GeometricField<Type, fvPatchField, volMesh>
-            (
-                IOobject
-                (
-                    fieldIoObject.name(),
-                    procMeshes_[procI].time().timeName(),
-                    procMeshes_[procI],
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                ),
-                procMeshes_[procI]
-            )
-        );
-    }
-
-
-    // Create the internalField
-    Field<Type> internalField(mesh_.nCells());
-
-    // Create the patch fields
-    PtrList<fvPatchField<Type> > patchFields(mesh_.boundary().size());
-
-
-    forAll (procMeshes_, procI)
-    {
-        const GeometricField<Type, fvPatchField, volMesh>& procField =
-            procFields[procI];
-
-        // Set the cell values in the reconstructed field
-        internalField.rmap
-        (
-            procField.internalField(),
-            cellProcAddressing_[procI]
-        );
-
-        // Set the boundary patch values in the reconstructed field
-        forAll (boundaryProcAddressing_[procI], patchI)
-        {
-            // Get patch index of the original patch
-            const label curBPatch = boundaryProcAddressing_[procI][patchI];
-
-            // Get addressing slice for this patch
-            const labelList::subList cp =
-                procMeshes_[procI].boundary()[patchI].patchSlice
-                (
-                    faceProcAddressing_[procI]
-                );
-
-            // check if the boundary patch is not a processor patch
-            if (curBPatch >= 0)
-            {
-                // Regular patch. Fast looping
-
-                if (!patchFields(curBPatch))
-                {
-                    patchFields.set
-                    (
-                        curBPatch,
-                        fvPatchField<Type>::New
-                        (
-                            procField.boundaryField()[patchI],
-                            mesh_.boundary()[curBPatch],
-                            DimensionedField<Type, volMesh>::null(),
-                            fvPatchFieldReconstructor
-                            (
-                                mesh_.boundary()[curBPatch].size(),
-                                procField.boundaryField()[patchI].size()
-                            )
-                        )
-                    );
-                }
-
-                const label curPatchStart =
-                    mesh_.boundaryMesh()[curBPatch].start();
-
-                labelList reverseAddressing(cp.size());
-
-                forAll (cp, faceI)
-                {
-                    // Subtract one to take into account offsets for
-                    // face direction.
-                    reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart;
-                }
-
-                patchFields[curBPatch].rmap
-                (
-                    procField.boundaryField()[patchI],
-                    reverseAddressing
-                );
-            }
-            else
-            {
-                const Field<Type>& curProcPatch =
-                    procField.boundaryField()[patchI];
-
-                // In processor patches, there's a mix of internal faces (some
-                // of them turned) and possible cyclics. Slow loop
-                forAll (cp, faceI)
-                {
-                    // Subtract one to take into account offsets for
-                    // face direction.
-                    label curF = cp[faceI] - 1;
-
-                    // Is the face on the boundary?
-                    if (curF >= mesh_.nInternalFaces())
-                    {
-                        label curBPatch =
-                            mesh_.boundaryMesh().whichPatch(curF);
-
-                        if (!patchFields(curBPatch))
-                        {
-                            patchFields.set
-                            (
-                                curBPatch,
-                                fvPatchField<Type>::New
-                                (
-                                    mesh_.boundary()[curBPatch].type(),
-                                    mesh_.boundary()[curBPatch],
-                                    DimensionedField<Type, volMesh>::null()
-                                )
-                            );
-                        }
-
-                        // add the face
-                        label curPatchFace =
-                            mesh_.boundaryMesh()
-                                [curBPatch].whichFace(curF);
-
-                        patchFields[curBPatch][curPatchFace] =
-                            curProcPatch[faceI];
-                    }
-                }
-            }
-        }
-    }
-
-    forAll (mesh_.boundary(), patchI)
-    {
-        // add empty patches
-        if
-        (
-            isType<emptyFvPatch>(mesh_.boundary()[patchI])
-         && !patchFields(patchI)
-        )
-        {
-            patchFields.set
-            (
-                patchI,
-                fvPatchField<Type>::New
-                (
-                    emptyFvPatchField<Type>::typeName,
-                    mesh_.boundary()[patchI],
-                    DimensionedField<Type, volMesh>::null()
-                )
-            );
-        }
-    }
-
-
-    // Now construct and write the field
-    // setting the internalField and patchFields
-    return tmp<GeometricField<Type, fvPatchField, volMesh> >
-    (
-        new GeometricField<Type, fvPatchField, volMesh>
-        (
-            IOobject
-            (
-                fieldIoObject.name(),
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            procFields[0].dimensions(),
-            internalField,
-            patchFields
-        )
-    );
-}
-
-
-template<class Type>
-Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
-Foam::fvFieldReconstructor::reconstructFvSurfaceField
-(
-    const IOobject& fieldIoObject
-)
-{
-    // Read the field for all the processors
-    PtrList<GeometricField<Type, fvsPatchField, surfaceMesh> > procFields
-    (
-        procMeshes_.size()
-    );
-
-    forAll (procMeshes_, procI)
-    {
-        procFields.set
-        (
-            procI,
-            new GeometricField<Type, fvsPatchField, surfaceMesh>
-            (
-                IOobject
-                (
-                    fieldIoObject.name(),
-                    procMeshes_[procI].time().timeName(),
-                    procMeshes_[procI],
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                ),
-                procMeshes_[procI]
-            )
-        );
-    }
-
-
-    // Create the internalField
-    Field<Type> internalField(mesh_.nInternalFaces());
-
-    // Create the patch fields
-    PtrList<fvsPatchField<Type> > patchFields(mesh_.boundary().size());
-
-    forAll (procMeshes_, procI)
-    {
-        const GeometricField<Type, fvsPatchField, surfaceMesh>& procField =
-            procFields[procI];
-
-        // Set the face values in the reconstructed field
-
-        // It is necessary to create a copy of the addressing array to
-        // take care of the face direction offset trick.
-        //
-        {
-            labelList curAddr(faceProcAddressing_[procI]);
-
-            forAll (curAddr, addrI)
-            {
-                curAddr[addrI] -= 1;
-            }
-
-            internalField.rmap
-            (
-                procField.internalField(),
-                curAddr
-            );
-        }
-
-        // Set the boundary patch values in the reconstructed field
-        forAll (boundaryProcAddressing_[procI], patchI)
-        {
-            // Get patch index of the original patch
-            const label curBPatch = boundaryProcAddressing_[procI][patchI];
-
-            // Get addressing slice for this patch
-            const labelList::subList cp =
-                procMeshes_[procI].boundary()[patchI].patchSlice
-                (
-                    faceProcAddressing_[procI]
-                );
-
-            // Check if the boundary patch is not a processor patch
-            if (curBPatch >= 0)
-            {
-                // Regular patch. Fast looping
-                if (!patchFields(curBPatch))
-                {
-                    patchFields.set
-                    (
-                        curBPatch,
-                        fvsPatchField<Type>::New
-                        (
-                            procField.boundaryField()[patchI],
-                            mesh_.boundary()[curBPatch],
-                            DimensionedField<Type, surfaceMesh>::null(),
-                            fvPatchFieldReconstructor
-                            (
-                                mesh_.boundary()[curBPatch].size(),
-                                procField.boundaryField()[patchI].size()
-                            )
-                        )
-                    );
-                }
-
-                const label curPatchStart =
-                    mesh_.boundaryMesh()[curBPatch].start();
-
-                labelList reverseAddressing(cp.size());
-
-                forAll (cp, faceI)
-                {
-                    // Subtract one to take into account offsets for
-                    // face direction.
-                    reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart;
-                }
-
-                patchFields[curBPatch].rmap
-                (
-                    procField.boundaryField()[patchI],
-                    reverseAddressing
-                );
-            }
-            else
-            {
-                const Field<Type>& curProcPatch =
-                    procField.boundaryField()[patchI];
-
-                // In processor patches, there's a mix of internal faces (some
-                // of them turned) and possible cyclics. Slow loop
-                forAll (cp, faceI)
-                {
-                    label curF = cp[faceI] - 1;
-
-                    // Is the face turned the right side round
-                    if (curF >= 0)
-                    {
-                        // Is the face on the boundary?
-                        if (curF >= mesh_.nInternalFaces())
-                        {
-                            label curBPatch =
-                                mesh_.boundaryMesh().whichPatch(curF);
-
-                            if (!patchFields(curBPatch))
-                            {
-                                patchFields.set
-                                (
-                                    curBPatch,
-                                    fvsPatchField<Type>::New
-                                    (
-                                        mesh_.boundary()[curBPatch].type(),
-                                        mesh_.boundary()[curBPatch],
-                                        DimensionedField<Type, surfaceMesh>
-                                           ::null()
-                                    )
-                                );
-                            }
-
-                            // add the face
-                            label curPatchFace =
-                                mesh_.boundaryMesh()
-                                [curBPatch].whichFace(curF);
-
-                            patchFields[curBPatch][curPatchFace] =
-                                curProcPatch[faceI];
-                        }
-                        else
-                        {
-                            // Internal face
-                            internalField[curF] = curProcPatch[faceI];
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-    forAll (mesh_.boundary(), patchI)
-    {
-        // add empty patches
-        if
-        (
-            isType<emptyFvPatch>(mesh_.boundary()[patchI])
-         && !patchFields(patchI)
-        )
-        {
-            patchFields.set
-            (
-                patchI,
-                fvsPatchField<Type>::New
-                (
-                    emptyFvsPatchField<Type>::typeName,
-                    mesh_.boundary()[patchI],
-                    DimensionedField<Type, surfaceMesh>::null()
-                )
-            );
-        }
-    }
-
-
-    // Now construct and write the field
-    // setting the internalField and patchFields
-    return tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
-    (
-        new GeometricField<Type, fvsPatchField, surfaceMesh>
-        (
-            IOobject
-            (
-                fieldIoObject.name(),
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            procFields[0].dimensions(),
-            internalField,
-            patchFields
-        )
-    );
-}
-
-
-// Reconstruct and write all/selected volume fields
-template<class Type>
-void Foam::fvFieldReconstructor::reconstructFvVolumeFields
-(
-    const IOobjectList& objects,
-    const HashSet<word>& selectedFields
-)
-{
-    const word& fieldClassName =
-        GeometricField<Type, fvPatchField, volMesh>::typeName;
-
-    IOobjectList fields = objects.lookupClass(fieldClassName);
-
-    if (fields.size())
-    {
-        Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
-
-        forAllConstIter(IOobjectList, fields, fieldIter)
-        {
-            if
-            (
-                !selectedFields.size()
-             || selectedFields.found(fieldIter()->name())
-            )
-            {
-                Info<< "        " << fieldIter()->name() << endl;
-
-                reconstructFvVolumeField<Type>(*fieldIter())().write();
-            }
-        }
-        Info<< endl;
-    }
-}
-
-
-// Reconstruct and write all/selected surface fields
-template<class Type>
-void Foam::fvFieldReconstructor::reconstructFvSurfaceFields
-(
-    const IOobjectList& objects,
-    const HashSet<word>& selectedFields
-)
-{
-    const word& fieldClassName =
-        GeometricField<Type, fvsPatchField, surfaceMesh>::typeName;
-
-    IOobjectList fields = objects.lookupClass(fieldClassName);
-
-    if (fields.size())
-    {
-        Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
-
-        forAllConstIter(IOobjectList, fields, fieldIter)
-        {
-            if
-            (
-                !selectedFields.size()
-             || selectedFields.found(fieldIter()->name())
-            )
-            {
-                Info<< "        " << fieldIter()->name() << endl;
-
-                reconstructFvSurfaceField<Type>(*fieldIter())().write();
-            }
-        }
-
-        Info<< endl;
-    }
-}
-
-
-// ************************************************************************* //