BUGFIX: Several minor bugfixes and corrected code formatting. Author: Hrvoje Jasak. Merge: Dominik Christ.

ThirdParty/AllMake.stage4

@@ -56,10 +56,10 @@ echo Starting ThirdParty AllMake: Stage4
 echo ========================================
 echo
 
-# qt-everywhere-opensource-src-4.8.5
+# qt-everywhere-opensource-src-4.8.6
 if [ ! -z "$QT_THIRD_PARTY" ]
 then
-    ( rpm_make -p qt-everywhere-opensource-src-4.8.5 -s qt-everywhere-opensource-src-4.8.5.spec -u http://download.qt-project.org/official_releases/qt/4.8/4.8.5/qt-everywhere-opensource-src-4.8.5.tar.gz )
+    ( rpm_make -p qt-everywhere-opensource-src-4.8.6 -s qt-everywhere-opensource-src-4.8.6.spec -u http://download.qt-project.org/official_releases/qt/4.8/4.8.6/qt-everywhere-opensource-src-4.8.6.tar.gz )
 else
     echo "Using system installed QT"
     echo ""
@@ -75,10 +75,14 @@ then
     then
         ( rpm_make -p ParaView-4.0.1 -s ParaView-4.0.1.spec -u http://downloads.sourceforge.net/project/openfoam-extend/foam-extend-3.1/ThirdParty/ParaView-v4.0.1-source.tgz \
             -f --define='_qmakePath $QT_BIN_DIR/qmake'
+#        ( rpm_make -p ParaView-4.2.0 -s ParaView-4.2.0.spec -u http://www.paraview.org/files/v4.2/ParaView-v4.2.0-source.tar.gz \
+#            -f --define='_qmakePath $QT_BIN_DIR/qmake'
+#        ( rpm_make -p ParaView-4.2.0 -s ParaView-4.2.0.spec -u http://downloads.sourceforge.net/project/openfoam-extend/foam-extend-3.1/ThirdParty/ParaView-v4.2.0-source.tgz \
+#            -f --define='_qmakePath $QT_BIN_DIR/qmake'
         )
     else
         echo "WARNING: "
-        echo "WARNING: Skipping the installation of ParaView-4.0.1."
+        echo "WARNING: Skipping the installation of ParaView-4.2.0."
         echo "WARNING: Please make sure the QT_BIN_DIR environment variable properly"
         echo "WARNING: initialized in the file prefs.sh or prefs.csh"
         echo "WARNING: The command \$QT_BIN_DIR/qmake needs to be valid"

applications/solvers/basic/potentialDyMFoam/potentialDyMFoam.C

@@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    potentialDyMFoam
+
+Description
+    Transient solver for potential flow with dynamic mesh.
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::validOptions.insert("resetU", "");
+    argList::validOptions.insert("writep", "");
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createDynamicFvMesh.H"
+#   include "createFields.H"
+#   include "initTotalVolume.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+#       include "readPISOControls.H"
+#       include "checkTotalVolume.H"
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        bool meshChanged = mesh.update();
+        reduce(meshChanged, orOp<bool>());
+
+        p.internalField() = 0;
+
+        if (args.optionFound("resetU"))
+        {
+            U.internalField() = vector::zero;
+        }
+
+#       include "volContinuity.H"
+#       include "meshCourantNo.H"
+
+        // Solve potential flow equations
+        adjustPhi(phi, U, p);
+
+        for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
+        {
+            p.storePrevIter();
+
+            fvScalarMatrix pEqn
+            (
+                fvm::laplacian
+                (
+                    dimensionedScalar
+                    (
+                        "1",
+                        dimTime/p.dimensions()*dimensionSet(0, 2, -2, 0, 0),
+                        1
+                    ),
+                    p
+                )
+             ==
+                fvc::div(phi)
+            );
+
+            pEqn.setReference(pRefCell, pRefValue);
+            pEqn.solve();
+
+            if (nonOrth == nNonOrthCorr)
+            {
+                phi -= pEqn.flux();
+            }
+            else
+            {
+                p.relax();
+            }
+        }
+
+        Info<< "continuity error = "
+            << mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
+            << endl;
+
+        U = fvc::reconstruct(phi);
+        U.correctBoundaryConditions();
+
+        Info<< "Interpolated U error = "
+            << (sqrt(sum(sqr((fvc::interpolate(U) & mesh.Sf()) - phi)))
+            /sum(mesh.magSf())).value()
+            << endl;
+
+        // Calculate velocity magnitude
+        {
+            volScalarField magU = mag(U);
+
+            Info<< "mag(U): max: " << gMax(magU.internalField())
+                << " min: " << gMin(magU.internalField()) << endl;
+        }
+
+        if (args.optionFound("writep"))
+        {
+            // Find reference patch
+            label refPatch = -1;
+            scalar maxMagU = 0;
+
+            // Go through all velocity patches and find the one that fixes
+            // velocity to the largest value
+
+            forAll (U.boundaryField(), patchI)
+            {
+                const fvPatchVectorField& Upatch = U.boundaryField()[patchI];
+
+                if (Upatch.fixesValue())
+                {
+                    // Calculate mean velocity
+                    scalar u = sum(mag(Upatch));
+                    label patchSize = Upatch.size();
+
+                    reduce(u, sumOp<scalar>());
+                    reduce(patchSize, sumOp<label>());
+
+                    if (patchSize > 0)
+                    {
+                        scalar curMag = u/patchSize;
+
+                        if (curMag > maxMagU)
+                        {
+                            refPatch = patchI;
+
+                            maxMagU = curMag;
+                        }
+                    }
+                }
+            }
+
+            if (refPatch > -1)
+            {
+                // Calculate reference pressure
+                const fvPatchVectorField& Upatch = U.boundaryField()[refPatch];
+                const fvPatchScalarField& pPatch = p.boundaryField()[refPatch];
+
+                scalar patchE = sum(mag(pPatch + 0.5*magSqr(Upatch)));
+                label patchSize = Upatch.size();
+
+                reduce(patchE, sumOp<scalar>());
+                reduce(patchSize, sumOp<label>());
+
+                scalar e = patchE/patchSize;
+
+                Info<< "Using reference patch " << refPatch
+                    << " with mag(U) = " << maxMagU
+                    << " p + 0.5*U^2 = " << e << endl;
+
+                p.internalField() = e - 0.5*magSqr(U.internalField());
+                p.correctBoundaryConditions();
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/compressible/sonicDyMFoam/Make/files

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

applications/solvers/compressible/sonicDyMFoam/Make/options

@@ -0,0 +1,21 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lfiniteVolume \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -ldynamicMesh \
+    -lmeshTools \
+    -llduSolvers

applications/solvers/compressible/sonicDyMFoam/UEqn.H

@@ -0,0 +1,24 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+    );
+
+    if (oCorr == nOuterCorr - 1)
+    {
+        if (mesh.solutionDict().relax("UFinal"))
+        {
+            UEqn.relax(mesh.solutionDict().relaxationFactor("UFinal"));
+        }
+        else
+        {
+            UEqn.relax(1);
+        }
+    }
+    else
+    {
+        UEqn.relax();
+    }
+
+    solve(UEqn == -fvc::grad(p));

applications/solvers/compressible/sonicDyMFoam/createFields.H

@@ -0,0 +1,65 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<basicPsiThermo> pThermo
+    (
+        basicPsiThermo::New(mesh)
+    );
+    basicPsiThermo& thermo = pThermo();
+
+    volScalarField& T = const_cast<volScalarField&>(thermo.T());
+    volScalarField& p = thermo.p();
+    volScalarField& e = thermo.e();
+    const volScalarField& psi = thermo.psi();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh
+        ),
+        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"
+
+    // Store energy source term for under-relaxation
+    volScalarField pDivU
+    (
+        IOobject
+        (
+            "pDivU",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        p*fvc::div(phi/fvc::interpolate(rho))
+    );
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );

applications/solvers/compressible/sonicDyMFoam/eEqn.H

@@ -0,0 +1,43 @@
+{
+    fvScalarMatrix eEqn
+    (
+        fvm::ddt(rho, e)
+      + fvm::div(phi, e)
+      - fvm::laplacian(turbulence->alphaEff(), e)
+     ==
+      - fvm::SuSp(pDivU/e, e)
+        // viscous heating?
+    );
+
+    if (oCorr == nOuterCorr - 1)
+    {
+        if (mesh.solutionDict().relax("eFinal"))
+        {
+            eEqn.relax(mesh.solutionDict().relaxationFactor("eFinal"));
+        }
+        else
+        {
+            eEqn.relax(1);
+        }
+    }
+    else
+    {
+        eEqn.relax();
+    }
+
+    eEqn.solve();
+
+    // Bound the energy using TMin and TMax
+    {
+        dimensionedScalar Tstd("Tstd", dimTemperature, specie::Tstd);
+
+        volScalarField Cv = thermo.Cv();
+        volScalarField R = thermo.Cp() - Cv;
+
+        e = Foam::min(e, TMax*Cv + R*Tstd);
+        e = Foam::max(e, TMin*Cv + R*Tstd);
+        e.correctBoundaryConditions();
+    }
+
+    thermo.correct();
+}

applications/solvers/compressible/sonicDyMFoam/limitU.H

@@ -0,0 +1,17 @@
+{
+    // Bound the velocity
+    volScalarField magU = mag(U);
+
+    if (max(magU) > UMax)
+    {
+        volScalarField Ulimiter =
+            pos(magU - UMax)*UMax/(magU + smallU)
+          + neg(magU - UMax);
+
+        Ulimiter.max(scalar(0));
+        Ulimiter.min(scalar(1));
+
+        U *= Ulimiter;
+        U.correctBoundaryConditions();
+    }
+}

applications/solvers/compressible/sonicDyMFoam/pEqn.H

@@ -0,0 +1,69 @@
+{
+    U = UEqn.H()/UEqn.A();
+
+#   include "limitU.H"
+
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)*
+        (
+            (fvc::interpolate(U) & mesh.Sf())
+          - fvc::meshPhi(rho, U)
+        )
+    );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        // Store pressure for under-relaxation
+        p.storePrevIter();
+
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rho/UEqn.A(), p)
+        );
+
+        if
+        (
+//              oCorr == nOuterCorr - 1
+            corr == nCorr - 1
+         && nonOrth == nNonOrthCorr
+        )
+        {
+            pEqn.solve
+            (
+                mesh.solutionDict().solver(p.name() + "Final")
+            );
+        }
+        else
+        {
+            pEqn.solve(mesh.solutionDict().solver(p.name()));
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi = pEqn.flux();
+        }
+
+
+        // Bound the pressure
+        if (min(p) < pMin || max(p) > pMax)
+        {
+            p.max(pMin);
+            p.min(pMax);
+            p.correctBoundaryConditions();
+        }
+
+        // Relax the pressure
+        p.relax();
+    }
+
+#   include "compressibleContinuityErrs.H"
+
+    U -= fvc::grad(p)/UEqn.A();
+    U.correctBoundaryConditions();
+
+#   include "limitU.H"
+}

