Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev

src/OpenFOAM/db/dictionary/dictionary.C

@@ -465,7 +465,8 @@ const Foam::entry* Foam::dictionary::lookupScopedEntryPtr
                         *this
                     )   << "keyword " << keyword
                         << " is undefined in dictionary "
-                        << name()
+                        << name() << endl
+                        << "Valid keywords are " << keys()
                         << exit(FatalIOError);
                 }
                 if (!entPtr->isDict())

src/finiteVolume/fields/fvPatchFields/derived/uniformTotalPressure/uniformTotalPressureFvPatchScalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -44,7 +44,6 @@ uniformTotalPressureFvPatchScalarField
     rhoName_("none"),
     psiName_("none"),
     gamma_(0.0),
-    p0_(0.0),
     pressure_()
 {}
 
@@ -63,7 +62,6 @@ uniformTotalPressureFvPatchScalarField
     rhoName_(dict.lookupOrDefault<word>("rho", "none")),
     psiName_(dict.lookupOrDefault<word>("psi", "none")),
     gamma_(readScalar(dict.lookup("gamma"))),
-    p0_(readScalar(dict.lookup("p0"))),
     pressure_(DataEntry<scalar>::New("pressure", dict))
 {
     if (dict.found("value"))
@@ -75,7 +73,8 @@ uniformTotalPressureFvPatchScalarField
     }
     else
     {
-        fvPatchField<scalar>::operator=(p0_);
+        scalar p0 = pressure_->value(this->db().time().timeOutputValue());
+        fvPatchField<scalar>::operator=(p0);
     }
 }
 
@@ -95,7 +94,6 @@ uniformTotalPressureFvPatchScalarField
     rhoName_(ptf.rhoName_),
     psiName_(ptf.psiName_),
     gamma_(ptf.gamma_),
-    p0_(ptf.p0_),
     pressure_(ptf.pressure_().clone().ptr())
 {}
 
@@ -112,7 +110,6 @@ uniformTotalPressureFvPatchScalarField
     rhoName_(tppsf.rhoName_),
     psiName_(tppsf.psiName_),
     gamma_(tppsf.gamma_),
-    p0_(tppsf.p0_),
     pressure_(tppsf.pressure_().clone().ptr())
 {}
 
@@ -130,7 +127,6 @@ uniformTotalPressureFvPatchScalarField
     rhoName_(tppsf.rhoName_),
     psiName_(tppsf.psiName_),
     gamma_(tppsf.gamma_),
-    p0_(tppsf.p0_),
     pressure_(tppsf.pressure_().clone().ptr())
 {}
 
@@ -147,14 +143,14 @@ void Foam::uniformTotalPressureFvPatchScalarField::updateCoeffs
         return;
     }
 
-    p0_ = pressure_->value(this->db().time().timeOutputValue());
+    scalar p0 = pressure_->value(this->db().time().timeOutputValue());
 
     const fvsPatchField<scalar>& phip =
         patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);
 
     if (psiName_ == "none" && rhoName_ == "none")
     {
-        operator==(p0_ - 0.5*(1.0 - pos(phip))*magSqr(Up));
+        operator==(p0 - 0.5*(1.0 - pos(phip))*magSqr(Up));
     }
     else if (rhoName_ == "none")
     {
@@ -167,7 +163,7 @@ void Foam::uniformTotalPressureFvPatchScalarField::updateCoeffs
 
             operator==
             (
-                p0_
+                p0
                /pow
                 (
                     (1.0 + 0.5*psip*gM1ByG*(1.0 - pos(phip))*magSqr(Up)),
@@ -177,7 +173,7 @@ void Foam::uniformTotalPressureFvPatchScalarField::updateCoeffs
         }
         else
         {
-            operator==(p0_/(1.0 + 0.5*psip*(1.0 - pos(phip))*magSqr(Up)));
+            operator==(p0/(1.0 + 0.5*psip*(1.0 - pos(phip))*magSqr(Up)));
         }
     }
     else if (psiName_ == "none")
@@ -185,7 +181,7 @@ void Foam::uniformTotalPressureFvPatchScalarField::updateCoeffs
         const fvPatchField<scalar>& rho =
             patch().lookupPatchField<volScalarField, scalar>(rhoName_);
 
-        operator==(p0_ - 0.5*rho*(1.0 - pos(phip))*magSqr(Up));
+        operator==(p0 - 0.5*rho*(1.0 - pos(phip))*magSqr(Up));
     }
     else
     {
@@ -219,7 +215,6 @@ void Foam::uniformTotalPressureFvPatchScalarField::write(Ostream& os) const
     os.writeKeyword("rho") << rhoName_ << token::END_STATEMENT << nl;
     os.writeKeyword("psi") << psiName_ << token::END_STATEMENT << nl;
     os.writeKeyword("gamma") << gamma_ << token::END_STATEMENT << nl;
-    os.writeKeyword("p0") << p0_ << token::END_STATEMENT << nl;
     pressure_->writeData(os);
     writeEntry("value", os);
 }

src/finiteVolume/fields/fvPatchFields/derived/uniformTotalPressure/uniformTotalPressureFvPatchScalarField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -75,9 +75,6 @@ class uniformTotalPressureFvPatchScalarField
         //- Heat capacity ratio
         scalar gamma_;
 
-        //- Total pressure
-        scalar p0_;
-
         //- Table of time vs total pressure, including the bounding treatment
         autoPtr<DataEntry<scalar> > pressure_;
 
@@ -178,18 +175,6 @@ public:
                 return gamma_;
             }
 
-            //- Return the total pressure
-            scalar p0() const
-            {
-                return p0_;
-            }
-
-            //- Return reference to the total pressure to allow adjustment
-            scalar p0()
-            {
-                return p0_;
-            }
-
 
         // Evaluation functions
 

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

@@ -470,8 +470,8 @@ public:
             //- Total rotational kinetic energy in the system
             inline scalar rotationalKineticEnergyOfSystem() const;
 
-            //- Penetration for percentage of the current total mass
+            //- Penetration for fraction [0-1] of the current total mass
-            inline scalar penetration(const scalar& prc) const;
+            inline scalar penetration(const scalar& fraction) const;
 
             //- Mean diameter Dij
             inline scalar Dij(const label i, const label j) const;

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

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "fvmSup.H"
+#include "SortableList.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -336,106 +337,130 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::Dmax() const
 
     reduce(d, maxOp<scalar>());
 
