Merge commit 'OpenCFD/master' into olesenm

Conflicts:

	src/OpenFOAM/db/IOstreams/Pstreams/IPstream.C
	src/OpenFOAM/db/IOstreams/Pstreams/OPstream.C

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files

@@ -5,7 +5,7 @@ coupleManager/coupleManager.C
 derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
 derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C
 derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C
-
+derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
 
 chtMultiRegionFoam.C
 

applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.C

@@ -26,6 +26,7 @@ License
 
 #include "coupleManager.H"
 #include "OFstream.H"
+#include "regionProperties.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -74,6 +75,51 @@ Foam::coupleManager::~coupleManager()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::coupleManager::regionOwner() const
+{
+    const fvMesh& nbrRegion = neighbourRegion();
+
+    const regionProperties& props =
+        localRegion_.objectRegistry::parent().lookupObject<regionProperties>
+        (
+            "regionProperties"
+        );
+
+    label myIndex = findIndex(props.fluidRegionNames(), localRegion_.name());
+    if (myIndex == -1)
+    {
+        label i = findIndex(props.solidRegionNames(), localRegion_.name());
+
+        if (i == -1)
+        {
+            FatalErrorIn("coupleManager::regionOwner() const")
+                << "Cannot find region " << localRegion_.name()
+                << " neither in fluids " << props.fluidRegionNames()
+                << " nor in solids " << props.solidRegionNames()
+                << exit(FatalError);
+        }
+        myIndex = props.fluidRegionNames().size() + i;
+    }
+    label nbrIndex = findIndex(props.fluidRegionNames(), nbrRegion.name());
+    if (nbrIndex == -1)
+    {
+        label i = findIndex(props.solidRegionNames(), nbrRegion.name());
+
+        if (i == -1)
+        {
+            FatalErrorIn("coupleManager::regionOwner() const")
+                << "Cannot find region " << nbrRegion.name()
+                << " neither in fluids " << props.fluidRegionNames()
+                << " nor in solids " << props.solidRegionNames()
+                << exit(FatalError);
+        }
+        nbrIndex = props.fluidRegionNames().size() + i;
+    }
+
+    return myIndex < nbrIndex;
+}
+
+
 void Foam::coupleManager::checkCouple() const
 {
     Info<< "neighbourRegionName_ = " << neighbourRegionName_ << endl;

applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.H

@@ -74,9 +74,6 @@ class coupleManager
 
     // Private Member Functions
 
-        //- Disallow default bitwise copy construct
-//        coupleManager(const coupleManager&);
-
         //- Disallow default bitwise assignment
         void operator=(const coupleManager&);
 
@@ -129,6 +126,9 @@ public:
             template<class Type>
             inline const fvPatchField<Type>& neighbourPatchField() const;
 
+            //- Am I owner (= first to evaluate) of this region interface? 
+            bool regionOwner() const;
+
             //- Check that the couple is valid
             void checkCouple() const;
 

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C

@@ -115,7 +115,7 @@ void Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::updateCoeffs()
     const fvPatchField<scalar>& K =
         patch().lookupPatchField<volScalarField, scalar>(KName_);
 
-    gradient() = refCast<const solidWallTemperatureCoupledFvPatchScalarField>
+    gradient() = -refCast<const solidWallTemperatureCoupledFvPatchScalarField>
         (neighbourField).flux()/K;
 
     fixedGradientFvPatchScalarField::updateCoeffs();

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C

@@ -0,0 +1,252 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "solidWallMixedTemperatureCoupledFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "regionProperties.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
+solidWallMixedTemperatureCoupledFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    coupleManager_(p),
+    KName_("undefined-K")
+{
+    this->refValue() = 0.0;
+    this->refGrad() = 0.0;
+    this->valueFraction() = 1.0;
+    this->fixesValue_ = true;
+}
+
+
+Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
+solidWallMixedTemperatureCoupledFvPatchScalarField
+(
+    const solidWallMixedTemperatureCoupledFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchScalarField(ptf, p, iF, mapper),
+    coupleManager_(ptf.coupleManager_),
+    KName_(ptf.KName_),
+    fixesValue_(ptf.fixesValue_)
+{}
+
+
+Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
+solidWallMixedTemperatureCoupledFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    coupleManager_(p, dict),
+    KName_(dict.lookup("K"))
+{
+    fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
+    refValue() = static_cast<scalarField>(*this);
+    refGrad() = 0.0;
+    valueFraction() = 1.0;
+    fixesValue_ = true;
+}
+
+
+Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
+solidWallMixedTemperatureCoupledFvPatchScalarField
+(
+    const solidWallMixedTemperatureCoupledFvPatchScalarField& wtcsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(wtcsf, iF),
+    coupleManager_(wtcsf.coupleManager_),
+    KName_(wtcsf.KName_),
+    fixesValue_(wtcsf.fixesValue_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::fvPatchScalarField&
+Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::K() const
+{
+    return this->patch().lookupPatchField<volScalarField, scalar>(KName_);
+}
+
+
+void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    tmp<scalarField> intFld = patchInternalField();
+
+    label nFixed = 0;
+
+    // Like snGrad but bypass switching on refValue/refGrad.
+    tmp<scalarField> normalGradient =
+        (*this-intFld())
+      * patch().deltaCoeffs();
+
+    if (debug)
+    {
+        Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
+            << "updateCoeffs() :"
+            << " walltemperature "
+            << " min:" << gMin(*this)
+            << " max:" << gMax(*this)
+            << " avg:" << gAverage(*this)
+            << endl;
+    }
+
+    forAll(*this, i)
+    {
+        // if outgoing flux use fixed value.
+        if (intFld()[i] > operator[](i))
+        {
+            this->refValue()[i] = operator[](i);
+            this->refGrad()[i] = 0.0;   // not used
+            this->valueFraction()[i] = 1.0;
+            nFixed++;
+        }
+        else
+        {
+            this->refValue()[i] = 0.0;  // not used
+            this->refGrad()[i] = normalGradient()[i];
+            this->valueFraction()[i] = 0.0;
+        }
+    }
+
+    reduce(nFixed, sumOp<label>());
+
+    fixesValue_ = (nFixed > 0);
+
+    if (debug)
+    {
+        label nTotSize = returnReduce(this->size(), sumOp<label>());
+
+        Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
+            << "updateCoeffs() : Out of " << nTotSize
+            << " fixedBC:" << nFixed
+            << " gradient:" << nTotSize-nFixed << endl;
+    }
+
+    mixedFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::evaluate
+(
+    const Pstream::commsTypes
+)
+{
+    if (!this->updated())
+    {
+        this->updateCoeffs();
+    }
+
+    if (!coupleManager_.regionOwner())
+    {
+        // I am the last one to evaluate.
+
+        tmp<scalarField> intFld = patchInternalField();
+
+        const fvPatch& nbrPatch = coupleManager_.neighbourPatch();
+
+        solidWallMixedTemperatureCoupledFvPatchScalarField& nbrField =
+        refCast<solidWallMixedTemperatureCoupledFvPatchScalarField>
+        (
+            const_cast<fvPatchField<scalar>&>
+            (
+                coupleManager_.neighbourPatchField<scalar>()
+            )
+        );
+        tmp<scalarField> nbrIntFld = nbrField.patchInternalField();
+        tmp<scalarField> myKDelta = K()*patch().deltaCoeffs();
+        tmp<scalarField> nbrKDelta = nbrField.K()*nbrPatch.deltaCoeffs();
+
+        // Calculate common wall temperature and assign to both sides
+        scalarField::operator=
+        (
+            (myKDelta()*intFld + nbrKDelta()*nbrIntFld)
+          / (myKDelta() + nbrKDelta())
+        );
+
+        nbrField.scalarField::operator=(*this);
+
+        if (debug)
+        {
+            Info<< "Setting master and slave to wall temperature "
+                << " min:" << gMin(*this)
+                << " max:" << gMax(*this)
+                << " avg:" << gAverage(*this)
+                << endl;
+        }
+    }
+
+    fvPatchScalarField::evaluate();
+}
+
+
+void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    mixedFvPatchScalarField::write(os);
+    coupleManager_.writeEntries(os);
+    os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    solidWallMixedTemperatureCoupledFvPatchScalarField
+);
+
+} // End namespace Foam
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H

@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    solidWallMixedTemperatureCoupledFvPatchScalarField
+
+Description
+    Mixed boundary condition for temperature, to be used by the
+    conjugate heat transfer solver.
+    If my temperature is T1, neighbour is T2:
+
+    T1 > T2: my side becomes fixedValue T2 bc, other side becomes fixedGradient.
+    
+
+    Example usage:
+        myInterfacePatchName
+        {
+            type                solidWallMixedTemperatureCoupled;
+            neighbourRegionName fluid;
+            neighbourPatchName  fluidSolidInterface;
+            neighbourFieldName  T;
+            K                   K;
+            value               uniform 300;
+        }
+
+SourceFiles
+    solidWallMixedTemperatureCoupledFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef solidWallMixedTemperatureCoupledFvPatchScalarField_H
+#define solidWallMixedTemperatureCoupledFvPatchScalarField_H
+
+#include "fvPatchFields.H"
+#include "mixedFvPatchFields.H"
+#include "coupleManager.H"
+#include "fvPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+        Class solidWallMixedTemperatureCoupledFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class solidWallMixedTemperatureCoupledFvPatchScalarField
+:
+    public mixedFvPatchScalarField
+{
+    // Private data
+
+        //- Couple manager object
+        coupleManager coupleManager_;
+
+        //- Name of thermal conductivity field
+        word KName_;
+
+        bool fixesValue_;
+
+public:
+
+    //- Runtime type information
+    TypeName("solidWallMixedTemperatureCoupled");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        solidWallMixedTemperatureCoupledFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        solidWallMixedTemperatureCoupledFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given solidWallMixedTemperatureCoupledFvPatchScalarField
+        //  onto a new patch
+        solidWallMixedTemperatureCoupledFvPatchScalarField
+        (
+            const solidWallMixedTemperatureCoupledFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new solidWallMixedTemperatureCoupledFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        solidWallMixedTemperatureCoupledFvPatchScalarField
+        (
+            const solidWallMixedTemperatureCoupledFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new solidWallMixedTemperatureCoupledFvPatchScalarField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+
+    // Member functions
+
+        //- Get corresponding K field
+        const fvPatchScalarField& K() const;
+
+        //- Return true if this patch field fixes a value.
+        //  Needed to check if a level has to be specified while solving
+        //  Poissons equations.
+        virtual bool fixesValue() const
+        {
+            return fixesValue_;
+        }
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Evaluate the patch field
+        virtual void evaluate
+        (
+            const Pstream::commsTypes commsType=Pstream::blocking
+        );
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C

@@ -136,7 +136,7 @@ Foam::solidWallTemperatureCoupledFvPatchScalarField::flux() const
 
     const fvPatchScalarField& Tw = *this;
 
-    return Tw.snGrad()*patch().magSf()*Kw;
+    return Tw.snGrad()*Kw;
 }
 
 

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.H

@@ -138,7 +138,7 @@ public:
 
     // Member functions
 
-        //- Flux
+        //- (intensive) flux
         tmp<scalarField> flux() const;
 
         //- Update the coefficients associated with the patch field

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

@@ -82,7 +82,7 @@
                     IOobject::NO_READ,
                     IOobject::NO_WRITE
                 ),
-                thermof[i].rho()*thermof[i].Cp()*thermof[i].alpha()
+                thermof[i].Cp()*thermof[i].alpha()
             )
         );
 

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

@@ -56,5 +56,5 @@
     }
 
     // Update thermal conductivity
-    Kf[i] = rhof[i]*thermof[i].Cp()*turb[i].alphaEff();
+    Kf[i] = thermof[i].Cp()*turb[i].alphaEff();
 }

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

@@ -6,4 +6,7 @@
             fvm::ddt(rhosCps[i], Ts[i]) - fvm::laplacian(Ks[i], Ts[i])
         );
     }
+
+    Info<< "Min/max T:" << min(Ts[i]) << ' ' << max(Ts[i])
+        << endl;
 }

applications/solvers/incompressible/pisoFoam/pisoFoam.C

@@ -70,6 +70,8 @@ int main(int argc, char *argv[])
               + turbulence->divDevReff(U)
             );
 