applications/solvers/compressible/sonicDyMFoam/readControls.H

@@ -0,0 +1,14 @@
+#   include "readTimeControls.H"
+#   include "readPIMPLEControls.H"
+
+    bool correctPhi = false;
+    if (pimple.found("correctPhi"))
+    {
+        correctPhi = Switch(pimple.lookup("correctPhi"));
+    }
+
+    bool checkMeshCourantNo = false;
+    if (pimple.found("checkMeshCourantNo"))
+    {
+        checkMeshCourantNo = Switch(pimple.lookup("checkMeshCourantNo"));
+    }

applications/solvers/compressible/sonicDyMFoam/readFieldBounds.H

@@ -0,0 +1,20 @@
+    // Read field bounds
+    dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds");
+
+    // Pressure bounds
+    dimensionedScalar pMin("pMin", dimPressure, 0);
+    dimensionedScalar pMax("pMax", dimPressure, GREAT);
+
+    fieldBounds.lookup("p") >> pMin.value() >> pMax.value();
+
+    // Temperature bounds
+    dimensionedScalar TMin("TMin", dimTemperature, 0);
+    dimensionedScalar TMax("TMax", dimTemperature, GREAT);
+
+    fieldBounds.lookup("T") >> TMin.value() >> TMax.value();
+
+    // Velocity bound
+    dimensionedScalar UMax("UMax", dimVelocity, GREAT);
+
+    fieldBounds.lookup(U.name()) >> UMax.value();
+    dimensionedScalar smallU("smallU", dimVelocity, 1e-10);

applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C

@@ -22,27 +22,36 @@ License
     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    sonicFoamAutoMotion
+    sonicDyMFoam
 
 Description
-    Transient solver for trans-sonic/supersonic, laminar flow of a
-    compressible gas with mesh motion..
+    Transient solver for trans-sonic/supersonic for laminar or flow of a
+    compressible gas with support for mesh motion and topological changes
+
+    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
+    pseudo-transient simulations.
+
+    Updated from sonicFoamAutoMotion by Hrvoje Jasak
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "motionSolver.H"
+#include "dynamicFvMesh.H"
+#include "specie.H"
+#include "basicPsiThermo.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
-
 #   include "setRootCase.H"
+
 #   include "createTime.H"
-#   include "createMesh.H"
-#   include "readThermodynamicProperties.H"
-#   include "readTransportProperties.H"
+#   include "createDynamicFvMesh.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
 