-    return d;
+    return max(0.0, d);
 }
 
 
 template<class CloudType>
 inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
 (
-    const scalar& prc
+    const scalar& fraction
 ) const
 {
-    scalar distance = 0.0;
+    if ((fraction < 0) || (fraction > 1))
-    scalar mTot = 0.0;
-
-    label np = this->size();
-
-    // arrays containing the parcels mass and
-    // distance from injector in ascending order
-    scalarField mass(np);
-    scalarField dist(np);
-
-    if (np > 0)
     {
-        label n = 0;
+        FatalErrorIn
-
+        (
-        // first arrange the parcels in ascending order
+            "inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration"
-        // the first parcel is closest to its injection position
+            "("
-        // and the last one is most far away.
+                "const scalar&"
-        forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
+            ") const"
-        {
+        )
-            const parcelType& p = iter();
+            << "fraction should be in the range 0 < fraction < 1"
-            scalar mi = p.nParticle()*p.mass();
+            << exit(FatalError);
-            scalar di = mag(p.position() - p.position0());
-            mTot += mi;
-
-            // insert at the last place
-            mass[n] = mi;
-            dist[n] = di;
-
-            label i = 0;
-            bool found = false;
-
-            // insert the parcel in the correct place
-            // and move the others
-            while ((i < n) && (!found))
-            {
-                if (di < dist[i])
-                {
-                    found = true;
-                    for (label j=n; j>i; j--)
-                    {
-                        mass[j] = mass[j-1];
-                        dist[j] = dist[j-1];
-                    }
-                    mass[i] = mi;
-                    dist[i] = di;
-                }
-                i++;
-            }
-            n++;
-        }
     }
 
-    reduce(mTot, sumOp<scalar>());
+    scalar distance = 0.0;
 
-    if (np > 0)
+    const label nParcel = this->size();
+    globalIndex globalParcels(nParcel);
+    const label nParcelSum = globalParcels.size();
+
+    if (nParcelSum == 0)
     {
-        scalar mLimit = prc*mTot;
+        return distance;
-        scalar mOff = (1.0 - prc)*mTot;
+    }
 
-        if (np > 1)
+    // lists of parcels mass and distance from initial injection point
+    List<scalar> mass(nParcel, 0.0);
+    List<scalar> dist(nParcel, 0.0);
+
+    label i = 0;
+    scalar mSum = 0.0;
+    forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
+    {
+        const parcelType& p = iter();
+        scalar m = p.nParticle()*p.mass();
+        scalar d = mag(p.position() - p.position0());
+        mSum += m;
+
+        mass[i] = m;
+        dist[i] = d;
+
+        i++;
+    }
+
+    // calculate total mass across all processors
+    reduce(mSum, sumOp<scalar>());
+
+    // flatten the mass list
+    List<scalar> allMass(nParcelSum, 0.0);
+    SubList<scalar>
+    (
+        allMass,
+        globalParcels.localSize(Pstream::myProcNo()),
+        globalParcels.offset(Pstream::myProcNo())
+    ).assign(mass);
+
+    // flatten the distance list
+    SortableList<scalar> allDist(nParcelSum, 0.0);
+    SubList<scalar>
+    (
+        allDist,
+        globalParcels.localSize(Pstream::myProcNo()),
+        globalParcels.offset(Pstream::myProcNo())
+    ).assign(dist);
+
+    // sort allDist distances into ascending order
+    // note: allMass masses are left unsorted
+    allDist.sort();
+
+    if (nParcelSum > 1)
+    {
+        const scalar mLimit = fraction*mSum;
+        const labelList& indices = allDist.indices();
+
+        if (mLimit > (mSum - allMass[indices.last()]))
         {
-            // 'prc' is large enough that the parcel most far
+            distance = allDist.last();
-            // away will be used, no need to loop...
-            if (mLimit > mTot - mass[np-1])
-            {
-                distance = dist[np-1];
-            }
-            else
-            {
-                scalar mOffSum = 0.0;
-                label i = np;
-
-                while ((mOffSum < mOff) && (i>0))
-                {
-                    i--;
-                    mOffSum += mass[i];
-                }
-                distance =
-                    dist[i+1]
-                  + (dist[i] - dist[i+1])*(mOffSum - mOff)
-                   /mass[i+1] ;
-            }
         }
         else
         {
-            distance = dist[0];
+            // assuming that 'fraction' is generally closer to 1 than 0, loop
+            // through in reverse distance order
+            const scalar mThreshold = (1.0 - fraction)*mSum;
+            scalar mCurrent = 0.0;
+            label i0 = 0;
+
+            forAllReverse(indices, i)
+            {
+                label indI = indices[i];
+
+                mCurrent += allMass[indI];
+
+                if (mCurrent > mThreshold)
+                {
+                    i0 = i;
+                    break;
+                }
+            }
+
+            if (i0 == indices.size() - 1)
+            {
+                distance = allDist.last();
+            }
+            else
+            {
+                // linearly interpolate to determine distance
+                scalar alpha = (mCurrent - mThreshold)/allMass[indices[i0]];
+                distance = allDist[i0] + alpha*(allDist[i0+1] - allDist[i0]);
+            }
         }
     }
-
+    else
-    reduce(distance, maxOp<scalar>());
+    {
+        distance = allDist.first();
+    }
 
     return distance;
 }

src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H

@@ -89,7 +89,7 @@ public:
             virtual scalar rotationalKineticEnergyOfSystem() const = 0;
 
             //- Penetration for percentage of the current total mass
-//            virtual scalar penetration(const scalar& prc) const = 0;
+//            virtual scalar penetration(const scalar& fraction) const = 0;
 
             //- Mean diameter Dij
             virtual scalar Dij(const label i, const label j) const = 0;

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C

@@ -636,7 +636,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
        *(sum(dMassSRGas) + sum(dMassSRLiquid) + sum(dMassSRSolid))
     );
 
-    const scalar xsi = min(T/5000.0, 1.0);
+    const scalar xsi = min(T/td.cloud().constProps().TMax(), 1.0);
     const scalar coeff =
         (1.0 - xsi*xsi)*td.cloud().constProps().hRetentionCoeff();
 

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "PressureGradientForce.H"
+#include "fvcDdt.H"
 #include "fvcGrad.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -33,12 +34,13 @@ Foam::PressureGradientForce<CloudType>::PressureGradientForce
 (
     CloudType& owner,
     const fvMesh& mesh,
-    const dictionary& dict
+    const dictionary& dict,
+    const word& forceType
 )
 :
-    ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
+    ParticleForce<CloudType>(owner, mesh, dict, forceType, true),
-    UName_(this->coeffs().lookup("U")),
+    UName_(this->coeffs().template lookupOrDefault<word>("U", "U")),
-    gradUPtr_(NULL)
+    DUcDtInterpPtr_(NULL)
 {}
 
 
@@ -50,7 +52,7 @@ Foam::PressureGradientForce<CloudType>::PressureGradientForce
 :
     ParticleForce<CloudType>(pgf),
     UName_(pgf.UName_),
-    gradUPtr_(NULL)
+    DUcDtInterpPtr_(NULL)
 {}
 
 