+            UEqn.relax();
+
             if (momentumPredictor)
             {
                 solve(UEqn == -fvc::grad(p));

applications/solvers/multiphase/cavitatingFoam/pEqn.H

@@ -37,7 +37,14 @@
           - fvm::laplacian(rUAf, p)
         );
 
-        pEqn.solve();
+        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+        {
+            pEqn.solve(mesh.solver(p.name() + "Final"));
+        }
+        else
+        {
+            pEqn.solve(mesh.solver(p.name()));
+        }
 
         if (nonOrth == nNonOrthCorr)
         {

applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H

@@ -1,13 +1,16 @@
-    surfaceScalarField muf =
+    surfaceScalarField muEff
+    (
+        "muEff",
         twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut());
+      + fvc::interpolate(rho*turbulence->nut())
+    );
 
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muf, U)
+      - fvm::laplacian(muEff, U)
-      - (fvc::grad(U) & fvc::grad(muf))
+      - (fvc::grad(U) & fvc::grad(muEff))
     //- fvc::div(muf*(mesh.Sf() & fvc::interpolate(fvc::grad(U)().T())))
     );
 

applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqnsSubCycle.H

@@ -12,7 +12,9 @@
     surfaceScalarField phic = mag(phi/mesh.magSf());
     phic = min(interface.cAlpha()*phic, max(phic));
 
+    fvc::makeAbsolute(phi, U);
     volScalarField divU = fvc::div(phi);
+    fvc::makeRelative(phi, U);
 
     if (nAlphaSubCycles > 1)
     {

applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H

@@ -38,6 +38,26 @@
           - fvm::laplacian(rUAf, pd)
         );
 
+        if
+        (
+            oCorr == nOuterCorr-1
+            && corr == nCorr-1
+            && nonOrth == nNonOrthCorr
+        )
+        {
+            solve
+            (
+                (
+                    max(alpha1, scalar(0))*(psi1/rho1)
+                  + max(alpha2, scalar(0))*(psi2/rho2)
+                )
+               *pdEqnComp()
+              + pdEqnIncomp,
+                mesh.solver(pd.name() + "Final")
+            );
+        }
+        else
+        {
             solve
             (
                 (
@@ -47,6 +67,7 @@
                *pdEqnComp()
               + pdEqnIncomp
             );
+        }
 
         if (nonOrth == nNonOrthCorr)
         {
@@ -62,7 +83,8 @@
 
     p = max
         (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
+            (pd + gh*(alpha1*rho10 + alpha2*rho20))
+           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
             pMin
         );
 

applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H

@@ -1,15 +1,18 @@
-    surfaceScalarField muf =
+    surfaceScalarField muEff
+    (
+        "muEff",
         twoPhaseProperties->muf()
-      + fvc::interpolate(rho*turbulence->nut());
+      + fvc::interpolate(rho*turbulence->nut())
+    );
 
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
       - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
-      - fvm::laplacian(muf, U)
+      - fvm::laplacian(muEff, U)
-      - (fvc::grad(U) & fvc::grad(muf))
+      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muf*(fvc::interpolate(dev2(fvc::grad(U))) & mesh.Sf()))
+    //- fvc::div(muEff*(fvc::interpolate(dev2(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();

applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C

@@ -167,7 +167,7 @@ Foam::kineticTheoryModel::kineticTheoryModel
             IOobject::NO_WRITE
         ),
         Ua_.mesh(),
-        dimensionedScalar("zero", dimensionSet(0, 0, 0, 0, 0), 0.0)
+        dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
     ),
     gs0_
     (

applications/utilities/mesh/manipulation/setSet/writeFuns.C

@@ -28,7 +28,7 @@ Description
 
 #include "writeFuns.H"
 
-#ifdef __mips
+#if defined(__mips) && !defined(__SICORTEX__)
 #include <standards.h>
 #include <sys/endian.h>
 #endif

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

@@ -27,7 +27,7 @@ License
 #include "writeFuns.H"
 #include "vtkTopo.H"
 
-#ifdef __mips
+#if defined(__mips) && !defined(__SICORTEX__)
 #include <standards.h>
 #include <sys/endian.h>
 #endif

applications/utilities/postProcessing/sampling/probeLocations/probeLocations.C

@@ -44,7 +44,14 @@ int main(int argc, char *argv[])
     instantList timeDirs = timeSelector::select0(runTime, args);
 #   include "createMesh.H"
 
-    IOprobes sniff(mesh, "probesDict", IOobject::MUST_READ, true);
+    IOprobes sniff
+    (
+        probes::typeName,
+        mesh,
+        "probesDict",
+        IOobject::MUST_READ,
+        true
+    );
 
     forAll(timeDirs, timeI)
     {

applications/utilities/postProcessing/sampling/sample/sample.C

@@ -101,8 +101,23 @@ int main(int argc, char *argv[])
     instantList timeDirs = timeSelector::select0(runTime, args);
 #   include "createMesh.H"
 
-    IOsampledSets sSets(mesh, "sampleDict", IOobject::MUST_READ, true);
+    IOsampledSets sSets
-    IOsampledSurfaces sSurfs(mesh, "sampleDict", IOobject::MUST_READ, true);
+    (
+        sampledSets::typeName,
+        mesh,
+        "sampleDict",
+        IOobject::MUST_READ,
+        true
+    );
+
+    IOsampledSurfaces sSurfs
+    (
+        sampledSurfaces::typeName,
+        mesh,
+        "sampleDict",
+        IOobject::MUST_READ,
+        true
+    );
 
     forAll(timeDirs, timeI)
     {

bin/tools/CleanFunctions

@@ -88,7 +88,7 @@ removeCase ()
 
 cleanSamples ()
 {
-    rm -rf {samples,sampleSurfaces} > /dev/null 2>&1
+    rm -rf {sets,samples,sampleSurfaces} > /dev/null 2>&1
 }
 
 cleanUcomponents ()

etc/bashrc

@@ -142,6 +142,16 @@ Linux)
         WM_ARCH=linuxIA64
         export WM_COMPILER=I64
         ;;
+    mips64)
+        WM_ARCH=SiCortex64
+        export WM_COMPILER_LIB_ARCH=64
+        export WM_CC='gcc'
+        export WM_CXX='g++'
+        export WM_CFLAGS='-mabi=64 -fPIC'
+        export WM_CXXFLAGS='-mabi=64 -fPIC'
+        export WM_LDFLAGS='-mabi=64 -G0'
+        export WM_MPLIB=MPI
+        ;;
     *)
         echo Unknown processor type `uname -m` for Linux
         ;;

etc/cshrc

@@ -139,6 +139,16 @@ case Linux:
         setenv WM_ARCH linuxIA64
         setenv WM_COMPILER I64
         breaksw
+    mips64)
+        setenv WM_ARCH SiCortex64
+        setenv WM_COMPILER_LIB_ARCH 64
+        setenv WM_CC 'gcc'
+        setenv WM_CXX 'g++'
+        setenv WM_CFLAGS '-mabi=64 -fPIC'
+        setenv WM_CXXFLAGS '-mabi=64 -fPIC'
+        setenv WM_LDFLAGS '-mabi=64 -G0'
+        setenv WM_MPLIB MPI
+        ;;
     default:
         echo Unknown processor type `uname -m` for Linux
         breaksw

src/OSspecific/Unix/printStack.C

@@ -149,7 +149,7 @@ void getSymbolForRaw
     const word& address
 )
 {
-    if (filename[0] == '/')
+    if (filename.size() > 0 && filename[0] == '/')
     {
         string fcnt = pOpen
         (
@@ -220,6 +220,7 @@ void error::printStack(Ostream& os)
             if (lPos != string::npos && rPos != string::npos && lPos<rPos)
             {
                 address = msg.substr(lPos+1, rPos-lPos-1);
+                msg = msg.substr(0, lPos);
             }
 
             string::size_type bracketPos = msg.find('(');

src/OpenFOAM/db/IOobject/IOobjectI.H

@@ -65,8 +65,8 @@ inline void Foam::IOobject::writeBanner(Stream& os, bool noHint)
         "| =========                 |                                                 |\n"
         "| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\n"
         "|  \\\\    /   O peration     | Version:  " << FOAMversion << spaces << "|\n"
-        "|   \\\\  /    A nd           |                                                 |\n"
+        "|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      |\n"
-        "|    \\\\/     M anipulation  |                                www.OpenFOAM.org |\n"
+        "|    \\\\/     M anipulation  |                                                 |\n"
         "\\*---------------------------------------------------------------------------*/\n";
 }
 

src/OpenFOAM/db/IOstreams/Pstreams/IPstream.C

@@ -45,15 +45,27 @@ inline void Foam::IPstream::checkEof()
 template<class T>
 inline void Foam::IPstream::readFromBuffer(T& t)
 {
+    const size_t align = sizeof(T);
+    bufPosition_ = align + ((bufPosition_ - 1) & ~(align - 1));
+
     t = reinterpret_cast<T&>(buf_[bufPosition_]);
     bufPosition_ += sizeof(T);
     checkEof();
-    // readFromBuffer(&t, sizeof(T));
 }
 
 
-inline void Foam::IPstream::readFromBuffer(void* data, size_t count)
+inline void Foam::IPstream::readFromBuffer
+(
+    void* data,
+    size_t count,
+    size_t align
+)
 {
+    if (align > 1)
+    {
+        bufPosition_ = align + ((bufPosition_ - 1) & ~(align - 1));
+    }
+
     register const char* bufPtr = &buf_[bufPosition_];
     register char* dataPtr = reinterpret_cast<char*>(data);
     register size_t i = count;
@@ -279,7 +291,7 @@ Foam::Istream& Foam::IPstream::read(char* data, std::streamsize count)
             << Foam::abort(FatalError);
     }
 
-    readFromBuffer(data, count);
+    readFromBuffer(data, count, 8);
     return *this;
 }
 

src/OpenFOAM/db/IOstreams/Pstreams/IPstream.H

@@ -70,7 +70,7 @@ class IPstream
         inline void readFromBuffer(T&);
 
         //- Read data from the transfer buffer
-        inline void readFromBuffer(void* data, size_t count);
+        inline void readFromBuffer(void* data, size_t count, size_t align);
 
 
 public:

src/OpenFOAM/db/IOstreams/Pstreams/OPstream.C

@@ -40,10 +40,7 @@ Description
 template<class T>
 inline void Foam::OPstream::writeToBuffer(const T& t)
 {
-    // (T&)(buf_[bufPosition_]) = t;
+    writeToBuffer(&t, sizeof(T), sizeof(T));
-    // bufPosition_ += sizeof(T);
-
-    writeToBuffer(&t, sizeof(T));
 }
 
 
@@ -59,11 +56,24 @@ inline void Foam::OPstream::writeToBuffer(const char& c)
 }
 
 
-inline void Foam::OPstream::writeToBuffer(const void* data, size_t count)
+inline void Foam::OPstream::writeToBuffer
+(
+    const void* data,
+    size_t count,
+    size_t align
+)
 {
+    label oldPos = bufPosition_;
+
+    if (align > 1)
+    {
+        // Align bufPosition. Pads bufPosition_ - oldPos characters.
+        bufPosition_ = align + ((bufPosition_ - 1) & ~(align - 1));
+    }
+
     if (size_t(buf_.size()) < bufPosition_ + count)
     {
-        enlargeBuffer(count);
+        enlargeBuffer(bufPosition_ - oldPos + count);
     }
 
     register char* bufPtr = &buf_[bufPosition_];
@@ -75,6 +85,7 @@ inline void Foam::OPstream::writeToBuffer(const void* data, size_t count)
 }
 
 
+
 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
 
 Foam::OPstream::OPstream
@@ -146,7 +157,7 @@ Foam::Ostream& Foam::OPstream::write(const word& str)
 
     size_t len = str.size();
     writeToBuffer(len);
-    writeToBuffer(str.c_str(), len + 1);
+    writeToBuffer(str.c_str(), len + 1, 1);
 
     return *this;
 }
@@ -158,7 +169,7 @@ Foam::Ostream& Foam::OPstream::write(const string& str)
 
     size_t len = str.size();
     writeToBuffer(len);
-    writeToBuffer(str.c_str(), len + 1);
+    writeToBuffer(str.c_str(), len + 1, 1);
 
     return *this;
 }
@@ -197,7 +208,8 @@ Foam::Ostream& Foam::OPstream::write(const char* data, std::streamsize count)
             << Foam::abort(FatalError);
     }
 