@@ -50,77 +59,77 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    autoPtr<Foam::motionSolver> motionPtr = motionSolver::New(mesh);
-
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.run())
     {
+#       include "readControls.H"
+#       include "readFieldBounds.H"
+#       include "compressibleCourantNo.H"
+#       include "setDeltaT.H"
+
+        runTime++;
+
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
+        bool meshChanged = mesh.update();
 
-        mesh.movePoints(motionPtr->newPoints());
+#       include "volContinuity.H"
+
+        if (checkMeshCourantNo)
+        {
+#           include "meshCourantNo.H"
+        }
+
+        // Mesh motion update
+        if (meshChanged)
+        {
+            T.correctBoundaryConditions();
+            p.correctBoundaryConditions();
+            e.correctBoundaryConditions();
+            thermo.correct();
+            rho = thermo.rho();
+
+            if (correctPhi)
+            {
+// #           include "correctPhi.H"
+            }
+        }
+
+        if (meshChanged)
+        {
+#           include "compressibleCourantNo.H"
+        }
+
+        // --- PIMPLE loop
+        label oCorr = 0;
+        do
+        {
+            // Under-relax pDivU term
+            pDivU.storePrevIter();
+
+            pDivU =
+                p*fvc::div
+                (
+                    phi/fvc::interpolate(rho)
+                  + fvc::meshPhi(rho, U)
+                );
+
+            pDivU.relax();
 
 #           include "rhoEqn.H"
-
-        fvVectorMatrix UEqn
-        (
-            fvm::ddt(rho, U)
-          + fvm::div(phi, U)
-          - fvm::laplacian(mu, U)
-        );
-
-        solve(UEqn == -fvc::grad(p));
-
-        solve
-        (
-            fvm::ddt(rho, e)
-          + fvm::div(phi, e)
-          - fvm::laplacian(mu, e)
-          ==
-          - p*fvc::div(phi/fvc::interpolate(rho) + fvc::meshPhi(rho, U))
-          + mu*magSqr(symm(fvc::grad(U)))
-        );
-
-        T = e/Cv;
-        psi = 1.0/(R*T);
+#           include "eEqn.H"
+#           include "UEqn.H"
 
             // --- PISO loop
-
             for (int corr = 0; corr < nCorr; corr++)
             {
-            U = UEqn.H()/UEqn.A();
-
-            surfaceScalarField phid
-            (
-                "phid",
-                fvc::interpolate(psi)*
-                (
-                    (fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)
-                )
-            );
-
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-            {
-                fvScalarMatrix pEqn
-                (
-                    fvm::ddt(psi, p)
-                  + fvm::div(phid, p)
-                  - fvm::laplacian(rho/UEqn.A(), p)
-                );
-
-                pEqn.solve();
-
-                phi = pEqn.flux();
+#               include "pEqn.H"
             }
 
-#           include "compressibleContinuityErrs.H"
+            // Recalculate density
+            rho = thermo.rho();
 
-            U -= fvc::grad(p)/UEqn.A();
-            U.correctBoundaryConditions();
-        }
-
-        rho = psi*p;
+            turbulence->correct();
+        } while (++oCorr < nOuterCorr);
 
         runTime.write();
 

applications/solvers/compressible/sonicFoamAutoMotion/Make/files

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

applications/solvers/compressible/sonicFoamAutoMotion/Make/options

@@ -1,10 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
-
-EXE_LIBS = \
-    -lfiniteVolume \
-    -ldynamicMesh \
-    -lmeshTools \
-    -llduSolvers

applications/solvers/compressible/sonicFoamAutoMotion/createFields.H

@@ -1,103 +0,0 @@
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-
-    Info<< "Reading field T\n" << endl;
-    volScalarField T
-    (
-        IOobject
-        (
-            "T",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-
-    Info<< "Calculating field e from T\n" << endl;
-    volScalarField e
-    (
-        IOobject
-        (
-            "e",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        Cv*T,
-        T.boundaryField().types()
-    );
-
-
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-
-    volScalarField psi
-    (
-        IOobject
-        (
-            "psi",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        1.0/(R*T)
-    );
-    psi.oldTime();
-
-    volScalarField rho
-    (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh
-        ),
-        psi*p
-    );
-
-#   include "compressibleCreatePhi.H"
-
-
-    Info<< "Creating field phid\n" << endl;
-    surfaceScalarField phid
-    (
-        IOobject
-        (
-            "phid",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        phi/fvc::interpolate(p),
-        phi.boundaryField().types()
-    );

applications/solvers/compressible/sonicFoamAutoMotion/readThermodynamicProperties.H

@@ -1,23 +0,0 @@
-    Info<< "Reading thermodynamicProperties\n" << endl;
-
-    IOdictionary thermodynamicProperties
-    (
-        IOobject
-        (
-            "thermodynamicProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    dimensionedScalar R
-    (
-        thermodynamicProperties.lookup("R")
-    );
-
-    dimensionedScalar Cv
-    (
-        thermodynamicProperties.lookup("Cv")
-    );

applications/solvers/compressible/steadyCompressibleFoam/readFieldBounds.H

@@ -2,19 +2,19 @@
     dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds");
 
     // Pressure bounds
-    dimensionedScalar pMin("pMin", p.dimensions(), 0);
-    dimensionedScalar pMax("pMax", p.dimensions(), 0);
+    dimensionedScalar pMin("pMin", dimPressure, 0);
+    dimensionedScalar pMax("pMax", dimPressure, GREAT);
 
-    fieldBounds.lookup(p.name()) >> pMin.value() >> pMax.value();
+    fieldBounds.lookup("p") >> pMin.value() >> pMax.value();
 
     // Temperature bounds
-    dimensionedScalar TMin("TMin", T.dimensions(), 0);
-    dimensionedScalar TMax("TMax", T.dimensions(), 0);
+    dimensionedScalar TMin("TMin", dimTemperature, 0);
+    dimensionedScalar TMax("TMax", dimTemperature, GREAT);
 
-    fieldBounds.lookup(T.name()) >> TMin.value() >> TMax.value();
+    fieldBounds.lookup("T") >> TMin.value() >> TMax.value();
 
     // Velocity bound
-    dimensionedScalar UMax("UMax", U.dimensions(), 0);
+    dimensionedScalar UMax("UMax", dimVelocity, GREAT);
 
     fieldBounds.lookup(U.name()) >> UMax.value();
     dimensionedScalar smallU("smallU", dimVelocity, 1e-10);

applications/solvers/compressible/steadyCompressibleMRFFoam/readFieldBounds.H

@@ -2,19 +2,19 @@
     dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds");
 
     // Pressure bounds
-    dimensionedScalar pMin("pMin", p.dimensions(), 0);
-    dimensionedScalar pMax("pMax", p.dimensions(), 0);
+    dimensionedScalar pMin("pMin", dimPressure, 0);
+    dimensionedScalar pMax("pMax", dimPressure, GREAT);
 
-    fieldBounds.lookup(p.name()) >> pMin.value() >> pMax.value();
+    fieldBounds.lookup("p") >> pMin.value() >> pMax.value();
 
     // Temperature bounds
-    dimensionedScalar TMin("TMin", T.dimensions(), 0);
-    dimensionedScalar TMax("TMax", T.dimensions(), 0);
+    dimensionedScalar TMin("TMin", dimTemperature, 0);
+    dimensionedScalar TMax("TMax", dimTemperature, GREAT);
 
-    fieldBounds.lookup(T.name()) >> TMin.value() >> TMax.value();
+    fieldBounds.lookup("T") >> TMin.value() >> TMax.value();
 
     // Velocity bound
-    dimensionedScalar UMax("UMax", U.dimensions(), 0);
+    dimensionedScalar UMax("UMax", dimVelocity, GREAT);
 
     fieldBounds.lookup(U.name()) >> UMax.value();
     dimensionedScalar smallU("smallU", dimVelocity, 1e-10);

applications/solvers/compressible/steadyCompressibleSRFFoam/readFieldBounds.H

@@ -1,19 +1,20 @@
+    // Read field bounds
     dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds");
 
     // Pressure bounds
-    dimensionedScalar pMin("pMin", p.dimensions(), 0);
-    dimensionedScalar pMax("pMax", p.dimensions(), 0);
+    dimensionedScalar pMin("pMin", dimPressure, 0);
+    dimensionedScalar pMax("pMax", dimPressure, GREAT);
 
-    fieldBounds.lookup(p.name()) >> pMin.value() >> pMax.value();
+    fieldBounds.lookup("p") >> pMin.value() >> pMax.value();
 
     // Temperature bounds
-    dimensionedScalar TMin("TMin", T.dimensions(), 0);
-    dimensionedScalar TMax("TMax", T.dimensions(), 0);
+    dimensionedScalar TMin("TMin", dimTemperature, 0);
+    dimensionedScalar TMax("TMax", dimTemperature, GREAT);
 
-    fieldBounds.lookup(T.name()) >> TMin.value() >> TMax.value();
+    fieldBounds.lookup("T") >> TMin.value() >> TMax.value();
 
     // Velocity bound
-    dimensionedScalar UrelMax("UrelMax", Urel.dimensions(), 0);
+    dimensionedScalar UrelMax("UrelMax", dimVelocity, GREAT);
 
     fieldBounds.lookup(Urel.name()) >> UrelMax.value();
     dimensionedScalar smallUrel("smallUrel", dimVelocity, 1e-10);

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

@@ -1616,7 +1616,7 @@ int main(int argc, char *argv[])
     // interior boundaries are handled via faceSets
     // cell zones will only be written if there is more than one
 
-    if (writeZones && cellGroupZoneID.size()>1)
+    if (writeZones && cellGroupZoneID.size() > 1)
     {
         Info<< "Adding Zones" << endl;
         List<pointZone*> pz(0);

applications/utilities/mesh/manipulation/createPatch/createPatchDict

@@ -38,7 +38,7 @@ FoamFile
 // Tolerance used in matching faces. Absolute tolerance is span of
 // face times this factor. To load incorrectly matches meshes set this
 // to a higher value.
-matchTolerance 1E-3;
+matchTolerance 1e-3;
 
 // Do a synchronisation of coupled points after creation of any patches.
 // Note: this does not work with points that are on multiple coupled patches

applications/utilities/mesh/manipulation/moveDyMEngineMesh/moveDyMEngineMesh.C

@@ -69,6 +69,7 @@ int main(int argc, char *argv[])
         mesh.update();
 
 #       include "checkVolContinuity.H"
+#       include "meshCourantNo.H"
 
         if
         (

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

@@ -62,6 +62,7 @@ int main(int argc, char *argv[])
         mesh.update();
 
 #       include "checkVolContinuity.H"
+#       include "meshCourantNo.H"
 
         if (runTime.timeIndex() % checkFrequency == 0)
         {

applications/utilities/mesh/manipulation/refineMesh/refineMesh.C

@@ -78,7 +78,7 @@ void printEdgeStats(const primitiveMesh& mesh)
 
     const edgeList& edges = mesh.edges();
 
-    forAll(edges, edgeI)
+    forAll (edges, edgeI)
     {
         const edge& e = edges[edgeI];
 
@@ -88,19 +88,19 @@ void printEdgeStats(const primitiveMesh& mesh)
 
         eVec /= eMag;
 
-        if (mag(eVec & x) > 1-edgeTol)
+        if (mag(eVec & x) > 1 - edgeTol)
         {
             minX = min(minX, eMag);
             maxX = max(maxX, eMag);
             nX++;
         }
-        else if (mag(eVec & y) > 1-edgeTol)
+        else if (mag(eVec & y) > 1 - edgeTol)
         {
             minY = min(minY, eMag);
             maxY = max(maxY, eMag);
             nY++;
         }
-        else if (mag(eVec & z) > 1-edgeTol)
+        else if (mag(eVec & z) > 1 - edgeTol)
         {
             minZ = min(minZ, eMag);
             maxZ = max(maxZ, eMag);
@@ -131,15 +131,15 @@ label axis(const vector& normal)
 {
     label axisIndex = -1;
 
-    if (mag(normal & point(1, 0, 0)) > (1-edgeTol))
+    if (mag(normal & point(1, 0, 0)) > (1 - edgeTol))
     {
         axisIndex = 0;
     }
-    else if (mag(normal & point(0, 1, 0)) > (1-edgeTol))
+    else if (mag(normal & point(0, 1, 0)) > (1 - edgeTol))
     {
         axisIndex = 1;
     }
-    else if (mag(normal & point(0, 0, 1)) > (1-edgeTol))
+    else if (mag(normal & point(0, 0, 1)) > (1 - edgeTol))
     {
         axisIndex = 2;
     }
@@ -194,13 +194,13 @@ label twoDNess(const polyMesh& mesh)
     plane cellPlane(ctrs[0], ctrs[1], ctrs[otherCellI]);
 
 
-    forAll(ctrs, cellI)
+    forAll (ctrs, cellI)
     {
         const labelList& cEdges = mesh.cellEdges()[cellI];
 
         scalar minLen = GREAT;
 
-        forAll(cEdges, i)
+        forAll (cEdges, i)
         {
             minLen = min(minLen, mesh.edges()[cEdges[i]].mag(mesh.points()));
         }
@@ -230,15 +230,15 @@ label twoDNess(const polyMesh& mesh)
     // Mark boundary points
     boolList boundaryPoint(mesh.allPoints().size(), false);
 
-    forAll(patches, patchI)
+    forAll (patches, patchI)
     {
         const polyPatch& patch = patches[patchI];
 
-        forAll(patch, patchFaceI)
+        forAll (patch, patchFaceI)
         {
             const face& f = patch[patchFaceI];
 
-            forAll(f, fp)
+            forAll (f, fp)
             {
                 boundaryPoint[f[fp]] = true;
             }
@@ -248,7 +248,7 @@ label twoDNess(const polyMesh& mesh)
 
     const edgeList& edges = mesh.edges();
 
-    forAll(edges, edgeI)
+    forAll (edges, edgeI)
     {
         const edge& e = edges[edgeI];
 
@@ -263,7 +263,7 @@ label twoDNess(const polyMesh& mesh)
     // 3. For all non-wedge patches: all faces either perp or aligned with
     //    cell-plane normal. (wedge patches already checked upon construction)
 
-    forAll(patches, patchI)
+    forAll (patches, patchI)
     {
         const polyPatch& patch = patches[patchI];
 
@@ -349,7 +349,9 @@ int main(int argc, char *argv[])
         // Select all cells
         refCells.setSize(mesh.nCells());
 
-        forAll(mesh.cells(), cellI)
+        const cellList& c = mesh.cells();
+
+        forAll (c, cellI)
         {
             refCells[cellI] = cellI;
         }
@@ -440,11 +442,11 @@ int main(int argc, char *argv[])
     // Create cellSet with added cells for easy inspection
     cellSet newCells(mesh, "refinedCells", refCells.size());
 
-    forAll(oldToNew, oldCellI)
+    forAll (oldToNew, oldCellI)
     {
         const labelList& added = oldToNew[oldCellI];
 
-        forAll(added, i)
+        forAll (added, i)
         {
             newCells.insert(added[i]);
         }
@@ -482,14 +484,13 @@ int main(int argc, char *argv[])
       + " to cells in mesh at "
       + oldTimeName;
 
-
-    forAll(oldToNew, oldCellI)
+    forAll (oldToNew, oldCellI)
     {
         const labelList& added = oldToNew[oldCellI];
 
         if (added.size())
         {
-            forAll(added, i)
+            forAll (added, i)
             {
                 newToOld[added[i]] = oldCellI;
             }
@@ -506,7 +507,6 @@ int main(int argc, char *argv[])
 
     newToOld.write();
 
-
     // Some statistics.
 
     printEdgeStats(mesh);

applications/utilities/postProcessing/stressField/solidStress/solidStress.C

@@ -36,7 +36,6 @@ Description
 
 int main(int argc, char *argv[])
 {
-
 #   include "addTimeOptions.H"
 #   include "setRootCase.H"
 
@@ -53,7 +52,7 @@ int main(int argc, char *argv[])
 #   include "createMesh.H"
 #   include "readMechanicalProperties.H"
 
-    for (label i=startTime; i<endTime; i++)
+    for (label i = startTime; i < endTime; i++)
     {
         runTime.setTime(Times[i], i);
 

applications/utilities/postProcessing/velocityField/Mach/Mach.C

@@ -36,6 +36,7 @@ Description
 #include "basicPsiThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
 {
     bool writeResults = !args.optionFound("noWrite");
@@ -134,8 +135,6 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
     {
         Info<< "    Missing U or T" << endl;
     }
-
-    Info<< "\nEnd\n" << endl;
 }
 
 

src/ODE/sixDOF/finiteRotation/finiteRotation.C

@@ -41,7 +41,18 @@ Foam::vector Foam::finiteRotation::rotVector(const tensor& rotT)
 
     // Scaling to a unit vector.  HJ, problems with round-off
     // HJ, 4/Aug/2008
+
+    if (mag(ur) > SMALL)
+    {
         return ur/(mag(ur) + SMALL);
+    }
+    else
+    {
+        // Rotation vector is undertermined at zero rotation
+        // Returning arbitrary unit vector
+        // HJ, 4/Mar/2015
+        return vector(0, 0, 1);
+    }
 }
 
 

src/ODE/sixDOF/finiteRotation/finiteRotation.H

@@ -52,7 +52,7 @@ class finiteRotation
     // Private data
 
         //- Initial rotation
-        const HamiltonRodriguezRot eInitial_;
+        HamiltonRodriguezRot eInitial_;
 
         //- Rotational tensor
         tensor rotTensor_;
@@ -63,13 +63,6 @@ class finiteRotation
 
     // Private Member Functions
 
-        //- Disallow default bitwise copy construct
-        finiteRotation(const finiteRotation&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const finiteRotation&);
-
-
         //- Calculate unit rotation vector from given rotation tensor
         static vector rotVector(const tensor& rotT);
 

src/ODE/sixDOF/quaternions/HamiltonRodriguezRot.H

@@ -74,17 +74,17 @@ class HamiltonRodriguezRot
         //- Calculate R_ - inertial to body coord. sys. rotation
         inline void calculateR() const
         {
-            R_.xx() = 2 * (e1_*e1_ + e0_*e0_ - 0.5);
-            R_.xy() = 2 * (e1_*e2_ + e0_*e3_);
-            R_.xz() = 2 * (e1_*e3_ - e0_*e2_);
+            R_.xx() = 2*(e1_*e1_ + e0_*e0_ - 0.5);
+            R_.xy() = 2*(e1_*e2_ + e0_*e3_);
+            R_.xz() = 2*(e1_*e3_ - e0_*e2_);
 
-            R_.yx() = 2 * (e2_*e1_ - e0_*e3_);
-            R_.yy() = 2 * (e2_*e2_ + e0_*e0_ - 0.5);
-            R_.yz() = 2 * (e2_*e3_ + e0_*e1_);
+            R_.yx() = 2*(e2_*e1_ - e0_*e3_);
+            R_.yy() = 2*(e2_*e2_ + e0_*e0_ - 0.5);
+            R_.yz() = 2*(e2_*e3_ + e0_*e1_);
 
-            R_.zx() = 2 * (e3_*e1_ + e0_*e2_);
-            R_.zy() = 2 * (e3_*e2_ - e0_*e1_);
-            R_.zz() = 2 * (e3_*e3_ + e0_*e0_ - 0.5);
+            R_.zx() = 2*(e3_*e1_ + e0_*e2_);
+            R_.zy() = 2*(e3_*e2_ - e0_*e1_);
+            R_.zz() = 2*(e3_*e3_ + e0_*e0_ - 0.5);
         }
 
         //- Calculate Gt()

src/ODE/sixDOF/sixDOFbodies/sixDOFbodies.C

@@ -30,7 +30,7 @@ Description
 Author
     Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
 
-\*----------------------------------------------------------------------------*/
+\*---------------------------------------------------------------------------*/
 
 #include "objectRegistry.H"
 #include "sixDOFbodies.H"

src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C

@@ -35,17 +35,6 @@ Author
 #include "objectRegistry.H"
 #include "sixDOFqODE.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-// Runtime type information
-// Not possible because of I/O error: incorrect type, expecting dictionary
-// HJ, 11/Feb/2008
-// namespace Foam
-// {
-//     defineTypeNameAndDebug(sixDOFqODE, 0);
-// }
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::sixDOFqODE::setCoeffs()
@@ -155,6 +144,49 @@ Foam::sixDOFqODE::sixDOFqODE(const IOobject& io)
 }
 
 
+// Construct as copy
+Foam::sixDOFqODE::sixDOFqODE
+(
+    const word& name,
+    const sixDOFqODE& sd
+)
+:
+    IOdictionary
+    (
+        IOobject
+        (
+            name,
+            sd.instance(),
+            sd.local(),
+            sd.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    mass_(sd.mass_.name(), sd.mass_),
+    momentOfInertia_(sd.momentOfInertia_.name(), sd.momentOfInertia_),
+    Xequilibrium_(sd.Xequilibrium_.name(), sd.Xequilibrium_),
+    linSpringCoeffs_(sd.linSpringCoeffs_.name(), sd.linSpringCoeffs_),
+    linDampingCoeffs_(sd.linDampingCoeffs_.name(), sd.linDampingCoeffs_),
+    Xrel_(sd.Xrel_.name(), sd.Xrel_),
+    U_(sd.U_.name(), sd.U_),
+    Uaverage_(sd.Uaverage_.name(), sd.Uaverage_),
+    rotation_(sd.rotation_),
+    omega_(sd.omega_.name(), sd.omega_),
+    omegaAverage_(sd.omegaAverage_.name(), sd.omegaAverage_),
+    omegaAverageAbsolute_
+    (
+        sd.omegaAverageAbsolute_.name(),
+        sd.omegaAverageAbsolute_
+    ),
+    force_(sd.force_.name(), sd.force_),
+    moment_(sd.moment_.name(), sd.moment_),
+    forceRelative_(sd.forceRelative_.name(), sd.forceRelative_),
+    momentRelative_(sd.momentRelative_.name(), sd.momentRelative_),
+    coeffs_(sd.coeffs_)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::sixDOFqODE::~sixDOFqODE()
@@ -163,6 +195,31 @@ Foam::sixDOFqODE::~sixDOFqODE()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::sixDOFqODE::setState(const sixDOFqODE& sd)
+{
+    mass_ = sd.mass_;
+    momentOfInertia_ = sd.momentOfInertia_;
+    Xequilibrium_ = sd.Xequilibrium_;
+    linSpringCoeffs_ = sd.linSpringCoeffs_;
+    linDampingCoeffs_ = sd.linDampingCoeffs_;
+    Xrel_ = sd.Xrel_;
+    U_ = sd.U_;
+    Uaverage_ = sd.Uaverage_;
+    rotation_ = sd.rotation_;
+    omega_ = sd.omega_;
+    omegaAverage_ = sd.omegaAverage_;
+    omegaAverageAbsolute_ = sd.omegaAverageAbsolute_;
+    force_ = sd.force_;
+    moment_ = sd.moment_;
+    forceRelative_ = sd.forceRelative_;
+    momentRelative_ = sd.momentRelative_;
+
+    // Copy ODE coefficients: this carries actual ODE state
+    // HJ, 23/Mar/2015
+    coeffs_ = sd.coeffs_;
+}
+
+
 void Foam::sixDOFqODE::derivatives
 (
     const scalar x,

src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H

@@ -173,18 +173,18 @@ class sixDOFqODE
 
 public:
 
-
-//     //- Runtime type information
-//       Not possible because of I/O error: incorrect type, expecting dictionary
-//       HJ, 11/Feb/2008
-//     TypeName("sixDOFqODE");
-
-
     // Constructors
 
         //- Construct from dictionary
         sixDOFqODE(const IOobject& io);
 
+        //- Construct sixDOFqODE, changing name
+        sixDOFqODE
+        (
+            const word& name,
+            const sixDOFqODE& sd
+        );
+
 
     // Destructor
 
@@ -246,14 +246,21 @@ public:
             inline dimensionedScalar rotAngle() const;
 
 
-        // Non-constant access to velocities
+        // Non-constant access to 
 
-            //- Return access to velocity of origin
-            inline dimensionedVector& U();
+            //- Set position of origin
+            inline void setXrel(const vector& x);
 
-            //- Return access to rotational velocity in relative
+            //- Set velocity of origin
+            inline void setU(const vector& u);
+
+            //- Set rotational angle in relative
             //  coordinate system
-            inline dimensionedVector& omega();
+            inline void setRotation(const HamiltonRodriguezRot& rot);
+
+            //- Set rotational velocity in relative
+            //  coordinate system
+            inline void setOmega(const vector& omega);
 
 
         // Average motion per time-step
@@ -261,10 +268,12 @@ public:
             //- Return average rotational vector of body
             inline vector rotVectorAverage() const;
 
-            //- Return average rotational velocity in relative coordinate system
+            //- Return average rotational velocity in
+            //  relative coordinate system
             inline const dimensionedVector& omegaAverage() const;
 
-            //- Return average rotational velocity in absolute coordinate system
+            //- Return average rotational velocity in
+            //  absolute coordinate system
             inline const dimensionedVector& omegaAverageAbsolute() const;
 
 
@@ -313,6 +322,9 @@ public:
             //- Return transformation tensor between new and previous rotation
             inline const tensor& rotIncrementTensor() const;
 
+            //-  Set ODE parameters from another ODE
+            void setState(const sixDOFqODE&);
+
 
         // ODE parameters
 

src/ODE/sixDOF/sixDOFqODE/sixDOFqODEI.H

@@ -118,15 +118,51 @@ Foam::dimensionedScalar Foam::sixDOFqODE::rotAngle() const
 }
 
 
-Foam::dimensionedVector& Foam::sixDOFqODE::U()
+void Foam::sixDOFqODE::setXrel(const vector& x)
 {
-    return U_;
+    Xrel_.value() = x;
+
+    // Set X in coefficients
+    coeffs_[0] = x.x();
+    coeffs_[1] = x.y();
+    coeffs_[2] = x.z();
 }
 
 
-Foam::dimensionedVector& Foam::sixDOFqODE::omega()
+void Foam::sixDOFqODE::setU(const vector& U)
 {
-    return omega_;
+    U_.value() = U;
+    Uaverage_ = U_;
+
+    // Set U in coefficients
+    coeffs_[3] = U.x();
+    coeffs_[4] = U.y();
+    coeffs_[5] = U.z();
+}
+
+
+void Foam::sixDOFqODE::setRotation(const HamiltonRodriguezRot& rot)
+{
+    // Cannot call update rotation: simply set up coefficients
+    // HJ, 23/Mar/2015
+
+    coeffs_[9] = rot.e0();
+    coeffs_[10] = rot.e1();
+    coeffs_[11] = rot.e2();
+    coeffs_[12] = rot.e3();
+}
+
+
+void Foam::sixDOFqODE::setOmega(const vector& omega)
+{
+    omega_.value() = omega;
+    omegaAverage_ = omega_;
+    omegaAverageAbsolute_ = omega_;
+
+    // Set omega in coefficients
+    coeffs_[6] = omega.x();
+    coeffs_[7] = omega.y();
+    coeffs_[8] = omega.z();
 }
 
 

src/coupledMatrix/coupledLduMatrix/coupledLduMatrix.C

@@ -30,9 +30,10 @@ Description
 Author
     Hrvoje Jasak, Wikki Ltd.  All rights reserved
 
-\*----------------------------------------------------------------------------*/
+\*---------------------------------------------------------------------------*/
 
 #include "coupledLduMatrix.H"
+#include "processorLduInterfaceField.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -42,9 +43,6 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct given size
@@ -202,7 +200,54 @@ void Foam::coupledLduMatrix::initMatrixInterfaces
 {
     const PtrList<lduMatrix>& matrices = *this;
 
+    // Note.  The comms design requires all non-processor interfaces
+    // to be updated first, followed by the update of processor
+    // interfaces.  The reason is that non-processor coupled
+    // interfaces require a complex comms pattern involving more than
+    // pairwise communications.
+    // Under normal circumstances this is achieved naturall, since
+    // processor interfaces come last on the list and other coupled
+    // interfaces execute complex comms at init() level.
+    // For coupled matrices, the update loop needs to be split over
+    // all matrices by hand
+    // Bug fix: Zeljko Tukovic, 7/Apr/2015
+
+    // Init update all non-processor coupled interfaces
     forAll (matrices, rowI)
+    {
+        if
+        (
+            Pstream::defaultCommsType == Pstream::blocking
+         || Pstream::defaultCommsType == Pstream::nonBlocking
+        )
+        {
+            forAll (interfaces[rowI], interfaceI)
+            {
+                if (interfaces[rowI].set(interfaceI))
+                {
+                    if
+                    (
+                       !isA<processorLduInterfaceField>
+                        (
+                            interfaces[rowI][interfaceI]
+                        )
+                    )
+                    {
+                        interfaces[rowI][interfaceI].initInterfaceMatrixUpdate
+                        (
+                            x[rowI],
+                            result[rowI],
+                            matrices[rowI],
+                            coupleCoeffs[rowI][interfaceI],
+                            cmpt,
+                            Pstream::defaultCommsType,
+                            false
+                        );
+                    }
+                }
+            }
+        }
+        else
         {
             matrices[rowI].initMatrixInterfaces
             (
@@ -213,6 +258,44 @@ void Foam::coupledLduMatrix::initMatrixInterfaces
                 cmpt
             );
         }
+    }
+
+    // Init update for all processor interfaces
+    forAll (matrices, rowI)
+    {
+        if
+        (
+            Pstream::defaultCommsType == Pstream::blocking
+         || Pstream::defaultCommsType == Pstream::nonBlocking
+        )
+        {
+            forAll (interfaces[rowI], interfaceI)
+            {
+                if (interfaces[rowI].set(interfaceI))
+                {
+                    if
+                    (
+                        isA<processorLduInterfaceField>
+                        (
+                            interfaces[rowI][interfaceI]
+                        )
+                    )
+                    {
+                        interfaces[rowI][interfaceI].initInterfaceMatrixUpdate
+                        (
+                            x[rowI],
+                            result[rowI],
+                            matrices[rowI],
+                            coupleCoeffs[rowI][interfaceI],
+                            cmpt,
+                            Pstream::defaultCommsType,
+                            false
+                        );
+                    }
+                }
+            }
+        }
+    }
 }
 
 
@@ -241,13 +324,4 @@ void Foam::coupledLduMatrix::updateMatrixInterfaces
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //

src/dynamicMesh/dynamicFvMesh/dynamicBodyFvMesh/dynamicBodyFvMesh.C

@@ -98,7 +98,7 @@ Foam::dynamicBodyFvMesh::dynamicBodyFvMesh(const IOobject& io)
 {
     bodyPatchID_ = boundaryMesh().findPatchID(bodyPatchName_);
 
-    if(bodyPatchID_<0)
+    if (bodyPatchID_<0)
     {
         FatalErrorIn
         (

src/dynamicMesh/dynamicMesh/polyMeshModifiers/slidingInterface/slidingInterfaceProjectPoints.C

@@ -297,7 +297,10 @@ bool Foam::slidingInterface::projectPoints() const
     }
 
     // Projected slave points are stored for mesh motion correction
-    if (projectedSlavePointsPtr_) delete projectedSlavePointsPtr_;
+    if (projectedSlavePointsPtr_)
+    {
+        delete projectedSlavePointsPtr_;
+    }
 
     projectedSlavePointsPtr_ =
         new pointField(slavePointFaceHits.size(), vector::zero);

src/dynamicMesh/topoChangerFvMesh/mixerFvMesh/mixerFvMesh.C

@@ -467,7 +467,7 @@ bool Foam::mixerFvMesh::update()
 
     if (globalMorphing)
     {
-        Info << "Attaching rotors" << endl;
+        Info<< "Attaching rotors" << endl;
 
         deleteDemandDrivenData(movingPointsMaskPtr_);
 

src/finiteArea/faMesh/faPatches/faPatch/faPatch.C

@@ -309,14 +309,15 @@ Foam::labelList Foam::faPatch::ngbPolyPatchFaces() const
 
 Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchFaceNormals() const
 {
-    tmp<vectorField> fN(new vectorField());
+    tmp<vectorField> tfN(new vectorField());
+    vectorField& fN = tfN();
 
     if (ngbPolyPatchIndex() == -1)
     {
-        return fN;
+        return tfN;
     }
 
-    fN().setSize(faPatch::size());
+    fN.setSize(faPatch::size());
 
     labelList ngbFaces = ngbPolyPatchFaces();
 
@@ -325,13 +326,13 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchFaceNormals() const
     const faceList& faces = pMesh.faces();
     const pointField& points = pMesh.points();
 
-    forAll(fN(), faceI)
+    forAll(fN, faceI)
     {
-        fN() = faces[ngbFaces[faceI]].normal(points)
+        fN[faceI] = faces[ngbFaces[faceI]].normal(points)
             /faces[ngbFaces[faceI]].mag(points);
     }
 
-    return fN;
+    return tfN;
 }
 
 
@@ -344,21 +345,22 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchPointNormals() const
 
     labelListList pntEdges = pointEdges();
 
-    tmp<vectorField> pN(new vectorField(pntEdges.size(), vector::zero));
+    tmp<vectorField> tpN(new vectorField(pntEdges.size(), vector::zero));
+    vectorField& pN = tpN();
 
     vectorField faceNormals = ngbPolyPatchFaceNormals();
 
-    forAll(pN(), pointI)
+    forAll(pN, pointI)
     {
         forAll(pntEdges[pointI], edgeI)
         {
-            pN()[pointI] += faceNormals[pntEdges[pointI][edgeI]];
+            pN[pointI] += faceNormals[pntEdges[pointI][edgeI]];
         }
     }
 
-    pN() /= mag(pN());
+    pN /= mag(pN);
 
-    return pN;
+    return tpN;
 }
 
 

src/finiteArea/fields/faPatchFields/basic/fixedGradient/fixedGradientFaPatchField.H

@@ -166,7 +166,7 @@ public:
             //- Evaluate the patch field
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             );
 
 

src/finiteArea/fields/faPatchFields/basic/mixed/mixedFaPatchField.H

@@ -199,7 +199,7 @@ public:
             //- Evaluate the patch field
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             );
 
             //- Return the matrix diagonal coefficients corresponding to the

src/finiteArea/fields/faPatchFields/basic/zeroGradient/zeroGradientFaPatchField.H

@@ -135,7 +135,7 @@ public:
             //- Evaluate the patch field
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             );
 
             //- Return the matrix diagonal coefficients corresponding to the

src/finiteArea/fields/faPatchFields/faPatchField/faPatchField.H

@@ -351,14 +351,14 @@ public:
             //- Initialise the evaluation of the patch field
             virtual void initEvaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             )
             {}
 
             //- Evaluate the patch field, sets Updated to false
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             );
 
 

src/finiteVolume/fields/fvPatchFields/basic/directionMixed/directionMixedFvPatchField.H

@@ -195,7 +195,7 @@ public:
             //- Evaluate the patch field
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             );
 
             //- Return face-gradient transform diagonal

src/finiteVolume/fields/fvPatchFields/basic/fixedGradient/fixedGradientFvPatchField.H

@@ -167,7 +167,7 @@ public:
             //- Evaluate the patch field
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             );
 
             //- Return the matrix diagonal coefficients corresponding to the

src/finiteVolume/fields/fvPatchFields/basic/mixed/mixedFvPatchField.H

@@ -200,7 +200,7 @@ public:
             //- Evaluate the patch field
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             );
 
             //- Return the matrix diagonal coefficients corresponding to the

src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.H

@@ -162,14 +162,14 @@ public:
             //- Initialise the evaluation of the patch field
             virtual void initEvaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             )
             {}
 
             //- Evaluate the patch field, sets Updated to false
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             )
             {}
 

src/finiteVolume/fields/fvPatchFields/basic/zeroGradient/zeroGradientFvPatchField.H

@@ -136,7 +136,7 @@ public:
             //- Evaluate the patch field
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             );
 
             //- Return the matrix diagonal coefficients corresponding to the

src/finiteVolume/fields/fvPatchFields/constraint/wedge/wedgeFvPatchField.H

@@ -130,7 +130,7 @@ public:
             //- Evaluate the patch field
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             );
 
             //- Return face-gradient transform diagonal

src/finiteVolume/fields/fvPatchFields/derived/fixedNormalSlip/fixedNormalSlipFvPatchField.H

@@ -164,7 +164,7 @@ public:
             //- Evaluate the patch field
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             );
 
             //- Return face-gradient transform diagonal

src/finiteVolume/fields/fvPatchFields/derived/fluxCorrectedVelocity/fluxCorrectedVelocityFvPatchVectorField.H

@@ -132,7 +132,7 @@ public:
             //- Evaluate the patch field
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             );
 
 

src/finiteVolume/fields/fvPatchFields/derived/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.C

@@ -109,7 +109,8 @@ void translatingWallVelocityFvPatchVectorField::updateCoeffs()
         return;
     }
 
-    // Remove the component of U normal to the wall in case the wall is not flat
+    // Remove the component of U normal to the wall in case the wall
+    // is not flat
     vectorField n = patch().nf();
     vectorField::operator=(U_ - n*(n & U_));
 

src/finiteVolume/fields/fvPatchFields/derived/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.H

@@ -52,7 +52,7 @@ class translatingWallVelocityFvPatchVectorField
 {
     // Private data
 
-        //- Origin of the rotation
+        //- Translating velocity
         vector U_;
 
 

src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H

@@ -332,14 +332,14 @@ public:
             //- Initialise the evaluation of the patch field
             virtual void initEvaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             )
             {}
 
             //- Evaluate the patch field, sets Updated to false
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             )
             {}
 

src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C

@@ -88,6 +88,7 @@ void Foam::fvc::makeRelative
     }
 }
 
+
 void Foam::fvc::makeRelative
 (
     surfaceScalarField& phi,
@@ -101,6 +102,7 @@ void Foam::fvc::makeRelative
     }
 }
 
+
 void Foam::fvc::makeRelative
 (
     surfaceScalarField& phi,
@@ -127,6 +129,7 @@ void Foam::fvc::makeAbsolute
     }
 }
 
+
 void Foam::fvc::makeAbsolute
 (
     surfaceScalarField& phi,
@@ -140,6 +143,7 @@ void Foam::fvc::makeAbsolute
     }
 }
 
+
 void Foam::fvc::makeAbsolute
 (
     surfaceScalarField& phi,

src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C

@@ -101,19 +101,15 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
     // Set local references to mesh data
     const unallocLabelList& owner = mesh().owner();
     const unallocLabelList& neighbour = mesh().neighbour();
-
     const volVectorField& C = mesh().C();
-    const surfaceScalarField& w = mesh().weights();
-//     const surfaceScalarField& magSf = mesh().magSf();
-
 
     // Set up temporary storage for the dd tensor (before inversion)
     symmTensorField dd(mesh().nCells(), symmTensor::zero);
 
-    forAll(owner, facei)
+    forAll(owner, faceI)
     {
-        label own = owner[facei];
-        label nei = neighbour[facei];
+        label own = owner[faceI];
+        label nei = neighbour[faceI];
 
         vector d = C[nei] - C[own];
         symmTensor wdd = (1.0/magSqr(d))*sqr(d);
@@ -122,113 +118,92 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
         dd[nei] += wdd;
     }
 
-
-    forAll(lsP.boundaryField(), patchi)
+    forAll(lsP.boundaryField(), patchI)
     {
-        const fvsPatchScalarField& pw = w.boundaryField()[patchi];
-        // Note: least squares in 1.4.1 and other OpenCFD versions
-        // are wrong because of incorrect weighting.  HJ, 23/Oct/2008
-//         const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
-
-        const fvPatch& p = pw.patch();
+        const fvPatch& p = mesh().boundary()[patchI];
         const unallocLabelList& faceCells = p.patch().faceCells();
 
-        // Build the d-vectors
-
-        // Original version: closest distance to boundary
-        vectorField pd =
-            mesh().Sf().boundaryField()[patchi]
-           /(
-               mesh().magSf().boundaryField()[patchi]
-              *mesh().deltaCoeffs().boundaryField()[patchi]
-           );
-
-        if (!mesh().orthogonal())
-        {
-            pd -= mesh().correctionVectors().boundaryField()[patchi]
-                /mesh().deltaCoeffs().boundaryField()[patchi];
-        }
-
         // Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010
-        // Experimental: review fixed gradient condition.  HJ, 30/Sep/2010
-//         vectorField pd = p.delta();
+        const vectorField pd = p.delta();
 
-        if (p.coupled())
+        forAll(pd, patchFaceI)
         {
-            forAll(pd, patchFacei)
-            {
-                const vector& d = pd[patchFacei];
-
-                dd[faceCells[patchFacei]] += (1.0/magSqr(d))*sqr(d);
-            }
-        }
-        else
-        {
-            forAll(pd, patchFacei)
-            {
-                const vector& d = pd[patchFacei];
-
-                dd[faceCells[patchFacei]] += (1.0/magSqr(d))*sqr(d);
-            }
+            const vector& d = pd[patchFaceI];
+            dd[faceCells[patchFaceI]] += (1.0/magSqr(d))*sqr(d);
         }
     }
 
+    // For easy access of neighbour coupled patch field needed for
+    // lsN vectors on implicitly coupled boundaries. VV, 18/June/2014
+    volSymmTensorField volInvDd
+    (
+        IOobject
+        (
+            "volInvDd",
+            mesh().pointsInstance(),
+            mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh(),
+        dimensionedSymmTensor("zero", dimless, symmTensor::zero),
+        "zeroGradient"
+    );
+    symmTensorField& invDd = volInvDd.internalField();
+//    invDd = inv(dd);
+    invDd = hinv(dd);
 
-    // Invert the dd tensor
-//     symmTensorField invDd = inv(dd);
-    // Fix: householder inverse.  HJ, 3/Nov/2009
-    symmTensorField invDd = hinv(dd);
-
+    // Evaluate coupled to exchange coupled neighbour field data
+    // across coupled boundaries.  HJ, 18/Mar/2015
+    volInvDd.boundaryField().evaluateCoupled();
 
     // Revisit all faces and calculate the lsP and lsN vectors
-    forAll(owner, facei)
+    forAll(owner, faceI)
     {
-        label own = owner[facei];
-        label nei = neighbour[facei];
+        label own = owner[faceI];
+        label nei = neighbour[faceI];
 
         vector d = C[nei] - C[own];
         scalar magSfByMagSqrd = 1.0/magSqr(d);
 
-        lsP[facei] = magSfByMagSqrd*(invDd[own] & d);
-        lsN[facei] = -magSfByMagSqrd*(invDd[nei] & d);
+        lsP[faceI] = magSfByMagSqrd*(invDd[own] & d);
+        lsN[faceI] = -magSfByMagSqrd*(invDd[nei] & d);
     }
 
-    forAll(lsP.boundaryField(), patchi)
+    forAll(lsP.boundaryField(), patchI)
     {
-        fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchi];
-
-        const fvsPatchScalarField& pw = w.boundaryField()[patchi];
-        // Note: least squares in 1.4.1 and other OpenCFD versions
-        // are wrong because of incorrect weighting.  HJ, 23/Oct/2008
-//         const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
-
-        const fvPatch& p = pw.patch();
+        fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchI];
+        fvsPatchVectorField& patchLsN = lsN.boundaryField()[patchI];
+        const fvPatch& p = mesh().boundary()[patchI];
         const unallocLabelList& faceCells = p.faceCells();
 
-        // Build the d-vectors
         // Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010
-        vectorField pd = p.delta();
+        const vectorField pd = p.delta();
 
         if (p.coupled())
         {
-            forAll(pd, patchFacei)
-            {
-                const vector& d = pd[patchFacei];
+            const symmTensorField invDdNei =
+                volInvDd.boundaryField()[patchI].patchNeighbourField();
 
-                patchLsP[patchFacei] =
-                    (1.0/magSqr(d))
-                   *(invDd[faceCells[patchFacei]] & d);
+            forAll(pd, patchFaceI)
+            {
+                const vector& d = pd[patchFaceI];
+
+                patchLsP[patchFaceI] = (1.0/magSqr(d))
+                   *(invDd[faceCells[patchFaceI]] & d);
+
+                patchLsN[patchFaceI] = - (1.0/magSqr(d))
+                   *(invDdNei[patchFaceI] & d);
             }
         }
         else
         {
-            forAll(pd, patchFacei)
+            forAll(pd, patchFaceI)
             {
-                const vector& d = pd[patchFacei];
+                const vector& d = pd[patchFaceI];
 
-                patchLsP[patchFacei] =
-                    (1.0/magSqr(d))
-                   *(invDd[faceCells[patchFacei]] & d);
+                patchLsP[patchFaceI] = (1.0/magSqr(d))
+                   *(invDd[faceCells[patchFaceI]] & d);
             }
         }
     }
@@ -278,12 +253,15 @@ bool Foam::leastSquaresVectors::movePoints() const
     return true;
 }
 
+
 bool Foam::leastSquaresVectors::updateMesh(const mapPolyMesh& mpm) const
 {
     if (debug)
     {
-        InfoIn("bool leastSquaresVectors::updateMesh(const mapPolyMesh&) const")
-            << "Clearing least square data" << endl;
+        InfoIn
+        (
+            "bool leastSquaresVectors::updateMesh(const mapPolyMesh&) const"
+        )   << "Clearing least square data" << endl;
     }
 
     deleteDemandDrivenData(pVectorsPtr_);

src/foam/coordinateSystems/coordinateRotation/coordinateRotation.H

@@ -116,6 +116,7 @@ public:
     //- Runtime type information
     TypeName("coordinateRotation");
 
+
     // Constructors
 
         //- Construct null

src/foam/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.C

@@ -76,6 +76,7 @@ void Foam::OutputFilterFunctionObject<OutputFilter>::destroyFilter()
     ptr_.reset();
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class OutputFilter>

src/foam/fields/PointPatchFields/constraint/processor/ProcessorPointPatchField.C

@@ -476,7 +476,7 @@ initEvaluate
 }
 
 
-// Initialise field transfer
+// Field transfer
 template
 <
     template<class> class PatchField,
@@ -511,7 +511,18 @@ evaluate
             }
 
             // Average over two sides
-            tpn = 0.5*(this->patchInternalField(this->internalField()) + tpn);
+
+            // ZT, 22-07-2014 - point ordering is not same
+            // at master and slave side after topology change
+            const labelList& neiPoints =
+                procPatch_.procPolyPatch().neighbPoints();
+
+            tpn =
+                0.5*
+                (
+                   this->patchInternalField(this->internalField())
+                 + Field<Type>(tpn, neiPoints)
+                );
 
             // Get internal field to insert values into
             Field<Type>& iF = const_cast<Field<Type>&>(this->internalField());

src/foam/fields/PointPatchFields/derived/componentMixed/ComponentMixedPointPatchVectorField.H

@@ -196,7 +196,7 @@ public:
             //- Insert boundary value into the internal field
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType = Pstream::blocking
             );
 
             //- Set boundary condition to matrix

src/foam/interpolations/GGIInterpolation/GGIInterpolationPolygonIntersection.C

@@ -192,7 +192,7 @@ GGIInterpolation<MasterPatch, SlavePatch>::polygonIntersection
             (
                 "GGIInterpolation<MasterPatch, SlavePatch>::"
                 "polygonIntersection"
-            )   << "Intersection might be wrong wrong: clipping side "
+            )   << "Intersection might be wrong.  Clipping side "
                 << intersectionArea/clippingArea << " subject: "
                 << intersectionArea/subjectArea << endl;
         }