@@ -66,18 +68,48 @@ Foam::PressureGradientForce<CloudType>::~PressureGradientForce()
 template<class CloudType>
 void Foam::PressureGradientForce<CloudType>::cacheFields(const bool store)
 {
+    static word fName("DUcDt");
+
+    bool fieldExists = this->mesh().template foundObject<volVectorField>(fName);
+
     if (store)
     {
-        const volVectorField& U = this->mesh().template
+        if (!fieldExists)
-            lookupObject<volVectorField>(UName_);
+        {
-        gradUPtr_ = fvc::grad(U).ptr();
+            const volVectorField& Uc = this->mesh().template
+                lookupObject<volVectorField>(UName_);
+
+            volVectorField* DUcDtPtr = new volVectorField
+            (
+                fName,
+                fvc::ddt(Uc) + (Uc & fvc::grad(Uc))
+            );
+
+            DUcDtPtr->store();
+        }
+
+        const volVectorField& DUcDt = this->mesh().template
+            lookupObject<volVectorField>(fName);
+
+        DUcDtInterpPtr_.reset
+        (
+            interpolation<vector>::New
+            (
+                this->owner().solution().interpolationSchemes(),
+                DUcDt
+            ).ptr()
+        );
     }
     else
     {
-        if (gradUPtr_)
+        DUcDtInterpPtr_.clear();
+
+        if (fieldExists)
         {
-            delete gradUPtr_;
+            const volVectorField& DUcDt = this->mesh().template
-            gradUPtr_ = NULL;
+                lookupObject<volVectorField>(fName);
+
+            const_cast<volVectorField&>(DUcDt).checkOut();
         }
     }
 }
@@ -95,11 +127,24 @@ Foam::forceSuSp Foam::PressureGradientForce<CloudType>::calcCoupled
 {
     forceSuSp value(vector::zero, 0.0);
 
-    const volTensorField& gradU = *gradUPtr_;
+    vector DUcDt =
-    value.Su() = mass*p.rhoc()/p.rho()*(p.U() & gradU[p.cell()]);
+        DUcDtInterp().interpolate(p.position(), p.currentTetIndices());
+
+    value.Su() = mass*p.rhoc()/p.rho()*DUcDt;
 
     return value;
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::PressureGradientForce<CloudType>::massAdd
+(
+    const typename CloudType::parcelType&,
+    const scalar
+) const
+{
+    return 0.0;
+}
+
+
 // ************************************************************************* //

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,6 +38,7 @@ SourceFiles
 
 #include "ParticleForce.H"
 #include "volFields.H"
+#include "interpolation.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,13 +54,15 @@ class PressureGradientForce
 :
     public ParticleForce<CloudType>
 {
-    // Private data
+protected:
+
+    // Protected data
 
         //- Name of velocity field
         const word UName_;
 
-        //- Velocity gradient field
+        //- Rate of change of carrier phase velocity interpolator
-        const volTensorField* gradUPtr_;
+        autoPtr<interpolation<vector> > DUcDtInterpPtr_;
 
 
 public:
@@ -75,7 +78,8 @@ public:
         (
             CloudType& owner,
             const fvMesh& mesh,
-            const dictionary& dict
+            const dictionary& dict,
+            const word& forceType = typeName
         );
 
         //- Construct copy
@@ -99,8 +103,8 @@ public:
 
         // Access
 
-            //- Return const access to the velocity gradient field
+            //- Return the rate of change of carrier phase velocity interpolator
-            inline const volTensorField& gradU() const;
+            inline const interpolation<vector>& DUcDtInterp() const;
 
 
         // Evaluation
@@ -117,6 +121,13 @@ public:
                 const scalar Re,
                 const scalar muc
             ) const;
+
+            //- Return the added mass
+            virtual scalar massAdd
+            (
+                const typename CloudType::parcelType& p,
+                const scalar mass
+            ) const;
 };
 
 

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForceI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,23 +26,20 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class CloudType>
-const Foam::volTensorField& Foam::PressureGradientForce<CloudType>::gradU()
+inline const Foam::interpolation<Foam::vector>&
-const
+Foam::PressureGradientForce<CloudType>::DUcDtInterp() const
 {
-    if (gradUPtr_)
+    if (!DUcDtInterpPtr_.valid())
-    {
-        return *gradUPtr_;
-    }
-    else
     {
         FatalErrorIn
         (
-            "const volTensorField& PressureGradientForce<CloudType>::gradU()"
+            "inline const Foam::interpolation<Foam::vector>&"
-            "const"
+            "Foam::PressureGradientForce<CloudType>::DUcDtInterp() const"
-        )   << "gradU field not allocated" << abort(FatalError);
+        )   << "Carrier phase DUcDt interpolation object not set"
-
+            << abort(FatalError);
-        return *reinterpret_cast<const volTensorField*>(0);
     }
+
+    return DUcDtInterpPtr_();
 }
 
 

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C

@@ -24,8 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "VirtualMassForce.H"
-#include "fvcDdt.H"
-#include "fvcGrad.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -34,14 +32,12 @@ Foam::VirtualMassForce<CloudType>::VirtualMassForce
 (
     CloudType& owner,
     const fvMesh& mesh,
-    const dictionary& dict
+    const dictionary& dict,
+    const word& forceType
 )
 :
-    ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
+    PressureGradientForce<CloudType>(owner, mesh, dict, forceType),
-    UName_(this->coeffs().template lookupOrDefault<word>("U", "U")),
+    Cvm_(readScalar(this->coeffs().lookup("Cvm")))
-    Cvm_(readScalar(this->coeffs().lookup("Cvm"))),
-    DUcDtPtr_(NULL),
-    DUcDtInterpPtr_(NULL)
 {}
 
 
@@ -51,11 +47,8 @@ Foam::VirtualMassForce<CloudType>::VirtualMassForce
     const VirtualMassForce& vmf
 )
 :
-    ParticleForce<CloudType>(vmf),
+    PressureGradientForce<CloudType>(vmf),
-    UName_(vmf.UName_),
+    Cvm_(vmf.Cvm_)
-    Cvm_(vmf.Cvm_),
-    DUcDtPtr_(NULL),
-    DUcDtInterpPtr_(NULL)
 {}
 
 
@@ -71,36 +64,7 @@ Foam::VirtualMassForce<CloudType>::~VirtualMassForce()
 template<class CloudType>
 void Foam::VirtualMassForce<CloudType>::cacheFields(const bool store)
 {
-    if (store && !DUcDtPtr_)
+    PressureGradientForce<CloudType>::cacheFields(store);
-    {
-        const volVectorField& Uc = this->mesh().template
-            lookupObject<volVectorField>(UName_);
-
-        DUcDtPtr_ = new volVectorField
-        (
-            "DUcDt",
-            fvc::ddt(Uc) + (Uc & fvc::grad(Uc))
-        );
-
-        DUcDtInterpPtr_.reset
-        (
-            interpolation<vector>::New
-            (
-                this->owner().solution().interpolationSchemes(),
-                *DUcDtPtr_
-            ).ptr()
-        );
-    }
-    else
-    {
-        DUcDtInterpPtr_.clear();
-
-        if (DUcDtPtr_)
-        {
-            delete DUcDtPtr_;
-            DUcDtPtr_ = NULL;
-        }
-    }
 }
 
 