-    writeToBuffer(data, count);
+    writeToBuffer(data, count, 8);
+
     return *this;
 }
 

src/OpenFOAM/db/IOstreams/Pstreams/OPstream.H

@@ -69,7 +69,7 @@ class OPstream
         inline void writeToBuffer(const char&);
 
         //- Write data to the transfer buffer
-        inline void writeToBuffer(const void* data, size_t count);
+        inline void writeToBuffer(const void* data, size_t count, size_t align);
 
 
 public:

src/OpenFOAM/db/Time/TimeIO.C

@@ -67,7 +67,12 @@ void Foam::Time::readDict()
             case wcAdjustableRunTime:
                 // Recalculate outputTimeIndex_ to be in units of current
                 // writeInterval.
-                outputTimeIndex_ *= oldWriteInterval/writeInterval_;
+                outputTimeIndex_ = label
+                (
+                    outputTimeIndex_
+                  * oldWriteInterval
+                  / writeInterval_
+                );
             break;
 
             default:

src/OpenFOAM/global/global.Cver

@@ -37,6 +37,11 @@ Description
 
 const char* const Foam::FOAMversion = "WM_PROJECT_VERSION";
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Static initializers for string::null, word::null and fileName::null
+
+#include "stringsGlobals.C"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Setup an error handler for the global new operator
 

src/OpenFOAM/primitives/strings/fileName/fileName.C

@@ -33,7 +33,6 @@ License
 
 const char* const Foam::fileName::typeName = "fileName";
 int Foam::fileName::debug(debug::debugSwitch(fileName::typeName, 0));
-const Foam::fileName Foam::fileName::null;
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 

src/OpenFOAM/primitives/strings/string/string.C

@@ -32,7 +32,6 @@ License
 
 const char* const Foam::string::typeName = "string";
 int Foam::string::debug(debug::debugSwitch(string::typeName, 0));
-const Foam::string Foam::string::null;
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 

src/OpenFOAM/primitives/strings/stringsGlobals.C

@@ -0,0 +1,45 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Static initializers for
+        Foam::string::null
+        Foam::word::null
+        Foam::fileName::null.
+    This file is included in global.Cver since these members are required by
+    debug.C.
+
+\*---------------------------------------------------------------------------*/
+
+#include "string.H"
+#include "word.H"
+#include "fileName.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::string Foam::string::null;
+const Foam::word Foam::word::null;
+const Foam::fileName Foam::fileName::null;
+
+// ************************************************************************* //

src/OpenFOAM/primitives/strings/word/word.C

@@ -31,6 +31,5 @@ License
 
 const char* const Foam::word::typeName = "word";
 int Foam::word::debug(Foam::debug::debugSwitch(word::typeName, 0));
-const Foam::word Foam::word::null;
 
 // ************************************************************************* //

src/OpenFOAM/primitives/uint/uintIO.C