src/foam/interpolations/GGIInterpolation/GGIInterpolationQuickRejectTests.C

@@ -52,7 +52,7 @@ GGIInterpolation<MasterPatch, SlavePatch>::faceBoundBoxExtendSpanFraction_
 
 template<class MasterPatch, class SlavePatch>
 const label
- GGIInterpolation<MasterPatch, SlavePatch>::octreeSearchMinNLevel_
+GGIInterpolation<MasterPatch, SlavePatch>::octreeSearchMinNLevel_
 (
     debug::optimisationSwitch("GGIOctreeSearchMinNLevel", 3)
 );

src/foam/meshes/polyMesh/globalMeshData/globalMeshData.H

@@ -74,7 +74,7 @@ SourceFiles
 #include "polyMesh.H"
 #include "Switch.H"
 #include "processorTopology.H"
-#include "labelPair.H"
+#include "labelPairList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -119,9 +119,6 @@ class globalMeshData
         };
 
 
-        typedef List<labelPair> labelPairList;
-
-
     // Private data
 
         //- Reference to mesh

src/foam/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C

@@ -168,7 +168,7 @@ void Foam::polyBoundaryMesh::calcGeometry()
 }
 
 
-const Foam::List<Foam::labelPairList>&
+const Foam::labelPairListList&
 Foam::polyBoundaryMesh::neighbourEdges() const
 {
     if (Pstream::parRun())
@@ -180,8 +180,8 @@ Foam::polyBoundaryMesh::neighbourEdges() const
 
     if (!neighbourEdgesPtr_)
     {
-        neighbourEdgesPtr_ = new List<labelPairList>(size());
-        List<labelPairList>& neighbourEdges = *neighbourEdgesPtr_;
+        neighbourEdgesPtr_ = new labelPairListList(size());
+        labelPairListList& neighbourEdges = *neighbourEdgesPtr_;
 
         // Initialize.
         label nEdgePairs = 0;

src/foam/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H

@@ -37,14 +37,13 @@ SourceFiles
 
 #include "polyPatchList.H"
 #include "regIOobject.H"
-#include "labelPair.H"
+#include "labelPairList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-typedef List<labelPair> labelPairList;
 class polyMesh;
 
 // Forward declaration of friend functions and operators
@@ -69,7 +68,7 @@ class polyBoundaryMesh
         const polyMesh& mesh_;
 
         //- Edges of neighbouring patches
-        mutable List<labelPairList>* neighbourEdgesPtr_;
+        mutable labelPairListList* neighbourEdgesPtr_;
 
 
     // Private Member Functions
@@ -138,13 +137,14 @@ public:
             return mesh_;
         }
 
-        //- Per patch the edges on the neighbouring patch. Is for every external
-        //  edge the neighbouring patch and neighbouring (external) patch edge
-        //  label. Note that edge indices are offset by nInternalEdges to keep
-        //  it as much as possible consistent with coupled patch addressing
+        //- Per patch the edges on the neighbouring patch
+        //  For every external edge the neighbouring patch and
+        //  neighbouring (external) patch edge label.  Note that edge
+        //  indices are offset by nInternalEdges to keep it as much as
+        //  possible consistent with coupled patch addressing
         //  (where coupling is by local patch face index).
         //  Only valid for singly connected polyBoundaryMesh and not parallel
-        const List<labelPairList>& neighbourEdges() const;
+        const labelPairListList& neighbourEdges() const;
 
         //- Return a list of patch names
         wordList names() const;

src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C

@@ -962,6 +962,17 @@ void Foam::ggiPolyPatch::calcTransforms() const
                     << abort(FatalError);
             }
         }