@@ -114,12 +78,10 @@ Foam::forceSuSp Foam::VirtualMassForce<CloudType>::calcCoupled
     const scalar muc
 ) const
 {
-    forceSuSp value(vector::zero, 0.0);
+    forceSuSp value =
+        PressureGradientForce<CloudType>::calcCoupled(p, dt, mass, Re, muc);
 
-    vector DUcDt =
+    value.Su() *= Cvm_;
-        DUcDtInterp().interpolate(p.position(), p.currentTetIndices());
-
-    value.Su() = mass*p.rhoc()/p.rho()*Cvm_*DUcDt;
 
     return value;
 }

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H

@@ -28,7 +28,6 @@ Description
     Calculates particle virtual mass force
 
 SourceFiles
-    VirtualMassForceI.H
     VirtualMassForce.C
 
 \*---------------------------------------------------------------------------*/
@@ -36,9 +35,7 @@ SourceFiles
 #ifndef VirtualMassForce_H
 #define VirtualMassForce_H
 
-#include "ParticleForce.H"
+#include "PressureGradientForce.H"
-#include "volFields.H"
-#include "interpolation.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,22 +49,13 @@ namespace Foam
 template<class CloudType>
 class VirtualMassForce
 :
-    public ParticleForce<CloudType>
+    public PressureGradientForce<CloudType>
 {
     // Private data
 
-        //- Name of velocity field
-        const word UName_;
-
         //- Virtual mass coefficient - typically 0.5
         scalar Cvm_;
 
-        //- Rate of change of carrier phase velocity
-        volVectorField* DUcDtPtr_;
-
-        //- Rate of change of carrier phase velocity interpolator
-        autoPtr<interpolation<vector> > DUcDtInterpPtr_;
-
 
 public:
 
@@ -82,7 +70,8 @@ public:
         (
             CloudType& owner,
             const fvMesh& mesh,
-            const dictionary& dict
+            const dictionary& dict,
+            const word& forceType = typeName
         );
 
         //- Construct copy
@@ -104,15 +93,6 @@ public:
 
     // Member Functions
 
-        // Access
-
-            //- Return the rate of change of carrier phase velocity
-            inline const volVectorField& DUcDt() const;
-
-            //- Return the rate of change of carrier phase velocity interpolator
-            inline const interpolation<vector>& DUcDtInterp() const;
-
-
         // Evaluation
 
             //- Cache fields
@@ -143,8 +123,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "VirtualMassForceI.H"
-
 #ifdef NoRepository
     #include "VirtualMassForce.C"
 #endif

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForceI.H

@@ -1,66 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-const Foam::volVectorField& Foam::VirtualMassForce<CloudType>::DUcDt() const
-{
-    if (DUcDtPtr_)
-    {
-        return *DUcDtPtr_;
-    }
-    else
-    {
-        FatalErrorIn
-        (
-            "const volVectorField& VirtualMassForce<CloudType>::DUcDt()"
-            "const"
-        )   << "DUcDt field not allocated" << abort(FatalError);
-
-        return *reinterpret_cast<const volVectorField*>(0);
-    }
-}
-
-
-template<class CloudType>
-inline const Foam::interpolation<Foam::vector>&
-Foam::VirtualMassForce<CloudType>::DUcDtInterp() const
-{
-    if (!DUcDtInterpPtr_.valid())
-    {
-        FatalErrorIn
-        (
-            "inline const Foam::interpolation<Foam::vector>&"
-            "Foam::VirtualMassForce<CloudType>::DUcDtInterp() const"
-        )   << "Carrier pahase DUcDt interpolation object not set"
-            << abort(FatalError);
-    }
-
-    return DUcDtInterpPtr_();
-}
-
-
-// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -48,6 +48,10 @@ Description
 #include "globalIndex.H"
 #include "DynamicField.H"
 #include "PatchTools.H"
+#include "slipPointPatchFields.H"
+#include "fixedValuePointPatchFields.H"
+#include "calculatedPointPatchFields.H"
+#include "cyclicSlipPointPatchFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -71,7 +75,7 @@ void Foam::autoLayerDriver::dumpDisplacement
 )
 {
     OFstream dispStr(prefix + "_disp.obj");
-    Info<< "Writing all displacements to " << dispStr.name() << nl << endl;
+    Info<< "Writing all displacements to " << dispStr.name() << endl;
 
     label vertI = 0;
 
@@ -87,7 +91,7 @@ void Foam::autoLayerDriver::dumpDisplacement
 
 
     OFstream illStr(prefix + "_illegal.obj");
-    Info<< "Writing invalid displacements to " << illStr.name() << nl << endl;
+    Info<< "Writing invalid displacements to " << illStr.name() << endl;
 
     vertI = 0;
 
@@ -801,6 +805,81 @@ void Foam::autoLayerDriver::setNumLayers
 }
 
 
+// Construct pointVectorField with correct boundary conditions for adding
+// layers
+Foam::tmp<Foam::pointVectorField>
+Foam::autoLayerDriver::makeLayerDisplacementField
+(
+    const pointMesh& pMesh,
+    const labelList& numLayers
+)
+{
+    // Construct displacement field.
+    const pointBoundaryMesh& pointPatches = pMesh.boundary();
+
+    wordList patchFieldTypes
+    (
+        pointPatches.size(),
+        slipPointPatchVectorField::typeName
+    );
+
+    forAll(numLayers, patchI)
+    {
+        //  0 layers: do not allow lslip so fixedValue 0
+        // >0 layers: fixedValue which gets adapted
+        if (numLayers[patchI] >= 0)
+        {
+            patchFieldTypes[patchI] = fixedValuePointPatchVectorField::typeName;
+        }
+    }
+
+    forAll(pointPatches, patchI)
+    {
+        if (isA<processorPointPatch>(pointPatches[patchI]))
+        {
+            patchFieldTypes[patchI] = calculatedPointPatchVectorField::typeName;
+        }
+        else if (isA<cyclicPointPatch>(pointPatches[patchI]))
+        {
+            patchFieldTypes[patchI] = cyclicSlipPointPatchVectorField::typeName;
+        }
+    }
+
+
+//Pout<< "*** makeLayerDisplacementField : boundary conditions:" << endl;
+//forAll(patchFieldTypes, patchI)
+//{
+//    Pout<< "\t" << patchI << " name:" << pointPatches[patchI].name()
+//        << " type:" << patchFieldTypes[patchI]
+//        << " nLayers:" << numLayers[patchI]
+//        << endl;
+//}
+
+    const polyMesh& mesh = pMesh();
+
+    // Note: time().timeName() instead of meshRefinement::timeName() since
+    // postprocessable field.
+    tmp<pointVectorField> tfld
+    (
+        new pointVectorField
+        (
+            IOobject
+            (
+                "pointDisplacement",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            pMesh,
+            dimensionedVector("displacement", dimLength, vector::zero),
+            patchFieldTypes
+        )
+    );
+    return tfld;
+}
+
+
 void Foam::autoLayerDriver::growNoExtrusion
 (
     const indirectPrimitivePatch& pp,
@@ -2392,10 +2471,10 @@ void Foam::autoLayerDriver::addLayers
             mesh,
             pp(),
             patchIDs,
-            meshRefinement::makeDisplacementField
+            makeLayerDisplacementField
             (
                 pointMesh::New(mesh),
-                patchIDs
+                layerParams.numLayers()
             ),
             motionDict
         )
@@ -3186,7 +3265,7 @@ void Foam::autoLayerDriver::doLayers
     // Merge coplanar boundary faces
     mergePatchFacesUndo(layerParams, motionDict);
 
-    // Per patch the number of layers (0 if no layer)
+    // Per patch the number of layers (-1 or 0 if no layer)
     const labelList& numLayers = layerParams.numLayers();
 
     // Patches that need to get a layer

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H

@@ -195,6 +195,19 @@ class autoLayerDriver
                     List<extrudeMode>& extrudeStatus
                 ) const;
 
