Merge remote-tracking branch 'upstream/develop' into wp3-directional-refinement

.gitmodules

@@ -0,0 +1,3 @@
+[submodule "cfmesh"]
+	path = modules/cfmesh
+	url = https://develop.openfoam.com/Community/integration-cfmesh.git

Allwmake

@@ -43,6 +43,15 @@ echo "Compile OpenFOAM applications"
 echo
 applications/Allwmake $targetType $*
 
+# Additional components/modules
+if [ -d "$WM_PROJECT_DIR/modules" ]
+then
+    echo "========================================"
+    echo "Compile OpenFOAM modules"
+    echo
+    (cd $WM_PROJECT_DIR/modules 2>/dev/null && wmake -all)
+fi
+
 # Some summary information
 echo
 date "+%Y-%m-%d %H:%M:%S %z" 2>/dev/null || echo "date is unknown"

applications/solvers/basic/potentialFoam/overPotentialFoam/Make/files

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

applications/solvers/basic/potentialFoam/overPotentialFoam/Make/options

@@ -0,0 +1,12 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/overset/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -loverset

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

@@ -0,0 +1,9 @@
+const dictionary& potentialFlow
+(
+    mesh.solutionDict().subDict("potentialFlow")
+);
+
+const int nNonOrthCorr
+(
+    potentialFlow.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0)
+);

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

@@ -0,0 +1,146 @@
+Info<< "Reading velocity field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Initialise the velocity internal field to zero
+U = dimensionedVector("0", U.dimensions(), Zero);
+
+surfaceScalarField phi
+(
+    IOobject
+    (
+        "phi",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    fvc::flux(U)
+);
+
+if (args.optionFound("initialiseUBCs"))
+{
+    U.correctBoundaryConditions();
+    phi = fvc::flux(U);
+}
+
+
+// Construct a pressure field
+// If it is available read it otherwise construct from the velocity BCs
+// converting fixed-value BCs to zero-gradient and vice versa.
+word pName("p");
+
+// Update name of the pressure field from the command-line option
+args.optionReadIfPresent("pName", pName);
+
+// Infer the pressure BCs from the velocity
+wordList pBCTypes
+(
+    U.boundaryField().size(),
+    fixedValueFvPatchScalarField::typeName
+);
+
+forAll(U.boundaryField(), patchi)
+{
+    if (U.boundaryField()[patchi].fixesValue())
+    {
+        pBCTypes[patchi] = zeroGradientFvPatchScalarField::typeName;
+    }
+}
+
+// Note that registerObject is false for the pressure field. The pressure
+// field in this solver doesn't have a physical value during the solution.
+// It shouldn't be looked up and used by sub models or boundary conditions.
+Info<< "Constructing pressure field " << pName << nl << endl;
+volScalarField p
+(
+    IOobject
+    (
+        pName,
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::NO_WRITE,
+        false
+    ),
+    mesh,
+    dimensionedScalar(pName, sqr(dimVelocity), 0),
+    pBCTypes
+);
+
+// Infer the velocity potential BCs from the pressure
+wordList PhiBCTypes
+(
+    p.boundaryField().size(),
+    zeroGradientFvPatchScalarField::typeName
+);
+
+forAll(p.boundaryField(), patchi)
+{
+    if (p.boundaryField()[patchi].fixesValue())
+    {
+        PhiBCTypes[patchi] = fixedValueFvPatchScalarField::typeName;
+    }
+}
+
+Info<< "Constructing velocity potential field Phi\n" << endl;
+volScalarField Phi
+(
+    IOobject
+    (
+        "Phi",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Phi", dimLength*dimVelocity, 0),
+    PhiBCTypes
+);
+
+label PhiRefCell = 0;
+scalar PhiRefValue = 0;
+setRefCell
+(
+    Phi,
+    potentialFlow.dict(),
+    PhiRefCell,
+    PhiRefValue
+);
+mesh.setFluxRequired(Phi.name());
+
+#include "createMRF.H"
+
+// Add overset specific interpolations
+{
+    dictionary oversetDict;
+    oversetDict.add("Phi", true);
+    oversetDict.add("U", true);
+
+    const_cast<dictionary&>
+    (
+        mesh.schemesDict()
+    ).add
+    (
+        "oversetInterpolationRequired",
+        oversetDict,
+        true
+    );
+}
+
+// Mask field for zeroing out contributions on hole cells
+#include "createCellMask.H"
+
+// Create bool field with interpolated cells
+#include "createInterpolatedCells.H"

applications/solvers/basic/potentialFoam/overPotentialFoam/potentialFoam.C

@@ -0,0 +1,260 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    potentialFoam
+
+Group
+    grpBasicSolvers
+
+Description
+    Potential flow solver which solves for the velocity potential, to
+    calculate the flux-field, from which the velocity field is obtained by
+    reconstructing the flux.
+
+    \heading Solver details
+    The potential flow solution is typically employed to generate initial fields
+    for full Navier-Stokes codes.  The flow is evolved using the equation:
+
+    \f[
+        \laplacian \Phi = \div(\vec{U})
+    \f]
+
+    Where:
+    \vartable
+        \Phi      | Velocity potential [m2/s]
+        \vec{U}   | Velocity [m/s]
+    \endvartable
+
+    The corresponding pressure field could be calculated from the divergence
+    of the Euler equation:
+
+    \f[
+        \laplacian p + \div(\div(\vec{U}\otimes\vec{U})) = 0
+    \f]
+
+    but this generates excessive pressure variation in regions of large
+    velocity gradient normal to the flow direction.  A better option is to
+    calculate the pressure field corresponding to velocity variation along the
+    stream-lines:
+
+    \f[
+        \laplacian p + \div(\vec{F}\cdot\div(\vec{U}\otimes\vec{U})) = 0
+    \f]
+    where the flow direction tensor \f$\vec{F}\f$ is obtained from
+    \f[
+        \vec{F} = \hat{\vec{U}}\otimes\hat{\vec{U}}
+    \f]
+
+    \heading Required fields
+    \plaintable
+        U         | Velocity [m/s]
+    \endplaintable
+
+    \heading Optional fields
+    \plaintable
+        p         | Kinematic pressure [m2/s2]
+        Phi       | Velocity potential [m2/s]
+                  | Generated from p (if present) or U if not present
+    \endplaintable
+
+    \heading Options
+    \plaintable
+        -writep   | write the Euler pressure
+        -writePhi | Write the final velocity potential
+        -initialiseUBCs | Update the velocity boundaries before solving for Phi
+    \endplaintable
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "pisoControl.H"
+#include "dynamicFvMesh.H"
+#include "cellCellStencilObject.H"
+#include "localMin.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addOption
+    (
+        "pName",
+        "pName",
+        "Name of the pressure field"
+    );
+
+    argList::addBoolOption
+    (
+        "initialiseUBCs",
+        "Initialise U boundary conditions"
+    );
+
+    argList::addBoolOption
+    (
+        "writePhi",
+        "Write the final velocity potential field"
+    );
+
+    argList::addBoolOption
+    (
+        "writep",
+        "Calculate and write the Euler pressure field"
+    );
+
+    argList::addBoolOption
+    (
+        "withFunctionObjects",
+        "execute functionObjects"
+    );
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createNamedDynamicFvMesh.H"
+
+    pisoControl potentialFlow(mesh, "potentialFlow");
+
+    #include "createFields.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< nl << "Calculating potential flow" << endl;
+
+    mesh.update();
+
+    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
+
+    // Since solver contains no time loop it would never execute
+    // function objects so do it ourselves
+    runTime.functionObjects().start();
+
+    MRF.makeRelative(phi);
+    adjustPhi(phi, U, p);
+
+    // Non-orthogonal velocity potential corrector loop
+    while (potentialFlow.correct())
+    {
+        phi = fvc::flux(U);
+
+        while (potentialFlow.correctNonOrthogonal())
+        {
+            fvScalarMatrix PhiEqn
+            (
+                fvm::laplacian(faceMask, Phi)
+            ==
+                fvc::div(phi)
+            );
+
+            PhiEqn.setReference(PhiRefCell, PhiRefValue);
+            PhiEqn.solve();
+
+            if (potentialFlow.finalNonOrthogonalIter())
+            {
+                phi -= PhiEqn.flux();
+            }
+        }
+
+        MRF.makeAbsolute(phi);
+
+        Info<< "Continuity error = "
+            << mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
+            << endl;
+
+        U = fvc::reconstruct(phi);
+        U.correctBoundaryConditions();
+
+        Info<< "Interpolated velocity error = "
+            << (sqrt(sum(sqr(fvc::flux(U) - phi)))/sum(mesh.magSf())).value()
+            << endl;
+    }
+
+    // Write U and phi
+    U.write();
+    phi.write();
+
+    // Optionally write Phi
+    if (args.optionFound("writePhi"))
+    {
+        Phi.write();
+    }
+
+    // Calculate the pressure field from the Euler equation
+    if (args.optionFound("writep"))
+    {
+        Info<< nl << "Calculating approximate pressure field" << endl;
+
+        label pRefCell = 0;
+        scalar pRefValue = 0.0;
+        setRefCell
+        (
+            p,
+            potentialFlow.dict(),
+            pRefCell,
+            pRefValue
+        );
+
+        // Calculate the flow-direction filter tensor
+        volScalarField magSqrU(magSqr(U));
+        volSymmTensorField F(sqr(U)/(magSqrU + SMALL*average(magSqrU)));
+
+        // Calculate the divergence of the flow-direction filtered div(U*U)
+        // Filtering with the flow-direction generates a more reasonable
+        // pressure distribution in regions of high velocity gradient in the
+        // direction of the flow
+        volScalarField divDivUU
+        (
+            fvc::div
+            (
+                F & fvc::div(phi, U),
+                "div(div(phi,U))"
+            )
+        );
+
+        // Solve a Poisson equation for the approximate pressure
+        while (potentialFlow.correctNonOrthogonal())
+        {
+            fvScalarMatrix pEqn
+            (
+                fvm::laplacian(p) + divDivUU
+            );
+
+            pEqn.setReference(pRefCell, pRefValue);
+            pEqn.solve();
+        }
+
+        p.write();
+    }
+
+    runTime.functionObjects().end();
+
+    Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+        << nl << endl;
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

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

@@ -7,8 +7,8 @@ surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU));
-volVectorField HbyA("HbyA", constrainHbyA(rAU*UEqn.H(), U, p));
-//mesh.interpolate(HbyA);
+volVectorField HbyA("HbyA", U);
+HbyA = constrainHbyA(rAU*UEqn.H(), U, p);
 
 if (pimple.nCorrPISO() <= 1)
 {

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/Make/files

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

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/Make/options

@@ -0,0 +1,27 @@
+EXE_INC = \
+    -I.. \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/overset/lnInclude
+
+EXE_LIBS = \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lsampling \
+    -lmeshTools \
+    -lfvOptions \
+    -loverset \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh

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

@@ -0,0 +1,23 @@
+    // Solve the Momentum equation
+    MRF.correctBoundaryVelocity(U);
+
+    tmp<fvVectorMatrix> tUEqn
+    (
+        fvm::div(phi, U)
+      + MRF.DDt(rho, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        fvOptions(rho, U)
+    );
+    fvVectorMatrix& UEqn = tUEqn.ref();
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (simple.momentumPredictor())
+    {
+        solve(UEqn == -cellMask*fvc::grad(p));
+    }
+
+    fvOptions.correct(U);

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

@@ -0,0 +1 @@
+const volScalarField& psi = thermo.psi();

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

@@ -0,0 +1,91 @@
+Info<< "Reading thermophysical properties\n" << endl;
+
+autoPtr<fluidThermo> pThermo
+(
+    fluidThermo::New(mesh)
+);
+fluidThermo& thermo = pThermo();
+thermo.validate(args.executable(), "h", "e");
+
+volScalarField& p = thermo.p();
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.rho()
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "compressibleCreatePhi.H"
+
+pressureControl pressureControl(p, rho, simple.dict());
+
+mesh.setFluxRequired(p.name());
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+#include "createMRF.H"
+
+//- Overset specific
+
+// Add solver-specific interpolations
+{
+    dictionary oversetDict;
+    oversetDict.add("U", true);
+    oversetDict.add("p", true);
+    oversetDict.add("HbyA", true);
+    oversetDict.add("grad(p)", true);
+    oversetDict.add("rho", true);
+
+    const_cast<dictionary&>
+    (
+        mesh.schemesDict()
+    ).add
+    (
+        "oversetInterpolationRequired",
+        oversetDict,
+        true
+    );
+}
+
+// Mask field for zeroing out contributions on hole cells
+#include "createCellMask.H"
+
+#include "createInterpolatedCells.H"
+
+bool adjustFringe
+(
+    simple.dict().lookupOrDefault("oversetAdjustPhi", false)
+);

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

@@ -0,0 +1,22 @@
+    Info<< "Create dynamic mesh for time = "
+        << runTime.timeName() << nl << endl;
+
+    autoPtr<dynamicFvMesh> meshPtr
+    (
+        dynamicFvMesh::New
+        (
+            IOobject
+            (
+                dynamicFvMesh::defaultRegion,
+                runTime.timeName(),
+                runTime,
+                IOobject::MUST_READ
+            )
+        )
+    );
+
+    dynamicFvMesh& mesh = meshPtr();
+
+    // Calculate initial mesh-to-mesh mapping. Note that this should be
+    // done under the hood, e.g. as a MeshObject
+    mesh.update();

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

@@ -0,0 +1,123 @@
+{
+    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
+
+    volScalarField rAU(1.0/UEqn.A());
+    surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU));
+    volVectorField HbyA("HbyA", U);
+    HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
+    tUEqn.clear();
+
+    bool closedVolume = false;
+
+    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
+    if (simple.transonic())
+    {
+        surfaceScalarField phid
+        (
+            "phid",
+            (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
+        );
+
+        phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+
+        while (simple.correctNonOrthogonal())
+        {
+            fvScalarMatrix pEqn
+            (
+                fvc::div(phiHbyA)
+              + fvm::div(phid, p)
+              - fvm::laplacian(rhorAUf, p)
+            ==
+                fvOptions(psi, p, rho.name())
+            );
+
+            // Relax the pressure equation to ensure diagonal-dominance
+            pEqn.relax();
+
+            pEqn.setReference
+            (
+                pressureControl.refCell(),
+                pressureControl.refValue()
+            );
+
+            pEqn.solve();
+
+            if (simple.finalNonOrthogonalIter())
+            {
+                phi = phiHbyA + pEqn.flux();
+            }
+        }
+    }
+    else
+    {
+        closedVolume = adjustPhi(phiHbyA, U, p);
+
+        if (adjustFringe)
+        {
+            oversetAdjustPhi(phiHbyA, U);
+        }
+
+        while (simple.correctNonOrthogonal())
+        {
+            fvScalarMatrix pEqn
+            (
+                fvc::div(phiHbyA)
+              - fvm::laplacian(rhorAUf, p)
+            ==
+                fvOptions(psi, p, rho.name())
+            );
+
+            pEqn.setReference
+            (
+                pressureControl.refCell(),
+                pressureControl.refValue()
+            );
+
+            pEqn.solve();
+
+            if (simple.finalNonOrthogonalIter())
+            {
+                phi = phiHbyA + pEqn.flux();
+            }
+        }
+    }
+
+    #include "incompressible/continuityErrs.H"
+
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    volVectorField gradP(fvc::grad(p));
+
+    U = HbyA - rAU*cellMask*gradP;
+    U.correctBoundaryConditions();
+    fvOptions.correct(U);
+
+
+    bool pLimited = pressureControl.limit(p);
+
+    // For closed-volume cases adjust the pressure and density levels
+    // to obey overall mass continuity
+    if (closedVolume)
+    {
+        p += (initialMass - fvc::domainIntegrate(psi*p))
+            /fvc::domainIntegrate(psi);
+    }
+
+    if (pLimited || closedVolume)
+    {
+        p.correctBoundaryConditions();
+    }
+
+    rho = thermo.rho();
+
+    if (!simple.transonic())
+    {
+        rho.relax();
+    }
+}

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/rhoSimpleFoam.C

@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    overRhoSimpleFoam
+
+Group
+    grpCompressibleSolvers
+
+Description
+    Overset steady-state solver for turbulent flow of compressible fluids.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "fluidThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "simpleControl.H"
+#include "pressureControl.H"
+#include "fvOptions.H"
+#include "cellCellStencilObject.H"
+#include "localMin.H"
+#include "oversetAdjustPhi.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #define CREATE_MESH createUpdatedDynamicFvMesh.H
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createUpdatedDynamicFvMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
+
+    turbulence->validate();
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // Pressure-velocity SIMPLE corrector
+        #include "UEqn.H"
+        #include "EEqn.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/doc/solver.dox

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,6 +39,7 @@ The available solvers are grouped into the following categories:
   - \ref grpLagrangianSolvers
   - \ref grpMultiphaseSolvers
   - \ref grpStressAnalysisSolvers
+  - \ref grpFiniteAreaSolvers
 
 
 \*---------------------------------------------------------------------------*/

applications/solvers/doc/solversDoc.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,4 +34,10 @@ License
     This group contains moving mesh solvers solvers
 @}
 
+\defgroup grpFiniteAreaSolvers Finite area solvers
+@{
+    \ingroup grpSolvers
+    This group contains finite area solvers
+@}
+
 \*---------------------------------------------------------------------------*/

applications/solvers/finiteArea/liquidFilmFoam/Make/files

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

applications/solvers/finiteArea/liquidFilmFoam/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/cfdTools/general/lnInclude
+
+EXE_LIBS = \
+    -lfiniteArea \
+    -lfiniteVolume \
+    -lmeshTools

applications/solvers/finiteArea/liquidFilmFoam/calcFrictionFactor.H

@@ -0,0 +1,7 @@
+{
+    // Stabilisation of friction factor calculation
+    // Friction factor is defined with standard gravity
+    frictionFactor.primitiveFieldRef() =
+        mag(2*9.81*sqr(manningField.primitiveField())/
+        pow(mag(h.primitiveField()) + 1e-7, 1.0/3.0));
+}

applications/solvers/finiteArea/liquidFilmFoam/capillaryCourantNo.H

@@ -0,0 +1,13 @@
+{
+    scalar CoNumSigma = max
+    (
+        sqrt
+        (
+            2*M_PI*sigma*sqr(aMesh.edgeInterpolation::deltaCoeffs())
+            *aMesh.edgeInterpolation::deltaCoeffs()
+            /rhol
+        )
+    ).value()*runTime.deltaT().value();
+
+    Info<< "Max Capillary Courant Number = " << CoNumSigma << '\n' << endl;
+}

applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H