+        else
+        {
+            InfoIn("label ggiPolyPatch::shadowIndex() const")
+                << "ggi patch " << name() << " with shadow "
+                << shadowName() << " has "
+                << patchToPatch().uncoveredMasterFaces().size()
+                << " uncovered master faces and "
+                << patchToPatch().uncoveredSlaveFaces().size()
+                << " uncovered slave faces.  Bridging is switched on. "
+                << abort(FatalError);
+        }
     }
 }
 

src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H

@@ -279,8 +279,7 @@ public:
         }
 
 
-    // Destructor
-
+    //- Destructor
     virtual ~ggiPolyPatch();
 
 

src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatchTemplates.C

@@ -31,7 +31,6 @@ Author
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-// New.  HJ, 12/Jun/2011
 template<class Type>
 Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::fastExpand
 (
@@ -40,22 +39,6 @@ Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::fastExpand
 {
     // Check and expand the field from patch size to zone size
     // with communication
-
-    // Algorithm:
-    // 1) Master processor holds maps of all zone addressing (data provided)
-    // and all remote zone addressing (data required)
-    // 2) Each processor will send the locally active data to the master
-    // 3) Master assembles all the data
-    // 4) Master sends to all processors the data they need to receive
-    //
-    // Notes:
-    // A) If the size of zone addressing is zero, data is not sent
-    // B) Communicated data on each processor has the size of live faces
-    // C) Expanded data will be equal to actual data from other processors
-    //    only for the faces marked in remote; for other faces, it will be
-    //    equal to zero
-    // D) On processor zero, complete data is available
-    // HJ, 4/Jun/2011
     if (ff.size() != size())
     {
         FatalErrorIn
@@ -81,6 +64,22 @@ Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::fastExpand
             << abort(FatalError);
     }
 
+    // Algorithm:
+    // 1) Master processor holds maps of all zone addressing (data provided)
+    // and all remote zone addressing (data required)
+    // 2) Each processor will send the locally active data to the master
+    // 3) Master assembles all the data
+    // 4) Master sends to all processors the data they need to receive
+    //
+    // Notes:
+    // A) If the size of zone addressing is zero, data is not sent
+    // B) Communicated data on each processor has the size of live faces
+    // C) Expanded data will be equal to actual data from other processors
+    //    only for the faces marked in remote; for other faces, it will be
+    //    equal to zero
+    // D) On processor zero, complete data is available
+    // HJ, 4/Jun/2011
+
     // HR, 10/Jul/2013
     // This function requires send-receive-addressing, but usage is not
     // symmetric across processors. Hence trigger re-calculate at this point

src/foam/primitives/Lists/labelPairList.H

@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Typedef
+    Foam::labelPairList
+
+Description
+    List of labelPairs.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef labelPairList_H
+#define labelPairList_H
+
+#include "labelPair.H"
+#include "List.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef List<labelPair> labelPairList;
+    typedef List<labelPairList> labelPairListList;
+    typedef List<labelPairListList> labelPairListListList;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/primitives/Tuple2/Tuple2.H

@@ -25,7 +25,7 @@ Class
     Foam::Tuple2
 
 Description
-    A 2-tuple.
+    A 2-tuple: storing two objects of different types.
 
 SeeAlso
     Foam::Pair for storing two objects of identical types.

src/foam/primitives/ints/ulong/ulongIO.C

@@ -42,6 +42,7 @@ Foam::word Foam::name(const unsigned long val)
     return buf.str();
 }
 
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 Foam::Istream& Foam::operator>>(Istream& is, unsigned long& val)

src/mesh/blockMesh/blockMesh/blockMesh.H

@@ -202,7 +202,7 @@ public:
             void clearGeom();
 
             //- Enable/disable verbose information about the progress
-            static void verbose(const bool on=true);
+            static void verbose(const bool on = true);
 
         // Write
 

src/mesh/blockMesh/blockMesh/blockMeshCreate.C

@@ -54,23 +54,23 @@ void Foam::blockMesh::createPoints() const
 
             label v0 = blocks[blockI].vtxLabel(0, 0, 0);
             label vi1 = blocks[blockI].vtxLabel(1, 0, 0);
-            scalar diStart = mag(blockPoints[vi1]-blockPoints[v0]);
+            scalar diStart = mag(blockPoints[vi1] - blockPoints[v0]);
 
-            label vinM1 = blocks[blockI].vtxLabel(density.x()-1, 0, 0);
+            label vinM1 = blocks[blockI].vtxLabel(density.x() - 1, 0, 0);
             label vin = blocks[blockI].vtxLabel(density.x(), 0, 0);
-            scalar diFinal = mag(blockPoints[vin]-blockPoints[vinM1]);
+            scalar diFinal = mag(blockPoints[vin] - blockPoints[vinM1]);
 
             label vj1 = blocks[blockI].vtxLabel(0, 1, 0);
-            scalar djStart = mag(blockPoints[vj1]-blockPoints[v0]);
-            label vjnM1 = blocks[blockI].vtxLabel(0, density.y()-1, 0);
+            scalar djStart = mag(blockPoints[vj1] - blockPoints[v0]);
+            label vjnM1 = blocks[blockI].vtxLabel(0, density.y() - 1, 0);
             label vjn = blocks[blockI].vtxLabel(0, density.y(), 0);
-            scalar djFinal = mag(blockPoints[vjn]-blockPoints[vjnM1]);
+            scalar djFinal = mag(blockPoints[vjn] - blockPoints[vjnM1]);
 
             label vk1 = blocks[blockI].vtxLabel(0, 0, 1);
-            scalar dkStart = mag(blockPoints[vk1]-blockPoints[v0]);
-            label vknM1 = blocks[blockI].vtxLabel(0, 0, density.z()-1);
+            scalar dkStart = mag(blockPoints[vk1] - blockPoints[v0]);
+            label vknM1 = blocks[blockI].vtxLabel(0, 0, density.z() - 1);
             label vkn = blocks[blockI].vtxLabel(0, 0, density.z());
-            scalar dkFinal = mag(blockPoints[vkn]-blockPoints[vknM1]);
+            scalar dkFinal = mag(blockPoints[vkn] - blockPoints[vknM1]);
 
             Info<< "    Block " << blockI << " cell size :" << nl
                 << "        i : " << scaleFactor_*diStart << " .. "
@@ -222,13 +222,13 @@ Foam::faceList Foam::blockMesh::createPatchFaces
                               + blockOffsets_[blockI]
                             ];
 