+                //- Helper function to make a pointVectorField with correct
+                //  bcs for layer addition:
+                //  - numLayers > 0         : fixedValue
+                //  - numLayers == 0        : fixedValue (always zero)
+                //  - processor             : calculated (so free to move)
+                //  - cyclic/wedge/symmetry : slip
+                //  - other                 : slip
+                static tmp<pointVectorField> makeLayerDisplacementField
+                (
+                    const pointMesh& pMesh,
+                    const labelList& numLayers
+                );
+
                 //- Grow no-extrusion layer.
                 void growNoExtrusion
                 (
@@ -444,7 +457,6 @@ class autoLayerDriver
                     const PackedBoolList& isMasterEdge,
                     const labelList& meshEdges,
                     const scalar minCosLayerTermination,
-                    scalarField& field,
                     List<extrudeMode>& extrudeStatus,
                     pointField& patchDisp,
                     labelList& patchNLayers

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C

@@ -99,6 +99,76 @@ void Foam::autoLayerDriver::sumWeights
 
 
 // Smooth field on moving patch
+//void Foam::autoLayerDriver::smoothField
+//(
+//    const motionSmoother& meshMover,
+//    const PackedBoolList& isMasterEdge,
+//    const labelList& meshEdges,
+//    const scalarField& fieldMin,
+//    const label nSmoothDisp,
+//    scalarField& field
+//) const
+//{
+//    const indirectPrimitivePatch& pp = meshMover.patch();
+//    const edgeList& edges = pp.edges();
+//    const labelList& meshPoints = pp.meshPoints();
+//
+//    scalarField invSumWeight(pp.nPoints());
+//    sumWeights
+//    (
+//        isMasterEdge,
+//        meshEdges,
+//        meshPoints,
+//        edges,
+//        invSumWeight
+//    );
+//
+//    // Get smoothly varying patch field.
+//    Info<< "shrinkMeshDistance : Smoothing field ..." << endl;
+//
+//    for (label iter = 0; iter < nSmoothDisp; iter++)
+//    {
+//        scalarField average(pp.nPoints());
+//        averageNeighbours
+//        (
+//            meshMover.mesh(),
+//            isMasterEdge,
+//            meshEdges,
+//            meshPoints,
+//            pp.edges(),
+//            invSumWeight,
+//            field,
+//            average
+//        );
+//
+//        // Transfer to field
+//        forAll(field, pointI)
+//        {
+//            //full smoothing neighbours + point value
+//            average[pointI] = 0.5*(field[pointI]+average[pointI]);
+//
+//            // perform monotonic smoothing
+//            if
+//            (
+//                average[pointI] < field[pointI]
+//             && average[pointI] >= fieldMin[pointI]
+//            )
+//            {
+//                field[pointI] = average[pointI];
+//            }
+//        }
+//
+//        // Do residual calculation every so often.
+//        if ((iter % 10) == 0)
+//        {
+//            Info<< "    Iteration " << iter << "   residual "
+//                <<  gSum(mag(field-average))
+//                   /returnReduce(average.size(), sumOp<label>())
+//                << endl;
+//        }
+//    }
+//}
+//XXXXXXXXX
 void Foam::autoLayerDriver::smoothField
 (
     const motionSmoother& meshMover,
@@ -126,9 +196,15 @@ void Foam::autoLayerDriver::smoothField
     // Get smoothly varying patch field.
     Info<< "shrinkMeshDistance : Smoothing field ..." << endl;
 
-    for (label iter = 0; iter < nSmoothDisp; iter++)
+
+    const scalar lambda = 0.33;
+    const scalar mu = -0.34;
+
+    for (label iter = 0; iter < 90; iter++)
     {
         scalarField average(pp.nPoints());
+
+        // Calculate average of field
         averageNeighbours
         (
             meshMover.mesh(),
@@ -141,23 +217,37 @@ void Foam::autoLayerDriver::smoothField
             average
         );
 
-        // Transfer to field
+        forAll(field, i)
-        forAll(field, pointI)
         {
-            //full smoothing neighbours + point value
+            if (field[i] >= fieldMin[i])
-            average[pointI] = 0.5*(field[pointI]+average[pointI]);
-
-            // perform monotonic smoothing
-            if
-            (
-                average[pointI] < field[pointI]
-             && average[pointI] >= fieldMin[pointI]
-            )
             {
-                field[pointI] = average[pointI];
+                field[i] = (1-lambda)*field[i]+lambda*average[i];
             }
         }
 
+
+        // Calculate average of field
+        averageNeighbours
+        (
+            meshMover.mesh(),
+            isMasterEdge,
+            meshEdges,
+            meshPoints,
+            pp.edges(),
+            invSumWeight,
+            field,
+            average
+        );
+
+        forAll(field, i)
+        {
+            if (field[i] >= fieldMin[i])
+            {
+                field[i] = (1-mu)*field[i]+mu*average[i];
+            }
+        }
+
+
         // Do residual calculation every so often.
         if ((iter % 10) == 0)
         {
@@ -168,7 +258,7 @@ void Foam::autoLayerDriver::smoothField
         }
     }
 }
-
+//XXXXXXXXX
 
 // Smooth normals on moving patch.
 void Foam::autoLayerDriver::smoothPatchNormals
@@ -480,7 +570,6 @@ void Foam::autoLayerDriver::findIsolatedRegions
     const PackedBoolList& isMasterEdge,
     const labelList& meshEdges,
     const scalar minCosLayerTermination,
-    scalarField& field,
     List<extrudeMode>& extrudeStatus,
     pointField& patchDisp,
     labelList& patchNLayers
@@ -661,7 +750,6 @@ void Foam::autoLayerDriver::findIsolatedRegions
                         )
                     )
                     {
-                        field[f[fp]] = 0.0;
                         nPointCounter++;
                     }
                 }