@@ -44,7 +44,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Return a string representation of an uint
 word name(const unsigned int i)
 {
     std::ostringstream osBuffer;
@@ -66,7 +65,7 @@ Istream& operator>>(Istream& is, unsigned int& i)
 
     if (t.isLabel())
     {
-        i = uint(t.labelToken());
+        i = unsigned(t.labelToken());
     }
     else
     {

src/autoMesh/Make/files

@@ -15,6 +15,7 @@ $(autoHexMeshDriver)/pointData/pointData.C
 $(autoHexMesh)/meshRefinement/meshRefinementBaffles.C
 $(autoHexMesh)/meshRefinement/meshRefinement.C
 $(autoHexMesh)/meshRefinement/meshRefinementMerge.C
+$(autoHexMesh)/meshRefinement/meshRefinementProblemCells.C
 $(autoHexMesh)/meshRefinement/meshRefinementRefine.C
 $(autoHexMesh)/refinementSurfaces/refinementSurfaces.C
 $(autoHexMesh)/shellSurfaces/shellSurfaces.C

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C

@@ -518,7 +518,9 @@ void Foam::autoRefineDriver::baffleAndSplitMesh
     // be like boundary face from now on so not coupled anymore.
     meshRefiner_.baffleAndSplitMesh
     (
-        handleSnapProblems,
+        handleSnapProblems,             // detect&remove potential snap problem
+        false,                          // perpendicular edge connected cells
+        scalarField(0),                 // per region perpendicular angle
         !handleSnapProblems,            // merge free standing baffles?
         motionDict,
         const_cast<Time&>(mesh.time()),
@@ -592,9 +594,13 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
         const_cast<Time&>(mesh.time())++;
     }
 
+    const scalarField& perpAngle = meshRefiner_.surfaces().perpendicularAngle();
+
     meshRefiner_.baffleAndSplitMesh
     (
         handleSnapProblems,
+        handleSnapProblems,                 // remove perp edge connected cells
+        perpAngle,                          // perp angle
         false,                              // merge free standing baffles?
         motionDict,
         const_cast<Time&>(mesh.time()),

src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -36,6 +36,7 @@ SourceFiles
     meshRefinement.C
     meshRefinementBaffles.C
     meshRefinementMerge.C
+    meshRefinementProblemCells.C
     meshRefinementRefine.C
 
 \*---------------------------------------------------------------------------*/
@@ -51,6 +52,7 @@ SourceFiles
 #include "indirectPrimitivePatch.H"
 #include "pointFieldsFwd.H"
 #include "Tuple2.H"
+#include "pointIndexHit.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -345,6 +347,8 @@ private:
                 polyTopoChange& meshMod
             ) const;
 
+        // Problem cell handling
+
             //- Helper function to mark face as being on 'boundary'. Used by
             //  markFacesOnProblemCells
             void markBoundaryFace
@@ -355,15 +359,32 @@ private:
                 boolList& isBoundaryPoint
             ) const;
 
-            //- Returns list with for every internal face -1 or the patch
+            void findNearest
-            //  they should be baffled into.
-            labelList markFacesOnProblemCells
             (
-                const labelList& globalToPatch
+                const labelList& meshFaces,
+                List<pointIndexHit>& nearestInfo,
+                labelList& nearestSurface,
+                labelList& nearestRegion,
+                vectorField& nearestNormal
+            ) const;
+
+            Map<label> findEdgeConnectedProblemCells
+            (
+                const scalarField& perpendicularAngle,
+                const labelList&
             ) const;
 
             //- Returns list with for every internal face -1 or the patch
-            //  they should be baffled into.
+            //  they should be baffled into. If removeEdgeConnectedCells is set
+            //  removes cells based on perpendicularAngle.
+            labelList markFacesOnProblemCells
+            (
+                const bool removeEdgeConnectedCells,
+                const scalarField& perpendicularAngle,
+                const labelList& globalToPatch
+            ) const;
+
+            //- Initial test of marking faces using geometric information.
             labelList markFacesOnProblemCellsGeometric
             (
                 const dictionary& motionDict,
@@ -589,6 +610,8 @@ public:
             void baffleAndSplitMesh
             (
                 const bool handleSnapProblems,
+                const bool removeEdgeConnectedCells,
+                const scalarField& perpendicularAngle,
                 const bool mergeFreeStanding,
                 const dictionary& motionDict,
                 Time& runTime,

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -44,10 +44,6 @@ License
 #include "OFstream.H"
 #include "regionSplit.H"
 #include "removeCells.H"
-#include "motionSmoother.H"
-#include "polyMeshGeometry.H"
-#include "IOmanip.H"
-#include "cellSet.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -230,7 +226,7 @@ void Foam::meshRefinement::getBafflePatches
     label vertI = 0;
     if (debug&OBJINTERSECTIONS)
     {
-        str.reset(new OFstream(mesh_.time().path()/"intersections.obj"));
+        str.reset(new OFstream(mesh_.time().timePath()/"intersections.obj"));
 
         Pout<< "getBafflePatches : Writing surface intersections to file "
             << str().name() << nl << endl;
@@ -510,642 +506,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
 }
 
 
-void Foam::meshRefinement::markBoundaryFace
-(
-    const label faceI,
-    boolList& isBoundaryFace,
-    boolList& isBoundaryEdge,
-    boolList& isBoundaryPoint
-) const
-{
-    isBoundaryFace[faceI] = true;
-
-    const labelList& fEdges = mesh_.faceEdges(faceI);
-
-    forAll(fEdges, fp)
-    {
-        isBoundaryEdge[fEdges[fp]] = true;
-    }
-
-    const face& f = mesh_.faces()[faceI];
-
-    forAll(f, fp)
-    {
-        isBoundaryPoint[f[fp]] = true;
-    }
-}
-
-
-// Returns list with for every internal face -1 or the patch they should
-// be baffled into. Gets run after createBaffles so all the surface
-// intersections have already been turned into baffles. Used to remove cells
-// by baffling all their faces and have the splitMeshRegions chuck away non
-// used regions.
-Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
-(
-    const labelList& globalToPatch
-) const
-{
-    const labelList& cellLevel = meshCutter_.cellLevel();
-    const labelList& pointLevel = meshCutter_.pointLevel();
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-
-    // Per internal face (boundary faces not used) the patch that the
-    // baffle should get (or -1)
-    labelList facePatch(mesh_.nFaces(), -1);
-
-    // Mark all points and edges on baffle patches (so not on any inlets,
-    // outlets etc.)
-    boolList isBoundaryPoint(mesh_.nPoints(), false);
-    boolList isBoundaryEdge(mesh_.nEdges(), false);
-    boolList isBoundaryFace(mesh_.nFaces(), false);
-
-    // Fill boundary data. All elements on meshed patches get marked.
-    // Get the labels of added patches.
-    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
-
-    forAll(adaptPatchIDs, i)
-    {
-        label patchI = adaptPatchIDs[i];
-
-        const polyPatch& pp = patches[patchI];
-
-        label faceI = pp.start();
-
-        forAll(pp, j)
-        {
-            markBoundaryFace
-            (
-                faceI,
-                isBoundaryFace,
-                isBoundaryEdge,
-                isBoundaryPoint
-            );
-
-            faceI++;
-        }
-    }
-
-    syncTools::syncPointList
-    (
-        mesh_,
-        isBoundaryPoint,
-        orEqOp<bool>(),
-        false,              // null value
-        false               // no separation
-    );
-
-    syncTools::syncEdgeList
-    (
-        mesh_,
-        isBoundaryEdge,
-        orEqOp<bool>(),
-        false,              // null value
-        false               // no separation
-    );
-
-    syncTools::syncFaceList
-    (
-        mesh_,
-        isBoundaryFace,
-        orEqOp<bool>(),
-        false               // no separation
-    );
-
-
-    // For each cell count the number of anchor points that are on
-    // the boundary:
-    // 8 : check the number of (baffle) boundary faces. If 3 or more block
-    //     off the cell since the cell would get squeezed down to a diamond
-    //     (probably; if the 3 or more faces are unrefined (only use the
-    //      anchor points))
-    // 7 : store. Used to check later on whether there are points with
-    //     3 or more of these cells. (note that on a flat surface a boundary
-    //     point will only have 4 cells connected to it)
-
-    // Does cell have exactly 7 of its 8 anchor points on the boundary?
-    PackedList<1> hasSevenBoundaryAnchorPoints(mesh_.nCells(), 0u);
-    // If so what is the remaining non-boundary anchor point?
-    labelHashSet nonBoundaryAnchors(mesh_.nCells()/10000);
-
-    // On-the-fly addressing storage.
-    DynamicList<label> dynFEdges;
-    DynamicList<label> dynCPoints;
-
-    // Count of faces marked for baffling
-    label nBaffleFaces = 0;
-
-    forAll(cellLevel, cellI)
-    {
-        const labelList& cPoints = mesh_.cellPoints(cellI, dynCPoints);
-
-        // Get number of anchor points (pointLevel == cellLevel)
-
-        label nBoundaryAnchors = 0;
-        label nNonAnchorBoundary = 0;
-        label nonBoundaryAnchor = -1;
-
-        forAll(cPoints, i)
-        {
-            label pointI = cPoints[i];
-
-            if (pointLevel[pointI] <= cellLevel[cellI])
-            {
-                // Anchor point
-                if (isBoundaryPoint[pointI])
-                {
-                    nBoundaryAnchors++;
-                }
-                else
-                {
-                    // Anchor point which is not on the surface
-                    nonBoundaryAnchor = pointI;
-                }
-            }
-            else if (isBoundaryPoint[pointI])
-            {
-                nNonAnchorBoundary++;
-            }
-        }
-
-        if (nBoundaryAnchors == 8)
-        {
-            const cell& cFaces = mesh_.cells()[cellI];
-
-            // Count boundary faces.
-            label nBfaces = 0;
-
-            forAll(cFaces, cFaceI)
-            {
-                if (isBoundaryFace[cFaces[cFaceI]])
-                {
-                    nBfaces++;
-                }
-            }
-
-            // If nBfaces > 1 make all non-boundary non-baffle faces baffles.
-            // We assume that this situation is where there is a single
-            // cell sticking out which would get flattened.
-
-            // Eugene: delete cell no matter what.
-            //if (nBfaces > 1)
-            {
-                forAll(cFaces, cf)
-                {
-                    label faceI = cFaces[cf];
-
-                    if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
-                    {
-                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
-                        nBaffleFaces++;
-
-                        // Mark face as a 'boundary'
-                        markBoundaryFace
-                        (
-                            faceI,
-                            isBoundaryFace,
-                            isBoundaryEdge,
-                            isBoundaryPoint
-                        );
-                    }
-                }
-            }
-        }
-        else if (nBoundaryAnchors == 7)
-        {
-            // Mark the cell. Store the (single!) non-boundary anchor point.
-            hasSevenBoundaryAnchorPoints.set(cellI, 1u);
-            nonBoundaryAnchors.insert(nonBoundaryAnchor);
-        }
-    }
-
-
-    // Loop over all points. If a point is connected to 4 or more cells
-    // with 7 anchor points on the boundary set those cell's non-boundary faces
-    // to baffles
-
-    DynamicList<label> dynPCells;
-
-    forAllConstIter(labelHashSet, nonBoundaryAnchors, iter)
-    {
-        label pointI = iter.key();
-
-        const labelList& pCells = mesh_.pointCells(pointI, dynPCells);
-
-        // Count number of 'hasSevenBoundaryAnchorPoints' cells.
-        label n = 0;
-
-        forAll(pCells, i)
-        {
-            if (hasSevenBoundaryAnchorPoints.get(pCells[i]) == 1u)
-            {
-                n++;
-            }
-        }
-
-        if (n > 3)
-        {
-            // Point in danger of being what? Remove all 7-cells.
-            forAll(pCells, i)
-            {
-                label cellI = pCells[i];
-
-                if (hasSevenBoundaryAnchorPoints.get(cellI) == 1u)
-                {
-                    const cell& cFaces = mesh_.cells()[cellI];
-
-                    forAll(cFaces, cf)
-                    {
-                        label faceI = cFaces[cf];
-
-                        if
-                        (
-                            facePatch[faceI] == -1
-                         && mesh_.isInternalFace(faceI)
-                        )
-                        {
-                            facePatch[faceI] = getBafflePatch(facePatch, faceI);
-                            nBaffleFaces++;
-
-                            // Mark face as a 'boundary'
-                            markBoundaryFace
-                            (
-                                faceI,
-                                isBoundaryFace,
-                                isBoundaryEdge,
-                                isBoundaryPoint
-                            );
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-
-    // Sync all. (note that pointdata and facedata not used anymore but sync
-    // anyway)
-
-    syncTools::syncPointList
-    (
-        mesh_,
-        isBoundaryPoint,
-        orEqOp<bool>(),
-        false,              // null value
-        false               // no separation
-    );
-
-    syncTools::syncEdgeList
-    (
-        mesh_,
-        isBoundaryEdge,
-        orEqOp<bool>(),
-        false,              // null value
-        false               // no separation
-    );
-
-    syncTools::syncFaceList
-    (
-        mesh_,
-        isBoundaryFace,
-        orEqOp<bool>(),
-        false               // no separation
-    );
-
-
-    // Find faces with all edges on the boundary and make them baffles
-    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-    {
-        if (facePatch[faceI] == -1)
-        {
-            const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
-            label nFaceBoundaryEdges = 0;
-
-            forAll(fEdges, fe)
-            {
-                if (isBoundaryEdge[fEdges[fe]])
-                {
-                    nFaceBoundaryEdges++;
-                }
-            }
-
-            if (nFaceBoundaryEdges == fEdges.size())
-            {
-                facePatch[faceI] = getBafflePatch(facePatch, faceI);
-                nBaffleFaces++;
-
-                // Do NOT update boundary data since this would grow blocked
-                // faces across gaps.
-            }
-        }
-    }
-
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
-        {
-            label faceI = pp.start();
-
-            forAll(pp, i)
-            {
-                if (facePatch[faceI] == -1)
-                {
-                    const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
-                    label nFaceBoundaryEdges = 0;
-
-                    forAll(fEdges, fe)
-                    {
-                        if (isBoundaryEdge[fEdges[fe]])
-                        {
-                            nFaceBoundaryEdges++;
-                        }
-                    }
-
-                    if (nFaceBoundaryEdges == fEdges.size())
-                    {
-                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
-                        nBaffleFaces++;
-
-                        // Do NOT update boundary data since this would grow
-                        // blocked faces across gaps.
-                    }
-                }
-
-                faceI++;
-            }
-        }
-    }
-
-    Info<< "markFacesOnProblemCells : marked "
-        << returnReduce(nBaffleFaces, sumOp<label>())
-        << " additional internal faces to be converted into baffles."
-        << endl;
-
-    return facePatch;
-}
-//XXXXXXXXXXXXXX
-// Mark faces to be baffled to prevent snapping problems. Does
-// test to find nearest surface and checks which faces would get squashed.
-Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
-(
-    const dictionary& motionDict,
-    const labelList& globalToPatch
-) const
-{
-    // Get the labels of added patches.
-    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
-
-    // Construct addressing engine.
-    autoPtr<indirectPrimitivePatch> ppPtr
-    (
-        meshRefinement::makePatch
-        (
-            mesh_,
-            adaptPatchIDs
-        )
-    );
-    const indirectPrimitivePatch& pp = ppPtr();
-    const pointField& localPoints = pp.localPoints();
-    const labelList& meshPoints = pp.meshPoints();
-
-    // Find nearest (non-baffle) surface
-    pointField newPoints(mesh_.points());
-    {
-        List<pointIndexHit> hitInfo;
-        labelList hitSurface;
-        surfaces_.findNearest
-        (
-            surfaces_.getUnnamedSurfaces(),
-            localPoints,
-            scalarField(localPoints.size(), sqr(GREAT)),    // sqr of attraction
-            hitSurface,
-            hitInfo
-        );
-    
-        forAll(hitInfo, i)
-        {
-            if (hitInfo[i].hit())
-            {
-                //label pointI = meshPoints[i];
-                //Pout<< "   " << pointI << " moved from "
-                //    << mesh_.points()[pointI] << " by "
-                //    << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
-                //    << endl;
-                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
-            }
-        }
-    }
-
-    // Per face (internal or coupled!) the patch that the
-    // baffle should get (or -1).
-    labelList facePatch(mesh_.nFaces(), -1);
-    // Count of baffled faces
-    label nBaffleFaces = 0;
-
-
-//    // Sync position? Or not since same face on both side so just sync
-//    // result of baffle.
-//
-//    const scalar minArea(readScalar(motionDict.lookup("minArea")));
-//
-//    Pout<< "markFacesOnProblemCellsGeometric : Comparing to minArea:"
-//        << minArea << endl;
-//
-//    pointField facePoints;
-//    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
-//    {
-//        const face& f = mesh_.faces()[faceI];
-//
-//        bool usesPatchPoint = false;
-//
-//        facePoints.setSize(f.size());
-//        forAll(f, fp)
-//        {
-//            Map<label>::const_iterator iter = pp.meshPointMap().find(f[fp]);
-//
-//            if (iter != pp.meshPointMap().end())
-//            {
-//                facePoints[fp] = newPosition[iter()];
-//                usesPatchPoint = true;
-//            }
-//            else
-//            {
-//                facePoints[fp] = mesh_.points()[f[fp]];
-//            }
-//        }
-//
-//        if (usesPatchPoint)
-//        {
-//            // Check area of face wrt original area
-//            face identFace(identity(f.size()));
-//
-//            if (identFace.mag(facePoints) < minArea)
-//            {
-//                facePatch[faceI] = getBafflePatch(facePatch, faceI);
-//                nBaffleFaces++;
-//            }
-//        }
-//    }
-//
-//
-//    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-//    forAll(patches, patchI)
-//    {
-//        const polyPatch& pp = patches[patchI];
-//
-//        if (pp.coupled())
-//        {
-//            forAll(pp, i)
-//            {
-//                label faceI = pp.start()+i;
-//
-//                const face& f = mesh_.faces()[faceI];
-//
-//                bool usesPatchPoint = false;
-//
-//                facePoints.setSize(f.size());
-//                forAll(f, fp)
-//                {
-//                    Map<label>::const_iterator iter =
-//                        pp.meshPointMap().find(f[fp]);
-//
-//                    if (iter != pp.meshPointMap().end())
-//                    {
-//                        facePoints[fp] = newPosition[iter()];
-//                        usesPatchPoint = true;
-//                    }
-//                    else
-//                    {
-//                        facePoints[fp] = mesh_.points()[f[fp]];
-//                    }
-//                }
-//
-//                if (usesPatchPoint)
-//                {
-//                    // Check area of face wrt original area
-//                    face identFace(identity(f.size()));
-//
-//                    if (identFace.mag(facePoints) < minArea)
-//                    {
-//                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
-//                        nBaffleFaces++;
-//                    }
-//                }
-//            }
-//        }
-//    }
-
-    {
-        pointField oldPoints(mesh_.points());
-        mesh_.movePoints(newPoints);
-        faceSet wrongFaces(mesh_, "wrongFaces", 100);
-        {
-            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
-
-            // Just check the errors from squashing
-            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-            const labelList allFaces(identity(mesh_.nFaces()));
-            label nWrongFaces = 0;
-
-            scalar minArea(readScalar(motionDict.lookup("minArea")));
-            if (minArea > -SMALL)
-            {
-                polyMeshGeometry::checkFaceArea
-                (
-                    false,
-                    minArea,
-                    mesh_,
-                    mesh_.faceAreas(),
-                    allFaces,
-                    &wrongFaces
-                );
-
-                label nNewWrongFaces = returnReduce
-                (
-                    wrongFaces.size(),
-                    sumOp<label>()
-                );
-
-                Info<< "    faces with area < "
-                    << setw(5) << minArea
-                    << " m^2                            : "
-                    << nNewWrongFaces-nWrongFaces << endl;
-
-                nWrongFaces = nNewWrongFaces;
-            }
-
-//            scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
-            scalar minDet = 0.01;
-            if (minDet > -1)
-            {
-                polyMeshGeometry::checkCellDeterminant
-                (
-                    false,
-                    minDet,
-                    mesh_,
-                    mesh_.faceAreas(),
-                    allFaces,
-                    polyMeshGeometry::affectedCells(mesh_, allFaces),
-                    &wrongFaces
-                );
-
-                label nNewWrongFaces = returnReduce
-                (
-                    wrongFaces.size(),
-                    sumOp<label>()
-                );
-
-                Info<< "    faces on cells with determinant < "
-                    << setw(5) << minDet << "                : "
-                    << nNewWrongFaces-nWrongFaces << endl;
-
-                nWrongFaces = nNewWrongFaces;
-            }
-        }
-
-
-        forAllConstIter(faceSet, wrongFaces, iter)
-        {
-            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
-
-            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
-            {
-                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
-                nBaffleFaces++;
-
-                //Pout<< "    " << iter.key()
-                //    //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
-                //    << " is destined for patch " << facePatch[iter.key()]
-                //    << endl;
-            }
-        }
-        // Restore points.
-        mesh_.movePoints(oldPoints);
-    }
-
-
-    Info<< "markFacesOnProblemCellsGeometric : marked "
-        << returnReduce(nBaffleFaces, sumOp<label>())
-        << " additional internal and coupled faces"
-        << " to be converted into baffles." << endl;
-
-    syncTools::syncFaceList
-    (
-        mesh_,
-        facePatch,
-        maxEqOp<label>(),
-        false               // no separation
-    );
-
-    return facePatch;
-}
-//XXXXXXXX
-
-
 // Return a list of coupled face pairs, i.e. faces that use the same vertices.
 // (this information is recalculated instead of maintained since would be too
 // hard across splitMeshRegions).
@@ -1904,6 +1264,8 @@ void Foam::meshRefinement::findCellZoneTopo
 void Foam::meshRefinement::baffleAndSplitMesh
 (
     const bool handleSnapProblems,
+    const bool removeEdgeConnectedCells,
+    const scalarField& perpendicularAngle,
     const bool mergeFreeStanding,
     const dictionary& motionDict,
     Time& runTime,
@@ -1981,6 +1343,8 @@ void Foam::meshRefinement::baffleAndSplitMesh
         (
             markFacesOnProblemCells
             (
+                removeEdgeConnectedCells,
+                perpendicularAngle,
                 globalToPatch
             )
             //markFacesOnProblemCellsGeometric
@@ -2024,6 +1388,8 @@ void Foam::meshRefinement::baffleAndSplitMesh
             (    
                 markFacesOnProblemCells
                 (
+                    removeEdgeConnectedCells,
+                    perpendicularAngle,
                     globalToPatch
                 )
             );
@@ -2099,7 +1465,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
     {
         Pout<< "Writing subsetted mesh to time " << mesh_.time().timeName()
             << endl;
-        write(debug, runTime.path()/runTime.timeName());
+        write(debug, runTime.timePath());
         Pout<< "Dumped debug data in = "
             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
     }

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C

@@ -0,0 +1,952 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "meshRefinement.H"
+#include "fvMesh.H"
+#include "syncTools.H"
+#include "Time.H"
+#include "refinementSurfaces.H"
+#include "pointSet.H"
+#include "faceSet.H"
+#include "indirectPrimitivePatch.H"
+#include "OFstream.H"
+#include "cellSet.H"
+#include "searchableSurfaces.H"
+#include "polyMeshGeometry.H"
+#include "IOmanip.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::meshRefinement::markBoundaryFace
+(
+    const label faceI,
+    boolList& isBoundaryFace,
+    boolList& isBoundaryEdge,
+    boolList& isBoundaryPoint
+) const
+{
+    isBoundaryFace[faceI] = true;
+
+    const labelList& fEdges = mesh_.faceEdges(faceI);
+
+    forAll(fEdges, fp)
+    {
+        isBoundaryEdge[fEdges[fp]] = true;
+    }
+
+    const face& f = mesh_.faces()[faceI];
+
+    forAll(f, fp)
+    {
+        isBoundaryPoint[f[fp]] = true;
+    }
+}
+
+
+void Foam::meshRefinement::findNearest
+(
+    const labelList& meshFaces,
+    List<pointIndexHit>& nearestInfo,
+    labelList& nearestSurface,
+    labelList& nearestRegion,
+    vectorField& nearestNormal
+) const
+{
+    pointField fc(meshFaces.size());
+    forAll(meshFaces, i)
+    {
+        fc[i] = mesh_.faceCentres()[meshFaces[i]];
+    }
+
+    const labelList allSurfaces(identity(surfaces_.surfaces().size()));
+
+    surfaces_.findNearest
+    (
+        allSurfaces,
+        fc,
+        scalarField(fc.size(), sqr(GREAT)),    // sqr of attraction
+        nearestSurface,
+        nearestInfo
+    );
+
+    // Do normal testing per surface.
+    nearestNormal.setSize(nearestInfo.size());
+    nearestRegion.setSize(nearestInfo.size());
+
+    forAll(allSurfaces, surfI)
+    {
+        DynamicList<pointIndexHit> localHits;
+
+        forAll(nearestSurface, i)
+        {
+            if (nearestSurface[i] == surfI)
+            {
+                localHits.append(nearestInfo[i]);
+            }
+        }
+
+        label geomI = surfaces_.surfaces()[surfI];
+
+        pointField localNormals;
+        surfaces_.geometry()[geomI].getNormal(localHits, localNormals);
+
+        labelList localRegion;
+        surfaces_.geometry()[geomI].getRegion(localHits, localRegion);
+
+        label localI = 0;
+        forAll(nearestSurface, i)
+        {
+            if (nearestSurface[i] == surfI)
+            {
+                nearestNormal[i] = localNormals[localI];
+                nearestRegion[i] = localRegion[localI];
+                localI++;
+            }
+        }
+    }
+}
+
+
+Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
+(
+    const scalarField& perpendicularAngle,
+    const labelList& globalToPatch
+) const
+{
+    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
+
+    // Construct addressing engine.
+    autoPtr<indirectPrimitivePatch> ppPtr
+    (
+        meshRefinement::makePatch
+        (
+            mesh_,
+            adaptPatchIDs
+        )
+    );
+    const indirectPrimitivePatch& pp = ppPtr();
+
+
+    // 1. Collect faces to test
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    DynamicList<label> candidateFaces(pp.size()/20);
+
+    const labelListList& edgeFaces = pp.edgeFaces();
+
+    const labelList& cellLevel = meshCutter_.cellLevel();
+
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        if (eFaces.size() == 2)
+        {
+            label face0 = pp.addressing()[eFaces[0]];
+            label face1 = pp.addressing()[eFaces[1]];
+
+            label cell0 = mesh_.faceOwner()[face0];
+            label cell1 = mesh_.faceOwner()[face1];
+
+            if (cellLevel[cell0] > cellLevel[cell1])
+            {
+                // cell0 smaller.
+                const vector& n0 = pp.faceNormals()[eFaces[0]];
+                const vector& n1 = pp.faceNormals()[eFaces[1]];
+
+                if (mag(n0 & n1) < 0.1)
+                {
+                    candidateFaces.append(face0);
+                }
+            }
+            else if (cellLevel[cell1] > cellLevel[cell0])
+            {
+                // cell1 smaller.
+                const vector& n0 = pp.faceNormals()[eFaces[0]];
+                const vector& n1 = pp.faceNormals()[eFaces[1]];
+
+                if (mag(n0 & n1) < 0.1)
+                {
+                    candidateFaces.append(face1);
+                }
+            }
+        }
+    }
+    candidateFaces.shrink();
+
+    Info<< "Testing " << returnReduce(candidateFaces.size(), sumOp<label>())
+        << " faces on edge-connected cells of differing level."
+        << endl;
+
+    if (debug)
+    {
+        faceSet fSet(mesh_, "edgeConnectedFaces", candidateFaces);
+        Pout<< "Writing " << fSet.size()
+            << " with problematic topology to faceSet "
+            << fSet.objectPath() << endl;
+        fSet.write();
+    }
+
+
+    // 2. Find nearest surface on candidate faces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    List<pointIndexHit> nearestInfo;
+    labelList nearestSurface;
+    labelList nearestRegion;
+    vectorField nearestNormal;
+    findNearest
+    (
+        candidateFaces,
+        nearestInfo,
+        nearestSurface,
+        nearestRegion,
+        nearestNormal
+    );
+
+
+    // 3. Test angle to surface
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    Map<label> candidateCells(candidateFaces.size());
+
+    faceSet perpFaces(mesh_, "perpendicularFaces", pp.size()/100);
+
+    forAll(candidateFaces, i)
+    {
+        label faceI = candidateFaces[i];
+
+        vector n = mesh_.faceAreas()[faceI];
+        n /= mag(n);
+
+        label region = surfaces_.globalRegion
+        (
+            nearestSurface[i],
+            nearestRegion[i]
+        );
+
+        scalar angle =
+            perpendicularAngle[region]
+          / 180.0
+          * mathematicalConstant::pi;
+
+        if (angle >= 0)
+        {
+            if (mag(n & nearestNormal[i]) < Foam::sin(angle))
+            {
+                perpFaces.insert(faceI);
+                candidateCells.insert
+                (
+                    mesh_.faceOwner()[faceI],
+                    globalToPatch[region]
+                );
+            }
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "Writing " << perpFaces.size()
+            << " faces that are perpendicular to the surface to set "
+            << perpFaces.objectPath() << endl;
+        perpFaces.write();
+    }
+    return candidateCells;
+}
+
+
+// Returns list with for every internal face -1 or the patch they should
+// be baffled into. Gets run after createBaffles so all the surface
+// intersections have already been turned into baffles. Used to remove cells
+// by baffling all their faces and have the splitMeshRegions chuck away non
+// used regions.
+Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
+(
+    const bool removeEdgeConnectedCells,
+    const scalarField& perpendicularAngle,
+    const labelList& globalToPatch
+) const
+{
+    const labelList& cellLevel = meshCutter_.cellLevel();
+    const labelList& pointLevel = meshCutter_.pointLevel();
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+
+    // Per internal face (boundary faces not used) the patch that the
+    // baffle should get (or -1)
+    labelList facePatch(mesh_.nFaces(), -1);
+
+    // Mark all points and edges on baffle patches (so not on any inlets,
+    // outlets etc.)
+    boolList isBoundaryPoint(mesh_.nPoints(), false);
+    boolList isBoundaryEdge(mesh_.nEdges(), false);
+    boolList isBoundaryFace(mesh_.nFaces(), false);
+
+    // Fill boundary data. All elements on meshed patches get marked.
+    // Get the labels of added patches.
+    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
+
+    forAll(adaptPatchIDs, i)
+    {
+        label patchI = adaptPatchIDs[i];
+
+        const polyPatch& pp = patches[patchI];
+
+        label faceI = pp.start();
+
+        forAll(pp, j)
+        {
+            markBoundaryFace
+            (
+                faceI,
+                isBoundaryFace,
+                isBoundaryEdge,
+                isBoundaryPoint
+            );
+
+            faceI++;
+        }
+    }
+
+    // Count of faces marked for baffling
+    label nBaffleFaces = 0;
+
+    if (removeEdgeConnectedCells && max(perpendicularAngle) >= 0)
+    {
+        Info<< "markFacesOnProblemCells :"
+            << " Checking for edge-connected cells of highly differing sizes."
+            << endl;
+
+        // Pick up the cells that need to be removed and (a guess for)
+        // the patch they should be patched with.
+        Map<label> problemCells
+        (
+            findEdgeConnectedProblemCells
+            (
+                perpendicularAngle,
+                globalToPatch
+            )
+        );
+
+        // Baffle all faces of cells that need to be removed
+        forAllConstIter(Map<label>, problemCells, iter)
+        {
+            const cell& cFaces = mesh_.cells()[iter.key()];
+
+            forAll(cFaces, i)
+            {
+                label faceI = cFaces[i];
+
+                if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
+                {
+                    facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                    nBaffleFaces++;
+
+                    // Mark face as a 'boundary'
+                    markBoundaryFace
+                    (
+                        faceI,
+                        isBoundaryFace,
+                        isBoundaryEdge,
+                        isBoundaryPoint
+                    );
+                }
+            }
+        }
+        Info<< "markFacesOnProblemCells : Marked "
+            << returnReduce(nBaffleFaces, sumOp<label>())
+            << " additional internal faces to be converted into baffles"
+            << " due to "
+            << returnReduce(problemCells.size(), sumOp<label>())
+            << " cells edge-connected to lower level cells." << endl;
+
+        if (debug)
+        {
+            cellSet problemCellSet(mesh_, "problemCells", problemCells.toc());
+            Pout<< "Writing " << problemCellSet.size()
+                << " cells that are edge connected to coarser cell to set "
+                << problemCellSet.objectPath() << endl;
+            problemCellSet.write();
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh_,
+        isBoundaryPoint,
+        orEqOp<bool>(),
+        false,              // null value
+        false               // no separation
+    );
+
+    syncTools::syncEdgeList
+    (
+        mesh_,
+        isBoundaryEdge,
+        orEqOp<bool>(),
+        false,              // null value
+        false               // no separation
+    );
+
+    syncTools::syncFaceList
+    (
+        mesh_,
+        isBoundaryFace,
+        orEqOp<bool>(),
+        false               // no separation
+    );
+
+
+    // For each cell count the number of anchor points that are on
+    // the boundary:
+    // 8 : check the number of (baffle) boundary faces. If 3 or more block
+    //     off the cell since the cell would get squeezed down to a diamond
+    //     (probably; if the 3 or more faces are unrefined (only use the
+    //      anchor points))
+    // 7 : store. Used to check later on whether there are points with
+    //     3 or more of these cells. (note that on a flat surface a boundary
+    //     point will only have 4 cells connected to it)
+
+    // Does cell have exactly 7 of its 8 anchor points on the boundary?
+    PackedList<1> hasSevenBoundaryAnchorPoints(mesh_.nCells(), 0u);
+    // If so what is the remaining non-boundary anchor point?
+    labelHashSet nonBoundaryAnchors(mesh_.nCells()/10000);
+
+    // On-the-fly addressing storage.
+    DynamicList<label> dynFEdges;
+    DynamicList<label> dynCPoints;
+
+    forAll(cellLevel, cellI)
+    {
+        const labelList& cPoints = mesh_.cellPoints(cellI, dynCPoints);
+
+        // Get number of anchor points (pointLevel == cellLevel)
+
+        label nBoundaryAnchors = 0;
+        label nNonAnchorBoundary = 0;
+        label nonBoundaryAnchor = -1;
+
+        forAll(cPoints, i)
+        {
+            label pointI = cPoints[i];
+
+            if (pointLevel[pointI] <= cellLevel[cellI])
+            {
+                // Anchor point
+                if (isBoundaryPoint[pointI])
+                {
+                    nBoundaryAnchors++;
+                }
+                else
+                {
+                    // Anchor point which is not on the surface
+                    nonBoundaryAnchor = pointI;
+                }
+            }
+            else if (isBoundaryPoint[pointI])
+            {
+                nNonAnchorBoundary++;
+            }
+        }
+
+        if (nBoundaryAnchors == 8)
+        {
+            const cell& cFaces = mesh_.cells()[cellI];
+
+            // Count boundary faces.
+            label nBfaces = 0;
+
+            forAll(cFaces, cFaceI)
+            {
+                if (isBoundaryFace[cFaces[cFaceI]])
+                {
+                    nBfaces++;
+                }
+            }
+
+            // If nBfaces > 1 make all non-boundary non-baffle faces baffles.
+            // We assume that this situation is where there is a single
+            // cell sticking out which would get flattened.
+
+            // Eugene: delete cell no matter what.
+            //if (nBfaces > 1)
+            {
+                forAll(cFaces, cf)
+                {
+                    label faceI = cFaces[cf];
+
+                    if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
+                    {
+                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                        nBaffleFaces++;
+
+                        // Mark face as a 'boundary'
+                        markBoundaryFace
+                        (
+                            faceI,
+                            isBoundaryFace,
+                            isBoundaryEdge,
+                            isBoundaryPoint
+                        );
+                    }
+                }
+            }
+        }
+        else if (nBoundaryAnchors == 7)
+        {
+            // Mark the cell. Store the (single!) non-boundary anchor point.
+            hasSevenBoundaryAnchorPoints.set(cellI, 1u);
+            nonBoundaryAnchors.insert(nonBoundaryAnchor);
+        }
+    }
+
+
+    // Loop over all points. If a point is connected to 4 or more cells
+    // with 7 anchor points on the boundary set those cell's non-boundary faces
+    // to baffles
+
+    DynamicList<label> dynPCells;
+
+    forAllConstIter(labelHashSet, nonBoundaryAnchors, iter)
+    {
+        label pointI = iter.key();
+
+        const labelList& pCells = mesh_.pointCells(pointI, dynPCells);
+
+        // Count number of 'hasSevenBoundaryAnchorPoints' cells.
+        label n = 0;
+
+        forAll(pCells, i)
+        {
+            if (hasSevenBoundaryAnchorPoints.get(pCells[i]) == 1u)
+            {
+                n++;
+            }
+        }
+
+        if (n > 3)
+        {
+            // Point in danger of being what? Remove all 7-cells.
+            forAll(pCells, i)
+            {
+                label cellI = pCells[i];
+
+                if (hasSevenBoundaryAnchorPoints.get(cellI) == 1u)
+                {
+                    const cell& cFaces = mesh_.cells()[cellI];
+
+                    forAll(cFaces, cf)
+                    {
+                        label faceI = cFaces[cf];
+
+                        if
+                        (
+                            facePatch[faceI] == -1
+                         && mesh_.isInternalFace(faceI)
+                        )
+                        {
+                            facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                            nBaffleFaces++;
+
+                            // Mark face as a 'boundary'
+                            markBoundaryFace
+                            (
+                                faceI,
+                                isBoundaryFace,
+                                isBoundaryEdge,
+                                isBoundaryPoint
+                            );
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+
+    // Sync all. (note that pointdata and facedata not used anymore but sync
+    // anyway)
+
+    syncTools::syncPointList
+    (
+        mesh_,
+        isBoundaryPoint,
+        orEqOp<bool>(),
+        false,              // null value
+        false               // no separation
+    );
+
+    syncTools::syncEdgeList
+    (
+        mesh_,
+        isBoundaryEdge,
+        orEqOp<bool>(),
+        false,              // null value
+        false               // no separation
+    );
+
+    syncTools::syncFaceList
+    (
+        mesh_,
+        isBoundaryFace,
+        orEqOp<bool>(),
+        false               // no separation
+    );
+
+
+    // Find faces with all edges on the boundary and make them baffles
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        if (facePatch[faceI] == -1)
+        {
+            const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
+            label nFaceBoundaryEdges = 0;
+
+            forAll(fEdges, fe)
+            {
+                if (isBoundaryEdge[fEdges[fe]])
+                {
+                    nFaceBoundaryEdges++;
+                }
+            }
+
+            if (nFaceBoundaryEdges == fEdges.size())
+            {
+                facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                nBaffleFaces++;
+
+                // Do NOT update boundary data since this would grow blocked
+                // faces across gaps.
+            }
+        }
+    }
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                if (facePatch[faceI] == -1)
+                {
+                    const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
+                    label nFaceBoundaryEdges = 0;
+
+                    forAll(fEdges, fe)
+                    {
+                        if (isBoundaryEdge[fEdges[fe]])
+                        {
+                            nFaceBoundaryEdges++;
+                        }
+                    }
+
+                    if (nFaceBoundaryEdges == fEdges.size())
+                    {
+                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                        nBaffleFaces++;
+
+                        // Do NOT update boundary data since this would grow
+                        // blocked faces across gaps.
+                    }
+                }
+
+                faceI++;
+            }
+        }
+    }
+
+    Info<< "markFacesOnProblemCells : marked "
+        << returnReduce(nBaffleFaces, sumOp<label>())
+        << " additional internal faces to be converted into baffles."
+        << endl;
+
+    return facePatch;
+}
+//XXXXXXXXXXXXXX
+// Mark faces to be baffled to prevent snapping problems. Does
+// test to find nearest surface and checks which faces would get squashed.
+Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
+(
+    const dictionary& motionDict,
+    const labelList& globalToPatch
+) const
+{
+    // Get the labels of added patches.
+    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
+
+    // Construct addressing engine.
+    autoPtr<indirectPrimitivePatch> ppPtr
+    (
+        meshRefinement::makePatch
+        (
+            mesh_,
+            adaptPatchIDs
+        )
+    );
+    const indirectPrimitivePatch& pp = ppPtr();
+    const pointField& localPoints = pp.localPoints();
+    const labelList& meshPoints = pp.meshPoints();
+
+    // Find nearest (non-baffle) surface
+    pointField newPoints(mesh_.points());
+    {
+        List<pointIndexHit> hitInfo;
+        labelList hitSurface;
+        surfaces_.findNearest
+        (
+            surfaces_.getUnnamedSurfaces(),
+            localPoints,
+            scalarField(localPoints.size(), sqr(GREAT)),    // sqr of attraction
+            hitSurface,
+            hitInfo
+        );
+    
+        forAll(hitInfo, i)
+        {
+            if (hitInfo[i].hit())
+            {
+                //label pointI = meshPoints[i];
+                //Pout<< "   " << pointI << " moved from "
+                //    << mesh_.points()[pointI] << " by "
+                //    << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
+                //    << endl;
+                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
+            }
+        }
+    }
+
+    // Per face (internal or coupled!) the patch that the
+    // baffle should get (or -1).
+    labelList facePatch(mesh_.nFaces(), -1);
+    // Count of baffled faces
+    label nBaffleFaces = 0;
+
+
+//    // Sync position? Or not since same face on both side so just sync
+//    // result of baffle.
+//
+//    const scalar minArea(readScalar(motionDict.lookup("minArea")));
+//
+//    Pout<< "markFacesOnProblemCellsGeometric : Comparing to minArea:"
+//        << minArea << endl;
+//
+//    pointField facePoints;
+//    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
+//    {
+//        const face& f = mesh_.faces()[faceI];
+//
+//        bool usesPatchPoint = false;
+//
+//        facePoints.setSize(f.size());
+//        forAll(f, fp)
+//        {
+//            Map<label>::const_iterator iter = pp.meshPointMap().find(f[fp]);
+//
+//            if (iter != pp.meshPointMap().end())
+//            {
+//                facePoints[fp] = newPosition[iter()];
+//                usesPatchPoint = true;
+//            }
+//            else
+//            {
+//                facePoints[fp] = mesh_.points()[f[fp]];
+//            }
+//        }
+//
+//        if (usesPatchPoint)
+//        {
+//            // Check area of face wrt original area
+//            face identFace(identity(f.size()));
+//
+//            if (identFace.mag(facePoints) < minArea)
+//            {
+//                facePatch[faceI] = getBafflePatch(facePatch, faceI);
+//                nBaffleFaces++;
+//            }
+//        }
+//    }
+//
+//
+//    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+//    forAll(patches, patchI)
+//    {
+//        const polyPatch& pp = patches[patchI];
+//
+//        if (pp.coupled())
+//        {
+//            forAll(pp, i)
+//            {
+//                label faceI = pp.start()+i;
+//
+//                const face& f = mesh_.faces()[faceI];
+//
+//                bool usesPatchPoint = false;
+//
+//                facePoints.setSize(f.size());
+//                forAll(f, fp)
+//                {
+//                    Map<label>::const_iterator iter =
+//                        pp.meshPointMap().find(f[fp]);
+//
+//                    if (iter != pp.meshPointMap().end())
+//                    {
+//                        facePoints[fp] = newPosition[iter()];
+//                        usesPatchPoint = true;
+//                    }
+//                    else
+//                    {
+//                        facePoints[fp] = mesh_.points()[f[fp]];
+//                    }
+//                }
+//
+//                if (usesPatchPoint)
+//                {
+//                    // Check area of face wrt original area
+//                    face identFace(identity(f.size()));
+//
+//                    if (identFace.mag(facePoints) < minArea)
+//                    {
+//                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
+//                        nBaffleFaces++;
+//                    }
+//                }
+//            }
+//        }
+//    }
+
+    {
+        pointField oldPoints(mesh_.points());
+        mesh_.movePoints(newPoints);
+        faceSet wrongFaces(mesh_, "wrongFaces", 100);
+        {
+            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
+
+            // Just check the errors from squashing
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            const labelList allFaces(identity(mesh_.nFaces()));
+            label nWrongFaces = 0;
+
+            scalar minArea(readScalar(motionDict.lookup("minArea")));
+            if (minArea > -SMALL)
+            {
+                polyMeshGeometry::checkFaceArea
+                (
+                    false,
+                    minArea,
+                    mesh_,
+                    mesh_.faceAreas(),
+                    allFaces,
+                    &wrongFaces
+                );
+
+                label nNewWrongFaces = returnReduce
+                (
+                    wrongFaces.size(),
+                    sumOp<label>()
+                );
+
+                Info<< "    faces with area < "
+                    << setw(5) << minArea
+                    << " m^2                            : "
+                    << nNewWrongFaces-nWrongFaces << endl;
+
+                nWrongFaces = nNewWrongFaces;
+            }
+
+//            scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
+            scalar minDet = 0.01;
+            if (minDet > -1)
+            {
+                polyMeshGeometry::checkCellDeterminant
+                (
+                    false,
+                    minDet,
+                    mesh_,
+                    mesh_.faceAreas(),
+                    allFaces,
+                    polyMeshGeometry::affectedCells(mesh_, allFaces),
+                    &wrongFaces
+                );
+
+                label nNewWrongFaces = returnReduce
+                (
+                    wrongFaces.size(),
+                    sumOp<label>()
+                );
+
+                Info<< "    faces on cells with determinant < "
+                    << setw(5) << minDet << "                : "
+                    << nNewWrongFaces-nWrongFaces << endl;
+
+                nWrongFaces = nNewWrongFaces;
+            }
+        }
+
+
+        forAllConstIter(faceSet, wrongFaces, iter)
+        {
+            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
+
+            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
+            {
+                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
+                nBaffleFaces++;
+
+                //Pout<< "    " << iter.key()
+                //    //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
+                //    << " is destined for patch " << facePatch[iter.key()]
+                //    << endl;
+            }
+        }
+        // Restore points.
+        mesh_.movePoints(oldPoints);
+    }
+
+
+    Info<< "markFacesOnProblemCellsGeometric : marked "
+        << returnReduce(nBaffleFaces, sumOp<label>())
+        << " additional internal and coupled faces"
+        << " to be converted into baffles." << endl;
+
+    syncTools::syncFaceList
+    (
+        mesh_,
+        facePatch,
+        maxEqOp<label>(),
+        false               // no separation
+    );
+
+    return facePatch;
+}
+//XXXXXXXX
+
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -53,8 +53,10 @@ Foam::refinementSurfaces::refinementSurfaces
 {
     labelList globalMinLevel(surfaceDicts.size(), 0);
     labelList globalMaxLevel(surfaceDicts.size(), 0);
+    scalarField globalAngle(surfaceDicts.size(), -GREAT);
     List<Map<label> > regionMinLevel(surfaceDicts.size());
     List<Map<label> > regionMaxLevel(surfaceDicts.size());
+    List<Map<scalar> > regionAngle(surfaceDicts.size());
 
     //wordList globalPatchType(surfaceDicts.size());
     //List<HashTable<word> > regionPatchType(surfaceDicts.size());
@@ -80,6 +82,12 @@ Foam::refinementSurfaces::refinementSurfaces
             dict.lookup("zoneInside") >> zoneInside_[surfI];
         }
 
+        // Global perpendicular angle
+        if (dict.found("perpendicularAngle"))
+        {
+            globalAngle[surfI] = readScalar(dict.lookup("perpendicularAngle"));
+        }
+
         //// Global patch name per surface
         //if (dict.found("patchType"))
         //{
@@ -130,6 +138,15 @@ Foam::refinementSurfaces::refinementSurfaces
                         << exit(FatalError);
                 }
                 regionMaxLevel[surfI].insert(regionI, max);
+
+                if (regionDict.found("perpendicularAngle"))
+                {
+                    regionAngle[surfI].insert
+                    (
+                        regionI,
+                        readScalar(regionDict.lookup("perpendicularAngle"))
+                    );
+                }
             }
         }
     }
@@ -170,6 +187,8 @@ Foam::refinementSurfaces::refinementSurfaces
     minLevel_ = 0;
     maxLevel_.setSize(nRegions);
     maxLevel_ = 0;
+    perpendicularAngle_.setSize(nRegions);
+    perpendicularAngle_ = -GREAT;
     //patchName_.setSize(nRegions);
     //patchType_.setSize(nRegions);
 
@@ -182,6 +201,7 @@ Foam::refinementSurfaces::refinementSurfaces
         {
             minLevel_[regionOffset_[surfI] + i] = globalMinLevel[surfI];
             maxLevel_[regionOffset_[surfI] + i] = globalMaxLevel[surfI];
+            perpendicularAngle_[regionOffset_[surfI] + i] = globalAngle[surfI];
         }
 
         // Overwrite with region specific information
@@ -191,6 +211,7 @@ Foam::refinementSurfaces::refinementSurfaces
 
             minLevel_[globalRegionI] = iter();
             maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
+            perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
 
             // Check validity
             if
@@ -240,8 +261,10 @@ Foam::refinementSurfaces::refinementSurfaces
 {
     labelList globalMinLevel(surfacesDict.size(), 0);
     labelList globalMaxLevel(surfacesDict.size(), 0);
+    scalarField globalAngle(surfacesDict.size(), -GREAT);
     List<Map<label> > regionMinLevel(surfacesDict.size());
     List<Map<label> > regionMaxLevel(surfacesDict.size());
+    List<Map<scalar> > regionAngle(surfacesDict.size());
 
     label surfI = 0;
     forAllConstIter(dictionary, surfacesDict, iter)
@@ -274,6 +297,12 @@ Foam::refinementSurfaces::refinementSurfaces
             dict.lookup("zoneInside") >> zoneInside_[surfI];
         }
 
+        // Global perpendicular angle
+        if (dict.found("perpendicularAngle"))
+        {
+            globalAngle[surfI] = readScalar(dict.lookup("perpendicularAngle"));
+        }
+
         if (dict.found("regions"))
         {
             const dictionary& regionsDict = dict.subDict("regions");
@@ -306,6 +335,15 @@ Foam::refinementSurfaces::refinementSurfaces
 
                     regionMinLevel[surfI].insert(regionI, refLevel[0]);
                     regionMaxLevel[surfI].insert(regionI, refLevel[1]);
+
+                    if (regionDict.found("perpendicularAngle"))
+                    {
+                        regionAngle[surfI].insert
+                        (
+                            regionI,
+                            readScalar(regionDict.lookup("perpendicularAngle"))
+                        );
+                    }
                 }
             }
         }
@@ -326,6 +364,8 @@ Foam::refinementSurfaces::refinementSurfaces
     minLevel_ = 0;
     maxLevel_.setSize(nRegions);
     maxLevel_ = 0;
+    perpendicularAngle_.setSize(nRegions);
+    perpendicularAngle_ = -GREAT;
 
 
     forAll(globalMinLevel, surfI)
@@ -337,6 +377,7 @@ Foam::refinementSurfaces::refinementSurfaces
         {
             minLevel_[regionOffset_[surfI] + i] = globalMinLevel[surfI];
             maxLevel_[regionOffset_[surfI] + i] = globalMaxLevel[surfI];
+            perpendicularAngle_[regionOffset_[surfI] + i] = globalAngle[surfI];
         }
 
         // Overwrite with region specific information
@@ -346,6 +387,7 @@ Foam::refinementSurfaces::refinementSurfaces
 
             minLevel_[globalRegionI] = iter();
             maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
+            perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
 
             // Check validity
             if
@@ -454,8 +496,6 @@ void Foam::refinementSurfaces::setMinLevelFields
     const shellSurfaces& shells
 )
 {
-    //minLevelFields_.setSize(surfaces_.size());
-
     forAll(surfaces_, surfI)
     {
         const searchableSurface& geom = allGeometry_[surfaces_[surfI]];

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -90,6 +90,9 @@ class refinementSurfaces
         //- From global region number to refinement level
         labelList maxLevel_;
 
+        //- From global region number to perpendicular angle
+        scalarField perpendicularAngle_;
+
 
     // Private Member Functions
 
@@ -178,6 +181,12 @@ public:
                 return maxLevel_;
             }
 
+            //- From global region number to perpendicular angle
+            const scalarField& perpendicularAngle() const
+            {
+                return perpendicularAngle_;
+            }
+
 
         // Helper
 

src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C

@@ -59,7 +59,7 @@ const NamedEnum<shellSurfaces::refineMode, 3> shellSurfaces::refineModeNames_;
 
 void Foam::shellSurfaces::setAndCheckLevels
 (
-    const scalar shellI,
+    const label shellI,
     const List<Tuple2<scalar, label> >& distLevels
 )
 {

src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.H

@@ -97,7 +97,7 @@ private:
         //- Helper function for initialisation.
         void setAndCheckLevels
         (
-            const scalar shellI,
+            const label shellI,
             const List<Tuple2<scalar, label> >&
         );
 

src/dynamicMesh/polyMeshAdder/polyMeshAdder.C

@@ -61,6 +61,29 @@ void Foam::polyMeshAdder::append
 }
 
 
+//- Append all mapped elements of a list to a DynamicList
+void Foam::polyMeshAdder::append
+(
+    const labelList& map,
+    const labelList& lst,
+    const SortableList<label>& sortedLst,
+    DynamicList<label>& dynLst
+)
+{
+    dynLst.setSize(dynLst.size() + lst.size());
+
+    forAll(lst, i)
+    {
+        label newElem = map[lst[i]];
+
+        if (newElem != -1 && findSortedIndex(sortedLst, newElem) == -1)
+        {
+            dynLst.append(newElem);
+        }
+    }
+}
+
+
 // Get index of patch in new set of patchnames/types
 Foam::label Foam::polyMeshAdder::patchIndex
 (
@@ -919,22 +942,33 @@ void Foam::polyMeshAdder::mergePointZones
 
     forAll(pz0, zoneI)
     {
-        DynamicList<label>& newZone = pzPoints[zoneI];
+        append(from0ToAllPoints, pz0[zoneI], pzPoints[zoneI]);
+    }
 
-        newZone.setCapacity(pz0[zoneI].size());
+    // Get sorted zone contents for duplicate element recognition
-
+    PtrList<SortableList<label> > pzPointsSorted(pzPoints.size());
-        append(from0ToAllPoints, pz0[zoneI], newZone);
+    forAll(pzPoints, zoneI)
+    {
+        pzPointsSorted.set
+        (
+            zoneI,
+            new SortableList<label>(pzPoints[zoneI])
+        );
     }
 
     // Now we have full addressing for points so do the pointZones of mesh1.
     forAll(pz1, zoneI)
     {
         // Relabel all points of zone and add to correct pzPoints.
-        DynamicList<label>& newZone = pzPoints[from1ToAll[zoneI]];
+        label allZoneI = from1ToAll[zoneI];
 
-        newZone.setCapacity(newZone.size() + pz1[zoneI].size());
+        append
-
+        (
-        append(from1ToAllPoints, pz1[zoneI], newZone);
+            from1ToAllPoints,
+            pz1[zoneI],
+            pzPointsSorted[allZoneI],
+            pzPoints[allZoneI]
+        );
     }
 
     forAll(pzPoints, i)
@@ -996,11 +1030,25 @@ void Foam::polyMeshAdder::mergeFaceZones
         }
     }
 
+    // Get sorted zone contents for duplicate element recognition
+    PtrList<SortableList<label> > fzFacesSorted(fzFaces.size());
+    forAll(fzFaces, zoneI)
+    {
+        fzFacesSorted.set
+        (
+            zoneI,
+            new SortableList<label>(fzFaces[zoneI])
+        );
+    }
+
     // Now we have full addressing for faces so do the faceZones of mesh1.
     forAll(fz1, zoneI)
     {
-        DynamicList<label>& newZone = fzFaces[from1ToAll[zoneI]];
+        label allZoneI = from1ToAll[zoneI];
-        DynamicList<bool>& newFlip = fzFlips[from1ToAll[zoneI]];
+
+        DynamicList<label>& newZone = fzFaces[allZoneI];
+        const SortableList<label>& newZoneSorted = fzFacesSorted[allZoneI];
+        DynamicList<bool>& newFlip = fzFlips[allZoneI];
 
         newZone.setCapacity(newZone.size() + fz1[zoneI].size());
         newFlip.setCapacity(newZone.size());
@@ -1011,10 +1059,15 @@ void Foam::polyMeshAdder::mergeFaceZones
         forAll(addressing, i)
         {
             label faceI = addressing[i];
+            label allFaceI = from1ToAllFaces[faceI];
 
-            if (from1ToAllFaces[faceI] != -1)
+            if
+            (
+                allFaceI != -1
+             && findSortedIndex(newZoneSorted, allFaceI) == -1
+            )
             {
-                newZone.append(from1ToAllFaces[faceI]);
+                newZone.append(allFaceI);
                 newFlip.append(flipMap[i]);
             }
         }
@@ -1055,7 +1108,6 @@ void Foam::polyMeshAdder::mergeCellZones
     czCells.setSize(zoneNames.size());
     forAll(cz0, zoneI)
     {
-        czCells[zoneI].setCapacity(cz0[zoneI].size());
         // Insert mesh0 cells
         append(cz0[zoneI], czCells[zoneI]);
     }
@@ -1064,11 +1116,9 @@ void Foam::polyMeshAdder::mergeCellZones
     // Cell mapping is trivial.
     forAll(cz1, zoneI)
     {
-        DynamicList<label>& newZone = czCells[from1ToAll[zoneI]];
+        label allZoneI = from1ToAll[zoneI];
 
-        newZone.setCapacity(newZone.size() + cz1[zoneI].size());
+        append(from1ToAllCells, cz1[zoneI], czCells[allZoneI]);
-
-        append(from1ToAllCells, cz1[zoneI], newZone);
     }
 
     forAll(czCells, i)

src/dynamicMesh/polyMeshAdder/polyMeshAdder.H

@@ -29,7 +29,7 @@ Description
     Adds two meshes without using any polyMesh morphing.
 
     Gets faces to couple as faceCoupleInfo which is list of faces on both
-    meshes. Holds map from last mesh addition.
+    meshes. Returns map from last mesh addition.
 
 SourceFiles
     polyMeshAdder.C
@@ -43,6 +43,7 @@ SourceFiles
 #include "polyMesh.H"
 #include "mapAddedPolyMesh.H"
 #include "faceCoupleInfo.H"
+#include "SortableList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -78,6 +79,16 @@ private:
             DynamicList<label>&
         );
 
+        //- Append all mapped elements of a list to a DynamicList that are
+        //  not already present in the sorted list.
+        static void append
+        (
+            const labelList& map,
+            const labelList& lst,
+            const SortableList<label>& sortedLst,
+            DynamicList<label>&
+        );
+
         //- Index of patch in allPatches. Add if nonexisting.
         static label patchIndex
         (

src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.C

@@ -175,9 +175,8 @@ void turbulentInletFvPatchField<Type>::updateCoeffs()
             ranGen_.randomise(randomField[facei]);
         }
 
-        // Correction-factor proposed by Yi Wang to compensate for the loss
+        // Correction-factor to compensate for the loss of RMS fluctuation
-        // of RMS fluctuation due to the temporal correlation introduced by
+        // due to the temporal correlation introduced by the alpha parameter.
-        // the alpha parameter.
         scalar rmsCorr = sqrt(12*(2*alpha_ - sqr(alpha_)))/alpha_;
 
         patchField =
@@ -206,6 +205,7 @@ void turbulentInletFvPatchField<Type>::write(Ostream& os) const
     os.writeKeyword("fluctuationScale")
         << fluctuationScale_ << token::END_STATEMENT << nl;
     referenceField_.writeEntry("referenceField", os);
+    os.writeKeyword("alpha") << alpha_ << token::END_STATEMENT << nl;
     this->writeEntry("value", os);
 }
 

src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.H

@@ -63,7 +63,6 @@ class turbulentInletFvPatchField
 :
     public fixedValueFvPatchField<Type>
 {
-
     // Private data
 
         Random ranGen_;

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C

@@ -35,65 +35,6 @@ License
 
 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
 
-template<class ParcelType>
-void Foam::KinematicCloud<ParcelType>::setInjectorCellAndPosition
-(
-    label& pCell,
-    vector& pPosition
-)
-{
-    const vector originalPosition = pPosition;
-
-    bool foundCell = false;
-
-    pCell = mesh_.findCell(pPosition);
-
-    if (pCell >= 0)
-    {
-        const vector& C = mesh_.C()[pCell];
-        pPosition += 1.0e-6*(C - pPosition);
-
-        foundCell = mesh_.pointInCell
-        (
-            pPosition,
-            pCell
-        );
-    }
-    reduce(foundCell, orOp<bool>());
-
-    // Last chance - find nearest cell and try that one
-    // - the point is probably on an edge
-    if (!foundCell)
-    {
-        pCell =  mesh_.findNearestCell(pPosition);
-
-        if (pCell >= 0)
-        {
-            const vector& C = mesh_.C()[pCell];
-            pPosition += 1.0e-6*(C - pPosition);
-
-            foundCell = mesh_.pointInCell
-            (
-                pPosition,
-                pCell
-            );
-        }
-        reduce(foundCell, orOp<bool>());
-    }
-
-    if (!foundCell)
-    {
-        FatalErrorIn
-        (
-            "void KinematicCloud<ParcelType>::findInjectorCell"
-            "(label&, vector&)"
-        )<< "Cannot find parcel injection cell. "
-         << "Parcel position = " << originalPosition << nl
-         << abort(FatalError);
-    }
-}
-
-
 template<class ParcelType>
 Foam::scalar Foam::KinematicCloud<ParcelType>::setNumberOfParticles
 (
@@ -324,7 +265,7 @@ void Foam::KinematicCloud<ParcelType>::evolve()
         g_.value()
     );
 
-    inject(td);
+    inject();
 
     if (coupled_)
     {
@@ -336,15 +277,11 @@ void Foam::KinematicCloud<ParcelType>::evolve()
 
 
 template<class ParcelType>
-template<class TrackingData>
+void Foam::KinematicCloud<ParcelType>::inject()
-void Foam::KinematicCloud<ParcelType>::inject
-(
-    TrackingData& td
-)
 {
     scalar time = this->db().time().value();
 
-    scalar pRho = td.constProps().rho0();
+    scalar pRho = constProps_.rho0();
 
     this->injection().prepareForNextTimeStep(time0_, time);
 
@@ -419,21 +356,21 @@ void Foam::KinematicCloud<ParcelType>::inject
 
         // Determine the injection cell
         label pCell = -1;
-        setInjectorCellAndPosition(pCell, pPosition);
+        this->injection().findInjectorCellAndPosition(pCell, pPosition);
 
         if (pCell >= 0)
         {
             // construct the parcel that is to be injected
             ParcelType* pPtr = new ParcelType
             (
-                td.cloud(),
+                *this,
                 parcelTypeId_,
                 pPosition,
                 pCell,
                 pDiameter,
                 pU,
                 pNumberOfParticles,
-                td.constProps()
+                constProps_
             );
 
             scalar dt = time - timeInj;
@@ -441,7 +378,7 @@ void Foam::KinematicCloud<ParcelType>::inject
             pPtr->stepFraction() = (this->db().time().deltaT().value() - dt)
                 /this->time().deltaT().value();
 
-            this->injectParcel(td, pPtr);
+            this->injectParcel(pPtr);
          }
     }
 
@@ -455,12 +392,7 @@ void Foam::KinematicCloud<ParcelType>::inject
 
 
 template<class ParcelType>
-template<class TrackingData>
+void Foam::KinematicCloud<ParcelType>::injectParcel(ParcelType* p)
-void Foam::KinematicCloud<ParcelType>::injectParcel
-(
-    TrackingData& td,
-    ParcelType* p
-)
 {
     addParticle(p);
     nParcelsAdded_++;

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H

@@ -223,13 +223,6 @@ protected:
 
     // Protected member functions
 
-        //- Set parcel position and cell into which parcel is introduced
-        void setInjectorCellAndPosition
-        (
-            label& pCell,
-            vector& pPosition
-        );
-
         //- Set the number of particles per parcel
         scalar setNumberOfParticles
         (
@@ -241,16 +234,10 @@ protected:
         );
 
         //- Inject more parcels
-        template<class TrackingData>
+        void inject();
-        void inject(TrackingData& td);
 
         //- Inject parcel if it is valid - delete otherwise
-        template<class TrackingData>
+        void injectParcel(ParcelType* p);
-        void injectParcel
-        (
-            TrackingData& td,
-            ParcelType* p
-        );
 
         //- Post-injection checks
         void postInjectCheck();

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C

@@ -174,7 +174,7 @@ void Foam::ReactingCloud<ParcelType>::evolve()
         this->g().value()
     );
 
-    inject(td);
+    inject();
 
     if (this->coupled())
     {
@@ -186,15 +186,11 @@ void Foam::ReactingCloud<ParcelType>::evolve()
 
 
 template<class ParcelType>
-template<class TrackingData>
+void Foam::ReactingCloud<ParcelType>::inject()
-void Foam::ReactingCloud<ParcelType>::inject
-(
-    TrackingData& td
-)
 {
     scalar time = this->db().time().value();
 
-    scalar pRho = td.constProps().rho0();
+    scalar pRho = this->constProps().rho0();
 
     this->injection().prepareForNextTimeStep(this->time0(), time);
 
@@ -269,14 +265,14 @@ void Foam::ReactingCloud<ParcelType>::inject
 
         // Determine the injection cell
         label pCell = -1;
-        this->setInjectorCellAndPosition(pCell, pPosition);
+        this->injection().findInjectorCellAndPosition(pCell, pPosition);
 
         if (pCell >= 0)
         {
             // construct the parcel that is to be injected
             ParcelType* pPtr = new ParcelType
             (
-                td.cloud(),
+                *this,
                 this->parcelTypeId(),
                 pPosition,
                 pCell,
@@ -287,7 +283,7 @@ void Foam::ReactingCloud<ParcelType>::inject
                 composition().YLiquid0(),
                 composition().YSolid0(),
                 composition().YMixture0(),
-                td.constProps()
+                this->constProps()
             );
 
             scalar dt = time - timeInj;
@@ -295,7 +291,7 @@ void Foam::ReactingCloud<ParcelType>::inject
             pPtr->stepFraction() = (this->db().time().deltaT().value() - dt)
                 /this->db().time().deltaT().value();
 
-            this->injectParcel(td, pPtr);
+            this->injectParcel(pPtr);
          }
     }
 

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H

@@ -117,8 +117,7 @@ class ReactingCloud
 protected:
 
     //- Inject more parcels
-    template<class TrackingData>
+    void inject();
-    void inject(TrackingData& td);
 
 
 public:

src/lagrangian/intermediate/submodels/IO/DataEntry/Table/Table.C

@@ -78,11 +78,14 @@ Type Foam::Table<Type>::value(const scalar x) const
         i++;
     }
 
-    // Linear interpolation to find value
+    // Linear interpolation to find value. Note constructor needed for
-    return
+    // Table<label> to convert intermediate scalar back to label.
+    return Type
+    (
         (x - table_[i].first())/(table_[i+1].first() - table_[i].first())
       * (table_[i+1].second() - table_[i].second())
-      + table_[i].second();
+      + table_[i].second()
+    );
 }
 
 

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C

@@ -159,6 +159,57 @@ void Foam::InjectionModel<CloudType>::prepareForNextTimeStep
 }
 
 
+template<class CloudType>
+void Foam::InjectionModel<CloudType>::findInjectorCellAndPosition
+(
+    label& cellI,
+    vector& position
+)
+{
+    const vector p0 = position;
+
+    bool foundCell = false;
+
+    cellI = owner_.mesh().findCell(position);
+
+    if (cellI >= 0)
+    {
+        const vector& C = owner_.mesh().C()[cellI];
+        position += 1.0e-6*(C - position);
+
+        foundCell = owner_.mesh().pointInCell(position, cellI);
+    }
+    reduce(foundCell, orOp<bool>());
+
+    // Last chance - find nearest cell and try that one
+    // - the point is probably on an edge
+    if (!foundCell)
+    {
+        cellI = owner_.mesh().findNearestCell(position);
+
+        if (cellI >= 0)
+        {
+            const vector& C = owner_.mesh().C()[cellI];
+            position += 1.0e-6*(C - position);
+
+            foundCell = owner_.mesh().pointInCell(position, cellI);
+        }
+        reduce(foundCell, orOp<bool>());
+    }
+
+    if (!foundCell)
+    {
+        FatalErrorIn
+        (
+            "InjectionModel<CloudType>::setInjectorCellAndPosition"
+            "(label&, vector&)"
+        )<< "Cannot find parcel injection cell. "
+         << "Parcel position = " << p0 << nl
+         << abort(FatalError);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "NewInjectionModel.C"

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H

@@ -208,6 +208,14 @@ public:
 
         // Injection geometry
 
+            //- Find the cell that contains the injector position
+            //  Will modify position slightly towards the owner cell centroid
+            virtual void findInjectorCellAndPosition
+            (
+                label& cellI,
+                vector& position
+            );
+
             //- Return the injection position
             virtual vector position
             (

src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C

@@ -98,7 +98,7 @@ Foam::scalar Foam::pairPotential::forceLookup(const scalar r) const
 {
     scalar k_rIJ = (r - rMin_)/dr_;
 
-    label k(k_rIJ);
+    label k = label(k_rIJ);
 
     if (k < 0)
     {
@@ -135,7 +135,7 @@ Foam::scalar Foam::pairPotential::energyLookup(const scalar r) const
 {
     scalar k_rIJ = (r - rMin_)/dr_;
 
-    label k(k_rIJ);
+    label k = label(k_rIJ);
 
     if (k < 0)
     {

src/postProcessing/Allwmake

@@ -3,10 +3,8 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake libo postCalc
-wmake libso forces
-wmake libso fieldAverage
 wmake libso foamCalcFunctions
-wmake libso minMaxFields
+
-wmake libso systemCall
+(cd functionObjects && ./Allwmake)
 
 # ----------------------------------------------------------------- end-of-file

src/postProcessing/functionObjects/Allwmake

@@ -0,0 +1,10 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake libso fieldAverage
+wmake libso forces
+wmake libso minMaxFields
+wmake libso systemCall
+
+# ----------------------------------------------------------------- end-of-file