-                        if (quadFace[nUnique] != quadFace[nUnique-1])
+                        if (quadFace[nUnique] != quadFace[nUnique - 1])
                         {
                             nUnique++;
                         }
                     }
 
-                    if (quadFace[nUnique-1] == quadFace[0])
+                    if (quadFace[nUnique - 1] == quadFace[0])
                     {
                         nUnique--;
                     }
@@ -290,4 +290,5 @@ void Foam::blockMesh::clearGeom()
     }
 }
 
+
 // ************************************************************************* //

src/mesh/blockMesh/blockMesh/blockMeshMerge.C

@@ -579,7 +579,7 @@ void Foam::blockMesh::calcMergeInfo()
     }
 
     nPoints_ = newPointLabel;
-
 }
 
+
 // ************************************************************************* //

src/postProcessing/functionObjects/forces/forces/forces.H

@@ -28,7 +28,7 @@ Description
     Calculates the forces and moments by integrating the pressure and
     skin-friction forces over a given list of patches.
 
-    Member function calcForcesMoment()calculates and returns the forces and
+    Member function calcForcesMoment() calculates and returns the forces and
     moments.
 
     Member function forces::write() calls calcForcesMoment() and writes the

src/sampling/probes/probes.C

@@ -323,6 +323,7 @@ void Foam::probes::write()
 {
     // Check if the mesh is changing and if so, resample
     const fvMesh& mesh = refCast<const fvMesh>(obr_);
+
     if (mesh.moving() || mesh.changing())
     {
         cellList_.clear();

src/solidModels/fvPatchFields/fixedDisplacement/fixedDisplacementFvPatchVectorField.C

@@ -135,10 +135,11 @@ snGrad() const
     return
     (
         *this
-      - (patchInternalField() + (k&gradField.patchInternalField()))
+      - (patchInternalField() + (k & gradField.patchInternalField()))
       )*this->patch().deltaCoeffs();
 }
 
+
 tmp<Field<vector> > fixedDisplacementFvPatchVectorField::
 gradientBoundaryCoeffs() const
 {
@@ -152,12 +153,13 @@ gradientBoundaryCoeffs() const
     vectorField delta = this->patch().delta();
 
     //- correction vector
-    vectorField k = delta - n*(n&delta);
+    vectorField k = delta - n*(n & delta);
 
     return this->patch().deltaCoeffs()
        *(*this - (k&gradField.patchInternalField()));
 }
 
+
 void fixedDisplacementFvPatchVectorField::write(Ostream& os) const
 {
     fixedValueFvPatchVectorField::write(os);

src/solidModels/fvPatchFields/fixedDisplacement/fixedDisplacementFvPatchVectorField.H

@@ -56,12 +56,12 @@ class fixedDisplacementFvPatchVectorField
 :
     public fixedValueFvPatchVectorField
 {
-
     // Private Data
 
         //- Name of the displacement field
         const word fieldName_;
 
+
 public:
 
     //- Runtime type information
@@ -141,6 +141,7 @@ public:
                 return fieldName_;
             }
 
+
         // Evaluation functions
 
             //- Return patch-normal gradient
@@ -152,6 +153,7 @@ public:
             //  evaluation of the gradient of this patchField
             virtual tmp<Field<vector> > gradientBoundaryCoeffs() const;
 
+
         //- Write
         virtual void write(Ostream&) const;
 };

src/tetFiniteElement/tetFemMatrix/tetFemMatrix.C

@@ -176,14 +176,6 @@ void tetFemMatrix<Type>::addConstraint
     }
     else
     {
-        WarningIn
-        (
-            "void tetFemMatrix<Type>::addConstraint(const label vertex, "
-            "const Type& value)"
-        )   << "Adding constraint on an already constrained point."
-            << "  Point: " << vertex
-            << endl;
-
         fixedEqns_[vertex].combine(cp);
     }
 }