@@ -670,7 +758,8 @@ void Foam::autoLayerDriver::findIsolatedRegions
     }
 
     reduce(nPointCounter, sumOp<label>());
-    Info<< "Number isolated points extrusion stopped : "<< nPointCounter<< endl;
+    Info<< "Number isolated points extrusion stopped : "<< nPointCounter
+        << endl;
 }
 
 
@@ -875,11 +964,14 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
         forAll(patches, patchI)
         {
             const polyPatch& pp = patches[patchI];
+            //const pointPatchVectorField& pvf =
+            //    meshMover.displacement().boundaryField()[patchI];
 
             if
             (
                 !pp.coupled()
              && !isA<emptyPolyPatch>(pp)
+             //&&  pvf.constraintType() != word::null //exclude fixedValue
              && !adaptPatches.found(patchI)
             )
             {
@@ -1168,12 +1260,22 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
         isMasterEdge,
         meshEdges,
         minCosLayerTermination,
-        thickness,
+
         extrudeStatus,
         patchDisp,
         patchNLayers
     );
 
+    // Update thickess for changed extrusion
+    forAll(thickness, patchPointI)
+    {
+        if (extrudeStatus[patchPointI] == NOEXTRUDE)
+        {
+            thickness[patchPointI] = 0.0;
+        }
+    }
+
+
     // smooth layer thickness on moving patch
     smoothField
     (
@@ -1182,6 +1284,7 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
         meshEdges,
         minThickness,
         nSmoothThickness,
+
         thickness
     );
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C

@@ -506,7 +506,7 @@ void Foam::autoSnapDriver::dumpMove
 )
 {
     // Dump direction of growth into file
-    Pout<< nl << "Dumping move direction to " << fName << endl;
+    Info<< "Dumping move direction to " << fName << endl;
 
     OFstream nearestStream(fName);
 
@@ -721,14 +721,14 @@ void Foam::autoSnapDriver::preSmoothPatch
     if (debug)
     {
         const_cast<Time&>(mesh.time())++;
-        Pout<< "Writing patch smoothed mesh to time "
+        Info<< "Writing patch smoothed mesh to time "
             << meshRefiner_.timeName() << '.' << endl;
         meshRefiner_.write
         (
             debug,
             mesh.time().path()/meshRefiner_.timeName()
         );
-        Pout<< "Dumped mesh in = "
+        Info<< "Dumped mesh in = "
             << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
     }
 
@@ -998,7 +998,7 @@ void Foam::autoSnapDriver::smoothDisplacement
     if (debug)
     {
         const_cast<Time&>(mesh.time())++;
-        Pout<< "Writing smoothed mesh to time " << meshRefiner_.timeName()
+        Info<< "Writing smoothed mesh to time " << meshRefiner_.timeName()
             << endl;
 
         // Moving mesh creates meshPhi. Can be cleared out by a mesh.clearOut
@@ -1010,13 +1010,12 @@ void Foam::autoSnapDriver::smoothDisplacement
             debug,
             mesh.time().path()/meshRefiner_.timeName()
         );
-
+        Info<< "Writing displacement field ..." << endl;
-        Pout<< "Writing displacement field ..." << endl;
         disp.write();
         tmp<pointScalarField> magDisp(mag(disp));
         magDisp().write();
 
-        Pout<< "Writing actual patch displacement ..." << endl;
+        Info<< "Writing actual patch displacement ..." << endl;
         vectorField actualPatchDisp(disp, pp.meshPoints());
         dumpMove
         (
@@ -1068,11 +1067,11 @@ bool Foam::autoSnapDriver::scaleMesh
         if (debug)
         {
             const_cast<Time&>(mesh.time())++;
-            Pout<< "Writing scaled mesh to time " << meshRefiner_.timeName()
+            Info<< "Writing scaled mesh to time " << meshRefiner_.timeName()
                 << endl;
             mesh.write();
 
-            Pout<< "Writing displacement field ..." << endl;
+            Info<< "Writing displacement field ..." << endl;
             meshMover.displacement().write();
             tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
             magDisp().write();
@@ -1521,14 +1520,14 @@ void Foam::autoSnapDriver::doSnap
             if (debug)
             {
                 const_cast<Time&>(mesh.time())++;
-                Pout<< "Writing scaled mesh to time "
+                Info<< "Writing scaled mesh to time "
                     << meshRefiner_.timeName() << endl;
                 meshRefiner_.write
                 (
                     debug,
                     mesh.time().path()/meshRefiner_.timeName()
                 );
-                Pout<< "Writing displacement field ..." << endl;
+                Info<< "Writing displacement field ..." << endl;
                 meshMover.displacement().write();
                 tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
                 magDisp().write();
@@ -1564,7 +1563,7 @@ void Foam::autoSnapDriver::doSnap
     if (nChanged > 0 && debug)
     {
         const_cast<Time&>(mesh.time())++;
-        Pout<< "Writing patchFace merged mesh to time "
+        Info<< "Writing patchFace merged mesh to time "
             << meshRefiner_.timeName() << endl;
         meshRefiner_.write
         (

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C

@@ -2565,6 +2565,76 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         }
     }
 
+
+    // Collect additionally 'normal' boundary faces for boundaryPoints of pp
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // points on the boundary of pp should pick up non-pp normals
+    // as well for the feature-reconstruction to behave correctly.
+    // (the movement is already constrained outside correctly so it
+    //  is only that the unconstrained attraction vector is calculated
+    //  correctly)
+    {
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        labelList patchID(pbm.patchID());
+
+        // Unmark all non-coupled boundary faces
+        forAll(pbm, patchI)
+        {
+            const polyPatch& pp = pbm[patchI];
+
+            if (pp.coupled() || isA<emptyPolyPatch>(pp))
+            {
+                forAll(pp, i)
+                {
+                    label meshFaceI = pp.start()+i;
+                    patchID[meshFaceI-mesh.nInternalFaces()] = -1;
+                }
+            }
+        }
+
+        // Remove any meshed faces
+        PackedBoolList ppFaces(mesh.nFaces());
+        forAll(pp.addressing(), i)
+        {
+            label meshFaceI = pp.addressing()[i];
+            patchID[meshFaceI-mesh.nInternalFaces()] = -1;
+        }
+
+        // See if pp point uses any non-meshed boundary faces
+
+        const labelList& boundaryPoints = pp.boundaryPoints();
+        forAll(boundaryPoints, i)
+        {
+            label pointI = boundaryPoints[i];
+            label meshPointI = pp.meshPoints()[pointI];
+            const point& pt = mesh.points()[meshPointI];
+            const labelList& pFaces = mesh.pointFaces()[meshPointI];
+
+            List<point>& pNormals = pointFaceSurfNormals[pointI];
+            List<point>& pDisp = pointFaceDisp[pointI];
+            List<point>& pFc = pointFaceCentres[pointI];
+            labelList& pFid = pointFacePatchID[pointI];
+
+            forAll(pFaces, i)
+            {
+                label meshFaceI = pFaces[i];
+                if (!mesh.isInternalFace(meshFaceI))
+                {
+                    label patchI = patchID[meshFaceI-mesh.nInternalFaces()];
+
+                    if (patchI != -1)
+                    {
+                        vector fn = mesh.faceAreas()[meshFaceI];
+                        pNormals.append(fn/mag(fn));
+                        pDisp.append(mesh.faceCentres()[meshFaceI]-pt);
+                        pFc.append(mesh.faceCentres()[meshFaceI]);
+                        pFid.append(patchI);
+                    }
+                }
+            }
+        }
+    }
+
     syncTools::syncPointList
     (
         mesh,

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C

@@ -43,7 +43,7 @@ Foam::layerParameters::layerParameters
     const polyBoundaryMesh& boundaryMesh
 )
 :
-    numLayers_(boundaryMesh.size(), 0),
+    numLayers_(boundaryMesh.size(), -1),
     expansionRatio_
     (
         boundaryMesh.size(),

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H

@@ -131,6 +131,10 @@ public:
             // Per patch information
 
                 //- How many layers to add.
+                //  -1 : no specification. Assume 0 layers but allow sliding
+                //       to make layers
+                //   0 : specified to have 0 layers. No sliding allowed.
+                //  >0 : number of layers
                 const labelList& numLayers() const
                 {
                     return numLayers_;

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C

@@ -571,7 +571,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
             hitInfo
         );
 
-        // Start of from current points
+        // Start off from current points
         newPoints = mesh_.points();
 
         forAll(hitInfo, i)
@@ -581,6 +581,17 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
                 newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
             }
         }
+
+        if (debug)
+        {
+            const_cast<Time&>(mesh_.time())++;
+            pointField oldPoints(mesh_.points());
+            mesh_.movePoints(newPoints);
+            Pout<< "Writing newPoints mesh to time " << timeName()
+                << endl;
+            write(debug, mesh_.time().path()/"newPoints");
+            mesh_.movePoints(oldPoints);
+        }
     }
 
 

src/turbulenceModels/compressible/LES/GenEddyVisc/GenEddyVisc.H

@@ -51,7 +51,7 @@ namespace LESModels
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class GenEddyVisc Declaration
+                         Class GenEddyVisc Declaration
 \*---------------------------------------------------------------------------*/
 
 class GenEddyVisc
@@ -108,7 +108,21 @@ public:
         //- Return sub-grid disipation rate
         virtual tmp<volScalarField> epsilon() const
         {
-            return ce_*k()*sqrt(k())/delta();
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (   
+                        "epsilon",
+                        runTime_.timeName(),
+                        mesh_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    ce_*k()*sqrt(k())/delta()
+                )
+            );
         }
 
         //- Return viscosity

src/turbulenceModels/compressible/LES/GenSGSStress/GenSGSStress.H

@@ -52,7 +52,7 @@ namespace LESModels
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class GenSGSStress Declaration
+                        Class GenSGSStress Declaration
 \*---------------------------------------------------------------------------*/
 
 class GenSGSStress
@@ -109,14 +109,43 @@ public:
         //- Return the SGS turbulent kinetic energy
         virtual tmp<volScalarField> k() const
         {
-            return 0.5*tr(B_);
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (   
+                        "k",
+                        runTime_.timeName(),
+                        mesh_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    0.5*tr(B_)
+                )
+            );
         }
 
         //- Return the SGS turbulent dissipation
         virtual tmp<volScalarField> epsilon() const
         {
             const volScalarField K(k());
-            return ce_*K*sqrt(K)/delta();
+
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (   
+                        "epsilon",
+                        runTime_.timeName(),
+                        mesh_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    ce_*K*sqrt(K)/delta()
+                )
+            );
         }
 
         //- Return the SGS viscosity