@@ -0,0 +1,158 @@
+    Info<< "Reading field h" << endl;
+    areaScalarField h
+    (
+        IOobject
+        (
+            "h",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        aMesh
+    );
+
+
+    Info<< "Reading field Us" << endl;
+    areaVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        aMesh
+    );
+
+
+    edgeScalarField phis
+    (
+        IOobject
+        (
+            "phis",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        fac::interpolate(Us) & aMesh.Le()
+    );
+
+
+    edgeScalarField phi2s
+    (
+        IOobject
+        (
+            "phi2s",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        fac::interpolate(h*Us) & aMesh.Le()
+    );
+
+
+    const areaVectorField& Ns = aMesh.faceAreaNormals();
+    areaVectorField Gs(g - Ns*(Ns & g));
+    areaScalarField Gn(mag(g - Gs));
+
+    // Mass source
+    areaScalarField Sm
+    (
+        IOobject
+        (
+            "Sm",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        aMesh,
+        dimensionedScalar("Sm", dimLength/dimTime, 0)
+    );
+
+    // Mass sink
+    areaScalarField Sd
+    (
+        IOobject
+        (
+            "Sd",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        aMesh,
+        dimensionedScalar("Sd", dimLength/dimTime, 0)
+    );
+
+    areaVectorField Ug
+    (
+        IOobject
+        (
+            "Sg",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        aMesh,
+        dimensionedVector("Ug", dimVelocity, vector::zero)
+    );
+
+
+    // Surface pressure
+    areaScalarField ps
+    (
+        IOobject
+        (
+            "ps",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        rhol*Gn*h - sigma*fac::laplacian(h)
+    );
+
+    // Friction factor
+    areaScalarField dotProduct
+    (
+        aMesh.faceAreaNormals() & (g/mag(g))
+    );
+
+    Info<< "View factor: min = " << min(dotProduct.internalField())
+        << " max = " << max(dotProduct.internalField()) << endl;
+
+    areaScalarField manningField
+    (
+        IOobject
+        (
+            "manningField",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        aMesh
+    );
+
+    areaScalarField frictionFactor
+    (
+        IOobject
+        (
+            "frictionFactor",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        aMesh,
+        dimensionedScalar("one", dimless, 0.01)
+    );
+
+    aMesh.setFluxRequired("h");

applications/solvers/finiteArea/liquidFilmFoam/createFvFields.H

@@ -0,0 +1,31 @@
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedVector("0", dimVelocity, vector::zero)
+);
+
+
+volScalarField H
+(
+    IOobject
+    (
+        "H",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("0", dimLength, 0)
+);
+
+// Create volume-to surface mapping object
+volSurfaceMapping vsm(aMesh);

applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    liquidFilmFoam
+
+Group
+    grpFiniteAreaSolvers
+
+Description
+    Transient solver for incompressible, laminar flow of Newtonian fluids in
+    liquid film formulation.
+
+Author
+    Zeljko Tukovic, FMENA
+    Hrvoje Jasak, Wikki Ltd.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "faCFD.H"
+#include "loopControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFaMesh.H"
+    #include "readGravitationalAcceleration.H"
+    #include "readTransportProperties.H"
+    #include "createFaFields.H"
+    #include "createFvFields.H"
+    #include "createTimeControls.H"
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readSolutionControls.H"
+        #include "readTimeControls.H"
+        #include "surfaceCourantNo.H"
+        #include "capillaryCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        while (iters.loop())
+        {
+            phi2s = fac::interpolate(h)*phis;
+
+            #include "calcFrictionFactor.H"
+
+            faVectorMatrix UsEqn
+            (
+                fam::ddt(h, Us)
+              + fam::div(phi2s, Us)
+              + fam::Sp(0.0125*frictionFactor*mag(Us), Us)
+             ==
+                Gs*h
+              - fam::Sp(Sd, Us)
+            );
+
+            UsEqn.relax();
+            solve(UsEqn == - fac::grad(ps*h)/rhol + ps*fac::grad(h)/rhol);
+
+            areaScalarField UsA(UsEqn.A());
+
+            Us = UsEqn.H()/UsA;
+            Us.correctBoundaryConditions();
+
+            phis =
+                (fac::interpolate(Us) & aMesh.Le())
+              - fac::interpolate(1.0/(rhol*UsA))*fac::lnGrad(ps*h)*aMesh.magLe()
+              + fac::interpolate(ps/(rhol*UsA))*fac::lnGrad(h)*aMesh.magLe();
+
+            faScalarMatrix hEqn
+            (
+                fam::ddt(h)
+              + fam::div(phis, h)
+             ==
+                Sm
+              - fam::Sp
+                (
+                    Sd/(h + dimensionedScalar("small", dimLength, SMALL)),
+                    h
+                )
+            );
+
+            hEqn.relax();
+            hEqn.solve();
+
+            phi2s = hEqn.flux();
+
+            // Bound h
+            h.primitiveFieldRef() = max
+            (
+                max
+                (
+                    h.primitiveField(),
+                    fac::average(max(h, h0))().primitiveField()
+                   *pos(h0.value() - h.primitiveField())
+                ),
+                h0.value()
+            );
+
+            ps = rhol*Gn*h - sigma*fac::laplacian(h);
+            ps.correctBoundaryConditions();
+
+            Us -= (1.0/(rhol*UsA))*fac::grad(ps*h)
+                - (ps/(rhol*UsA))*fac::grad(h);
+            Us.correctBoundaryConditions();
+        }
+
+        if (runTime.outputTime())
+        {
+            vsm.mapToVolume(h, H.boundaryFieldRef());
+            vsm.mapToVolume(Us, U.boundaryFieldRef());
+
+            runTime.write();
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/finiteArea/liquidFilmFoam/readSolutionControls.H

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

applications/solvers/finiteArea/liquidFilmFoam/readTransportProperties.H

@@ -0,0 +1,41 @@
+IOdictionary transportProperties
+(
+    IOobject
+    (
+        "transportProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+dimensionedScalar mug
+(
+    transportProperties.lookup("mug")
+);
+
+dimensionedScalar mul
+(
+    transportProperties.lookup("mul")
+);
+
+dimensionedScalar sigma
+(
+    transportProperties.lookup("sigma")
+);
+
+dimensionedScalar rhol
+(
+    transportProperties.lookup("rhol")
+);
+
+dimensionedScalar rhog
+(
+    transportProperties.lookup("rhog")
+);
+
+dimensionedScalar h0
+(
+    transportProperties.lookup("h0")
+);

applications/solvers/finiteArea/liquidFilmFoam/surfaceCourantNo.H

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    surfaceCourantNo
+
+Author
+    Hrvoje Jasak, Wikki Ltd.
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers for the
+    Finite Area method.
+
+\*---------------------------------------------------------------------------*/
+
+scalar CoNum = 0.0;
+scalar meanCoNum = 0.0;
+scalar velMag = 0.0;
+
+if (aMesh.nInternalEdges())
+{
+    edgeScalarField SfUfbyDelta
+    (
+        aMesh.edgeInterpolation::deltaCoeffs()*mag(phis)
+    );
+
+    CoNum = max(SfUfbyDelta/aMesh.magLe())
+        .value()*runTime.deltaT().value();
+
+    meanCoNum = (sum(SfUfbyDelta)/sum(aMesh.magLe()))
+        .value()*runTime.deltaT().value();
+
+    velMag = max(mag(phis)/aMesh.magLe()).value();
+}
+
+Info<< "Courant Number mean: " << meanCoNum
+    << " max: " << CoNum
+    << " velocity magnitude: " << velMag << endl;
+
+
+// ************************************************************************* //

applications/solvers/finiteArea/sphereSurfactantFoam/Make/files

@@ -0,0 +1,3 @@
+surfactantFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/sphereSurfactantFoam

applications/solvers/finiteArea/sphereSurfactantFoam/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/cfdTools/general/lnInclude
+
+EXE_LIBS = \
+    -lfiniteArea \
+    -lfiniteVolume \
+    -lmeshTools

applications/solvers/finiteArea/sphereSurfactantFoam/createFaFields.H

@@ -0,0 +1,78 @@
+Info<< "Reading field Cs" << endl;
+areaScalarField Cs
+(
+    IOobject
+    (
+        "Cs",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    aMesh
+);
+
+dimensioned<scalar> Cs0
+(
+    "Cs0",
+    dimensionSet(1, -2, 0, 0, 0, 0, 0),
+    1.0
+);
+
+const areaVectorField& R = aMesh.areaCentres();
+
+Cs = Cs0*(1.0 + R.component(vector::X)/mag(R));
+
+
+dimensioned<scalar> Ds
+(
+    "Ds",
+    dimensionSet(0, 2, -1, 0, 0, 0, 0),
+    1.0
+);
+
+
+areaVectorField Us
+(
+    IOobject
+    (
+        "Us",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    aMesh,
+    dimensioned<vector>("Us", dimVelocity, vector::zero)
+);
+
+dimensioned<scalar> Uinf("Uinf", dimVelocity, 1.0);
+
+forAll (Us, faceI)
+{
+    Us[faceI].x() =
+        Uinf.value()*(0.25*(3.0 + sqr(R[faceI].x()/mag(R[faceI]))) - 1.0);
+
+    Us[faceI].y() =
+        Uinf.value()*0.25*R[faceI].x()*R[faceI].y()/sqr(mag(R[faceI]));
+
+    Us[faceI].z() =
+        Uinf.value()*0.25*R[faceI].x()*R[faceI].z()/sqr(mag(R[faceI]));
+}
+
+Us -= aMesh.faceAreaNormals()*(aMesh.faceAreaNormals() & Us);
+
+
+edgeScalarField phis
+(
+    IOobject
+    (
+        "phis",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    linearEdgeInterpolate(Us) & aMesh.Le()
+);
+

applications/solvers/finiteArea/sphereSurfactantFoam/createFaMesh.H

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

applications/solvers/finiteArea/sphereSurfactantFoam/createVolFields.H

@@ -0,0 +1,36 @@
+    // Create volume-to surface mapping object
+    volSurfaceMapping vsm(aMesh);
+
+    volScalarField Cvf
+    (
+        IOobject
+        (
+            "Cvf",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimless/dimLength, 0)
+    );
+
+    vsm.mapToVolume(Cs, Cvf.boundaryFieldRef());
+    Cvf.write();
+
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimVelocity, vector::zero)
+    );
+
+    vsm.mapToVolume(Us, U.boundaryFieldRef());
+    U.write();

applications/solvers/finiteArea/sphereSurfactantFoam/surfactantFoam.C

@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    surfactantFoam for sphere transport
+
+Group
+    grpFiniteAreaSolvers
+
+Description
+    Passive scalar transport equation solver on a sphere
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "faCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    #include "createFaMesh.H"
+    #include "createFaFields.H"
+    #include "createVolFields.H"
+
+    Info<< "Total mass of surfactant: "
+        << sum(Cs.internalField()*aMesh.S()) << endl;
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.value() << endl;
+
+        faScalarMatrix CsEqn
+        (
+            fam::ddt(Cs)
+          + fam::div(phis, Cs)
+          - fam::laplacian(Ds, Cs)
+        );
+
+        CsEqn.solve();
+
+        if (runTime.writeTime())
+        {
+            vsm.mapToVolume(Cs, Cvf.boundaryFieldRef());
+
+            runTime.write();
+        }
+
+        Info<< "Total mass of surfactant: "
+            << sum(Cs.internalField()*aMesh.S()) << endl;
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/finiteArea/surfactantFoam/Make/files

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

applications/solvers/finiteArea/surfactantFoam/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/cfdTools/general/lnInclude
+
+EXE_LIBS = \
+    -lfiniteArea \
+    -lfiniteVolume \
+    -lmeshTools

applications/solvers/finiteArea/surfactantFoam/createFaFields.H

@@ -0,0 +1,63 @@
+Info<< "Reading field Cs" << endl;
+areaScalarField Cs
+(
+    IOobject
+    (
+        "Cs",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    aMesh
+);
+
+Info<< "Reading transportProperties\n" << endl;
+
+IOdictionary transportProperties
+(
+    IOobject
+    (
+        "transportProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+
+Info<< "Reading diffusivity D\n" << endl;
+
+dimensionedScalar Ds
+(
+    transportProperties.lookup("Ds")
+);
+
+areaVectorField Us
+(
+    IOobject
+    (
+        "Us",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    ),
+    aMesh
+);
+
+
+edgeScalarField phis
+(
+    IOobject
+    (
+        "phis",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    linearEdgeInterpolate(Us) & aMesh.Le()
+);
+

applications/solvers/finiteArea/surfactantFoam/createFaMesh.H

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

applications/solvers/finiteArea/surfactantFoam/createVolFields.H

@@ -0,0 +1,36 @@
+    // Create volume-to surface mapping object
+    volSurfaceMapping vsm(aMesh);
+
+    volScalarField Cvf
+    (
+        IOobject
+        (
+            "Cvf",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimless/dimLength, 0)
+    );
+
+    vsm.mapToVolume(Cs, Cvf.boundaryFieldRef());
+    Cvf.write();
+
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimVelocity, vector::zero)
+    );
+
+    vsm.mapToVolume(Us, U.boundaryFieldRef());
+    U.write();

applications/solvers/finiteArea/surfactantFoam/surfactantFoam.C

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    surfactantFoam
+
+Group
+    grpFiniteAreaSolvers
+
+Description
+    Passive scalar transport equation solver.
+
+    \heading Solver details
+    The equation is given by:
+
+    \f[
+        \ddt{Cs} + \div \left(\vec{U} Cs\right) - \div \left(D_T \grad Cs \right)
+        = 0
+    \f]
+
+    Where:
+    \vartable
+        Cs      | Passive scalar
+        Ds      | Diffusion coefficient
+    \endvartable
+
+    \heading Required fields
+    \plaintable
+        Cs      | Passive scalar
+        Us      | Velocity [m/s]
+    \endplaintable
+
+Author
+    Zeljko Tukovic, FMENA
+    Hrvoje Jasak, Wikki Ltd.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "faCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFaMesh.H"
+    #include "createFaFields.H"
+    #include "createVolFields.H"
+
+    Info<< "Total mass of surfactant: "
+        << sum(Cs.internalField()*aMesh.S()) << endl;
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.value() << endl;
+
+        faScalarMatrix CsEqn
+        (
+            fam::ddt(Cs)
+          + fam::div(phis, Cs)
+          - fam::laplacian(Ds, Cs)
+        );
+
+        CsEqn.solve();
+
+        if (runTime.writeTime())
+        {
+            vsm.mapToVolume(Cs, Cvf.boundaryFieldRef());
+
+            runTime.write();
+        }
+
+        Info<< "Total mass of surfactant: "
+            << sum(Cs.internalField()*aMesh.S()) << endl;
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,6 +47,7 @@ Description
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "coordinateSystem.H"
+#include "loopControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -89,11 +90,10 @@ int main(int argc, char *argv[])
             }
         }
 
-
         // --- PIMPLE loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for (int oCorr=0; oCorr<nOuterCorr; ++oCorr)
         {
-            bool finalIter = oCorr == nOuterCorr-1;
+            const bool finalIter = (oCorr == nOuterCorr-1);
 
             forAll(fluidRegions, i)
             {
@@ -113,6 +113,35 @@ int main(int argc, char *argv[])
                 #include "solveSolid.H"
             }
 
+            // Additional loops for energy solution only
+            if (!oCorr && nOuterCorr > 1)
+            {
+                loopControl looping(runTime, pimple, "energyCoupling");
+
+                while (looping.loop())
+                {
+                    Info<< nl << looping << nl;
+
+                    forAll(fluidRegions, i)
+                    {
+                        Info<< "\nSolving for fluid region "
+                            << fluidRegions[i].name() << endl;
+                       #include "setRegionFluidFields.H"
+                       #include "readFluidMultiRegionPIMPLEControls.H"
+                       frozenFlow = true;
+                       #include "solveFluid.H"
+                    }
+
+                    forAll(solidRegions, i)
+                    {
+                        Info<< "\nSolving for solid region "
+                            << solidRegions[i].name() << endl;
+                        #include "setRegionSolidFields.H"
+                        #include "readSolidMultiRegionPIMPLEControls.H"
+                        #include "solveSolid.H"
+                    }
+                }
+            }
         }
 
         runTime.write();

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,6 +42,7 @@ Description
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "coordinateSystem.H"
+#include "loopControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -57,7 +58,6 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
-
     while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
@@ -80,6 +80,35 @@ int main(int argc, char *argv[])
             #include "solveSolid.H"
         }
 
+        // Additional loops for energy solution only
+        {
+            loopControl looping(runTime, "SIMPLE", "energyCoupling");
+
+            while (looping.loop())
+            {
+                Info<< nl << looping << nl;
+
+                forAll(fluidRegions, i)
+                {
+                    Info<< "\nSolving for fluid region "
+                        << fluidRegions[i].name() << endl;
+                   #include "setRegionFluidFields.H"
+                   #include "readFluidMultiRegionSIMPLEControls.H"
+                   frozenFlow = true;
+                   #include "solveFluid.H"
+                }
+
+                forAll(solidRegions, i)
+                {
+                    Info<< "\nSolving for solid region "
+                        << solidRegions[i].name() << endl;
+                    #include "setRegionSolidFields.H"
+                    #include "readSolidMultiRegionSIMPLEControls.H"
+                    #include "solveSolid.H"
+                }
+            }
+        }
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

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

@@ -12,6 +12,14 @@
     volScalarField& p = thermo.p();
     const volScalarField& psi = thermo.psi();
 
+    volScalarField& p_rgh = p_rghFluid[i];
+
+    const dimensionedVector& g = gFluid[i];
+    const volScalarField& gh = ghFluid[i];
+    const surfaceScalarField& ghf = ghfFluid[i];
+
+    radiation::radiationModel& rad = radiation[i];
+
     IOMRFZoneList& MRF = MRFfluid[i];
     fv::options& fvOptions = fluidFvOptions[i];
 
@@ -22,14 +30,7 @@
         initialMassFluid[i]
     );
 
-    radiation::radiationModel& rad = radiation[i];
+    bool frozenFlow = frozenFlowFluid[i];
 
     const label pRefCell = pRefCellFluid[i];
     const scalar pRefValue = pRefValueFluid[i];
-    const bool frozenFlow = frozenFlowFluid[i];
-
-    volScalarField& p_rgh = p_rghFluid[i];
-
-    const dimensionedVector& g = gFluid[i];
-    const volScalarField& gh = ghFluid[i];
-    const surfaceScalarField& ghf = ghfFluid[i];

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

@@ -1,5 +1,5 @@
 {
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; ++nonOrth)
     {
         fvScalarMatrix hEqn
         (
@@ -20,9 +20,9 @@
 
         fvOptions.correct(h);
     }
-}
 
     thermo.correct();
 
     Info<< "Min/max T:" << min(thermo.T()).value() << ' '
         << max(thermo.T()).value() << endl;
+}

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

@@ -19,8 +19,8 @@ List<bool> frozenFlowFluid(fluidRegions.size(), false);
 PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
 PtrList<fv::options> fluidFvOptions(fluidRegions.size());
 
-List<label> refCellFluid(fluidRegions.size());
-List<scalar> refValueFluid(fluidRegions.size());
+List<label> pRefCellFluid(fluidRegions.size());
+List<scalar> pRefValueFluid(fluidRegions.size());
 
 // Populate fluid field pointer lists
 forAll(fluidRegions, i)
@@ -252,8 +252,8 @@ forAll(fluidRegions, i)
 
     turbulence[i].validate();
 
-    refCellFluid[i] = 0;
-    refValueFluid[i] = 0.0;
+    pRefCellFluid[i] = 0;
+    pRefValueFluid[i] = 0.0;
 
     if (p_rghFluid[i].needReference())
     {
@@ -262,8 +262,8 @@ forAll(fluidRegions, i)
             thermoFluid[i].p(),
             p_rghFluid[i],
             pimpleDict,
-            refCellFluid[i],
-            refValueFluid[i]
+            pRefCellFluid[i],
+            pRefValueFluid[i]
         );
     }
 

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

@@ -32,7 +32,8 @@
         initialMassFluid[i]
     );
 
-    const bool frozenFlow = frozenFlowFluid[i];
+    bool frozenFlow = frozenFlowFluid[i];
+
+    const label pRefCell = pRefCellFluid[i];
+    const scalar pRefValue = pRefValueFluid[i];
 
-    const label pRefCell = refCellFluid[i];
-    const scalar pRefValue = refValueFluid[i];

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

@@ -1,10 +1,10 @@
-    const wordList solidsNames(rp["solid"]);
+    const wordList solidNames(rp["solid"]);
 
-    PtrList<fvMesh> solidRegions(solidsNames.size());
+    PtrList<fvMesh> solidRegions(solidNames.size());
 
-    forAll(solidsNames, i)
+    forAll(solidNames, i)
     {
-        Info<< "Create solid mesh for region " << solidsNames[i]
+        Info<< "Create solid mesh for region " << solidNames[i]
             << " for time = " << runTime.timeName() << nl << endl;
 
         solidRegions.set
@@ -14,7 +14,7 @@
             (
                 IOobject
                 (
-                    solidsNames[i],
+                    solidNames[i],
                     runTime.timeName(),
                     runTime,
                     IOobject::MUST_READ

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

@@ -4,7 +4,7 @@ if (finalIter)
 }
 
 {
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; ++nonOrth)
     {
         fvScalarMatrix hEqn
         (
@@ -26,12 +26,12 @@ if (finalIter)
 
         fvOptions.correct(h);
     }
-}
 
     thermo.correct();
 
     Info<< "Min/max T:" << min(thermo.T()).value() << ' '
         << max(thermo.T()).value() << endl;
+}
 
 if (finalIter)
 {

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pimpleDyMFoam.C

@@ -54,10 +54,8 @@ Description
 
 int main(int argc, char *argv[])
 {
-    argList::addNote
-    (
-        "Experimental version of pimpleDyMFoam with support for overset meshes"
-    );
+    #include "postProcess.H"
+
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"

applications/solvers/incompressible/simpleFoam/overSimpleFoam/UEqn.H

@@ -16,6 +16,9 @@
 
     fvOptions.constrain(UEqn);
 
+    if (simple.momentumPredictor())
+    {
         solve(UEqn == -cellMask*fvc::grad(p));
+    }
 
     fvOptions.correct(U);

applications/solvers/incompressible/simpleFoam/overSimpleFoam/pEqn.H

@@ -5,7 +5,7 @@
     surfaceScalarField rAUf("rAUf", faceMask*fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
-    HbyA = constrainHbyA(rAU*UEqn.H(), U, p);
+    HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
 
     //mesh.interpolate(HbyA);
     if (massFluxInterpolation)

applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H

@@ -138,3 +138,4 @@ volScalarField pDivU
     mesh,
     dimensionedScalar("pDivU", p.dimensions()/dimTime, 0)
 );
+

applications/solvers/multiphase/interFoam/overInterDyMFoam/readControls.H

@@ -1,6 +1,6 @@
 #include "readTimeControls.H"
 
-correctPhi = pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
+correctPhi = pimple.dict().lookupOrDefault<Switch>("correctPhi", false);
 
 checkMeshCourantNo =
     pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);

applications/test/argList/Test-argList.C

@@ -58,6 +58,9 @@ int main(int argc, char *argv[])
 
     argList args(argc, argv);
 
+    Info<<"have: "
+        <<args.optionCount({"label", "scalar"}) << " options" << nl;
+
     label ival;
     scalar sval;
 

applications/test/quaternion/Test-quaternion.C

@@ -29,9 +29,13 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
+#include "argList.H"
 #include "quaternion.H"
 #include "septernion.H"
+#include "mathematicalConstants.H"
+#include "Tuple2.H"
 #include "IOstreams.H"
+#include "transform.H"
 
 using namespace Foam;
 
@@ -39,6 +43,95 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
+    argList::addOption
+    (
+        "rollPitchYaw",
+        "vector",
+        "Rotate by '(roll pitch yaw)' in degrees"
+    );
+    argList::addOption
+    (
+        "yawPitchRoll",
+        "vector",
+        "Rotate by '(yaw pitch roll)' in degrees"
+    );
+
+    argList::addOption
+    (
+        "rotate-angle",
+        "(vector angle)",
+        "Rotate about the <vector> by <angle> degrees - eg, '((1 0 0) 45)'"
+    );
+    argList args(argc, argv);
+
+    vector rotVector;
+
+    if (args.optionReadIfPresent("rollPitchYaw", rotVector))
+    {
+        Info<< "Rotate by" << nl
+            << "    roll  " << rotVector.x() << nl
+            << "    pitch " << rotVector.y() << nl
+            << "    yaw   " << rotVector.z() << nl;
+
+        // degToRad
+        rotVector *= constant::mathematical::pi/180.0;
+
+        const quaternion quat(quaternion::rotationSequence::XYZ, rotVector);
+
+        Info<< "quaternion " << quat << endl;
+        Info<< "rotation   = " << quat.R() << endl;
+    }
+    if (args.optionReadIfPresent("yawPitchRoll", rotVector))
+    {
+        Info<< "Rotate by" << nl
+            << "    yaw   " << rotVector.x() << nl
+            << "    pitch " << rotVector.y() << nl
+            << "    roll  " << rotVector.z() << nl;
+
+        // degToRad
+        rotVector *= constant::mathematical::pi/180.0;
+
+        const quaternion quat(quaternion::rotationSequence::ZYX, rotVector);
+
+        Info<< "quaternion " << quat << endl;
+        Info<< "rotation   = " << quat.R() << endl;
+    }
+    if (args.optionFound("rotate-angle"))
+    {
+        const Tuple2<vector, scalar> axisAngle
+        (
+            args.optionLookup("rotate-angle")()
+        );
+
+        Info<< "Rotate" << nl
+            << "    about " << axisAngle.first() << nl
+            << "    angle " << axisAngle.second() << nl;
+
+        const quaternion quat
+        (
+            axisAngle.first(),
+            axisAngle.second() * constant::mathematical::pi/180.0  // degToRad
+        );
+
+        Info<< "quaternion " << quat << endl;
+        Info<< "rotation   = " << quat.R() << endl;
+
+        Info<< "transform Ra = "
+            << Ra
+               (
+                   axisAngle.first() / mag(axisAngle.first()),
+                   axisAngle.second() * constant::mathematical::pi/180.0
+               ) << endl;
+        Info<< "-ve Ra = "
+            << Ra
+               (
+                   axisAngle.first() / mag(axisAngle.first()),
+                   -axisAngle.second() * constant::mathematical::pi/180.0
+               ) << endl;
+    }
+
+    Info<< nl << nl;
+
     quaternion q(vector(1, 2, 3), 0.7853981);
     Info<< "q " << q << endl;
 

applications/test/rigidBodyDynamics/pendulum/pendulum

@@ -32,6 +32,7 @@ bodies
         type            sphere;
         mass            1;
         radius          0.05;
+        centreOfMass    (0 10 0);
         transform       (1 0 0 0 1 0 0 0 1) (0 -1 0);
         mergeWith       hinge;
     }

applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring

@@ -14,6 +14,7 @@ bodies
         parent          root;
         mass            6e4;
         radius          0.01;
+        centreOfMass    (0 0 0);
         transform       (1 0 0 0 1 0 0 0 1) (0 0 0);
         joint
         {
@@ -37,6 +38,7 @@ bodies
         parent          M;
         mass            6e3;
         radius          0.01;
+        centreOfMass    (0 0 0);
         transform       (1 0 0 0 1 0 0 0 1) (0 -5 0);
         mergeWith       M;
     }

applications/utilities/finiteArea/checkFaMesh/Make/files

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

applications/utilities/finiteArea/checkFaMesh/Make/options

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfiniteArea

applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C

@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    makeFaMesh
+
+Description
+    Check a Finite Area mesh
+
+Author
+    Zeljko Tukovic, FAMENA
+    Hrvoje Jasak, Wikki Ltd.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "faCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "addRegionOption.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createNamedMesh.H"
+    #include "createFaMesh.H"
+
+    Info<< "Time = " << runTime.timeName() << nl << endl;
+
+    // General mesh statistics
+    Info<< "Number of points: " << aMesh.nPoints() << nl
+        << "Number of internal edges: " << aMesh.nInternalEdges() << nl
+        << "Number of edges: " << aMesh.nEdges() << nl
+        << "Number of faces: " << aMesh.nFaces() << nl
+        << endl;
+
+    // Check geometry
+    Info<< "Face area:  min = " << min(aMesh.S().field())
+        << " max = "  << max(aMesh.S().field()) << nl
+        << "Internal edge length: min = "
+        << min(aMesh.magLe().internalField()) << nl
+        << " max = "  << max(aMesh.magLe().internalField()) << nl
+        << "Edge length: min = "
+        << min(aMesh.magLe()).value() << nl
+        << " max = "  << max(aMesh.magLe()).value() << nl
+        << "Face area normals:  min = " << min(aMesh.faceAreaNormals().field())
+        << " max = "  << max(aMesh.faceAreaNormals().field()) << nl
+        << endl;
+
+
+    Info << "\nEnd" << endl;
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/finiteArea/makeFaMesh/Make/files

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

applications/utilities/finiteArea/makeFaMesh/Make/options

@@ -0,0 +1,8 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/cfdTools/general/lnInclude
+
+EXE_LIBS = \
+    -lfiniteArea \
+    -lfiniteVolume

applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C

@@ -0,0 +1,355 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    makeFaMesh
+
+Description
+    A mesh generator for finite area mesh.
+
+Author
+    Zeljko Tukovic, FAMENA
+    Hrvoje Jasak, Wikki Ltd.
+
+\*---------------------------------------------------------------------------*/
+
+#include "objectRegistry.H"
+#include "Time.H"
+#include "argList.H"
+#include "OSspecific.H"
+#include "faMesh.H"
+#include "fvMesh.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+class faPatchData
+{
+public:
+    word name_;
+    word type_;
+    dictionary dict_;
+    label ownPolyPatchID_;
+    label ngbPolyPatchID_;
+    labelList edgeLabels_;
+    faPatchData()
+    :
+        name_(word::null),
+        type_(word::null),
+        ownPolyPatchID_(-1),
+        ngbPolyPatchID_(-1)
+    {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "addRegionOption.H"
+    argList::noParallel();
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createNamedMesh.H"
+
+    // Reading faMeshDefinition dictionary
+    IOdictionary faMeshDefinition
+    (
+        IOobject
+        (
+            "faMeshDefinition",
+            runTime.constant(),
+            "faMesh",
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    wordList polyMeshPatches
+    (
+        faMeshDefinition.lookup("polyMeshPatches")
+    );
+
+    const dictionary& bndDict = faMeshDefinition.subDict("boundary");
+
+    const wordList faPatchNames(bndDict.toc());
+
+    List<faPatchData> faPatches(faPatchNames.size()+1);
+
+    forAll(faPatchNames, patchI)
+    {
+        const dictionary& curPatchDict = bndDict.subDict(faPatchNames[patchI]);
+
+        faPatches[patchI].name_ = faPatchNames[patchI];
+
+        faPatches[patchI].type_ = word(curPatchDict.lookup("type"));
+
+        const word ownName = curPatchDict.lookup("ownerPolyPatch");
+
+        faPatches[patchI].ownPolyPatchID_ =
+            mesh.boundaryMesh().findPatchID(ownName);
+
+        if (faPatches[patchI].ownPolyPatchID_ < 0)
+        {
+            FatalErrorIn("makeFaMesh:")
+                << "neighbourPolyPatch " << ownName << " does not exist"
+                << exit(FatalError);
+        }
+
+        const word neiName = curPatchDict.lookup("neighbourPolyPatch");
+
+        faPatches[patchI].ngbPolyPatchID_ =
+            mesh.boundaryMesh().findPatchID(neiName);
+
+        if (faPatches[patchI].ngbPolyPatchID_ < 0)
+        {
+            FatalErrorIn("makeFaMesh:")
+                << "neighbourPolyPatch " << neiName << " does not exist"
+                << exit(FatalError);
+        }
+    }
+
+    // Setting faceLabels list size
+    label size = 0;
+
+    labelList patchIDs(polyMeshPatches.size(), -1);
+
+    forAll(polyMeshPatches, patchI)
+    {
+        patchIDs[patchI] =
+            mesh.boundaryMesh().findPatchID(polyMeshPatches[patchI]);
+
+        if (patchIDs[patchI] < 0)
+        {
+            FatalErrorIn("makeFaMesh:")
+                << "Patch " << polyMeshPatches[patchI] << " does not exist"
+                << exit(FatalError);
+        }
+
+        size += mesh.boundaryMesh()[patchIDs[patchI]].size();
+    }
+
+    labelList faceLabels(size, -1);
+
+    sort(patchIDs);
+
+
+    // Filling of faceLabels list
+    label faceI = -1;
+
+    forAll(polyMeshPatches, patchI)
+    {
+        label start = mesh.boundaryMesh()[patchIDs[patchI]].start();
+
+        label size  = mesh.boundaryMesh()[patchIDs[patchI]].size();
+
+        for (label i = 0; i < size; ++i)
+        {
+            faceLabels[++faceI] = start + i;
+        }
+    }
+
+    // Creating faMesh
+    Info << "Create faMesh ... ";
+
+    faMesh areaMesh
+    (
+        mesh,
+        faceLabels
+    );
+    Info << "Done" << endl;
+
+
+    // Determination of faPatch ID for each boundary edge.
+    // Result is in the bndEdgeFaPatchIDs list
+    const indirectPrimitivePatch& patch = areaMesh.patch();
+
+    labelList faceCells(faceLabels.size(), -1);
+
+    forAll(faceCells, faceI)
+    {
+        label faceID = faceLabels[faceI];
+
+        faceCells[faceI] = mesh.faceOwner()[faceID];
+    }
+
+    labelList meshEdges =
+        patch.meshEdges
+        (
+            mesh.edges(),
+            mesh.cellEdges(),
+            faceCells
+        );
+
+    const labelListList& edgeFaces = mesh.edgeFaces();
+
+    const label nTotalEdges = patch.nEdges();
+    const label nInternalEdges = patch.nInternalEdges();
+
+    labelList bndEdgeFaPatchIDs(nTotalEdges - nInternalEdges, -1);
+
+    for (label edgeI = nInternalEdges; edgeI < nTotalEdges; ++edgeI)
+    {
+        label curMeshEdge = meshEdges[edgeI];
+
+        labelList curEdgePatchIDs(2, label(-1));
+
+        label patchI = -1;
+
+        forAll(edgeFaces[curMeshEdge], faceI)
+        {
+            label curFace = edgeFaces[curMeshEdge][faceI];
+
+            label curPatchID = mesh.boundaryMesh().whichPatch(curFace);
+
+            if (curPatchID != -1)
+            {
+                curEdgePatchIDs[++patchI] = curPatchID;
+            }
+        }
+
+        for (label pI = 0; pI < faPatches.size() - 1; ++pI)
+        {
+            if
+            (
+                (
+                    curEdgePatchIDs[0] == faPatches[pI].ownPolyPatchID_
+                 && curEdgePatchIDs[1] == faPatches[pI].ngbPolyPatchID_
+                )
+                ||
+                (
+                    curEdgePatchIDs[1] == faPatches[pI].ownPolyPatchID_
+                 && curEdgePatchIDs[0] == faPatches[pI].ngbPolyPatchID_
+                )
+            )
+            {
+                bndEdgeFaPatchIDs[edgeI - nInternalEdges] = pI;
+                break;
+            }
+        }
+    }
+
+
+    // Set edgeLabels for each faPatch
+    for (label pI=0; pI<(faPatches.size()-1); ++pI)
+    {
+        SLList<label> tmpList;
+
+        forAll(bndEdgeFaPatchIDs, eI)
+        {
+            if (bndEdgeFaPatchIDs[eI] == pI)
+            {
+                tmpList.append(nInternalEdges + eI);
+            }
+        }
+
+        faPatches[pI].edgeLabels_ = tmpList;
+    }
+
+    // Check for undefined edges
+    SLList<label> tmpList;
+
+    forAll(bndEdgeFaPatchIDs, eI)
+    {
+        if (bndEdgeFaPatchIDs[eI] == -1)
+        {
+            tmpList.append(nInternalEdges + eI);
+        }
+    }
+
+    if (tmpList.size() > 0)
+    {
+        label pI = faPatches.size()-1;
+
+        faPatches[pI].name_ = "undefined";
+        faPatches[pI].type_ = "patch";
+        faPatches[pI].edgeLabels_ = tmpList;
+    }
+
+    // Add good patches to faMesh
+    SLList<faPatch*> faPatchLst;
+
+    for (label pI = 0; pI < faPatches.size(); ++pI)
+    {
+        faPatches[pI].dict_.add("type", faPatches[pI].type_);
+        faPatches[pI].dict_.add("edgeLabels", faPatches[pI].edgeLabels_);
+        faPatches[pI].dict_.add
+        (
+            "ngbPolyPatchIndex",
+            faPatches[pI].ngbPolyPatchID_
+        );
+
+        if(faPatches[pI].edgeLabels_.size() > 0)
+        {
+            faPatchLst.append
+            (
+                faPatch::New
+                (
+                    faPatches[pI].name_,
+                    faPatches[pI].dict_,
+                    pI,
+                    areaMesh.boundary()
+                ).ptr()
+            );
+        }
+    }
+
+    if (args.optionFound("addEmptyPatch"))
+    {
+        word emptyPatchName(args.optionLookup("addEmptyPatch")());
+        dictionary emptyPatchDict;
+        emptyPatchDict.add("type", "empty");
+        emptyPatchDict.add("edgeLabels", labelList());
+        emptyPatchDict.add("ngbPolyPatchIndex", -1);
+
+        faPatchLst.append
+        (
+            faPatch::New
+            (
+                emptyPatchName,
+                emptyPatchDict,
+                faPatchLst.size(),
+                areaMesh.boundary()
+            ).ptr()
+        );
+    }
+
+    Info << "Add faPatches ... ";
+    areaMesh.addFaPatches(List<faPatch*>(faPatchLst));
+    Info << "Done" << endl;
+
+    // Writing faMesh
+    Info << "Write finite area mesh ... ";
+    areaMesh.write();
+
+    Info << "\nEnd" << endl;
+
+    return 0;
+}
+
+// ************************************************************************* //

applications/utilities/mesh/conversion/ansysToFoam/ansysToFoam.L

@@ -75,11 +75,10 @@ PtrList<SLList<label>> slPatchCellFaces;
 Map<word> cellTypes;
 label currentTypei = -1;
 
-
 // Dummy yywrap to keep yylex happy at compile time.
 // It is called by yylex but is not used as the mechanism to change file.
 // See <<EOF>>
-#if YY_FLEX_MINOR_VERSION < 6 && YY_FLEX_SUBMINOR_VERSION < 34
+#if YY_FLEX_MAJOR_VERSION <= 2 && YY_FLEX_MINOR_VERSION <= 5 && YY_FLEX_SUBMINOR_VERSION < 34
 extern "C" int yywrap()
 #else
 int yyFlexLexer::yywrap()

applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L

@@ -126,7 +126,7 @@ void uniquify(word& name, HashSet<word>& patchNames)
 // Dummy yywrap to keep yylex happy at compile time.
 // It is called by yylex but is not used as the mechanism to change file.
 // See <<EOF>>
-#if YY_FLEX_MINOR_VERSION < 6 && YY_FLEX_SUBMINOR_VERSION < 34
+#if YY_FLEX_MAJOR_VERSION <= 2 && YY_FLEX_MINOR_VERSION <= 5 && YY_FLEX_SUBMINOR_VERSION < 34
 extern "C" int yywrap()
 #else
 int yyFlexLexer::yywrap()

applications/utilities/mesh/conversion/fluentMeshToFoam/fluentMeshToFoam.L

@@ -103,7 +103,7 @@ wordList patchNameIDs(maxZoneID);
 // Dummy yywrap to keep yylex happy at compile time.
 // It is called by yylex but is not used as the mechanism to change file.
 // See <<EOF>>
-#if YY_FLEX_MINOR_VERSION < 6 && YY_FLEX_SUBMINOR_VERSION < 34
+#if YY_FLEX_MAJOR_VERSION <= 2 && YY_FLEX_MINOR_VERSION <= 5 && YY_FLEX_SUBMINOR_VERSION < 34
 extern "C" int yywrap()
 #else
 int yyFlexLexer::yywrap()

applications/utilities/mesh/conversion/gambitToFoam/gambitToFoam.L

@@ -82,7 +82,7 @@ label nValuesForPatchFaces = 0;
 // Dummy yywrap to keep yylex happy at compile time.
 // It is called by yylex but is not used as the mechanism to change file.
 // See <<EOF>>
-#if YY_FLEX_MINOR_VERSION < 6 && YY_FLEX_SUBMINOR_VERSION < 34
+#if YY_FLEX_MAJOR_VERSION <= 2 && YY_FLEX_MINOR_VERSION <= 5 && YY_FLEX_SUBMINOR_VERSION < 34
 extern "C" int yywrap()
 #else
 int yyFlexLexer::yywrap()

applications/utilities/mesh/manipulation/checkMesh/checkTopology.C

@@ -178,16 +178,16 @@ Foam::label Foam::checkTopology
             {
                 cells.insert(celli);
             }
-            forAll(cFaces, i)
+            for (const label facei : cFaces)
             {
-                if (cFaces[i] < 0 || cFaces[i] >= mesh.nFaces())
+                if (facei < 0 || facei >= mesh.nFaces())
                 {
                     cells.insert(celli);
                     break;
                 }
             }
         }
-        label nCells = returnReduce(cells.size(), sumOp<label>());
+        const label nCells = returnReduce(cells.size(), sumOp<label>());
 
         if (nCells > 0)
         {
@@ -203,7 +203,6 @@ Foam::label Foam::checkTopology
             {
                 mergeAndWrite(surfWriter(), cells);
             }
-
         }
         else
         {
@@ -238,7 +237,7 @@ Foam::label Foam::checkTopology
             noFailedChecks++;
         }
 
-        label nFaces = returnReduce(faces.size(), sumOp<label>());
+        const label nFaces = returnReduce(faces.size(), sumOp<label>());
 
         if (nFaces > 0)
         {
@@ -259,7 +258,7 @@ Foam::label Foam::checkTopology
         {
             noFailedChecks++;
 
-            label nFaces = returnReduce(faces.size(), sumOp<label>());
+            const label nFaces = returnReduce(faces.size(), sumOp<label>());
 
             Info<< "  <<Writing " << nFaces
                 << " faces with out-of-range or duplicate vertices to set "
@@ -280,7 +279,7 @@ Foam::label Foam::checkTopology
         {
             noFailedChecks++;
 
-            label nCells = returnReduce(cells.size(), sumOp<label>());
+            const label nCells = returnReduce(cells.size(), sumOp<label>());
 
             Info<< "  <<Writing " << nCells
                 << " cells with over used edges to set " << cells.name()
@@ -303,7 +302,7 @@ Foam::label Foam::checkTopology
             noFailedChecks++;
         }
 
-        label nFaces = returnReduce(faces.size(), sumOp<label>());
+        const label nFaces = returnReduce(faces.size(), sumOp<label>());
         if (nFaces > 0)
         {
             Info<< "  <<Writing " << nFaces
@@ -334,9 +333,9 @@ Foam::label Foam::checkTopology
             {
                 const labelUList& owners = patches[patchi].faceCells();
 
-                forAll(owners, i)
+                for (const label facei : owners)
                 {
-                    nInternalFaces[owners[i]]++;
+                    nInternalFaces[facei]++;
                 }
             }
         }
@@ -444,14 +443,14 @@ Foam::label Foam::checkTopology
                 (
                     label facei = mesh.nInternalFaces();
                     facei < mesh.nFaces();
-                    facei++
+                    ++facei
                 )
                 {
-                    label regioni = rs[mesh.faceOwner()[facei]];
+                    const label regioni = rs[mesh.faceOwner()[facei]];
                     const face& f = mesh.faces()[facei];
-                    forAll(f, fp)
+                    for (const label verti : f)
                     {
-                        label& pRegion = pointToRegion[f[fp]];
+                        label& pRegion = pointToRegion[verti];
                         if (pRegion == -1)
                         {
                             pRegion = regioni;
@@ -467,7 +466,7 @@ Foam::label Foam::checkTopology
                             regionDisconnected[regioni] = false;
                             regionDisconnected[pRegion] = false;
                             pRegion = -2;
-                            points.insert(f[fp]);
+                            points.insert(verti);
                         }
                     }
                 }
@@ -618,9 +617,8 @@ Foam::label Foam::checkTopology
             }
             Info<< endl;
 
-            forAll(faceZones, zoneI)
+            for (const faceZone& fz : faceZones)
             {
-                const faceZone& fz = faceZones[zoneI];
                 checkPatch(allGeometry, fz.name(), fz(), points);
                 Info<< endl;
             }
@@ -631,7 +629,7 @@ Foam::label Foam::checkTopology
         }
     }
 
-    label nPoints = returnReduce(points.size(), sumOp<label>());
+    const label nPoints = returnReduce(points.size(), sumOp<label>());
 
     if (nPoints)
     {
@@ -654,45 +652,44 @@ Foam::label Foam::checkTopology
 
         if (cellZones.size())
         {
-
             Info<< "    "
                 << setw(20) << "CellZone"
                 << setw(9) << "Cells"
                 << setw(9) << "Points"
-                << setw(13) << "BoundingBox" <<endl;
+                << ' ' << "BoundingBox" <<endl;
 
             const cellList& cells = mesh.cells();
             const faceList& faces = mesh.faces();
-            treeBoundBox bb(boundBox::invertedBox);
             PackedBoolList isZonePoint(mesh.nPoints());
 
-            forAll(cellZones, zoneI)
+            for (const cellZone& cZone : cellZones)
             {
-                const cellZone& cZone = cellZones[zoneI];
+                boundBox bb;
+                isZonePoint.reset();  // clears all bits (reset count)
 
-                forAll(cZone, i)
+                for (const label celli : cZone)
                 {
-                    const label cellI = cZone[i];
-                    const cell& cFaces = cells[cellI];
-                    forAll(cFaces, cFacei)
+                    for (const label facei : cells[celli])
                     {
-                        const face& f = faces[cFaces[cFacei]];
-                        forAll(f, fp)
+                        const face& f = faces[facei];
+                        for (const label verti : f)
                         {
-                            if (isZonePoint.set(f[fp]))
+                            if (isZonePoint.set(verti))
                             {
-                                bb.add(mesh.points()[f[fp]]);
+                                bb.add(mesh.points()[verti]);
                             }
                          }
                     }
                 }
 
+                bb.reduce();  // Global min/max
+
                 Info<< "    "
                     << setw(20) << cZone.name()
                     << setw(9) << returnReduce(cZone.size(), sumOp<label>())
                     << setw(9)
                     << returnReduce(isZonePoint.count(), sumOp<label>())
-                    << setw(3) << bb << endl;
+                    << ' ' << bb << endl;
             }
         }
         else

applications/utilities/mesh/manipulation/rotateMesh/rotateMesh.C

@@ -48,19 +48,19 @@ void RotateFields
 (
     const fvMesh& mesh,
     const IOobjectList& objects,
-    const tensor& T
+    const tensor& rotT
 )
 {
-    // Search list of objects for volScalarFields
+    // Objects of field type
     IOobjectList fields(objects.lookupClass(GeometricField::typeName));
 
     forAllIter(IOobjectList, fields, fieldIter)
     {
         Info<< "    Rotating " << fieldIter()->name() << endl;
 
-        GeometricField theta(*fieldIter(), mesh);
-        transform(theta, dimensionedTensor(T), theta);
-        theta.write();
+        GeometricField fld(*fieldIter(), mesh);
+        transform(fld, dimensionedTensor(rotT), fld);
+        fld.write();
     }
 }
 
@@ -69,6 +69,12 @@ void RotateFields
 
 int main(int argc, char *argv[])
 {
+    argList::addNote
+    (
+        "Rotate mesh points and vector/tensor fields\n"
+        "Rotation from the <n1> vector to the <n2> vector"
+    );
+
     timeSelector::addOptions();
 
     argList::addArgument("n1");
@@ -83,7 +89,7 @@ int main(int argc, char *argv[])
     vector n2(args.argRead<vector>(2));
     n2 /= mag(n2);
 
-    tensor T(rotationTensor(n1, n2));
+    const tensor rotT(rotationTensor(n1, n2));
 
     {
         pointIOField points
@@ -100,7 +106,7 @@ int main(int argc, char *argv[])
             )
         );
 
-        points = transform(T, points);
+        points = transform(rotT, points);
 
         // Set the precision of the points data to 10
         IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
@@ -123,15 +129,15 @@ int main(int argc, char *argv[])
         // Search for list of objects for this time
         IOobjectList objects(mesh, runTime.timeName());
 
-        RotateFields<volVectorField>(mesh, objects, T);
-        RotateFields<volSphericalTensorField>(mesh, objects, T);
-        RotateFields<volSymmTensorField>(mesh, objects, T);
-        RotateFields<volTensorField>(mesh, objects, T);
+        RotateFields<volVectorField>(mesh, objects, rotT);
+        RotateFields<volSphericalTensorField>(mesh, objects, rotT);
+        RotateFields<volSymmTensorField>(mesh, objects, rotT);
+        RotateFields<volTensorField>(mesh, objects, rotT);
 
-        RotateFields<surfaceVectorField>(mesh, objects, T);
-        RotateFields<surfaceSphericalTensorField>(mesh, objects, T);
-        RotateFields<surfaceSymmTensorField>(mesh, objects, T);
-        RotateFields<surfaceTensorField>(mesh, objects, T);
+        RotateFields<surfaceVectorField>(mesh, objects, rotT);
+        RotateFields<surfaceSphericalTensorField>(mesh, objects, rotT);
+        RotateFields<surfaceSymmTensorField>(mesh, objects, rotT);
+        RotateFields<surfaceTensorField>(mesh, objects, rotT);
     }
 
     Info<< "End\n" << endl;

applications/utilities/mesh/manipulation/transformPoints/transformPoints.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,6 +40,9 @@ Usage
     -rotate (vector vector)
         Rotates the points from the first vector to the second,
 
+    -rotate-angle (vector angle)
+        Rotate angle degrees about vector axis.
+
      or -yawPitchRoll (yawdegrees pitchdegrees rolldegrees)
      or -rollPitchYaw (rolldegrees pitchdegrees yawdegrees)
 
@@ -68,7 +71,6 @@ Usage
 #include "pointFields.H"
 #include "transformField.H"
 #include "transformGeometricField.H"
-#include "StringStream.H"
 #include "mathematicalConstants.H"
 
 using namespace Foam;
@@ -81,7 +83,7 @@ void readAndRotateFields
 (
     PtrList<GeoField>& flds,
     const fvMesh& mesh,
-    const tensor& T,
+    const tensor& rotT,
     const IOobjectList& objects
 )
 {
@@ -89,7 +91,7 @@ void readAndRotateFields
     forAll(flds, i)
     {
         Info<< "Transforming " << flds[i].name() << endl;
-        dimensionedTensor dimT("t", flds[i].dimensions(), T);
+        const dimensionedTensor dimT("t", flds[i].dimensions(), rotT);
         transform(flds[i], dimT, flds[i]);
     }
 }
@@ -155,23 +157,34 @@ int main(int argc, char *argv[])
         "Translate by specified <vector> - eg, '(1 0 0)' before rotations"
     );
     argList::addOption
+    (
+        "origin",
+        "point",
+        "Use specified <point> as origin for rotations"
+    );
+    argList::addOption
     (
         "rotate",
         "(vectorA vectorB)",
-        "Transform as a rotation between <vectorA> and <vectorB> "
-        "- eg, '( (1 0 0) (0 0 1) )'"
+        "Rotate from <vectorA> to <vectorB> - eg, '((1 0 0) (0 0 1))'"
+    );
+    argList::addOption
+    (
+        "rotate-angle",
+        "(vector scalar)",
+        "Rotate <angle> degrees about <vector> - eg, '((1 0 0) 45)'"
     );
     argList::addOption
     (
         "rollPitchYaw",
         "vector",
-        "Rotate by '(roll pitch yaw)' in degrees"
+        "Rotate by '(roll pitch yaw)' degrees"
     );
     argList::addOption
     (
         "yawPitchRoll",
         "vector",
-        "Rotate by '(yaw pitch roll)' in degrees"
+        "Rotate by '(yaw pitch roll)' degrees"
     );
     argList::addBoolOption
     (
@@ -182,25 +195,54 @@ int main(int argc, char *argv[])
     (
         "scale",
         "scalar | vector",
-        "Scale by the specified amount - eg, for a uniform [mm] to [m] scaling "
-        "use either (0.001 0.001 0.001)' or simply '0.001'"
+        "Scale by the specified amount - Eg, for uniform [mm] to [m] scaling "
+        "use either '(0.001 0.001 0.001)' or simply '0.001'"
     );
 
     #include "addRegionOption.H"
     #include "setRootCase.H"
+
+    const bool doRotateFields = args.optionFound("rotateFields");
+
+    // Verify that an operation has been specified
+    {
+        const List<word> operationNames
+        {
+            "translate",
+            "rotate",
+            "rotate-angle",
+            "rollPitchYaw",
+            "yawPitchRoll",
+            "scale"
+        };
+
+        if (!args.optionCount(operationNames))
+        {
+            FatalError
+                << "No operation supplied, "
+                << "use least one of the following:" << nl
+                << "   ";
+
+            for (const auto& opName : operationNames)
+            {
+                FatalError
+                    << " -" << opName;
+            }
+
+            FatalError
+                << nl << exit(FatalError);
+        }
+    }
+
     #include "createTime.H"
 
     word regionName = polyMesh::defaultRegion;
-    fileName meshDir;
+    fileName meshDir = polyMesh::meshSubDir;
 
     if (args.optionReadIfPresent("region", regionName))
     {
         meshDir = regionName/polyMesh::meshSubDir;
     }
-    else
-    {
-        meshDir = polyMesh::meshSubDir;
-    }
 
     pointIOField points
     (
@@ -216,17 +258,6 @@ int main(int argc, char *argv[])
         )
     );
 
-    const bool doRotateFields = args.optionFound("rotateFields");
-
-    // this is not actually stringent enough:
-    if (args.options().empty())
-    {
-        FatalErrorInFunction
-            << "No options supplied, please use one or more of "
-               "-translate, -rotate or -scale options."
-            << exit(FatalError);
-    }
-
     vector v;
     if (args.optionReadIfPresent("translate", v))
     {
@@ -235,6 +266,14 @@ int main(int argc, char *argv[])
         points += v;
     }
 
+    vector origin;
+    const bool useOrigin = args.optionReadIfPresent("origin", origin);
+    if (useOrigin)
+    {
+        Info<< "Set origin for rotations to " << origin << endl;
+        points -= origin;
+    }
+
     if (args.optionFound("rotate"))
     {
         Pair<vector> n1n2
@@ -243,15 +282,41 @@ int main(int argc, char *argv[])
         );
         n1n2[0] /= mag(n1n2[0]);
         n1n2[1] /= mag(n1n2[1]);
-        tensor T = rotationTensor(n1n2[0], n1n2[1]);
 
-        Info<< "Rotating points by " << T << endl;
+        const tensor rotT = rotationTensor(n1n2[0], n1n2[1]);
 
-        points = transform(T, points);
+        Info<< "Rotating points by " << rotT << endl;
+
+        points = transform(rotT, points);
 
         if (doRotateFields)
         {
-            rotateFields(args, runTime, T);
+            rotateFields(args, runTime, rotT);
+        }
+    }
+    else if (args.optionFound("rotate-angle"))
+    {
+        const Tuple2<vector, scalar> axisAngle
+        (
+            args.optionLookup("rotate-angle")()
+        );
+
+        Info<< "Rotating points " << nl
+            << "    about " << axisAngle.first() << nl
+            << "    angle " << axisAngle.second() << nl;
+
+        const quaternion quat
+        (
+            axisAngle.first(),
+            axisAngle.second() * pi/180.0  // degToRad
+        );
+
+        Info<< "Rotating points by quaternion " << quat << endl;
+        points = transform(quat, points);
+
+        if (doRotateFields)
+        {
+            rotateFields(args, runTime, quat.R());
         }
     }
     else if (args.optionReadIfPresent("rollPitchYaw", v))
@@ -264,14 +329,14 @@ int main(int argc, char *argv[])
         // degToRad
         v *= pi/180.0;
 
-        const quaternion R(quaternion::rotationSequence::XYZ, v);
+        const quaternion quat(quaternion::rotationSequence::XYZ, v);
 
-        Info<< "Rotating points by quaternion " << R << endl;
-        points = transform(R, points);
+        Info<< "Rotating points by quaternion " << quat << endl;
+        points = transform(quat, points);
 
         if (doRotateFields)
         {
-            rotateFields(args, runTime, R.R());
+            rotateFields(args, runTime, quat.R());
         }
     }
     else if (args.optionReadIfPresent("yawPitchRoll", v))
@@ -284,14 +349,14 @@ int main(int argc, char *argv[])
         // degToRad
         v *= pi/180.0;
 
-        const quaternion R(quaternion::rotationSequence::ZYX, v);
+        const quaternion quat(quaternion::rotationSequence::ZYX, v);
 
-        Info<< "Rotating points by quaternion " << R << endl;
-        points = transform(R, points);
+        Info<< "Rotating points by quaternion " << quat << endl;
+        points = transform(quat, points);
 
         if (doRotateFields)
         {
-            rotateFields(args, runTime, R.R());
+            rotateFields(args, runTime, quat.R());
         }
     }
 
@@ -325,6 +390,13 @@ int main(int argc, char *argv[])
         }
     }
 
+    if (useOrigin)
+    {
+        Info<< "Unset origin for rotations from " << origin << endl;
+        points += origin;
+    }
+
+
     // Set the precision of the points data to 10
     IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
 

applications/utilities/parallelProcessing/decomposePar/Make/files

@@ -4,6 +4,8 @@ domainDecompositionMesh.C
 domainDecompositionDistribute.C
 dimFieldDecomposer.C
 pointFieldDecomposer.C
+faMeshDecomposition.C
+faFieldDecomposer.C
 lagrangianFieldDecomposer.C
 
 EXE = $(FOAM_APPBIN)/decomposePar

applications/utilities/parallelProcessing/decomposePar/Make/options

@@ -2,6 +2,7 @@ EXE_INC = \
     -I$(LIB_SRC)/parallel/decompose/decompose/lnInclude \
     -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -106,6 +106,11 @@ Usage
 #include "decompositionModel.H"
 #include "collatedFileOperation.H"
 
+#include "faCFD.H"
+#include "emptyFaPatch.H"
+#include "faMeshDecomposition.H"
+#include "faFieldDecomposer.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -1239,6 +1244,178 @@ int main(int argc, char *argv[])
                         processorDbList.set(proci, nullptr);
                     }
                 }
+
+                // Finite area mesh and field decomposition
+
+                IOobject faMeshBoundaryIOobj
+                (
+                    "faBoundary",
+                    mesh.time().findInstance
+                    (
+                        mesh.dbDir()/polyMesh::meshSubDir,
+                        "boundary"
+                    ),
+                    faMesh::meshSubDir,
+                    mesh,
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::NO_WRITE
+                );
+
+
+                if (faMeshBoundaryIOobj.typeHeaderOk<faBoundaryMesh>(true))
+                {
+                    Info << "\nFinite area mesh decomposition" << endl;
+
+                    faMeshDecomposition aMesh(mesh);
+
+                    aMesh.decomposeMesh();
+
+                    aMesh.writeDecomposition();
+
+
+                    // Construct the area fields
+                    // ~~~~~~~~~~~~~~~~~~~~~~~~
+                    PtrList<areaScalarField> areaScalarFields;
+                    readFields(aMesh, objects, areaScalarFields);
+
+                    PtrList<areaVectorField> areaVectorFields;
+                    readFields(aMesh, objects, areaVectorFields);
+
+                    PtrList<areaSphericalTensorField> areaSphericalTensorFields;
+                    readFields(aMesh, objects, areaSphericalTensorFields);
+
+                    PtrList<areaSymmTensorField> areaSymmTensorFields;
+                    readFields(aMesh, objects, areaSymmTensorFields);
+
+                    PtrList<areaTensorField> areaTensorFields;
+                    readFields(aMesh, objects, areaTensorFields);
+
+
+                    // Construct the edge fields
+                    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+                    PtrList<edgeScalarField> edgeScalarFields;
+                    readFields(aMesh, objects, edgeScalarFields);
+
+                    Info << endl;
+
+                    // Split the fields over processors
+                    for (label procI = 0; procI < mesh.nProcs(); procI++)
+                    {
+                        Info<< "Processor " << procI
+                            << ": finite area field transfer" << endl;
+
+                        // open the database
+                        Time processorDb
+                        (
+                            Time::controlDictName,
+                            args.rootPath(),
+                            args.caseName()/
+                            fileName(word("processor") + name(procI))
+                        );
+
+                        processorDb.setTime(runTime);
+
+                        // Read the mesh
+                        fvMesh procFvMesh
+                        (
+                            IOobject
+                            (
+                                regionName,
+                                processorDb.timeName(),
+                                processorDb
+                            )
+                        );
+
+                        faMesh procMesh(procFvMesh);
+
+                        // // Does not work.  HJ, 15/Aug/2017
+                        // const labelIOList& faceProcAddressing =
+                        //     procAddressing
+                        //     (
+                        //         procMeshList,
+                        //         procI,
+                        //         "faceProcAddressing",
+                        //         faceProcAddressingList
+                        //     );
+
+                        // const labelIOList& boundaryProcAddressing =
+                        //     procAddressing
+                        //     (
+                        //         procMeshList,
+                        //         procI,
+                        //         "boundaryProcAddressing",
+                        //         boundaryProcAddressingList
+                        //     );
+
+                        labelIOList faceProcAddressing
+                        (
+                            IOobject
+                            (
+                                "faceProcAddressing",
+                                "constant",
+                                procMesh.meshSubDir,
+                                procFvMesh,
+                                IOobject::MUST_READ,
+                                IOobject::NO_WRITE
+                            )
+                        );
+
+                        labelIOList boundaryProcAddressing
+                        (
+                            IOobject
+                            (
+                                "boundaryProcAddressing",
+                                "constant",
+                                procMesh.meshSubDir,
+                                procFvMesh,
+                                IOobject::MUST_READ,
+                                IOobject::NO_WRITE
+                            )
+                        );
+
+                        // FA fields
+                        if
+                        (
+                            areaScalarFields.size()
+                         || areaVectorFields.size()
+                         || areaSphericalTensorFields.size()
+                         || areaSymmTensorFields.size()
+                         || areaTensorFields.size()
+                         || edgeScalarFields.size()
+                        )
+                        {
+                            labelIOList edgeProcAddressing
+                            (
+                                IOobject
+                                (
+                                    "edgeProcAddressing",
+                                    "constant",
+                                    procMesh.meshSubDir,
+                                    procFvMesh,
+                                    IOobject::MUST_READ,
+                                    IOobject::NO_WRITE
+                                )
+                            );
+
+                            faFieldDecomposer fieldDecomposer
+                            (
+                                aMesh,
+                                procMesh,
+                                edgeProcAddressing,
+                                faceProcAddressing,
+                                boundaryProcAddressing
+                            );
+
+                            fieldDecomposer.decomposeFields(areaScalarFields);
+                            fieldDecomposer.decomposeFields(areaVectorFields);
+                            fieldDecomposer.decomposeFields(areaSphericalTensorFields);
+                            fieldDecomposer.decomposeFields(areaSymmTensorFields);
+                            fieldDecomposer.decomposeFields(areaTensorFields);
+
+                            fieldDecomposer.decomposeFields(edgeScalarFields);
+                        }
+                    }
+                }
             }
         }
     }