src/thermophysicalModels/basic/psiThermo/ePsiThermo/ePsiThermo.C

@@ -98,7 +98,11 @@ void Foam::ePsiThermo<MixtureType>::calculate()
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class MixtureType>
-Foam::ePsiThermo<MixtureType>::ePsiThermo(const fvMesh& mesh, const objectRegistry& obj)
+Foam::ePsiThermo<MixtureType>::ePsiThermo
+(
+    const fvMesh& mesh,
+    const objectRegistry& obj
+)
 :
     basicPsiThermo(mesh, obj),
     MixtureType(*this, mesh, obj),

src/thermophysicalModels/basic/psiThermo/hPsiThermo/hPsiThermo.C

@@ -96,7 +96,11 @@ void Foam::hPsiThermo<MixtureType>::calculate()
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class MixtureType>
-Foam::hPsiThermo<MixtureType>::hPsiThermo(const fvMesh& mesh, const objectRegistry& obj)
+Foam::hPsiThermo<MixtureType>::hPsiThermo
+(
+    const fvMesh& mesh,
+    const objectRegistry& obj
+)
 :
     basicPsiThermo(mesh, obj),
     MixtureType(*this, mesh, obj),

src/thermophysicalModels/reactionThermo/mixtures/veryInhomogeneousMixture/veryInhomogeneousMixture.H