src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.C

@@ -79,7 +79,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField
     fixedGradientFvPatchScalarField(p, iF),
     temperatureCoupledBase(patch(), "undefined", "undefined-K"),
     heatSource_(hsPower),
-    q_(p.size(), 0.0)
+    q_(p.size(), 0.0),
+    QrName_("undefinedQr")
 {}
 
 
@@ -95,7 +96,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField
     fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
     temperatureCoupledBase(patch(), ptf.KMethod(), ptf.kappaName()),
     heatSource_(ptf.heatSource_),
-    q_(ptf.q_, mapper)
+    q_(ptf.q_, mapper),
+    QrName_(ptf.QrName_)
 {}
 
 
@@ -110,7 +112,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField
     fixedGradientFvPatchScalarField(p, iF),
     temperatureCoupledBase(patch(), dict),
     heatSource_(heatSourceTypeNames_.read(dict.lookup("heatSource"))),
-    q_("q", dict, p.size())
+    q_("q", dict, p.size()),
+    QrName_(dict.lookupOrDefault<word>("Qr", "none"))
 {
     fvPatchField<scalar>::operator=(patchInternalField());
     gradient() = 0.0;
@@ -126,7 +129,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField
     fixedGradientFvPatchScalarField(thftpsf),
     temperatureCoupledBase(patch(), thftpsf.KMethod(), thftpsf.kappaName()),
     heatSource_(thftpsf.heatSource_),
-    q_(thftpsf.q_)
+    q_(thftpsf.q_),
+    QrName_(thftpsf.QrName_)
 {}
 
 
@@ -140,7 +144,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField
     fixedGradientFvPatchScalarField(thftpsf, iF),
     temperatureCoupledBase(patch(), thftpsf.KMethod(), thftpsf.kappaName()),
     heatSource_(thftpsf.heatSource_),
-    q_(thftpsf.q_)
+    q_(thftpsf.q_),
+    QrName_(thftpsf.QrName_)
 {}
 
 
@@ -183,17 +188,25 @@ void turbulentHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
 
     const scalarField& Tp = *this;
 
+    scalarField qr(this->size(), 0.0);
+
+    //- qr is negative going into the domain
+    if (QrName_ != "none")
+    {
+        qr = patch().lookupPatchField<volScalarField, scalar>(QrName_);
+    }
+
     switch (heatSource_)
     {
         case hsPower:
         {
             const scalar Ap = gSum(patch().magSf());
-            gradient() = q_/(Ap*kappa(Tp));
+            gradient() = (q_/Ap + qr)/kappa(Tp);
             break;
         }
         case hsFlux:
         {
-            gradient() = q_/kappa(Tp);
+            gradient() = (q_ + qr)/kappa(Tp);
             break;
         }
         default:
@@ -220,6 +233,7 @@ void turbulentHeatFluxTemperatureFvPatchScalarField::write
         << token::END_STATEMENT << nl;
     temperatureCoupledBase::write(os);
     q_.writeEntry("q", os);
+    os.writeKeyword("Qr")<< QrName_ << token::END_STATEMENT << nl;
     writeEntry("value", os);
 }
 

src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.H

@@ -37,6 +37,7 @@ Description
             heatSource      flux;        // power [W]; flux [W/m2]
             q               uniform 10;  // heat power or flux
             kappa           fluidThermo; // calculate kappa=alphaEff*thermo.Cp
+            Qr              none;        // name of the radiative flux
             value           uniform 300; // initial temperature value
         }
 