applications/utilities/parallelProcessing/decomposePar/faFieldDecomposer.C

@@ -0,0 +1,238 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "faFieldDecomposer.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+faFieldDecomposer::patchFieldDecomposer::patchFieldDecomposer
+(
+    const label sizeBeforeMapping,
+    const labelUList& addressingSlice,
+    const label addressingOffset
+)
+:
+    sizeBeforeMapping_(sizeBeforeMapping),
+    directAddressing_(addressingSlice)
+{
+    forAll (directAddressing_, i)
+    {
+        // Subtract one to align addressing.
+        // directAddressing_[i] -= addressingOffset + 1;
+        // ZT, 12/Nov/2010
+        directAddressing_[i] -= addressingOffset;
+    }
+}
+
+
+faFieldDecomposer::processorAreaPatchFieldDecomposer::
+processorAreaPatchFieldDecomposer
+(
+    const faMesh& mesh,
+    const labelUList& addressingSlice
+)
+:
+    sizeBeforeMapping_(mesh.nFaces()),
+    addressing_(addressingSlice.size()),
+    weights_(addressingSlice.size())
+{
+    const scalarField& weights = mesh.weights().internalField();
+    const labelList& own = mesh.edgeOwner();
+    const labelList& neighb = mesh.edgeNeighbour();
+
+    forAll (addressing_, i)
+    {
+        // Subtract one to align addressing.
+        label ai = addressingSlice[i];
+//         label ai = mag(addressingSlice[i]) - 1;
+
+        if (ai < neighb.size())
+        {
+            // This is a regular edge. it has been an internal edge
+            // of the original mesh and now it has become a edge
+            // on the parallel boundary
+            addressing_[i].setSize(2);
+            weights_[i].setSize(2);
+
+            addressing_[i][0] = own[ai];
+            addressing_[i][1] = neighb[ai];
+
+            weights_[i][0] = weights[ai];
+            weights_[i][1] = 1.0 - weights[ai];
+        }
+        else
+        {
+            // This is a edge that used to be on a cyclic boundary
+            // but has now become a parallel patch edge. I cannot
+            // do the interpolation properly (I would need to look
+            // up the different (edge) list of data), so I will
+            // just grab the value from the owner face
+            //
+            addressing_[i].setSize(1);
+            weights_[i].setSize(1);
+
+            addressing_[i][0] = own[ai];
+
+            weights_[i][0] = 1.0;
+        }
+    }
+}
+
+
+faFieldDecomposer::processorEdgePatchFieldDecomposer::
+processorEdgePatchFieldDecomposer
+(
+    label sizeBeforeMapping,
+    const labelUList& addressingSlice
+)
+:
+    sizeBeforeMapping_(sizeBeforeMapping),
+    addressing_(addressingSlice.size()),
+    weights_(addressingSlice.size())
+{
+    forAll (addressing_, i)
+    {
+        addressing_[i].setSize(1);
+        weights_[i].setSize(1);
+
+        addressing_[i][0] = mag(addressingSlice[i]) - 1;
+        weights_[i][0] = sign(addressingSlice[i]);
+    }
+}
+
+
+faFieldDecomposer::faFieldDecomposer
+(
+    const faMesh& completeMesh,
+    const faMesh& procMesh,
+    const labelList& edgeAddressing,
+    const labelList& faceAddressing,
+    const labelList& boundaryAddressing
+)
+:
+    completeMesh_(completeMesh),
+    procMesh_(procMesh),
+    edgeAddressing_(edgeAddressing),
+    faceAddressing_(faceAddressing),
+    boundaryAddressing_(boundaryAddressing),
+    patchFieldDecomposerPtrs_
+    (
+        procMesh_.boundary().size(),
+        static_cast<patchFieldDecomposer*>(NULL)
+    ),
+    processorAreaPatchFieldDecomposerPtrs_
+    (
+        procMesh_.boundary().size(),
+        static_cast<processorAreaPatchFieldDecomposer*>(NULL)
+    ),
+    processorEdgePatchFieldDecomposerPtrs_
+    (
+        procMesh_.boundary().size(),
+        static_cast<processorEdgePatchFieldDecomposer*>(NULL)
+    )
+{
+    forAll (boundaryAddressing_, patchi)
+    {
+        if (boundaryAddressing_[patchi] >= 0)
+        {
+            patchFieldDecomposerPtrs_[patchi] = new patchFieldDecomposer
+            (
+                completeMesh_.boundary()[boundaryAddressing_[patchi]].size(),
+                procMesh_.boundary()[patchi].patchSlice(edgeAddressing_),
+//                 completeMesh_.boundaryMesh()
+                completeMesh_.boundary()
+                [
+                    boundaryAddressing_[patchi]
+                ].start()
+            );
+        }
+        else
+        {
+            processorAreaPatchFieldDecomposerPtrs_[patchi] =
+                new processorAreaPatchFieldDecomposer
+                (
+                    completeMesh_,
+                    procMesh_.boundary()[patchi].patchSlice(edgeAddressing_)
+                );
+
+            processorEdgePatchFieldDecomposerPtrs_[patchi] =
+                new processorEdgePatchFieldDecomposer
+                (
+                    procMesh_.boundary()[patchi].size(),
+                    static_cast<const labelUList&>
+                    (
+                        procMesh_.boundary()[patchi].patchSlice
+                        (
+                            edgeAddressing_
+                        )
+                    )
+                );
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+faFieldDecomposer::~faFieldDecomposer()
+{
+    forAll (patchFieldDecomposerPtrs_, patchi)
+    {
+        if (patchFieldDecomposerPtrs_[patchi])
+        {
+            delete patchFieldDecomposerPtrs_[patchi];
+        }
+    }
+
+    forAll (processorAreaPatchFieldDecomposerPtrs_, patchi)
+    {
+        if (processorAreaPatchFieldDecomposerPtrs_[patchi])
+        {
+            delete processorAreaPatchFieldDecomposerPtrs_[patchi];
+        }
+    }
+
+    forAll (processorEdgePatchFieldDecomposerPtrs_, patchi)
+    {
+        if (processorEdgePatchFieldDecomposerPtrs_[patchi])
+        {
+            delete processorEdgePatchFieldDecomposerPtrs_[patchi];
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/decomposePar/faFieldDecomposer.H

@@ -0,0 +1,319 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    faFieldDecomposer
+
+Description
+    Finite Area area and edge field decomposer.
+
+Author
+    Zeljko Tukovic, FSB Zagreb
+    Hrvoje Jasak, Wikki Ltd.
+
+SourceFiles
+    faFieldDecomposer.C
+    faFieldDecomposerDecomposeFields.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faFieldDecomposer_H
+#define faFieldDecomposer_H
+
+#include "faMesh.H"
+#include "faPatchFieldMapper.H"
+#include "edgeFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class IOobjectList;
+
+/*---------------------------------------------------------------------------*\
+                    Class faFieldDecomposer Declaration
+\*---------------------------------------------------------------------------*/
+
+class faFieldDecomposer
+{
+public:
+
+        //- Patch field decomposer class
+        class patchFieldDecomposer
+        :
+            public faPatchFieldMapper
+        {
+            // Private data
+
+                label sizeBeforeMapping_;
+                labelList directAddressing_;
+
+        public:
+
+            // Constructors
+
+                //- Construct given addressing
+                patchFieldDecomposer
+                (
+                    const label sizeBeforeMapping,
+                    const labelUList& addressingSlice,
+                    const label addressingOffset
+                );
+
+
+            // Member functions
+
+                label size() const
+                {
+                    return directAddressing_.size();
+                }
+
+                virtual label sizeBeforeMapping() const
+                {
+                    return sizeBeforeMapping_;
+                }
+
+                bool direct() const
+                {
+                    return true;
+                }
+
+                virtual bool hasUnmapped() const
+                {
+                    return false;
+                }
+
+                const labelUList& directAddressing() const
+                {
+                    return directAddressing_;
+                }
+        };
+
+
+        //- Processor patch field decomposer class
+        class processorAreaPatchFieldDecomposer
+        :
+            public faPatchFieldMapper
+        {
+            // Private data
+
+                label sizeBeforeMapping_;
+                labelListList addressing_;
+                scalarListList weights_;
+
+        public:
+
+            //- Construct given addressing
+            processorAreaPatchFieldDecomposer
+            (
+                const faMesh& mesh,
+                const labelUList& addressingSlice
+            );
+
+
+            // Member functions
+
+                label size() const
+                {
+                    return addressing_.size();
+                }
+
+                virtual label sizeBeforeMapping() const
+                {
+                    return sizeBeforeMapping_;
+                }
+
+                bool direct() const
+                {
+                    return false;
+                }
+
+                virtual bool hasUnmapped() const
+                {
+                    return false;
+                }
+
+                const labelListList& addressing() const
+                {
+                    return addressing_;
+                }
+
+                const scalarListList& weights() const
+                {
+                    return weights_;
+                }
+        };
+
+
+        //- Processor patch field decomposer class
+        class processorEdgePatchFieldDecomposer
+        :
+            public faPatchFieldMapper
+        {
+            label sizeBeforeMapping_;
+            labelListList addressing_;
+            scalarListList weights_;
+
+        public:
+
+            //- Construct given addressing
+            processorEdgePatchFieldDecomposer
+            (
+                label sizeBeforeMapping,
+                const labelUList& addressingSlice
+            );
+
+
+            // Member functions
+
+                label size() const
+                {
+                    return addressing_.size();
+                }
+
+                virtual label sizeBeforeMapping() const
+                {
+                    return sizeBeforeMapping_;
+                }
+
+                bool direct() const
+                {
+                    return false;
+                }
+
+                virtual bool hasUnmapped() const
+                {
+                    return false;
+                }
+
+                const labelListList& addressing() const
+                {
+                    return addressing_;
+                }
+
+                const scalarListList& weights() const
+                {
+                    return weights_;
+                }
+        };
+
+
+private:
+
+    // Private data
+
+        //- Reference to complete mesh
+        const faMesh& completeMesh_;
+
+        //- Reference to processor mesh
+        const faMesh& procMesh_;
+
+        //- Reference to edge addressing
+        const labelList& edgeAddressing_;
+
+        //- Reference to face addressing
+        const labelList& faceAddressing_;
+
+        //- Reference to boundary addressing
+        const labelList& boundaryAddressing_;
+
+        //- List of patch field decomposers
+        List<patchFieldDecomposer*> patchFieldDecomposerPtrs_;
+
+        List<processorAreaPatchFieldDecomposer*>
+            processorAreaPatchFieldDecomposerPtrs_;
+
+        List<processorEdgePatchFieldDecomposer*>
+            processorEdgePatchFieldDecomposerPtrs_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        faFieldDecomposer(const faFieldDecomposer&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const faFieldDecomposer&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        faFieldDecomposer
+        (
+            const faMesh& completeMesh,
+            const faMesh& procMesh,
+            const labelList& edgeAddressing,
+            const labelList& faceAddressing,
+            const labelList& boundaryAddressing
+        );
+
+
+    // Destructor
+
+        ~faFieldDecomposer();
+
+
+    // Member Functions
+
+        //- Decompose area field
+        template<class Type>
+        tmp<GeometricField<Type, faPatchField, areaMesh> >
+        decomposeField
+        (
+            const GeometricField<Type, faPatchField, areaMesh>& field
+        ) const;
+
+        //- Decompose surface field
+        template<class Type>
+        tmp<GeometricField<Type, faePatchField, edgeMesh> >
+        decomposeField
+        (
+            const GeometricField<Type, faePatchField, edgeMesh>& field
+        ) const;
+
+        template<class GeoField>
+        void decomposeFields(const PtrList<GeoField>& fields) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "faFieldDecomposerDecomposeFields.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/decomposePar/faFieldDecomposerDecomposeFields.C

@@ -0,0 +1,237 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "faFieldDecomposer.H"
+#include "processorFaPatchField.H"
+#include "processorFaePatchField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh> >
+faFieldDecomposer::decomposeField
+(
+    const GeometricField<Type, faPatchField, areaMesh>& field
+) const
+{
+    // Create and map the internal field values
+    Field<Type> internalField(field.internalField(), faceAddressing_);
+
+    // Create and map the patch field values
+    PtrList<faPatchField<Type> > patchFields(boundaryAddressing_.size());
+
+    forAll (boundaryAddressing_, patchi)
+    {
+        if (boundaryAddressing_[patchi] >= 0)
+        {
+            patchFields.set
+            (
+                patchi,
+                faPatchField<Type>::New
+                (
+                    field.boundaryField()[boundaryAddressing_[patchi]],
+                    procMesh_.boundary()[patchi],
+                    DimensionedField<Type, areaMesh>::null(),
+                    *patchFieldDecomposerPtrs_[patchi]
+                )
+            );
+        }
+        else
+        {
+            patchFields.set
+            (
+                patchi,
+                new processorFaPatchField<Type>
+                (
+                    procMesh_.boundary()[patchi],
+                    DimensionedField<Type, areaMesh>::null(),
+                    Field<Type>
+                    (
+                        field.internalField(),
+                        *processorAreaPatchFieldDecomposerPtrs_[patchi]
+                    )
+                )
+            );
+        }
+    }
+
+    // Create the field for the processor
+    return tmp<GeometricField<Type, faPatchField, areaMesh> >
+    (
+        new GeometricField<Type, faPatchField, areaMesh>
+        (
+            IOobject
+            (
+                field.name(),
+                procMesh_.time().timeName(),
+                procMesh_(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            procMesh_,
+            field.dimensions(),
+            internalField,
+            patchFields
+        )
+    );
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, faePatchField, edgeMesh> >
+faFieldDecomposer::decomposeField
+(
+    const GeometricField<Type, faePatchField, edgeMesh>& field
+) const
+{
+    labelList mapAddr
+    (
+        labelList::subList
+        (
+            edgeAddressing_,
+            procMesh_.nInternalEdges()
+        )
+    );
+    forAll (mapAddr, i)
+    {
+        mapAddr[i] -= 1;
+    }
+
+    // Create and map the internal field values
+    Field<Type> internalField
+    (
+        field.internalField(),
+        mapAddr
+    );
+
+    // Problem with addressing when a processor patch picks up both internal
+    // edges and edges from cyclic boundaries. This is a bit of a hack, but
+    // I cannot find a better solution without making the internal storage
+    // mechanism for edgeFields correspond to the one of edges in polyMesh
+    // (i.e. using slices)
+    Field<Type> allEdgeField(field.mesh().nEdges());
+
+    forAll (field.internalField(), i)
+    {
+        allEdgeField[i] = field.internalField()[i];
+    }
+
+    forAll (field.boundaryField(), patchi)
+    {
+        const Field<Type> & p = field.boundaryField()[patchi];
+
+        const label patchStart = field.mesh().boundary()[patchi].start();
+
+        forAll (p, i)
+        {
+            allEdgeField[patchStart + i] = p[i];
+        }
+    }
+
+    // Create and map the patch field values
+    PtrList<faePatchField<Type> > patchFields(boundaryAddressing_.size());
+
+    forAll (boundaryAddressing_, patchi)
+    {
+        if (boundaryAddressing_[patchi] >= 0)
+        {
+            patchFields.set
+            (
+                patchi,
+                faePatchField<Type>::New
+                (
+                    field.boundaryField()[boundaryAddressing_[patchi]],
+                    procMesh_.boundary()[patchi],
+                    DimensionedField<Type, edgeMesh>::null(),
+                    *patchFieldDecomposerPtrs_[patchi]
+                )
+            );
+        }
+        else
+        {
+            patchFields.set
+            (
+                patchi,
+                new processorFaePatchField<Type>
+                (
+                    procMesh_.boundary()[patchi],
+                    DimensionedField<Type, edgeMesh>::null(),
+                    Field<Type>
+                    (
+                        allEdgeField,
+                        *processorEdgePatchFieldDecomposerPtrs_[patchi]
+                    )
+                )
+            );
+        }
+    }
+
+    // Create the field for the processor
+    return tmp<GeometricField<Type, faePatchField, edgeMesh> >
+    (
+        new GeometricField<Type, faePatchField, edgeMesh>
+        (
+            IOobject
+            (
+                field.name(),
+                procMesh_.time().timeName(),
+                procMesh_(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            procMesh_,
+            field.dimensions(),
+            internalField,
+            patchFields
+        )
+    );
+}
+
+
+template<class GeoField>
+void faFieldDecomposer::decomposeFields
+(
+    const PtrList<GeoField>& fields
+) const
+{
+    forAll (fields, fieldI)
+    {
+        decomposeField(fields[fieldI])().write();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.H

@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    faMeshDecomposition
+
+Description
+    Automatic faMesh decomposition class
+
+Author
+    Zeljko Tukovic, FSB Zagreb
+    Hrvoje Jasak, Wikki Ltd.
+
+SourceFiles
+    faMeshDecomposition.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faMeshDecomposition_H
+#define faMeshDecomposition_H
+
+#include "fvMesh.H"
+#include "faMesh.H"
+#include "labelList.H"
+#include "SLList.H"
+#include "PtrList.H"
+#include "point.H"
+
+#ifndef namespaceFoam
+#define namespaceFoam
+    using namespace Foam;
+#endif
+
+/*---------------------------------------------------------------------------*\
+                           Class faMeshDecomposition Declaration
+\*---------------------------------------------------------------------------*/
+
+class faMeshDecomposition
+:
+    public faMesh
+{
+    // Private data
+
+        //- Mesh decomposition control dictionary
+        IOdictionary decompositionDict_;
+
+        //- Number of processors in decomposition
+        label nProcs_;
+
+        //- Is the decomposition data to be distributed for each processor
+        bool distributed_;
+
+        //- Processor label for each cell
+        labelList faceToProc_;
+
+        //- Face labels for each processor mesh
+        labelListList procFaceLabels_;
+
+        //-
+        List<Map<label> > procMeshEdgesMap_;
+
+        //- Number of internal edges for each processor mesh
+        labelList procNInternalEdges_;
+
+        //- Edge labels for patches of processor meshes
+        List<List<List<label> > > procPatchEdgeLabels_;
+
+        //- Labels of points for each processor
+        labelListList procPatchPointAddressing_;
+
+        //- Labels of edges for each processor
+        labelListList procPatchEdgeAddressing_;
+
+        //- Labels of edges for each processor
+        labelListList procEdgeAddressing_;
+
+        //- Labels of faces for each processor
+        labelListList procFaceAddressing_;
+
+        //- Original patch index for every processor patch
+        labelListList procBoundaryAddressing_;
+
+        //- Sizes for processor mesh patches
+        // Excludes inter-processor boundaries
+        labelListList procPatchSize_;
+
+        //- Start indices for processor patches
+        // Excludes inter-processor boundaries
+        labelListList procPatchStartIndex_;
+
+        //- Neighbour processor ID for inter-processor boundaries
+        labelListList procNeighbourProcessors_;
+
+        //- Sizes for inter-processor patches
+        labelListList procProcessorPatchSize_;
+
+        //- Start indices for inter-processor patches
+        labelListList procProcessorPatchStartIndex_;
+
+        //- List of globally shared point labels
+        labelList globallySharedPoints_;
+
+        //- Are there cyclic-parallel faces
+        bool cyclicParallel_;
+
+
+    // Private Member Functions
+
+        void distributeFaces();
+
+public:
+
+    // Constructors
+
+        //- Construct from fvMesh
+        faMeshDecomposition(const fvMesh& mesh);
+
+
+    // Destructor
+
+        ~faMeshDecomposition();
+
+
+    // Member Functions
+
+        //- Number of processor in decomposition
+        label nProcs() const
+        {
+            return nProcs_;
+        }
+
+        //- Is the decomposition data to be distributed for each processor
+        bool distributed() const
+        {
+            return distributed_;
+        }
+
+        //- Decompose mesh
+        void decomposeMesh();
+
+        //- Write decomposition
+        bool writeDecomposition();
+
+        //- Cell-processor decomposition labels
+        const labelList& faceToProc() const
+        {
+            return faceToProc_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/reconstructPar/Make/files

@@ -1,3 +1,5 @@
+processorFaMeshes.C
+faFieldReconstructor.C
 reconstructPar.C
 
 EXE = $(FOAM_APPBIN)/reconstructPar

applications/utilities/parallelProcessing/reconstructPar/Make/options

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \

applications/utilities/parallelProcessing/reconstructPar/faFieldReconstructor.C

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "faFieldReconstructor.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::faFieldReconstructor::faFieldReconstructor
+(
+    faMesh& mesh,
+    const PtrList<faMesh>& procMeshes,
+    const PtrList<labelIOList>& edgeProcAddressing,
+    const PtrList<labelIOList>& faceProcAddressing,
+    const PtrList<labelIOList>& boundaryProcAddressing
+)
+:
+    mesh_(mesh),
+    procMeshes_(procMeshes),
+    edgeProcAddressing_(edgeProcAddressing),
+    faceProcAddressing_(faceProcAddressing),
+    boundaryProcAddressing_(boundaryProcAddressing)
+{}
+
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/reconstructPar/faFieldReconstructor.H

@@ -0,0 +1,205 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    faFieldReconstructor
+
+Description
+    FA area and edge field reconstructor.
+
+Author
+    Zeljko Tukovic, FSB Zagreb
+    Hrvoje Jasak, Wikki Ltd.
+
+SourceFiles
+    faFieldReconstructor.C
+    faFieldReconstructorReconstructFields.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faFieldReconstructor_H
+#define faFieldReconstructor_H
+
+#include "PtrList.H"
+#include "faMesh.H"
+#include "IOobjectList.H"
+#include "faPatchFieldMapper.H"
+#include "labelIOList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+
+/*---------------------------------------------------------------------------*\
+                     Class faFieldReconstructor Declaration
+\*---------------------------------------------------------------------------*/
+
+class faFieldReconstructor
+{
+    // Private data
+
+        //- Reconstructed mesh reference
+        faMesh& mesh_;
+
+        //- List of processor meshes
+        const PtrList<faMesh>& procMeshes_;
+
+        //- List of processor edge addressing lists
+        const PtrList<labelIOList>& edgeProcAddressing_;
+
+        //- List of processor face addressing lists
+        const PtrList<labelIOList>& faceProcAddressing_;
+
+        //- List of processor boundary addressing lists
+        const PtrList<labelIOList>& boundaryProcAddressing_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        faFieldReconstructor(const faFieldReconstructor&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const faFieldReconstructor&);
+
+
+public:
+
+        class faPatchFieldReconstructor
+        :
+            public faPatchFieldMapper
+        {
+            label size_;
+            label sizeBeforeMapping_;
+
+        public:
+
+            // Constructors
+
+                //- Construct given size
+                faPatchFieldReconstructor
+                (
+                    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 bool hasUnmapped() const
+                {
+                    return false;
+                }
+
+                virtual const labelUList& directAddressing() const
+                {
+                    return labelUList::null();
+                }
+        };
+
+
+    // Constructors
+
+        //- Construct from components
+        faFieldReconstructor
+        (
+            faMesh& mesh,
+            const PtrList<faMesh>& procMeshes,
+            const PtrList<labelIOList>& edgeProcAddressing,
+            const PtrList<labelIOList>& faceProcAddressing,
+            const PtrList<labelIOList>& boundaryProcAddressing
+        );
+
+
+    // Member Functions
+
+        //- Reconstruct area field
+        template<class Type>
+        tmp<GeometricField<Type, faPatchField, areaMesh> >
+        reconstructFaAreaField
+        (
+            const IOobject& fieldIoObject
+        );
+
+        //- Reconstruct edge field
+        template<class Type>
+        tmp<GeometricField<Type, faePatchField, edgeMesh> >
+        reconstructFaEdgeField
+        (
+            const IOobject& fieldIoObject
+        );
+
+        //- Reconstruct and write all area fields
+        template<class Type>
+        void reconstructFaAreaFields
+        (
+            const IOobjectList& objects
+        );
+
+        //- Reconstruct and write all area fields
+        template<class Type>
+        void reconstructFaEdgeFields
+        (
+            const IOobjectList& objects
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "faFieldReconstructorReconstructFields.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/reconstructPar/faFieldReconstructorReconstructFields.C

@@ -0,0 +1,642 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "faFieldReconstructor.H"
+#include "Time.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
+        (
+            isA<emptyFaPatch>(mesh_.boundary()[patchI])
+         && !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
+        (
+            isA<emptyFaPatch>(mesh_.boundary()[patchI])
+         && !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/processorFaMeshes.C

@@ -0,0 +1,260 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "processorFaMeshes.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::processorFaMeshes::read()
+{
+    forAll (fvMeshes_, procI)
+    {
+        meshes_.set
+        (
+            procI,
+            new faMesh(fvMeshes_[procI])
+        );
+
+        pointProcAddressing_.set
+        (
+            procI,
+            new labelIOList
+            (
+                IOobject
+                (
+                    "pointProcAddressing",
+                    meshes_[procI].time().findInstance
+                    (
+                        meshes_[procI].meshDir(),
+                        "pointProcAddressing"
+                    ),
+                    meshes_[procI].meshSubDir,
+                    fvMeshes_[procI],
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                )
+            )
+        );
+
+        edgeProcAddressing_.set
+        (
+            procI,
+            new labelIOList
+            (
+                IOobject
+                (
+                    "edgeProcAddressing",
+                    meshes_[procI].time().findInstance
+                    (
+                        meshes_[procI].meshDir(),
+                        "edgeProcAddressing"
+                    ),
+                    meshes_[procI].meshSubDir,
+                    fvMeshes_[procI],
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                )
+            )
+        );
+
+        faceProcAddressing_.set
+        (
+            procI,
+            new labelIOList
+            (
+                IOobject
+                (
+                    "faceProcAddressing",
+                    meshes_[procI].time().findInstance
+                    (
+                        meshes_[procI].meshDir(),
+                        "faceProcAddressing"
+                    ),
+                    meshes_[procI].meshSubDir,
+                    fvMeshes_[procI],
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                )
+            )
+        );
+
+        boundaryProcAddressing_.set
+        (
+            procI,
+            new labelIOList
+            (
+                IOobject
+                (
+                    "boundaryProcAddressing",
+                    meshes_[procI].time().findInstance
+                    (
+                        meshes_[procI].meshDir(),
+                        "faceProcAddressing"
+                    ),
+                    meshes_[procI].meshSubDir,
+                    fvMeshes_[procI],
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                )
+            )
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::processorFaMeshes::processorFaMeshes
+(
+    const PtrList<fvMesh>& processorFvMeshes
+)
+:
+    fvMeshes_(processorFvMeshes),
+    meshes_(processorFvMeshes.size()),
+    pointProcAddressing_(processorFvMeshes.size()),
+    edgeProcAddressing_(processorFvMeshes.size()),
+    faceProcAddressing_(processorFvMeshes.size()),
+    boundaryProcAddressing_(processorFvMeshes.size())
+{
+    read();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+// Foam::fvMesh::readUpdateState Foam::processorFaMeshes::readUpdate()
+// {
+//     fvMesh::readUpdateState stat = fvMesh::UNCHANGED;
+
+//     forAll (databases_, procI)
+//     {
+//         // Check if any new meshes need to be read.
+//         fvMesh::readUpdateState procStat = meshes_[procI].readUpdate();
+
+//         /*
+//         if (procStat != fvMesh::UNCHANGED)
+//         {
+//             Info<< "Processor " << procI
+//                 << " at time " << databases_[procI].timeName()
+//                 << " detected mesh change " << procStat
+//                 << endl;
+//         }
+//         */
+
+//         // Combine into overall mesh change status
+//         if (stat == fvMesh::UNCHANGED)
+//         {
+//             stat = procStat;
+//         }
+//         else
+//         {
+//             if (stat != procStat)
+//             {
+//                 FatalErrorIn("processorFaMeshes::readUpdate()")
+//                     << "Processor " << procI
+//                     << " has a different polyMesh at time "
+//                     << databases_[procI].timeName()
+//                     << " compared to any previous processors." << nl
+//                     << "Please check time " << databases_[procI].timeName()
+//                     << " directories on all processors for consistent"
+//                     << " mesh files."
+//                     << exit(FatalError);
+//             }
+//         }
+//     }
+
+//     if
+//     (
+//         stat == fvMesh::TOPO_CHANGE
+//      || stat == fvMesh::TOPO_PATCH_CHANGE
+//     )
+//     {
+//         // Reread all meshes and addresssing
+//         read();
+//     }
+//     return stat;
+// }
+
+
+// void Foam::processorFaMeshes::reconstructPoints(fvMesh& mesh)
+// {
+//     // Read the field for all the processors
+//     PtrList<pointIOField> procsPoints(meshes_.size());
+
+//     forAll (meshes_, procI)
+//     {
+//         procsPoints.set
+//         (
+//             procI,
+//             new pointIOField
+//             (
+//                 IOobject
+//                 (
+//                     "points",
+//                     meshes_[procI].time().timeName(),
+//                     polyMesh::meshSubDir,
+//                     meshes_[procI],
+//                     IOobject::MUST_READ,
+//                     IOobject::NO_WRITE
+//                 )
+//             )
+//         );
+//     }
+
+//     // Create the new points
+//     vectorField newPoints(mesh.nPoints());
+
+//     forAll (meshes_, procI)
+//     {
+//         const vectorField& procPoints = procsPoints[procI];
+
+//         // Set the cell values in the reconstructed field
+
+//         const labelList& pointProcAddressingI = pointProcAddressing_[procI];
+
+//         if (pointProcAddressingI.size() != procPoints.size())
+//         {
+//             FatalErrorIn("processorFaMeshes")
+//                 << "problem :"
+//                 << " pointProcAddressingI:" << pointProcAddressingI.size()
+//                 << " procPoints:" << procPoints.size()
+//                 << abort(FatalError);
+//         }
+
+//         forAll(pointProcAddressingI, pointI)
+//         {
+//             newPoints[pointProcAddressingI[pointI]] = procPoints[pointI];
+//         }
+//     }
+
+//     mesh.movePoints(newPoints);
+//     mesh.write();
+// }
+
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/reconstructPar/processorFaMeshes.H

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    processorFaMeshes
+
+Description
+    Container for processor mesh addressing.
+
+Author
+    Zeljko Tukovic, FSB Zagreb
+
+SourceFiles
+    processorFaMeshes.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef processorFaMeshes_H
+#define processorFaMeshes_H
+
+#include "PtrList.H"
+#include "fvMesh.H"
+#include "faMesh.H"
+#include "IOobjectList.H"
+#include "labelIOList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+
+/*---------------------------------------------------------------------------*\
+                     Class processorFaMeshes Declaration
+\*---------------------------------------------------------------------------*/
+
+class processorFaMeshes
+{
+    // Private data
+
+        //- List of processor finite volume meshes
+        const PtrList<fvMesh>& fvMeshes_;
+
+        //- List of processor finite area meshes
+        PtrList<faMesh> meshes_;
+
+        //- List of processor point addressing lists
+        PtrList<labelIOList> pointProcAddressing_;
+
+        //- List of processor face addressing lists
+        PtrList<labelIOList> edgeProcAddressing_;
+
+        //- List of processor cell addressing lists
+        PtrList<labelIOList> faceProcAddressing_;
+
+        //- List of processor boundary addressing lists
+        PtrList<labelIOList> boundaryProcAddressing_;
+
+
+    // Private Member Functions
+
+        //- Read all meshes
+        void read();
+
+        //- Disallow default bitwise copy construct
+        processorFaMeshes(const processorFaMeshes&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const processorFaMeshes&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        processorFaMeshes(const PtrList<fvMesh>& processorFvMeshes);
+
+
+    // Member Functions
+
+        //- Update the meshes based on the mesh files saved in
+        //  time directories
+//         fvMesh::readUpdateState readUpdate();
+
+        //- Reconstruct point position after motion in parallel
+//         void reconstructPoints(faMesh& mesh);
+
+        PtrList<faMesh>& meshes()
+        {
+            return meshes_;
+        }
+        const PtrList<labelIOList>& pointProcAddressing() const
+        {
+            return pointProcAddressing_;
+        }
+        PtrList<labelIOList>& edgeProcAddressing()
+        {
+            return edgeProcAddressing_;
+        }
+        const PtrList<labelIOList>& faceProcAddressing() const
+        {
+            return faceProcAddressing_;
+        }
+        const PtrList<labelIOList>& boundaryProcAddressing() const
+        {
+            return boundaryProcAddressing_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C

@@ -44,6 +44,11 @@ Description
 #include "pointFieldReconstructor.H"
 #include "reconstructLagrangian.H"
 
+#include "faCFD.H"
+#include "faMesh.H"
+#include "processorFaMeshes.H"
+#include "faFieldReconstructor.H"
+
 #include "cellSet.H"
 #include "faceSet.H"
 #include "pointSet.H"
@@ -710,6 +715,47 @@ int main(int argc, char *argv[])
             }
 
 
+            // If there are any FA fields, reconstruct them
+
+            if
+            (
+                objects.lookupClass(areaScalarField::typeName).size()
+             || objects.lookupClass(areaVectorField::typeName).size()
+             || objects.lookupClass(areaSphericalTensorField::typeName).size()
+             || objects.lookupClass(areaSymmTensorField::typeName).size()
+             || objects.lookupClass(areaTensorField::typeName).size()
+             || objects.lookupClass(edgeScalarField::typeName).size()
+            )
+            {
+                Info << "Reconstructing FA fields" << nl << endl;
+                
+                faMesh aMesh(mesh);
+
+                processorFaMeshes procFaMeshes(procMeshes.meshes());
+
+                faFieldReconstructor faReconstructor
+                (
+                    aMesh,
+                    procFaMeshes.meshes(),
+                    procFaMeshes.edgeProcAddressing(),
+                    procFaMeshes.faceProcAddressing(),
+                    procFaMeshes.boundaryProcAddressing()
+                );
+
+                faReconstructor.reconstructFaAreaFields<scalar>(objects);
+                faReconstructor.reconstructFaAreaFields<vector>(objects);
+                faReconstructor
+                    .reconstructFaAreaFields<sphericalTensor>(objects);
+                faReconstructor.reconstructFaAreaFields<symmTensor>(objects);
+                faReconstructor.reconstructFaAreaFields<tensor>(objects);
+                
+                faReconstructor.reconstructFaEdgeFields<scalar>(objects);
+            }
+            else
+            {
+                Info << "No FA fields" << nl << endl;
+            }
+
             if (!noReconstructSets)
             {
                 // Scan to find all sets

applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C

@@ -97,25 +97,6 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// file-scope helper
-static bool inFileNameList
-(
-    const fileNameList& nameList,
-    const word& name
-)
-{
-    forAll(nameList, i)
-    {
-        if (nameList[i] == name)
-        {
-            return true;
-        }
-    }
-
-    return false;
-}
-
-
 int main(int argc, char *argv[])
 {
     timeSelector::addOptions();
@@ -641,7 +622,7 @@ int main(int argc, char *argv[])
 
             Info<< "Write " << cloudName << " (";
 
-            bool cloudExists = inFileNameList(currentCloudDirs, cloudName);
+            bool cloudExists = currentCloudDirs.found(cloudName);
             reduce(cloudExists, orOp<bool>());
 
             {

applications/utilities/postProcessing/dataConversion/foamToVTK/Make/options

@@ -3,11 +3,13 @@ EXE_INC = \
     -I$(LIB_SRC)/conversion/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
     -lconversion \
     -ldynamicMesh \
+    -lfiniteArea \
     -llagrangian \
     -lgenericPatchFields

applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C

@@ -155,6 +155,7 @@ Note
 #include "Cloud.H"
 #include "passiveParticle.H"
 #include "stringOps.H"
+#include "areaFields.H"
 
 #include "meshSubsetHelper.H"
 #include "readFields.H"
@@ -862,6 +863,119 @@ int main(int argc, char *argv[])
         );
 
 
+        // Finite-area mesh and fields - need not exist
+
+        autoPtr<faMesh> aMeshPtr;
+        {
+            const bool throwing = FatalError.throwExceptions();
+            try
+            {
+                aMeshPtr.reset(new faMesh(meshRef.baseMesh()));
+            }
+            catch (Foam::error& err)
+            {
+                aMeshPtr.clear();
+            }
+            FatalError.throwExceptions(throwing);
+        }
+
+        if (aMeshPtr.valid())
+        {
+            // Construct the area fields
+
+            PtrList<const areaScalarField> aScalarFld;
+            PtrList<const areaVectorField> aVectorFld;
+            PtrList<const areaSphericalTensorField> aSphTensorf;
+            PtrList<const areaSymmTensorField> aSymTensorFld;
+            PtrList<const areaTensorField> aTensorFld;
+
+            const faMesh& aMesh = aMeshPtr();
+
+            if (!specifiedFields || selectedFields.size())
+            {
+                readFields
+                (
+                    aMesh,
+                    objects,
+                    selectedFields,
+                    aScalarFld
+                );
+                print("    areaScalar           :", Info, aScalarFld);
+
+                readFields
+                (
+                    aMesh,
+                    objects,
+                    selectedFields,
+                    aVectorFld
+                );
+                print("    areaVector           :", Info, aVectorFld);
+
+                readFields
+                (
+                    aMesh,
+                    objects,
+                    selectedFields,
+                    aSphTensorf
+                );
+                print("    areaSphericalTensor   :", Info, aSphTensorf);
+
+                readFields
+                (
+                    aMesh,
+                    objects,
+                    selectedFields,
+                    aSymTensorFld
+                );
+                print("    areaSymmTensor        :", Info, aSymTensorFld);
+
+                readFields
+                (
+                    aMesh,
+                    objects,
+                    selectedFields,
+                    aTensorFld
+                );
+                print("    areaTensor            :", Info, aTensorFld);
+            }
+
+            const label nAreaFields =
+            (
+                aScalarFld.size()
+              + aVectorFld.size()
+              + aSphTensorf.size()
+              + aSymTensorFld.size()
+              + aTensorFld.size()
+            );
+
+            fileName outputName(fvPath/"finiteArea");
+
+            mkDir(outputName);
+
+            const auto& pp = aMesh.patch();
+
+            vtk::surfaceMeshWriter writer
+            (
+                pp,
+                aMesh.name(),
+                outputName/"finiteArea" + "_" + timeDesc,
+                fmtType
+            );
+
+            // Number of fields
+            writer.beginCellData(nAreaFields);
+
+            writer.write(aScalarFld);
+            writer.write(aVectorFld);
+            writer.write(aSphTensorf);
+            writer.write(aSymTensorFld);
+            writer.write(aTensorFld);
+
+            writer.endCellData();
+
+            writer.writeFooter();
+        }
+
         PtrList<const pointScalarField> pScalarFld;
         PtrList<const pointVectorField> pVectorFld;
         PtrList<const pointSphericalTensorField> pSphTensorFld;

applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.C

@@ -72,6 +72,38 @@ label readFields
 }
 
 
+template<class GeoField>
+void readFields
+(
+    const typename GeoField::Mesh& mesh,
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields,
+    PtrList<const GeoField>& fields
+)
+{
+    // Search list of objects for fields of type GeomField
+    IOobjectList fieldObjects(objects.lookupClass(GeoField::typeName));
+
+    // Construct the fields
+    fields.setSize(fieldObjects.size());
+    label nFields = 0;
+
+    forAllIters(fieldObjects, iter)
+    {
+        if (selectedFields.empty() || selectedFields.found(iter()->name()))
+        {
+            fields.set
+            (
+                nFields++,
+                new GeoField(*iter(), mesh)
+            );
+        }
+    }
+
+    fields.setSize(nFields);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/CMakeLists.txt

@@ -13,6 +13,7 @@ include_directories(
     $ENV{WM_PROJECT_DIR}/src/OpenFOAM/lnInclude
     $ENV{WM_PROJECT_DIR}/src/OSspecific/$ENV{WM_OSTYPE}/lnInclude
     $ENV{WM_PROJECT_DIR}/src/conversion/lnInclude
+    $ENV{WM_PROJECT_DIR}/src/finiteArea/lnInclude
     $ENV{WM_PROJECT_DIR}/src/finiteVolume/lnInclude
     ${PROJECT_SOURCE_DIR}/../foamPv
     ${PROJECT_SOURCE_DIR}/../vtkPVFoam
@@ -66,6 +67,7 @@ target_link_libraries(
     vtkPVFoam-pv${PARAVIEW_VERSION_MAJOR}.${PARAVIEW_VERSION_MINOR}
     foamPv-pv${PARAVIEW_VERSION_MAJOR}.${PARAVIEW_VERSION_MINOR}
     conversion
+    finiteArea
     finiteVolume
     OpenFOAM
 )

applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoamReader.cxx

@@ -309,16 +309,16 @@ int vtkPVFoamReader::RequestData
     // Get the requested time step.
     // We only support requests for a single time step
 
-    int nRequestTime = 0;
-    double requestTime[nInfo];
+    std::vector<double> requestTime;
+    requestTime.reserve(nInfo);
 
     // taking port0 as the lead for other outputs would be nice, but fails when
     // a filter is added - we need to check everything
     // but since PREVIOUS_UPDATE_TIME_STEPS() is protected, relay the logic
     // to the vtkPVFoam::setTime() method
-    for (int infoI = 0; infoI < nInfo; ++infoI)
+    for (int infoi = 0; infoi < nInfo; ++infoi)
     {
-        vtkInformation *outInfo = outputVector->GetInformationObject(infoI);
+        vtkInformation *outInfo = outputVector->GetInformationObject(infoi);
 
         const int nsteps =
             outInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
@@ -329,7 +329,7 @@ int vtkPVFoamReader::RequestData
          && nsteps > 0
         )
         {
-            requestTime[nRequestTime] =
+            const double timeValue =
             (
                 1 == nsteps
                 // Only one time-step available, UPDATE_TIME_STEP is unreliable
@@ -337,16 +337,13 @@ int vtkPVFoamReader::RequestData
               : outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP())
             );
 
-            // outInfo->Set(vtkDataObject::DATA_TIME_STEP(), requestTime[nRequestTime]);
-            // this->SetTimeValue(requestedTimeValue);
-            ++nRequestTime;
+            // outInfo->Set(vtkDataObject::DATA_TIME_STEP(), timeValue);
+            // this->SetTimeValue(timeValue);
+            requestTime.push_back(timeValue);
         }
     }
 
-    if (nRequestTime)
-    {
-        backend_->setTime(nRequestTime, requestTime);
-    }
+    backend_->setTime(requestTime);
 
     vtkMultiBlockDataSet* output = vtkMultiBlockDataSet::SafeDownCast
     (

applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoamReader.h

@@ -25,9 +25,9 @@ Class
     vtkPVFoamReader
 
 Description
-    reads a dataset in OpenFOAM format
+    Reads a dataset in OpenFOAM format
 
-    vtkPVblockMeshReader creates an multiblock dataset.
+    vtkPVFoamReader creates an multiblock dataset.
     It uses the OpenFOAM infrastructure (fvMesh, etc) to handle mesh and
     field data.
 

applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCore.C

@@ -43,6 +43,35 @@ namespace Foam
     defineTypeNameAndDebug(foamPvCore, 0);
 }
 
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::foamPvCore::printDataArraySelection
+(
+    Ostream& os,
+    vtkDataArraySelection* select
+)
+{
+    if (!select)
+    {
+        return os;
+    }
+
+    const int n = select->GetNumberOfArrays();
+
+    os << n << '(';
+    for (int i=0; i < n; ++i)
+    {
+        if (i) os << ' ';
+        os  << select->GetArrayName(i) << '='
+            << (select->GetArraySetting(i) ? 1 : 0);
+    }
+    os << ')';
+
+    return os;
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::foamPvCore::addToBlock
@@ -229,11 +258,9 @@ Foam::word Foam::foamPvCore::getFoamName(const std::string& str)
         // Already checked for valid/invalid chars
         return word(str.substr(beg, beg+end), false);
     }
-    else
-    {
+
     return word::null;
 }
-}
 
 
 void Foam::foamPvCore::printMemory()

applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCore.H

@@ -119,13 +119,51 @@ public:
             setStart(startAt);
         }
 
-
         //- Increment the size
         void operator+=(label n)
         {
             setSize(size() + n);
         }
 
+        //- True if the labelRange intersects any key in the Map
+        template<class T>
+        bool intersects(const Map<T>& map) const
+        {
+            for (const label idx : *this)
+            {
+                if (map.found(idx))
+                {
+                    return true;
+                }
+            }
+
+            return false;
+        }
+
+        //- The intersection ids with keys in the Map
+        template<class T>
+        List<label> intersection(const Map<T>& map) const
+        {
+            List<label> indices(Foam::min(map.size(), this->size()));
+
+            if (indices.size())
+            {
+                label n = 0;
+
+                for (const label idx : *this)
+                {
+                    if (map.found(idx))
+                    {
+                        indices[n++] = idx;
+                    }
+                }
+
+                indices.setSize(n);
+            }
+
+            return indices;
+        }
+
     }; // End class arrayRange
 
 
@@ -139,6 +177,7 @@ private:
         //- Disallow default bitwise assignment
         void operator=(const foamPvCore&) = delete;
 
+
 public:
 
     //- Static data members
@@ -146,8 +185,15 @@ public:
 
 
     //- Construct null
-    foamPvCore()
-    {}
+    foamPvCore() = default;
+
+
+    //- Print information about vtkDataArraySelection
+    static Ostream& printDataArraySelection
+    (
+        Ostream& os,
+        vtkDataArraySelection* select
+    );
 
 
     //- Convenience method for the VTK multiblock API

applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCoreTemplates.C

@@ -152,7 +152,8 @@ void Foam::foamPvCore::setSelectedArrayEntries
 )
 {
     const int n = select->GetNumberOfArrays();
-    // disable everything not explicitly enabled
+
+    // Disable everything not explicitly enabled
     select->DisableAllArrays();
 
     // Loop through entries, enabling as required

applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvFields.H

@@ -29,7 +29,7 @@ Description
 #ifndef foamPvFields_H
 #define foamPvFields_H
 
-#include "volFields.H"
+#include "symmTensor.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -40,28 +40,18 @@ namespace Foam
                         Class foamPvFields Declaration
 \*---------------------------------------------------------------------------*/
 
-class foamPvFields
+struct foamPvFields
 {
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        foamPvFields(const foamPvFields&) = delete;
-
-        //- Disallow default bitwise assignment
-        void operator=(const foamPvFields&) = delete;
-
-public:
-
     //- Remapping for some data types
     template<class Type>
     inline static void remapTuple(float vec[])
     {}
 
+};
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-}; // End class foamPvFields
-
-
 // Template specialization for symmTensor
 template<>
 inline void Foam::foamPvFields::remapTuple<Foam::symmTensor>(float vec[])

applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/Make/options

@@ -4,6 +4,7 @@ EXE_INC = \
     ${c++LESSWARN} \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
@@ -14,6 +15,7 @@ EXE_INC = \
     -I../PVFoamReader
 
 LIB_LIBS = \
+    -lfiniteArea \
     -ldynamicMesh \
     -lconversion \
     -lgenericPatchFields \

applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoam.C

@@ -27,9 +27,12 @@ License
 #include "vtkPVFoamReader.h"
 
 // OpenFOAM includes
+#include "areaFaMesh.H"
+#include "faMesh.H"
 #include "fvMesh.H"
 #include "Time.H"
 #include "patchZones.H"
+#include "IOobjectList.H"
 
 // VTK includes
 #include "vtkDataArraySelection.h"
@@ -108,7 +111,7 @@ void Foam::vtkPVFoam::resetCounters()
     // Reset array range information (ids and sizes)
     rangeVolume_.reset();
     rangePatches_.reset();
-    rangeLagrangian_.reset();
+    rangeClouds_.reset();
     rangeCellZones_.reset();
     rangeFaceZones_.reset();
     rangePointZones_.reset();
@@ -131,9 +134,9 @@ bool Foam::vtkPVFoam::addOutputBlock
     vtkSmartPointer<vtkMultiBlockDataSet> block;
     int datasetNo = 0;
 
-    for (auto partId : selector)
-    {
-        if (selectedPartIds_.found(partId))
+    const List<label> partIds = selector.intersection(selectedPartIds_);
+
+    for (const auto partId : partIds)
     {
         const auto& longName = selectedPartIds_[partId];
         const word shortName = getFoamName(longName);
@@ -176,23 +179,27 @@ bool Foam::vtkPVFoam::addOutputBlock
             ++datasetNo;
         }
     }
-    }
 
     return datasetNo;
 }
 
 
-int Foam::vtkPVFoam::setTime(int nRequest, const double requestTimes[])
+int Foam::vtkPVFoam::setTime(const std::vector<double>& requestTimes)
 {
+    if (requestTimes.empty())
+    {
+        return 0;
+    }
+
     Time& runTime = dbPtr_();
 
     // Get times list
     instantList Times = runTime.times();
 
     int nearestIndex = timeIndex_;
-    for (int requestI = 0; requestI < nRequest; ++requestI)
+    for (const double& timeValue : requestTimes)
     {
-        int index = Time::findClosestTimeIndex(Times, requestTimes[requestI]);
+        const int index = Time::findClosestTimeIndex(Times, timeValue);
         if (index >= 0 && index != timeIndex_)
         {
             nearestIndex = index;
@@ -208,13 +215,15 @@ int Foam::vtkPVFoam::setTime(int nRequest, const double requestTimes[])
     if (debug)
     {
         Info<< "<beg> setTime(";
-        for (int requestI = 0; requestI < nRequest; ++requestI)
+        unsigned reqi = 0;
+        for (const double& timeValue : requestTimes)
         {
-            if (requestI) Info<< ", ";
-            Info<< requestTimes[requestI];
+            if (reqi) Info<< ", ";
+            Info<< timeValue;
+            ++reqi;
         }
         Info<< ") - previousIndex = " << timeIndex_
-            << ", nearestIndex = " << nearestIndex << endl;
+            << ", nearestIndex = " << nearestIndex << nl;
     }
 
     // See what has changed
@@ -223,8 +232,13 @@ int Foam::vtkPVFoam::setTime(int nRequest, const double requestTimes[])
         timeIndex_ = nearestIndex;
         runTime.setTime(Times[nearestIndex], nearestIndex);
 
-        // When the changes, so do the fields
-        meshState_ = meshPtr_ ? meshPtr_->readUpdate() : polyMesh::TOPO_CHANGE;
+        // When mesh changes, so do fields
+        meshState_ =
+        (
+            volMeshPtr_
+          ? volMeshPtr_->readUpdate()
+          : polyMesh::TOPO_CHANGE
+        );
 
         reader_->UpdateProgress(0.05);
 
@@ -237,7 +251,7 @@ int Foam::vtkPVFoam::setTime(int nRequest, const double requestTimes[])
         Info<< "<end> setTime() - selectedTime="
             << Times[nearestIndex].name() << " index=" << timeIndex_
             << "/" << Times.size()
-            << " meshUpdateState=" << updateStateName(meshState_) << endl;
+            << " meshUpdateState=" << updateStateName(meshState_) << nl;
     }
 
     return nearestIndex;
@@ -260,15 +274,17 @@ Foam::vtkPVFoam::vtkPVFoam
 :
     reader_(reader),
     dbPtr_(nullptr),
-    meshPtr_(nullptr),
+    volMeshPtr_(nullptr),
+    areaMeshPtr_(nullptr),
     meshRegion_(polyMesh::defaultRegion),
     meshDir_(polyMesh::meshSubDir),
     timeIndex_(-1),
     decomposePoly_(false),
     meshState_(polyMesh::TOPO_CHANGE),
-    rangeVolume_("unzoned"),
+    rangeVolume_("volMesh"),
+    rangeArea_("areaMesh"),
     rangePatches_("patch"),
-    rangeLagrangian_("lagrangian"),
+    rangeClouds_("lagrangian"),
     rangeCellZones_("cellZone"),
     rangeFaceZones_("faceZone"),
     rangePointZones_("pointZone"),
@@ -278,7 +294,7 @@ Foam::vtkPVFoam::vtkPVFoam
 {
     if (debug)
     {
-        Info<< "vtkPVFoam - " << FileName << endl;
+        Info<< "vtkPVFoam - " << FileName << nl;
         printMemory();
     }
 
@@ -315,13 +331,13 @@ Foam::vtkPVFoam::vtkPVFoam
     // could be stringent and insist the prefix match the directory name...
     // Note: cannot use fileName::name() due to the embedded '{}'
     string caseName(fileName(FileName).lessExt());
-    string::size_type beg = caseName.find_last_of("/{");
-    string::size_type end = caseName.find('}', beg);
+    const auto beg = caseName.find_last_of("/{");
+    const auto end = caseName.find('}', beg);
 
     if
     (
-        beg != string::npos && caseName[beg] == '{'
-     && end != string::npos && end == caseName.size()-1
+        beg != std::string::npos && caseName[beg] == '{'
+     && end != std::string::npos && end == caseName.size()-1
     )
     {
         meshRegion_ = caseName.substr(beg+1, end-beg-1);
@@ -343,7 +359,7 @@ Foam::vtkPVFoam::vtkPVFoam
         Info<< "fullCasePath=" << fullCasePath << nl
             << "FOAM_CASE=" << getEnv("FOAM_CASE") << nl
             << "FOAM_CASENAME=" << getEnv("FOAM_CASENAME") << nl
-            << "region=" << meshRegion_ << endl;
+            << "region=" << meshRegion_ << nl;
     }
 
     // Create time object
@@ -369,10 +385,11 @@ Foam::vtkPVFoam::~vtkPVFoam()
 {
     if (debug)
     {
-        Info<< "~vtkPVFoam" << endl;
+        Info<< "~vtkPVFoam" << nl;
     }
 
-    delete meshPtr_;
+    delete volMeshPtr_;
+    delete areaMeshPtr_;
 }
 
 
@@ -383,15 +400,18 @@ void Foam::vtkPVFoam::updateInfo()
     if (debug)
     {
         Info<< "<beg> updateInfo"
-            << " [meshPtr=" << (meshPtr_ ? "set" : "nullptr") << "] timeIndex="
-            << timeIndex_ << endl;
+            << " [volMeshPtr=" << (volMeshPtr_ ? "set" : "nullptr")
+            << "] timeIndex="
+            << timeIndex_ << nl;
     }
 
     resetCounters();
 
-    vtkDataArraySelection* partSelection = reader_->GetPartSelection();
+    // Part selection
+    {
+        vtkDataArraySelection* select = reader_->GetPartSelection();
 
-    // there are two ways to ensure we have the correct list of parts:
+        // There are two ways to ensure we have the correct list of parts:
         // 1. remove everything and then set particular entries 'on'
         // 2. build a 'char **' list and call SetArraysWithDefault()
         //
@@ -400,47 +420,42 @@ void Foam::vtkPVFoam::updateInfo()
         // layer doesn't care about that anyhow.
 
         HashSet<string> enabled;
-    if (!partSelection->GetNumberOfArrays() && !meshPtr_)
+        if (!select->GetNumberOfArrays() && !volMeshPtr_)
         {
             // Fake enable 'internalMesh' on the first call
             enabled = { "internalMesh" };
         }
         else
         {
-        // preserve the enabled selections
-        enabled = getSelectedArraySet(partSelection);
+            // Preserve the enabled selections
+            enabled = getSelectedArraySet(select);
         }
 
-    // Clear current mesh parts list
-    partSelection->RemoveAllArrays();
+        select->RemoveAllArrays();   // Clear existing list
 
         // Update mesh parts list - add Lagrangian at the bottom
-    updateInfoInternalMesh(partSelection);
-    updateInfoPatches(partSelection, enabled);
-    updateInfoSets(partSelection);
-    updateInfoZones(partSelection);
-    updateInfoLagrangian(partSelection);
+        updateInfoInternalMesh(select);
+        updateInfoPatches(select, enabled);
+        updateInfoSets(select);
+        updateInfoZones(select);
+        updateInfoAreaMesh(select);
+        updateInfoLagrangian(select);
 
-    // Adjust/restore the enabled selections
-    setSelectedArrayEntries(partSelection, enabled);
+        setSelectedArrayEntries(select, enabled);  // Adjust/restore selected
+    }
 
-    // Update volume, point and lagrangian fields
-    updateInfoFields<fvPatchField, volMesh>
-    (
-        reader_->GetVolFieldSelection()
-    );
-    updateInfoFields<pointPatchField, pointMesh>
-    (
-        reader_->GetPointFieldSelection()
-    );
-    updateInfoLagrangianFields
-    (
-        reader_->GetLagrangianFieldSelection()
-    );
+    // Volume and area fields - includes save/restore of selected
+    updateInfoContinuumFields(reader_->GetVolFieldSelection());
+
+    // Point fields - includes save/restore of selected
+    updateInfoPointFields(reader_->GetPointFieldSelection());
+
+    // Lagrangian fields - includes save/restore of selected
+    updateInfoLagrangianFields(reader_->GetLagrangianFieldSelection());
 
     if (debug)
     {
-        Info<< "<end> updateInfo" << endl;
+        Info<< "<end> updateInfo" << nl;
     }
 }
 
@@ -531,26 +546,29 @@ void Foam::vtkPVFoam::Update
     {
         if (debug)
         {
-            Info<< "<beg> updateFoamMesh" << endl;
+            Info<< "<beg> updateFoamMesh" << nl;
             printMemory();
         }
 
         if (!caching)
         {
-            delete meshPtr_;
-            meshPtr_ = nullptr;
+            delete volMeshPtr_;
+            delete areaMeshPtr_;
+
+            volMeshPtr_ = nullptr;
+            areaMeshPtr_ = nullptr;
         }
 
         // Check to see if the OpenFOAM mesh has been created
-        if (!meshPtr_)
+        if (!volMeshPtr_)
         {
             if (debug)
             {
                 Info<< "Creating OpenFOAM mesh for region " << meshRegion_
-                    << " at time=" << dbPtr_().timeName() << endl;
+                    << " at time=" << dbPtr_().timeName() << nl;
             }
 
-            meshPtr_ = new fvMesh
+            volMeshPtr_ = new fvMesh
             (
                 IOobject
                 (
@@ -567,13 +585,27 @@ void Foam::vtkPVFoam::Update
         {
             if (debug)
             {
-                Info<< "Using existing OpenFOAM mesh" << endl;
+                Info<< "Using existing OpenFOAM mesh" << nl;
             }
         }
 
+        if (rangeArea_.intersects(selectedPartIds_))
+        {
+            if (!areaMeshPtr_)
+            {
+                areaMeshPtr_ = new faMesh(*volMeshPtr_);
+            }
+        }
+        else
+        {
+            delete areaMeshPtr_;
+
+            areaMeshPtr_ = nullptr;
+        }
+
         if (debug)
         {
-            Info<< "<end> updateFoamMesh" << endl;
+            Info<< "<end> updateFoamMesh" << nl;
             printMemory();
         }
     }
@@ -600,6 +632,8 @@ void Foam::vtkPVFoam::Update
         reader_->UpdateProgress(0.7);
     }
 
+    convertMeshArea();
+
     convertMeshLagrangian();
 
     reader_->UpdateProgress(0.8);
@@ -607,8 +641,9 @@ void Foam::vtkPVFoam::Update
     // Update fields
     convertVolFields();
     convertPointFields();
-    convertLagrangianFields();
+    convertAreaFields();
 
+    convertLagrangianFields();
 
     // Assemble multiblock output
     addOutputBlock(output, cachedVtu_, rangeVolume_, true); // One dataset
@@ -619,16 +654,17 @@ void Foam::vtkPVFoam::Update
     addOutputBlock(output, cachedVtu_, rangeCellSets_);
     addOutputBlock(output, cachedVtp_, rangeFaceSets_);
     addOutputBlock(output, cachedVtp_, rangePointSets_);
+    addOutputBlock(output, cachedVtp_, rangeArea_);
     addOutputBlock
     (
         (outputLagrangian ? outputLagrangian : output),
         cachedVtp_,
-        rangeLagrangian_
+        rangeClouds_
     );
 
     if (debug)
     {
-        Info<< "done reader part" << nl << endl;
+        Info<< "done reader part" << nl << nl;
     }
     reader_->UpdateProgress(0.95);
 
@@ -655,8 +691,11 @@ void Foam::vtkPVFoam::UpdateFinalize()
 {
     if (!reader_->GetMeshCaching())
     {
-        delete meshPtr_;
-        meshPtr_ = nullptr;
+        delete volMeshPtr_;
+        delete areaMeshPtr_;
+
+        volMeshPtr_ = nullptr;
+        areaMeshPtr_ = nullptr;
     }
 
     reader_->UpdateProgress(1.0);
@@ -674,20 +713,20 @@ std::vector<double> Foam::vtkPVFoam::findTimes(const bool skipZero) const
 
         // find the first time for which this mesh appears to exist
         label begIndex = timeLst.size();
-        forAll(timeLst, timeI)
+        forAll(timeLst, timei)
         {
             if
             (
                 IOobject
                 (
                     "points",
-                    timeLst[timeI].name(),
+                    timeLst[timei].name(),
                     meshDir_,
                     runTime
                 ).typeHeaderOk<pointIOField>(false, false)
             )
             {
-                begIndex = timeI;
+                begIndex = timei;
                 break;
             }
         }
@@ -735,15 +774,15 @@ void Foam::vtkPVFoam::renderPatchNames
     const bool show
 )
 {
-    // always remove old actors first
+    // Always remove old actors first
 
-    forAll(patchTextActors_, patchi)
+    for (auto& actor : patchTextActors_)
     {
-        renderer->RemoveViewProp(patchTextActors_[patchi]);
+        renderer->RemoveViewProp(actor);
     }
     patchTextActors_.clear();
 
-    if (show && meshPtr_)
+    if (show && volMeshPtr_)
     {
         // get the display patches, strip off any prefix/suffix
         hashedWordList selectedPatches = getSelected
@@ -757,7 +796,7 @@ void Foam::vtkPVFoam::renderPatchNames
             return;
         }
 
-        const polyBoundaryMesh& pbMesh = meshPtr_->boundaryMesh();
+        const polyBoundaryMesh& pbMesh = volMeshPtr_->boundaryMesh();
 
         // Find the total number of zones
         // Each zone will take the patch name
@@ -784,19 +823,19 @@ void Foam::vtkPVFoam::renderPatchNames
 
             boolList featEdge(pp.nEdges(), false);
 
-            forAll(edgeFaces, edgeI)
+            forAll(edgeFaces, edgei)
             {
-                const labelList& eFaces = edgeFaces[edgeI];
+                const labelList& eFaces = edgeFaces[edgei];
 
                 if (eFaces.size() == 1)
                 {
                     // Note: could also do ones with > 2 faces but this gives
                     // too many zones for baffles
-                    featEdge[edgeI] = true;
+                    featEdge[edgei] = true;
                 }
                 else if (mag(n[eFaces[0]] & n[eFaces[1]]) < 0.5)
                 {
-                    featEdge[edgeI] = true;
+                    featEdge[edgei] = true;
                 }
             }
 
@@ -808,7 +847,7 @@ void Foam::vtkPVFoam::renderPatchNames
             labelList zoneNFaces(pZones.nZones(), 0);
 
             // Create storage for additional zone centres
-            forAll(zoneNFaces, zoneI)
+            forAll(zoneNFaces, zonei)
             {
                 zoneCentre[patchi].append(Zero);
             }
@@ -816,14 +855,14 @@ void Foam::vtkPVFoam::renderPatchNames
             // Do averaging per individual zone
             forAll(pp, facei)
             {
-                label zoneI = pZones[facei];
-                zoneCentre[patchi][zoneI] += pp[facei].centre(pp.points());
-                zoneNFaces[zoneI]++;
+                const label zonei = pZones[facei];
+                zoneCentre[patchi][zonei] += pp[facei].centre(pp.points());
+                zoneNFaces[zonei]++;
             }
 
-            forAll(zoneCentre[patchi], zoneI)
+            forAll(zoneCentre[patchi], zonei)
             {
-                zoneCentre[patchi][zoneI] /= zoneNFaces[zoneI];
+                zoneCentre[patchi][zonei] /= zoneNFaces[zonei];
             }
         }
 
@@ -842,7 +881,7 @@ void Foam::vtkPVFoam::renderPatchNames
         if (debug)
         {
             Info<< "displayed zone centres = " << displayZoneI << nl
-                << "zones per patch = " << nZones << endl;
+                << "zones per patch = " << nZones << nl;
         }
 
         // Set the size of the patch labels to max number of zones
@@ -850,7 +889,7 @@ void Foam::vtkPVFoam::renderPatchNames
 
         if (debug)
         {
-            Info<< "constructing patch labels" << endl;
+            Info<< "constructing patch labels" << nl;
         }
 
         // Actor index
@@ -875,7 +914,7 @@ void Foam::vtkPVFoam::renderPatchNames
                 {
                     Info<< "patch name = " << pp.name() << nl
                         << "anchor = " << zoneCentre[patchi][globalZoneI] << nl
-                        << "globalZoneI = " << globalZoneI << endl;
+                        << "globalZoneI = " << globalZoneI << nl;
                 }
 
                 // Into a list for later removal
@@ -892,9 +931,9 @@ void Foam::vtkPVFoam::renderPatchNames
     }
 
     // Add text to each renderer
-    forAll(patchTextActors_, actori)
+    for (auto& actor : patchTextActors_)
     {
-        renderer->AddViewProp(patchTextActors_[actori]);
+        renderer->AddViewProp(actor);
     }
 }
 
@@ -902,10 +941,10 @@ void Foam::vtkPVFoam::renderPatchNames
 void Foam::vtkPVFoam::PrintSelf(ostream& os, vtkIndent indent) const
 {
     os  << indent << "Number of nodes: "
-        << (meshPtr_ ? meshPtr_->nPoints() : 0) << "\n";
+        << (volMeshPtr_ ? volMeshPtr_->nPoints() : 0) << "\n";
 
     os  << indent << "Number of cells: "
-        << (meshPtr_ ? meshPtr_->nCells() : 0) << "\n";
+        << (volMeshPtr_ ? volMeshPtr_->nCells() : 0) << "\n";
 
     os  << indent << "Number of available time steps: "
         << (dbPtr_.valid() ? dbPtr_().times().size() : 0) << "\n";
@@ -918,8 +957,8 @@ void Foam::vtkPVFoam::printInfo() const
 {
     std::cout
         << "Region: " << meshRegion_ << "\n"
-        << "nPoints: " << (meshPtr_ ? meshPtr_->nPoints() : 0) << "\n"
-        << "nCells: "  << (meshPtr_ ? meshPtr_->nCells() : 0) << "\n"
+        << "nPoints: " << (volMeshPtr_ ? volMeshPtr_->nPoints() : 0) << "\n"
+        << "nCells: "  << (volMeshPtr_ ? volMeshPtr_->nCells() : 0) << "\n"
         << "nTimes: "
         << (dbPtr_.valid() ? dbPtr_().times().size() : 0) << "\n";
 

applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoam.H

@@ -99,6 +99,7 @@ namespace Foam
 // OpenFOAM class forward declarations
 class argList;
 class Time;
+class faMesh;
 class fvMesh;
 class IOobjectList;
 class polyPatch;
@@ -127,9 +128,8 @@ class vtkPVFoam
         //  with the output fields.
         //  The original copy is reused for different timestep
         template<class DataType>
-        class foamVtkCaching
+        struct foamVtkCaching
         {
-        public:
             typedef DataType dataType;
 
             //- The geometry, without any cell/point data
@@ -208,7 +208,7 @@ class vtkPVFoam
 
 
         //- Bookkeeping for vtkPolyData
-        class foamVtpData
+        struct foamVtpData
         :
             public foamVtkCaching<vtkPolyData>,
             public foamVtkMeshMaps
@@ -216,7 +216,7 @@ class vtkPVFoam
 
 
         //- Bookkeeping for vtkUnstructuredGrid
-        class foamVtuData
+        struct foamVtuData
         :
             public foamVtkCaching<vtkUnstructuredGrid>,
             public foamVtkMeshMaps
@@ -231,8 +231,11 @@ class vtkPVFoam
         //- OpenFOAM time control
         autoPtr<Time> dbPtr_;
 
-        //- OpenFOAM mesh
-        fvMesh* meshPtr_;
+        //- OpenFOAM finite volume mesh
+        fvMesh* volMeshPtr_;
+
+        //- OpenFOAM finite area mesh
+        faMesh* areaMeshPtr_;
 
         //- The mesh region
         word meshRegion_;
@@ -262,8 +265,9 @@ class vtkPVFoam
         //  used to index into selectedPartIds and thus indirectly into
         //  cachedVtp, cachedVtu
         arrayRange rangeVolume_;
+        arrayRange rangeArea_;
         arrayRange rangePatches_;
-        arrayRange rangeLagrangian_;
+        arrayRange rangeClouds_;
         arrayRange rangeCellZones_;
         arrayRange rangeFaceZones_;
         arrayRange rangePointZones_;
@@ -295,6 +299,9 @@ class vtkPVFoam
         //- Internal mesh info
         void updateInfoInternalMesh(vtkDataArraySelection* select);
 
+        //- Finite area mesh info
+        void updateInfoAreaMesh(vtkDataArraySelection* select);
+
         //- Lagrangian info
         void updateInfoLagrangian(vtkDataArraySelection* select);
 
@@ -326,18 +333,28 @@ class vtkPVFoam
         template<template<class> class patchType, class meshType>
         void updateInfoFields
         (
-            vtkDataArraySelection* select
+            vtkDataArraySelection* select,
+            const IOobjectList& objects
         );
 
+        //- Volume/Area field info
+        void updateInfoContinuumFields(vtkDataArraySelection* select);
+
+        //- Point field info
+        void updateInfoPointFields(vtkDataArraySelection* select);
+
         //- Lagrangian field info
         void updateInfoLagrangianFields(vtkDataArraySelection* select);
 
 
       // Mesh conversion functions
 
-        //- Convert InternalMesh
+        //- Convert internalMesh
         void convertMeshVolume();
 
+        //- Convert areaMesh
+        void convertMeshArea();
+
         //- Convert Lagrangian points
         void convertMeshLagrangian();
 
@@ -490,9 +507,12 @@ class vtkPVFoam
         ) const;
 
 
-        //- Convert volume fields
+        //- Convert finite volume fields
         void convertVolFields();
 
+        //- Convert finite area fields
+        void convertAreaFields();
+
         //- Convert point fields
         void convertPointFields();
 
@@ -528,6 +548,14 @@ class vtkPVFoam
             const IOobjectList& objects
         );
 
+        //- Area fields - all types
+        template<class Type>
+        void convertAreaFields
+        (
+            const faMesh& mesh,
+            const IOobjectList& objects
+        );
+
         //- Volume field - all selected parts
         template<class Type>
         void convertVolFieldBlock
@@ -631,7 +659,7 @@ public:
         //- Set the runTime to the first plausible request time,
         //  returns the timeIndex
         //  sets to "constant" on error
-        int setTime(int count, const double requestTimes[]);
+        int setTime(const std::vector<double>& requestTimes);
 
 
         //- The current time index