@@ -86,7 +86,12 @@ public:
     // Constructors
 
         //- Construct from dictionary and mesh
-        veryInhomogeneousMixture(const dictionary&, const fvMesh&, const objectRegistry&);
+        veryInhomogeneousMixture
+        (
+            const dictionary&,
+            const fvMesh&,
+            const objectRegistry&
+        );
 
 
     //- Destructor

src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H

@@ -69,7 +69,9 @@ inline Foam::scalar Foam::specieThermo<thermo>::T
                 "scalar (specieThermo<thermo>::*F)(const scalar) const, "
                 "scalar (specieThermo<thermo>::*dFdT)(const scalar) const"
                 ") const"
-            )   << "Maximum number of iterations exceeded.  Rescue by HJ"
+            )   << "Maximum number of iterations exceeded.  T0 = "
+                << T0 << " F = " << (this->*F)(Test)
+                << " dFdT = " << (this->*dFdT)(Test)
                 << endl;
 
             // Use value where dFdT is calculated using T0.  HJ, 11/Oct/2010

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C

@@ -119,7 +119,8 @@ turbulentMixingLengthDissipationRateInletFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void turbulentMixingLengthDissipationRateInletFvPatchScalarField::updateCoeffs()
+void
+turbulentMixingLengthDissipationRateInletFvPatchScalarField::updateCoeffs()
 {
     if (updated())
     {

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C

@@ -177,7 +177,7 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
     if (!db().foundObject<volScalarField>(GName_))
     {
         InfoIn("void epsilonWallFunctionFvPatchScalarField::updateCoeffs()")
-            << "Cannot access " << GName_ << " field.  for patch "
+            << "Cannot access " << GName_ << " field for patch "
             << patch().name() << ".  Evaluating as zeroGradient"
             << endl;
 

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C

@@ -182,7 +182,7 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
     if (!db().foundObject<volScalarField>(GName_))
     {
         InfoIn("void omegaWallFunctionFvPatchScalarField::updateCoeffs()")
-            << "Cannot access " << GName_ << " field.  for patch "
+            << "Cannot access " << GName_ << " field for patch "
             << patch().name() << ".  Evaluating as zeroGradient"
             << endl;
 

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C

@@ -172,7 +172,7 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
     if (!db().foundObject<volScalarField>(GName_))
     {
         InfoIn("void epsilonWallFunctionFvPatchScalarField::updateCoeffs()")
-            << "Cannot access " << GName_ << " field.  for patch "
+            << "Cannot access " << GName_ << " field for patch "
             << patch().name() << ".  Evaluating as zeroGradient"
             << endl;
 

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C

@@ -177,7 +177,7 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
     if (!db().foundObject<volScalarField>(GName_))
     {
         InfoIn("void omegaWallFunctionFvPatchScalarField::updateCoeffs()")
-            << "Cannot access " << GName_ << " field.  for patch "
+            << "Cannot access " << GName_ << " field for patch "
             << patch().name() << ".  Evaluating as zeroGradient"
             << endl;
 

tutorials/compressible/dbnsTurbFoam/naca0012/Allclean

@@ -5,3 +5,4 @@
 
 cleanCase
 \rm -rf constant/polyMesh/*
+\rm -f naca0012.cas

tutorials/compressible/dbnsTurbFoam/naca0012/Allrun

@@ -7,4 +7,5 @@ application="dbnsTurbFoam"
 
 gunzip naca0012.cas.gz
 runApplication fluentMeshToFoam naca0012.cas
+gzip naca0012.cas
 runApplication $application

tutorials/compressible/dbnsTurbFoam/naca0012/system/fvSchemes

@@ -32,11 +32,14 @@ divSchemes
     div(phi,epsilon) Gauss upwind;
     div(devRhoReff) Gauss linear;
     div((devRhoReff&U)) Gauss linear;
+
+    div((muEff*dev2(grad(U).T()))) Gauss linear;
 }
 
 laplacianSchemes
 {
     default         none;
+    laplacian(muEff,U) Gauss linear limited 0.5;
     laplacian(DkEff,k) Gauss linear limited 0.5;
     laplacian(DepsilonEff,epsilon) Gauss linear limited 0.5;
     laplacian(alphaEff,e) Gauss linear limited 0.5;

tutorials/incompressible/pimpleFoam/t-junction/constant/polyMesh/boundary

@@ -37,7 +37,7 @@ FoamFile
     }
     defaultFaces
     {
-        type            wall;
+        type            empty;
         nFaces          3075;
         startFace       10125;
     }