@@ -93,6 +94,9 @@ private:
         //- Heat power [W] or flux [W/m2]
         scalarField q_;
 
+        //- Name of radiative in flux field
+        word QrName_;
+
 
 public:
 

src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H

@@ -51,7 +51,7 @@ namespace LESModels
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class GenEddyVisc Declaration
+                         Class GenEddyVisc Declaration
 \*---------------------------------------------------------------------------*/
 
 class GenEddyVisc
@@ -104,7 +104,21 @@ public:
         //- Return sub-grid disipation rate
         virtual tmp<volScalarField> epsilon() const
         {
-            return ce_*k()*sqrt(k())/delta();
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (   
+                        "epsilon",
+                        runTime_.timeName(),
+                        mesh_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    ce_*k()*sqrt(k())/delta()
+                )
+            );
         }
 
         //- Return the SGS viscosity.

src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H

@@ -52,7 +52,7 @@ namespace LESModels
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class GenSGSStress Declaration
+                        Class GenSGSStress Declaration
 \*---------------------------------------------------------------------------*/
 
 class GenSGSStress
@@ -104,14 +104,43 @@ public:
         //- Return the SGS turbulent kinetic energy.
         virtual tmp<volScalarField> k() const
         {
-            return 0.5*tr(B_);
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (   
+                        "k",
+                        runTime_.timeName(),
+                        mesh_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    0.5*tr(B_)
+                )
+            );
         }
 
         //- Return the SGS turbulent dissipation.
         virtual tmp<volScalarField> epsilon() const
         {
             const volScalarField K(k());
-            return ce_*K*sqrt(K)/delta();
+            
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (   
+                        "epsilon",
+                        runTime_.timeName(),
+                        mesh_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    ce_*K*sqrt(K)/delta()
+                )
+            );
         }
 
         //- Return the SGS viscosity.

tutorials/mesh/cvMesh/flange/system/cvMeshDict

@@ -217,7 +217,9 @@ polyMeshFiltering
 }
 
 
-#include "meshQualityControls"
+meshQualityControls
-
+{
+    #include "meshQualityDict"
+}
 
 // ************************************************************************* //

tutorials/mesh/cvMesh/flange/system/meshQualityControls

@@ -1,76 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-
-FoamFile
-{
-    version         2.0;
-    format          ascii;
-
-    root            "";
-    case            "";
-    instance        "";
-    local           "";
-
-    class           dictionary;
-    object          meshQualityControls;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-meshQualityControls
-{
-    //- Maximum non-orthogonality allowed. Set to 180 to disable.
-    maxNonOrtho         65;
-
-    //- Max skewness allowed. Set to <0 to disable.
-    maxBoundarySkewness 50;
-    maxInternalSkewness 10;
-
-    //- Max concaveness allowed. Is angle (in degrees) below which concavity
-    //  is allowed. 0 is straight face, <0 would be convex face.
-    //  Set to 180 to disable.
-    maxConcave          80;
-
-    //- Minimum quality of the tet formed by the face-centre
-    //  and variable base point minimum decomposition triangles and
-    //  the cell centre. This has to be a positive number for tracking
-    //  to work. Set to very negative number (e.g. -1E30) to
-    //  disable.
-    //     <0 = inside out tet,
-    //      0 = flat tet
-    //      1 = regular tet
-    minTetQuality       1e-30;
-
-    //- Minimum pyramid volume. Is absolute volume of cell pyramid.
-    //  Set to a sensible fraction of the smallest cell volume expected.
-    //  Set to very negative number (e.g. -1E30) to disable.
-    minVol              0;
-
-    //- Minimum face area. Set to <0 to disable.
-    minArea             -1;
-
-    //- Minimum face twist. Set to <-1 to disable. dot product of face normal
-    //- and face centre triangles normal
-    minTwist            0.001;
-
-    //- minimum normalised cell determinant
-    //- 1 = hex, <= 0 = folded or flattened illegal cell
-    minDeterminant      0.001;
-
-    //- minFaceWeight (0 -> 0.5)
-    minFaceWeight       0.02;
-
-    //- minVolRatio (0 -> 1)
-    minVolRatio         0.01;
-
-    //must be >0 for Fluent compatibility
-    minTriangleTwist    -1;
-}
-
-
-// ************************************************************************* //

tutorials/mesh/cvMesh/flange/system/meshQualityDict

@@ -0,0 +1,73 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          meshQualityDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Maximum non-orthogonality allowed. Set to 180 to disable.
+maxNonOrtho         65;
+
+//- Max skewness allowed. Set to <0 to disable.
+maxBoundarySkewness 50;
+maxInternalSkewness 10;
+
+//- Max concaveness allowed. Is angle (in degrees) below which concavity
+//  is allowed. 0 is straight face, <0 would be convex face.
+//  Set to 180 to disable.
+maxConcave          80;
+
+//- Minimum quality of the tet formed by the face-centre
+//  and variable base point minimum decomposition triangles and
+//  the cell centre. This has to be a positive number for tracking
+//  to work. Set to very negative number (e.g. -1E30) to
+//  disable.
+//     <0 = inside out tet,
+//      0 = flat tet
+//      1 = regular tet
+minTetQuality       1e-30;
+
+//- Minimum pyramid volume. Is absolute volume of cell pyramid.
+//  Set to a sensible fraction of the smallest cell volume expected.
+//  Set to very negative number (e.g. -1E30) to disable.
+minVol              0;
+
+//- Minimum face area. Set to <0 to disable.
+minArea             -1;
+
+//- Minimum face twist. Set to <-1 to disable. dot product of face normal
+//- and face centre triangles normal
+minTwist            0.001;
+
+//- minimum normalised cell determinant
+//- 1 = hex, <= 0 = folded or flattened illegal cell
+minDeterminant      0.001;
+
+//- minFaceWeight (0 -> 0.5)
+minFaceWeight       0.02;
+
+//- minVolRatio (0 -> 1)
+minVolRatio         0.01;
+
+//must be >0 for Fluent compatibility
+minTriangleTwist    -1;
+
+
+// ************************************************************************* //

tutorials/mesh/cvMesh/flange/system/snappyHexMeshDict

@@ -278,12 +278,12 @@ addLayersControls
 
 
 
-// Generic mesh quality settings. At any undoable phase these determine
-// where to undo.
-#include "meshQualityControls"
-
 meshQualityControls
 {
+    // Generic mesh quality settings. At any undoable phase these determine
+    // where to undo.
+    #include "meshQualityDict"
+
     //- Number of error distribution iterations
     nSmoothScale 4;
     //- amount to scale back displacement at error points