Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelFoam/createFields.H

@@ -1,4 +1,4 @@
-    Info<< "Reading transportProperties\n" << endl;
+    Info<< "\nReading transportProperties\n" << endl;
 
     IOdictionary transportProperties
     (
@@ -31,7 +31,7 @@
         rhoInfValue
     );
 
-    Info<< "\nReading field U\n" << endl;
+    Info<< "Reading field U\n" << endl;
     volVectorField U
     (
         IOobject
@@ -127,7 +127,7 @@
 
     if (HdotGradHheader.headerOk())
     {
-        Info<< "\nReading field HdotGradH\n" << endl;
+        Info<< "Reading field HdotGradH" << endl;
 
         HdotGradHPtr_.reset
         (

src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H

@@ -110,13 +110,11 @@ class timeVaryingMappedFixedValueFvPatchField
             label& hi
         ) const;
 
+
         //- Read boundary points and determine interpolation weights to patch
         //  faceCentres
         void readSamplePoints();
 
-        //- Find boundary data inbetween current time and interpolate
-        void checkTable();
-
         //- Do actual interpolation using current weights
         tmp<Field<Type> > interpolate(const Field<Type>&) const;
 
@@ -199,6 +197,12 @@ public:
                 return referenceCS_;
             }
 
+            //- Return startSampledValues
+            const Field<Type> startSampledValues()
+            {
+                 return startSampledValues_;
+            }
+
 
         // Mapping functions
 
@@ -216,6 +220,12 @@ public:
             );
 
 
+        // Utility functions
+
+            //- Find boundary data inbetween current time and interpolate
+            void checkTable();
+
+
         // Evaluation functions
 
             //- Update the coefficients associated with the patch field

src/lagrangian/basic/InteractionLists/InteractionLists.H

@@ -287,10 +287,13 @@ public:
             wallFaceIndexAndTransformToDistribute() const;
 
             //- Return access to the referred wall faces
-            const List<referredWallFace>& referredWallFaces() const;
+            inline const List<referredWallFace>& referredWallFaces() const;
+
+            //- Return the name of the velocity field
+            inline const word& UName() const;
 
             //- Return access to the referred wall data
-            const List<vector>& referredWallData() const;
+            inline const List<vector>& referredWallData() const;
 
             //- Return access to the referred particle container
             inline const List<IDLList<ParticleType> >&

src/lagrangian/basic/InteractionLists/InteractionListsI.H

@@ -128,6 +128,13 @@ Foam::InteractionLists<ParticleType>::referredWallFaces() const
 }
 
 
+template<class ParticleType>
+const Foam::word& Foam::InteractionLists<ParticleType>::UName() const
+{
+    return UName_;
+}
+
+
 template<class ParticleType>
 const Foam::List<Foam::vector>&
 Foam::InteractionLists<ParticleType>::referredWallData() const

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

@@ -42,7 +42,7 @@ template<class ParcelType>
 void Foam::KinematicCloud<ParcelType>::preEvolve()
 {
     this->dispersion().cacheFields(true);
-    forces_.cacheFields(true);
+    forces_.cacheFields(true, interpolationSchemes_);
 }
 
 
@@ -152,7 +152,7 @@ void Foam::KinematicCloud<ParcelType>::postEvolve()
     }
 
     this->dispersion().cacheFields(false);
-    forces_.cacheFields(false);
+    forces_.cacheFields(false, interpolationSchemes_);
 
     this->postProcessing().post();
 }

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

@@ -162,6 +162,7 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     // Momentum source due to particle forces
     const vector FCoupled = mass*td.cloud().forces().calcCoupled
     (
+        this->position(),
         cellI,
         dt,
         rhoc_,
@@ -173,6 +174,7 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 
     const vector FNonCoupled = mass*td.cloud().forces().calcNonCoupled
     (
+        this->position(),
         cellI,
         dt,
         rhoc_,

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H

@@ -182,8 +182,8 @@ public:
 
         // Constructors
 
-           //- Construct from components
+            //- Construct from components
-           inline trackData
+            inline trackData
             (
                 KinematicCloud<ParcelType>& cloud,
                 const constantProperties& constProps,

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H

@@ -439,14 +439,23 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::wallImpactDistance
     const vector&
 ) const
 {
-    // Do not use a wall impact distance to allow proper multiple face
+    const KinematicCloud<ParcelType>& c =
-    // collisions.  In a twisted mesh the particle can be within range
+        dynamic_cast<const KinematicCloud<ParcelType>&>(this->cloud());
-    // of a wall but not in the cell attached to a wall face, hence
-    // miss the interaction.
 
-    return 0.0;
+    if (c.collision().controlsWallInteraction())
+    {
+        // Do not use a wall impact distance if the collision model
+        // controls wall interactions to allow proper multiple face
+        // collisions.  In a twisted mesh the particle can be within
+        // range of a wall but not in the cell attached to a wall
+        // face, hence miss the interaction.
 
-    // return 0.5*d_;
+        return 0.0;
+    }
+    else
+    {
+        return 0.5*d_;
+    }
 }
 
 

src/lagrangian/intermediate/particleForces/particleForces.C

@@ -30,6 +30,24 @@ License
 #include "mathematicalConstants.H"
 #include "electromagneticConstants.H"
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::particleForces::deleteFields()
+{
+    if (gradUPtr_)
+    {
+        delete gradUPtr_;
+        gradUPtr_ = NULL;
+    }
+
+    if (HdotGradHInterPtr_)
+    {
+        delete HdotGradHInterPtr_;
+        HdotGradHInterPtr_ = NULL;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::particleForces::particleForces
@@ -85,7 +103,7 @@ Foam::particleForces::particleForces(const particleForces& f)
 
 Foam::particleForces::~particleForces()
 {
-    cacheFields(false);
+    deleteFields();
 }
 
 
@@ -151,26 +169,46 @@ const Foam::word& Foam::particleForces::HdotGradHName() const
 }
 
 
-void Foam::particleForces::cacheFields(const bool store)
+void Foam::particleForces::cacheFields
+(
+    const bool store,
+    const dictionary& interpolationSchemes
+)
 {
-    if (store && pressureGradient_)
+    if (store)
     {
-        const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
+        if (pressureGradient_)
-        gradUPtr_ = fvc::grad(U).ptr();
+        {
+            const volVectorField& U =
+                mesh_.lookupObject<volVectorField>(UName_);
+
+            gradUPtr_ = fvc::grad(U).ptr();
+        }
+
+        if (paramagnetic_)
+        {
+            const volVectorField& HdotGradH = mesh_.lookupObject<volVectorField>
+            (
+                HdotGradHName_
+            );
+
+            HdotGradHInterPtr_ = interpolation<vector>::New
+            (
+                interpolationSchemes,
+                HdotGradH
+            ).ptr();
+        }
     }
     else
     {
-        if (gradUPtr_)
+        deleteFields();
-        {
-            delete gradUPtr_;
-            gradUPtr_ = NULL;
-        }
     }
 }
 
 
 Foam::vector Foam::particleForces::calcCoupled
 (
+    const vector& position,
     const label cellI,
     const scalar dt,
     const scalar rhoc,
@@ -205,6 +243,7 @@ Foam::vector Foam::particleForces::calcCoupled
 
 Foam::vector Foam::particleForces::calcNonCoupled
 (
+    const vector& position,
     const label cellI,
     const scalar dt,
     const scalar rhoc,
@@ -226,15 +265,12 @@ Foam::vector Foam::particleForces::calcNonCoupled
 
     if (paramagnetic_)
     {
-        const volVectorField& HdotGradH = mesh_.lookupObject<volVectorField>
+        const interpolation<vector>& HdotGradHInter = *HdotGradHInterPtr_;
-        (
-            HdotGradHName_
-        );
 
         accelTot +=
             3.0*constant::electromagnetic::mu0.value()/rho
            *magneticSusceptibility_/(magneticSusceptibility_ + 3)
-           *HdotGradH[cellI];
+           *HdotGradHInter.interpolate(position, cellI);
 
         // force is:
 

src/lagrangian/intermediate/particleForces/particleForces.H

@@ -40,6 +40,7 @@ SourceFiles
 #include "Switch.H"
 #include "vector.H"
 #include "volFieldsFwd.H"
+#include "interpolation.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -69,6 +70,9 @@ class particleForces
         //- Velocity gradient field
         const volTensorField* gradUPtr_;
 
+        //- HdotGradH interpolator
+        const interpolation<vector>* HdotGradHInterPtr_;
+
 
         // Forces to include in particle motion evaluation
 
@@ -101,6 +105,12 @@ class particleForces
             const word HdotGradHName_;
 
 
+    // Private Member Functions
+
+        //- Delete cached carrier fields
+        void deleteFields();
+
+
 public:
 
     // Constructors
@@ -159,11 +169,16 @@ public:
        // Evaluation
 
             //- Cache carrier fields
-            void cacheFields(const bool store);
+            void cacheFields
+            (
+                const bool store,
+                const dictionary& interpolationSchemes
+            );
 
             //- Calculate action/reaction forces between carrier and particles
             vector calcCoupled
             (
+                const vector& position,
                 const label cellI,
                 const scalar dt,
                 const scalar rhoc,
@@ -176,6 +191,7 @@ public:
             //- Calculate external forces applied to the particles
             vector calcNonCoupled
             (
+                const vector& position,
                 const label cellI,
                 const scalar dt,
                 const scalar rhoc,

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

@@ -141,6 +141,10 @@ public:
         //- Flag to indicate whether model activates injection model
         virtual bool active() const = 0;
 
+        //- Indicates whether model determines wall collisions or not,
+        //  used to determine what value to use for wallImpactDistance
+        virtual bool controlsWallInteraction() const = 0;
+
         // Collision function
         virtual void collide() = 0;
 };

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/NoCollision/NoCollision.C

@@ -61,6 +61,13 @@ bool Foam::NoCollision<CloudType>::active() const
 }
 
 
+template<class CloudType>
+bool Foam::NoCollision<CloudType>::controlsWallInteraction() const
+{
+    return false;
+}
+
+
 template<class CloudType>
 void Foam::NoCollision<CloudType>::collide()
 {}

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/NoCollision/NoCollision.H

@@ -82,6 +82,10 @@ public:
         //- Flag to indicate whether model activates injection model
         virtual bool active() const;
 
+        //- Indicates whether model determines wall collisions or not,
+        //  used to determine what value to use for wallImpactDistance
+        virtual bool controlsWallInteraction() const;
+
         // Collision function
         virtual void collide();
 };

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C

@@ -175,13 +175,20 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
     const labelListList directWallFaces = il_.dwfil();
 
+    const labelList& patchID = mesh.boundaryMesh().patchID();
+
+    const volVectorField& U = mesh.lookupObject<volVectorField>(il_.UName());
+
     // Storage for the wall interaction sites
-    DynamicList<point> flatSites;
+    DynamicList<point> flatSitePoints;
     DynamicList<scalar> flatSiteExclusionDistancesSqr;
-    DynamicList<point> otherSites;
+    DynamicList<WallSiteData<vector> > flatSiteData;
+    DynamicList<point> otherSitePoints;
     DynamicList<scalar> otherSiteDistances;
-    DynamicList<point> sharpSites;
+    DynamicList<WallSiteData<vector> > otherSiteData;
+    DynamicList<point> sharpSitePoints;
     DynamicList<scalar> sharpSiteExclusionDistancesSqr;
+    DynamicList<WallSiteData<vector> > sharpSiteData;
 
     forAll(dil, realCellI)
     {
@@ -191,19 +198,22 @@ void Foam::PairCollision<CloudType>::wallInteraction()
         // Loop over all Parcels in cell
         forAll(cellOccupancy_[realCellI], cellParticleI)
         {
-            flatSites.clear();
+            flatSitePoints.clear();
             flatSiteExclusionDistancesSqr.clear();
-            otherSites.clear();
+            flatSiteData.clear();
+            otherSitePoints.clear();
             otherSiteDistances.clear();
-            sharpSites.clear();
+            otherSiteData.clear();
+            sharpSitePoints.clear();
             sharpSiteExclusionDistancesSqr.clear();
+            sharpSiteData.clear();
 
             typename CloudType::parcelType& p =
-            *cellOccupancy_[realCellI][cellParticleI];
+                *cellOccupancy_[realCellI][cellParticleI];
 
             const point& pos = p.position();
 
-            scalar r = p.d()/2;
+            scalar r = wallModel_->pREff(p);
 
             // real wallFace interactions
 
@@ -229,6 +239,18 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
                     scalar normalAlignment = normal & pW/mag(pW);
 
+                    // Find the patchIndex and wallData for WallSiteData object
+                    label patchI = patchID[realFaceI - mesh.nInternalFaces()];
+
+                    label patchFaceI =
+                        realFaceI - mesh.boundaryMesh()[patchI].start();
+
+                    WallSiteData<vector> wSD
+                    (
+                        patchI,
+                        U.boundaryField()[patchI][patchFaceI]
+                    );
+
                     if (normalAlignment > cosPhiMinFlatWall)
                     {
                         // Guard against a flat interaction being
@@ -239,25 +261,29 @@ void Foam::PairCollision<CloudType>::wallInteraction()
                         (
                             !duplicatePointInList
                             (
-                                flatSites,
+                                flatSitePoints,
                                 nearPt,
                                 sqr(r*flatWallDuplicateExclusion)
                             )
                         )
                         {
-                            flatSites.append(nearPt);
+                            flatSitePoints.append(nearPt);
 
                             flatSiteExclusionDistancesSqr.append
                             (
                                 sqr(r) - sqr(nearest.distance())
                             );
+
+                            flatSiteData.append(wSD);
                         }
                     }
                     else
                     {
-                        otherSites.append(nearPt);
+                        otherSitePoints.append(nearPt);
 
                         otherSiteDistances.append(nearest.distance());
+
+                        otherSiteData.append(wSD);
                     }
                 }
             }
@@ -269,8 +295,10 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
             forAll(cellRefWallFaces, rWFI)
             {
+                label refWallFaceI = cellRefWallFaces[rWFI];
+
                 const referredWallFace& rwf =
-                il_.referredWallFaces()[cellRefWallFaces[rWFI]];
+                    il_.referredWallFaces()[refWallFaceI];
 
                 const pointField& pts = rwf.points();
 
@@ -288,6 +316,14 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
                     scalar normalAlignment = normal & pW/mag(pW);
 
+                    // Find the patchIndex and wallData for WallSiteData object
+
+                    WallSiteData<vector> wSD
+                    (
+                        rwf.patchIndex(),
+                        il_.referredWallData()[refWallFaceI]
+                    );
+
                     if (normalAlignment > cosPhiMinFlatWall)
                     {
                         // Guard against a flat interaction being
@@ -298,25 +334,29 @@ void Foam::PairCollision<CloudType>::wallInteraction()
                         (
                             !duplicatePointInList
                             (
-                                flatSites,
+                                flatSitePoints,
                                 nearPt,
                                 sqr(r*flatWallDuplicateExclusion)
                             )
                         )
                         {
-                            flatSites.append(nearPt);
+                            flatSitePoints.append(nearPt);
 
                             flatSiteExclusionDistancesSqr.append
                             (
                                 sqr(r) - sqr(nearest.distance())
                             );
+
+                            flatSiteData.append(wSD);
                         }
                     }
                     else
                     {
-                        otherSites.append(nearPt);
+                        otherSitePoints.append(nearPt);
 
                         otherSiteDistances.append(nearest.distance());
+
+                        otherSiteData.append(wSD);
                     }
                 }
             }
@@ -338,13 +378,13 @@ void Foam::PairCollision<CloudType>::wallInteraction()
             {
                 label orderedIndex = sortedOtherSiteIndices[siteI];
 
-                const point& otherPt = otherSites[orderedIndex];
+                const point& otherPt = otherSitePoints[orderedIndex];
 
                 if
                 (
                     !duplicatePointInList
                     (
-                        flatSites,
+                        flatSitePoints,
                         otherPt,
                         flatSiteExclusionDistancesSqr
                     )
@@ -357,23 +397,32 @@ void Foam::PairCollision<CloudType>::wallInteraction()
                     (
                         !duplicatePointInList
                         (
-                            sharpSites,
+                            sharpSitePoints,
                             otherPt,
                             sharpSiteExclusionDistancesSqr
                         )
                     )
                     {
-                        sharpSites.append(otherPt);
+                        sharpSitePoints.append(otherPt);
 
                         sharpSiteExclusionDistancesSqr.append
                         (
                             sqr(r) - sqr(otherSiteDistances[orderedIndex])
                         );
+
+                        sharpSiteData.append(otherSiteData[orderedIndex]);
                     }
                 }
             }
 
-            evaluateWall(p, flatSites, sharpSites);
+            evaluateWall
+            (
+                p,
+                flatSitePoints,
+                flatSiteData,
+                sharpSitePoints,
+                sharpSiteData
+            );
         }
     }
 }
@@ -463,11 +512,20 @@ template<class CloudType>
 void Foam::PairCollision<CloudType>::evaluateWall
 (
     typename CloudType::parcelType& p,
-    const List<point>& flatSites,
+    const List<point>& flatSitePoints,
-    const List<point>& sharpSites
+    const List<WallSiteData<vector> >& flatSiteData,
+    const List<point>& sharpSitePoints,
+    const List<WallSiteData<vector> >& sharpSiteData
 ) const
 {
-    wallModel_->evaluateWall(p, flatSites, sharpSites);
+    wallModel_->evaluateWall
+    (
+        p,
+        flatSitePoints,
+        flatSiteData,
+        sharpSitePoints,
+        sharpSiteData
+    );
 }
 
 
@@ -502,7 +560,8 @@ Foam::PairCollision<CloudType>::PairCollision
     (
         owner.mesh(),
         readScalar(this->coeffDict().lookup("maxInteractionDistance")),
-        Switch(this->coeffDict().lookup("writeReferredParticleCloud"))
+        Switch(this->coeffDict().lookup("writeReferredParticleCloud")),
+        this->coeffDict().lookupOrDefault("UName", word("U"))
     )
 {}
 
@@ -552,6 +611,13 @@ bool Foam::PairCollision<CloudType>::active() const
 }
 
 
+template<class CloudType>
+bool Foam::PairCollision<CloudType>::controlsWallInteraction() const
+{
+    return true;
+}
+
+
 template<class CloudType>
 void Foam::PairCollision<CloudType>::collide()
 {

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.H

@@ -36,6 +36,7 @@ SourceFiles
 
 #include "CollisionModel.H"
 #include "InteractionLists.H"
+#include "WallSiteData.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -137,8 +138,10 @@ class PairCollision
         void evaluateWall
         (
             typename CloudType::parcelType& p,
-            const List<point>& flatSites,
+            const List<point>& flatSitePoints,
-            const List<point>& sharpSites
+            const List<WallSiteData<vector> >& flatSiteData,
+            const List<point>& sharpSitePoints,
+            const List<WallSiteData<vector> >& sharpSiteData
         ) const;
 
 
@@ -171,6 +174,10 @@ public:
         //- Flag to indicate whether model activates injection model
         virtual bool active() const;
 
+        //- Indicates whether model determines wall collisions or not,
+        //  used to determine what value to use for wallImpactDistance
+        virtual bool controlsWallInteraction() const;
+
         // Collision function
         virtual void collide();
 };

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairModel/PairModel.H

@@ -54,7 +54,7 @@ class PairModel
 {
     // Private data
 
-        //- The cloud dictionary
+        //- The CollisionModel dictionary
         const dictionary& dict_;
 
         //- Reference to the owner cloud class

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C

@@ -45,13 +45,15 @@ void Foam::PairSpringSliderDashpot<CloudType>::findMinMaxProperties
 
         // Finding minimum diameter to avoid excessive arithmetic
 
-        RMin = min(p.d(), RMin);
+        scalar dEff = p.d()*cbrt(p.nParticle()*volumeFactor_);
+
+        RMin = min(dEff, RMin);
 
         rhoMax = max(p.rho(), rhoMax);
 
         UMagMax = max
         (
-            mag(p.U()) + mag(p.omega())*p.d()/2,
+            mag(p.U()) + mag(p.omega())*dEff/2,
             UMagMax
         );
     }
@@ -91,7 +93,8 @@ Foam::PairSpringSliderDashpot<CloudType>::PairSpringSliderDashpot
         (
             this->coeffDict().lookup("collisionResolutionSteps")
         )
-    )
+    ),
+    volumeFactor_(this->dict().lookupOrDefault("volumeFactor", 1.0))
 {
     scalar nu = this->owner().constProps().poissonsRatio();
 
@@ -155,7 +158,11 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
 {
     vector r_AB = (pA.position() - pB.position());
 
-    scalar normalOverlapMag = 0.5*(pA.d() + pB.d()) - mag(r_AB);
+    scalar dAEff = pA.d()*cbrt(pA.nParticle()*volumeFactor_);
+
+    scalar dBEff = pB.d()*cbrt(pB.nParticle()*volumeFactor_);
+
+    scalar normalOverlapMag = 0.5*(dAEff + dBEff) - mag(r_AB);
 
     if (normalOverlapMag > 0)
     {
@@ -166,7 +173,7 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
         vector U_AB = pA.U() - pB.U();
 
         // Effective radius
-        scalar R = 0.5*pA.d()*pB.d()/(pA.d() + pB.d());
+        scalar R = 0.5*dAEff*dBEff/(dAEff + dBEff);
 
         // Effective mass
         scalar M = pA.mass()*pB.mass()/(pA.mass() + pB.mass());
@@ -185,8 +192,8 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
 
         vector USlip_AB =
             U_AB - (U_AB & rHat_AB)*rHat_AB
-          + (pA.omega() ^ (pA.d()/2*-rHat_AB))
+          + (pA.omega() ^ (dAEff/2*-rHat_AB))
-          - (pB.omega() ^ (pB.d()/2*rHat_AB));
+          - (pB.omega() ^ (dBEff/2*rHat_AB));
 
         scalar deltaT = this->owner().mesh().time().deltaTValue();
 
@@ -215,7 +222,7 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
         {
             scalar kT = 8.0*sqrt(R*normalOverlapMag)*Gstar_;
 
-            scalar& etaT = etaN;
+            scalar etaT = etaN;
 
             // Tangential force
             vector fT_AB;
@@ -233,16 +240,16 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
             else
             {
                 fT_AB =
-                -kT*tangentialOverlapMag
+                    -kT*tangentialOverlapMag
-               *tangentialOverlap_AB/tangentialOverlapMag
+                   *tangentialOverlap_AB/tangentialOverlapMag
-              - etaT*USlip_AB;
+                  - etaT*USlip_AB;
             }
 
             pA.f() += fT_AB;
             pB.f() += -fT_AB;
 
-            pA.torque() += (pA.d()/2*-rHat_AB) ^ fT_AB;
+            pA.torque() += (dAEff/2*-rHat_AB) ^ fT_AB;
-            pB.torque() += (pB.d()/2*rHat_AB) ^ -fT_AB;
+            pB.torque() += (dBEff/2*rHat_AB) ^ -fT_AB;
         }
     }
 }

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.H

@@ -74,6 +74,24 @@ class PairSpringSliderDashpot
         //  harmonic approximation of the collision period
         scalar collisionResolutionSteps_;
 
+        //- Volume factor for determining the equivalent size of a
+        //  parcel where nParticles is not 1.  The equivalent size of
+        //  the parcel is
+        //      parcelEquivVolume = volumeFactor*nParticles*p.volume()
+        //  so
+        //      parcelEquivD = cbrt(volumeFactor*nParticles)*p.d()
+        //  + When volumeFactor = 1, the particles are compressed
+        //    together so that the equivalent volume of the parcel is
+        //    the sum of the constituent particles
+        //  + When volumeFactor = 3*sqrt(2)/pi, the particles are
+        //    close packed, but uncompressed.
+        //  + When volumeFactor > 3*sqrt(2)/pi, the particles loosely
+        //    grouped.
+        // 3*sqrt(2)/pi = 1.350474 is the volume factor for close
+        // packing, i.e pi/(3*sqrt(2)) is the maximum close packing
+        // factor
+        scalar volumeFactor_;
+
 
     // Private Member Functions
 
@@ -104,6 +122,12 @@ public:
 
     // Member Functions
 
+        //- Return the volumeFactor
+        inline scalar volumeFactor() const
+        {
+            return volumeFactor_;
+        }
+
         //- Whether the PairModel has a timestep limit that will
         //  require subCycling
         virtual bool controlsTimestep() const;

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallModel/WallModel.H

@@ -54,7 +54,7 @@ class WallModel
 {
     // Private data
 
-        //- The cloud dictionary
+        //- The CollisionModel dictionary
         const dictionary& dict_;
 
         //- Reference to the owner cloud class
@@ -120,6 +120,9 @@ public:
 
     // Member Functions
 
+        //- Return the effective radius for a particle for the model
+        virtual scalar pREff(const typename CloudType::parcelType& p) const = 0;
+
         //- Whether the WallModel has a timestep limit that will
         //  require subCycling
         virtual bool controlsTimestep() const = 0;
@@ -133,8 +136,10 @@ public:
         virtual void evaluateWall
         (
             typename CloudType::parcelType& p,
-            const List<point>& flatSites,
+            const List<point>& flatSitePoints,
-            const List<point>& sharpSites
+            const List<WallSiteData<vector> >& flatSiteData,
+            const List<point>& sharpSitePoints,
+            const List<WallSiteData<vector> >& sharpSiteData
         ) const = 0;
 };
 

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C

@@ -44,13 +44,16 @@ void Foam::WallSpringSliderDashpot<CloudType>::findMinMaxProperties
         const typename CloudType::parcelType& p = iter();
 
         // Finding minimum diameter to avoid excessive arithmetic
-        rMin = min(p.d(), rMin);
+
+        scalar dEff = p.d()*cbrt(p.nParticle()*volumeFactor_);
+
+        rMin = min(dEff, rMin);
 
         rhoMax = max(p.rho(), rhoMax);
 
         UMagMax = max
         (
-            mag(p.U()) + mag(p.omega())*p.d()/2,
+            mag(p.U()) + mag(p.omega())*dEff/2,
             UMagMax
         );
     }
@@ -66,25 +69,75 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
 (
     typename CloudType::parcelType& p,
     const point& site,
+    const WallSiteData<vector>& data,
     scalar pNu,
     scalar pE,
+    scalar pREff,
     scalar Estar,
-    scalar kN
+    scalar kN,
+    scalar Gstar
 ) const
 {
     vector r_PW = p.position() - site;
 
-    scalar normalOverlapMag = p.d()/2 - mag(r_PW);
+    vector U_PW = p.U() - data.wallData();
+
+    scalar normalOverlapMag = pREff - mag(r_PW);
 
     vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
 
     scalar etaN = alpha_*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
 
     vector fN_PW =
-    rHat_PW
+        rHat_PW
-    *(kN*pow(normalOverlapMag, b_) - etaN*(p.U() & rHat_PW));
+       *(kN*pow(normalOverlapMag, b_) - etaN*(U_PW & rHat_PW));
 
     p.f() += fN_PW;
+
+    vector USlip_PW =
+        U_PW - (U_PW & rHat_PW)*rHat_PW
+      + (p.omega() ^ (pREff*-rHat_PW));
+
+    scalar deltaT = this->owner().mesh().time().deltaTValue();
+
+    // For remembering previous overlap
+    // vector deltaTangentialOverlap_PW = USlip_PW * deltaT;
+    // tangentialOverlap_PW += deltaTangentialOverlap_PW;
+
+    vector tangentialOverlap_PW = USlip_PW * deltaT;
+
+    scalar tangentialOverlapMag = mag(tangentialOverlap_PW);
+
+    if (tangentialOverlapMag > VSMALL)
+    {
+        scalar kT = 8.0*sqrt(pREff*normalOverlapMag)*Gstar;
+
+        scalar etaT = etaN;
+
+        // Tangential force
+        vector fT_PW;
+
+        if (kT*tangentialOverlapMag > mu_*mag(fN_PW))
+        {
+            // Tangential force greater than sliding friction,
+            // particle slips
+
+            fT_PW = -mu_*mag(fN_PW)*USlip_PW/mag(USlip_PW);
+
+            // tangentialOverlap_PW = vector::zero;
+        }
+        else
+        {
+            fT_PW =
+                -kT*tangentialOverlapMag
+               *tangentialOverlap_PW/tangentialOverlapMag
+              - etaT*USlip_PW;
+        }
+
+        p.f() += fT_PW;
+
+        p.torque() += (pREff*-rHat_PW) ^ fT_PW;
+    }
 }
 
 
@@ -109,7 +162,8 @@ Foam::WallSpringSliderDashpot<CloudType>::WallSpringSliderDashpot
         (
             this->coeffDict().lookup("collisionResolutionSteps")
         )
-    )
+    ),
+    volumeFactor_(this->dict().lookupOrDefault("volumeFactor", 1.0))
 {}
 
 
@@ -122,6 +176,15 @@ Foam::WallSpringSliderDashpot<CloudType>::~WallSpringSliderDashpot()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+Foam::scalar Foam::WallSpringSliderDashpot<CloudType>::pREff
+(
+    const typename CloudType::parcelType& p
+) const
+{
+    return p.d()/2*cbrt(p.nParticle()*volumeFactor_);
+}
+
 template<class CloudType>
 bool Foam::WallSpringSliderDashpot<CloudType>::controlsTimestep() const
 {
@@ -164,28 +227,56 @@ template<class CloudType>
 void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
 (
     typename CloudType::parcelType& p,
-    const List<point>& flatSites,
+    const List<point>& flatSitePoints,
-    const List<point>& sharpSites
+    const List<WallSiteData<vector> >& flatSiteData,
+    const List<point>& sharpSitePoints,
+    const List<WallSiteData<vector> >& sharpSiteData
 ) const
 {
     scalar pNu = this->owner().constProps().poissonsRatio();
 
     scalar pE = this->owner().constProps().youngsModulus();
 
+    scalar pREff = this->pREff(p);
+
     scalar Estar = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu_))/E_);
 
-    scalar kN = (4.0/3.0)*sqrt(p.d()/2)*Estar;
+    scalar kN = (4.0/3.0)*sqrt(pREff)*Estar;
 
-    forAll(flatSites, siteI)
+    scalar GStar = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu_ - sqr(nu_))/E_));
+
+    forAll(flatSitePoints, siteI)
     {
-        evaluateWall(p, flatSites[siteI], pNu, pE, Estar, kN);
+        evaluateWall
+        (
+            p,
+            flatSitePoints[siteI],
+            flatSiteData[siteI],
+            pNu,
+            pE,
+            pREff,
+            Estar,
+            kN,
+            GStar
+        );
     }
 
-    forAll(sharpSites, siteI)
+    forAll(sharpSitePoints, siteI)
     {
         // Treating sharp sites like flat sites
 
-        evaluateWall(p, sharpSites[siteI], pNu, pE, Estar, kN);
+        evaluateWall
+        (
+            p,
+            sharpSitePoints[siteI],
+            sharpSiteData[siteI],
+            pNu,
+            pE,
+            pREff,
+            Estar,
+            kN,
+            GStar
+        );
     }
 }
 

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H

@@ -69,6 +69,24 @@ class WallSpringSliderDashpot
         //  harmonic approximation of the collision period
         scalar collisionResolutionSteps_;
 
+        //- Volume factor for determining the equivalent size of a
+        //  parcel where nParticles is not 1.  The equivalent size of
+        //  the parcel is
+        //      parcelEquivVolume = volumeFactor*nParticles*p.volume()
+        //  so
+        //      parcelEquivD = cbrt(volumeFactor*nParticles)*p.d()
+        //  + When volumeFactor = 1, the particles are compressed
+        //    together so that the equivalent volume of the parcel is
+        //    the sum of the constituent particles
+        //  + When volumeFactor = 3*sqrt(2)/pi, the particles are
+        //    close packed, but uncompressed.
+        //  + When volumeFactor > 3*sqrt(2)/pi, the particles loosely
+        //    grouped.
+        // 3*sqrt(2)/pi = 1.350474 is the volume factor for close
+        // packing, i.e pi/(3*sqrt(2)) is the maximum close packing
+        // factor
+        scalar volumeFactor_;
+
 
     // Private Member Functions
 
@@ -86,10 +104,13 @@ class WallSpringSliderDashpot
         (
             typename CloudType::parcelType& p,
             const point& site,
+            const WallSiteData<vector>& data,
             scalar pNu,
             scalar pE,
+            scalar pREff,
             scalar Estar,
-            scalar kN
+            scalar kN,
+            scalar Gstar
         ) const;
 
 
@@ -111,6 +132,15 @@ public:
 
     // Member Functions
 
+        //- Return the volumeFactor
+        inline scalar volumeFactor() const
+        {
+            return volumeFactor_;
+        }
+
+        //- Return the effective radius for a particle for the model
+        virtual scalar pREff(const typename CloudType::parcelType& p) const;
+
         //- Whether the WallModel has a timestep limit that will
         //  require subCycling
         virtual bool controlsTimestep() const;
@@ -124,8 +154,10 @@ public:
         virtual void evaluateWall
         (
             typename CloudType::parcelType& p,
-            const List<point>& flatSites,
+            const List<point>& flatSitePoints,
-            const List<point>& sharpSites
+            const List<WallSiteData<vector> >& flatSiteData,
+            const List<point>& sharpSitePoints,
+            const List<WallSiteData<vector> >& sharpSiteData
         ) const;
 };
 

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallSiteData/WallSiteData.C

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "WallSiteData.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::WallSiteData<Type>::WallSiteData()
+:
+    patchI_(),
+    wallData_()
+{}
+
+
+template<class Type>
+Foam::WallSiteData<Type>::WallSiteData
+(
+    label patchI,
+    const Type& wallData
+)
+:
+    patchI_(patchI),
+    wallData_(wallData)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::WallSiteData<Type>::~WallSiteData()
+{}
+
+
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+template<class Type>
+bool Foam::WallSiteData<Type>::operator==
+(
+    const WallSiteData<Type>& rhs
+) const
+{
+    return patchI_ == rhs.patchI_ && wallData_ == rhs.wallData_;
+}
+
+
+template<class Type>
+bool Foam::WallSiteData<Type>::operator!=
+(
+    const WallSiteData<Type>& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Istream& Foam::operator>>
+(
+    Istream& is,
+    WallSiteData<Type>& wIS
+)
+{
+    is  >> wIS.patchI_ >> wIS.wallData_;
+
+    // Check state of Istream
+    is.check
+    (
+        "Foam::Istream& Foam::operator>>"
+        "(Foam::Istream&, Foam::WallSiteData<Type>&)"
+    );
+
+    return is;
+}
+
+
+template<class Type>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const WallSiteData<Type>& wIS
+)
+{
+    os  << wIS.patchI_ << token::SPACE << wIS.wallData_;
+
+    // Check state of Ostream
+    os.check
+    (
+        "Foam::Ostream& Foam::operator<<"
+        "(Ostream&, const WallSiteData<Type>&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallSiteData/WallSiteData.H

@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::WallSiteData
+
+Description
+    Stores the patch ID and templated data to represent a collision
+    with a wall to be passed to the wall model.
+
+SourceFiles
+    WallSiteDataI.H
+    WallSiteData.C
+    WallSiteDataIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef WallSiteData_H
+#define WallSiteData_H
+
+#include "label.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class Type>
+class WallSiteData;
+
+template<class Type>
+Istream& operator>>(Istream&, WallSiteData<Type>&);
+
+template<class Type>
+Ostream& operator<<(Ostream&, const WallSiteData<Type>&);
+
+
+/*---------------------------------------------------------------------------*\
+                         Class WallSiteData Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class WallSiteData
+{
+    // Private data
+
+        //- Index of originating patch
+        label patchI_;
+
+        //- Wall data
+        Type wallData_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        WallSiteData();
+
+        //- Construct from components
+        WallSiteData
+        (
+            label patchI,
+            const Type& wallData
+        );
+
+
+    //- Destructor
+    ~WallSiteData();
+
+
+    // Member Functions
+
+            //- Return access to the patch index
+            inline label patchIndex() const;
+
+            //- Return non-const access to the patch index
+            inline label& patchIndex();
+
+            //- Return access to wall data
+            inline const Type& wallData() const;
+
+            //- Return non-const access to wall data
+            inline Type& wallData();
+
+
+    // Member Operators
+
+        bool operator==(const WallSiteData<Type>&) const;
+        bool operator!=(const WallSiteData<Type>&) const;
+
+
+    // IOstream Operators
+
+        friend Istream& operator>> <Type>
+        (Istream&, WallSiteData<Type>&);
+
+        friend Ostream& operator<< <Type>
+        (Ostream&, const WallSiteData<Type>&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "WallSiteDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "WallSiteData.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallSiteData/WallSiteDataI.H

@@ -0,0 +1,60 @@
+#/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::label Foam::WallSiteData<Type>::patchIndex() const
+{
+    return patchI_;
+}
+
+
+template<class Type>
+Foam::label& Foam::WallSiteData<Type>::patchIndex()
+{
+    return patchI_;
+}
+
+
+template<class Type>
+const Type& Foam::WallSiteData<Type>::wallData() const
+{
+    return wallData_;
+}
+
+
+template<class Type>
+Type& Foam::WallSiteData<Type>::wallData()
+{
+    return wallData_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// ************************************************************************* //

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

@@ -129,10 +129,11 @@ void Foam::InjectionModel<CloudType>::prepareForNextTimeStep
 
 
 template<class CloudType>
-void Foam::InjectionModel<CloudType>::findCellAtPosition
+bool Foam::InjectionModel<CloudType>::findCellAtPosition
 (
     label& cellI,
-    vector& position
+    vector& position,
+    bool errorOnNotFound
 )
 {
     const volVectorField& cellCentres = owner_.mesh().C();
@@ -176,17 +177,26 @@ void Foam::InjectionModel<CloudType>::findCellAtPosition
 
     if (procI == -1)
     {
-        FatalErrorIn
+        if (errorOnNotFound)
-        (
+        {
-            "Foam::InjectionModel<CloudType>::findCellAtPosition"
+            FatalErrorIn
-            "("
+            (
-                "label&, "
+                "Foam::InjectionModel<CloudType>::findCellAtPosition"
-                "vector&"
+                "("
-            ")"
+                    "label&, "
-        )<< "Cannot find parcel injection cell. "
+                    "vector&"
-         << "Parcel position = " << p0 << nl
+                ")"
-         << abort(FatalError);
+            )   << "Cannot find parcel injection cell. "
+                << "Parcel position = " << p0 << nl
+                << abort(FatalError);
+        }
+        else
+        {
+            return false;
+        }
     }
+
+    return true;
 }
 
 
@@ -217,7 +227,7 @@ Foam::scalar Foam::InjectionModel<CloudType>::setNumberOfParticles
         }
         case pbFixed:
         {
-            nP = nParticlesFixed_;
+            nP = nParticleFixed_;
             break;
         }
         default:
@@ -290,7 +300,7 @@ Foam::InjectionModel<CloudType>::InjectionModel(CloudType& owner)
     nInjections_(0),
     parcelsAddedTotal_(0),
     parcelBasis_(pbNumber),
-    nParticlesFixed_(0.0),
+    nParticleFixed_(0.0),
     time0_(0.0),
     timeStep0_(0.0)
 {
@@ -316,7 +326,7 @@ Foam::InjectionModel<CloudType>::InjectionModel
     nInjections_(0),
     parcelsAddedTotal_(0),
     parcelBasis_(pbNumber),
-    nParticlesFixed_(0.0),
+    nParticleFixed_(0.0),
     time0_(owner.db().time().value()),
     timeStep0_(0.0)
 {
@@ -340,11 +350,11 @@ Foam::InjectionModel<CloudType>::InjectionModel
     {
         parcelBasis_ = pbFixed;
 
-        Info<< "    Choosing nParticles to be a fixed value, massTotal "
+        Info<< "    Choosing nParticle to be a fixed value, massTotal "
             << "variable now does not determine anything."
             << endl;
 
-        nParticlesFixed_ = readScalar(coeffDict_.lookup("nParticles"));
+        nParticleFixed_ = readScalar(coeffDict_.lookup("nParticle"));
     }
     else
     {

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

@@ -138,9 +138,9 @@ protected:
             //- Parcel basis enumeration
             parcelBasis parcelBasis_;
 
-            //- nParticles to assign to parcels when the 'fixed' basis
+            //- nParticle to assign to parcels when the 'fixed' basis
             //  is selected
-            scalar nParticlesFixed_;
+            scalar nParticleFixed_;
 
             //- Continuous phase time at start of injection time step [s]
             scalar time0_;
@@ -180,7 +180,12 @@ protected:
         //- Find the cell that contains the supplied position
         //  Will modify position slightly towards the owner cell centroid to
         //  ensure that it lies in a cell and not edge/face
-        virtual void findCellAtPosition(label& cellI, vector& position);
+        virtual bool findCellAtPosition
+        (
+            label& cellI,
+            vector& position,
+            bool errorOnNotFound = true
+        );
 
         //- Set number of particles to inject given parcel properties
         virtual scalar setNumberOfParticles

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

@@ -25,6 +25,8 @@ License
 
 #include "ManualInjection.H"
 #include "mathematicalConstants.H"
+#include "PackedBoolList.H"
+#include "Switch.H"
 
 using namespace Foam::constant::mathematical;
 
@@ -100,6 +102,40 @@ Foam::ManualInjection<CloudType>::ManualInjection
         )
     )
 {
+    Switch checkAndIgnoreOutOfBounds
+    (
+        this->coeffDict().lookupOrDefault("checkAndIgnoreOutOfBounds", false)
+    );
+
+    label nRejected = 0;
+
+    if (checkAndIgnoreOutOfBounds)
+    {
+        // Dummy cell
+        label cellI = -1;
+
+        PackedBoolList keep(positions_.size(), true);
+
+        forAll(positions_, pI)
+        {
+            if (!this->findCellAtPosition(cellI, positions_[pI], false))
+            {
+                keep[pI] = false;
+
+                nRejected++;
+            }
+        }
+
+        if (nRejected > 0)
+        {
+            inplaceSubset(keep, positions_);
+            inplaceSubset(keep, diameters_);
+
+            Info<< "    " << nRejected
+                << " particles ignored, out of bounds." << endl;
+        }
+    }
+
     // Construct parcel diameters
     forAll(diameters_, i)
     {

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H

@@ -36,7 +36,7 @@ Description
 
 SourceFiles
     sixDoFRigidBodyMotionConstraint.C
-    dynamicFvMeshNew.C
+    sixDoFRigidBodyMotionConstraintNew.C
 
 \*---------------------------------------------------------------------------*/
 

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H

@@ -36,7 +36,7 @@ Description
 
 SourceFiles
     sixDoFRigidBodyMotionRestraint.C
-    dynamicFvMeshNew.C
+    sixDoFRigidBodyMotionRestraintNew.C
 
 \*---------------------------------------------------------------------------*/
 

src/turbulenceModels/incompressible/RAS/Make/files

@@ -44,6 +44,8 @@ derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFv
 derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
 derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
 derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.C
+derivedFvPatchFields/ABLInletEpsilon/ABLInletEpsilonFvPatchScalarField.C
+derivedFvPatchFields/ABLInletVelocity/ABLInletVelocityFvPatchVectorField.C
 
 backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
 

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletEpsilon/ABLInletEpsilonFvPatchScalarField.C

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "ABLInletEpsilonFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+ABLInletEpsilonFvPatchScalarField::ABLInletEpsilonFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    Ustar_(0),
+    z_(pTraits<vector>::zero),
+    z0_(0),
+    kappa_(0.41)
+{}
+
+
+ABLInletEpsilonFvPatchScalarField::ABLInletEpsilonFvPatchScalarField
+(
+    const ABLInletEpsilonFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    Ustar_(ptf.Ustar_),
+    z_(ptf.z_),
+    z0_(ptf.z0_),
+    kappa_(ptf.kappa_)
+{}
+
+
+ABLInletEpsilonFvPatchScalarField::ABLInletEpsilonFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    Ustar_(readScalar(dict.lookup("Ustar"))),
+    z_(dict.lookup("z")),
+    z0_(readScalar(dict.lookup("z0"))),
+    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41))
+{
+    if (mag(z_) < SMALL)
+    {
+        FatalErrorIn("ABLInletEpsilonFvPatchScalarField(dict)")
+            << "z is not correct"
+            << abort(FatalError);
+    }
+
+    z_ /= mag(z_);
+
+    evaluate();
+}
+
+
+ABLInletEpsilonFvPatchScalarField::ABLInletEpsilonFvPatchScalarField
+(
+    const ABLInletEpsilonFvPatchScalarField& fcvpvf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(fcvpvf, iF),
+    Ustar_(fcvpvf.Ustar_),
+    z_(fcvpvf.z_),
+    z0_(fcvpvf.z0_),
+    kappa_(fcvpvf.kappa_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void ABLInletEpsilonFvPatchScalarField::updateCoeffs()
+{
+    const vectorField& c = patch().Cf();
+    scalarField coord = (c & z_);
+    scalarField::operator=(pow(Ustar_, 3.0)/(kappa_*(coord + z0_)));
+}
+
+
+// Write
+void ABLInletEpsilonFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("Ustar")
+        << Ustar_ << token::END_STATEMENT << nl;
+    os.writeKeyword("z")
+        << z_ << token::END_STATEMENT << nl;
+    os.writeKeyword("z0")
+        << z0_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa")
+        << kappa_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField(fvPatchScalarField, ABLInletEpsilonFvPatchScalarField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletEpsilon/ABLInletEpsilonFvPatchScalarField.H

@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    ABLInletEpsilonFvPatchScalarField
+
+Description
+    Boundary condition specifies a epsilon inlet for the atmospheric boundary
+    layer (ABL). This boundaty is to be used in conjunction with
+    ABLInletVelocity.
+
+    @verbatim
+        epsilon = Ustar^3 / (K(z + z0))
+
+    where:
+
+        Ustar is the frictional velocity
+        K is karman's constant
+        z is the verical coordinate
+        z0 is the surface roughness lenght
+
+    @endverbatim
+
+    Reference:
+    D.M. Hargreaves and N.G. Wright
+    "On the use of the k-epsilon model in commercial CFD software to model the
+     neutral atmospheric boundary layer"
+    Journal of Wind Engineering and Industrial Aerodynamics 95(2007) 355-369.
+
+SourceFiles
+    ABLInletEpsilonFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ABLInletEpsilonFvPatchScalarField_H
+#define ABLInletEpsilonFvPatchScalarField_H
+
+#include "fvPatchFields.H"
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+
+/*---------------------------------------------------------------------------*\
+              Class ABLInletEpsilonFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class ABLInletEpsilonFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+    // Private data
+
+        //- Peak velocity magnitude
+        scalar Ustar_;
+
+        //- Direction of the z-coordinate
+        vector z_;
+
+         //- Surface roughness lenght
+        scalar z0_;
+
+         //- Von Karman constant
+        scalar kappa_;
+
+public:
+
+    //- Runtime type information
+    TypeName("ABLInletEpsilon");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        ABLInletEpsilonFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        ABLInletEpsilonFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given ABLInletEpsilonFvPatchScalarField
+        //  onto a new patch
+        ABLInletEpsilonFvPatchScalarField
+        (
+            const ABLInletEpsilonFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new ABLInletEpsilonFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        ABLInletEpsilonFvPatchScalarField
+        (
+            const ABLInletEpsilonFvPatchScalarField&,
+            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 ABLInletEpsilonFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Return max value
+        scalar& Ustar()
+        {
+            return Ustar_;
+        }
+
+        //- Return z direction
+        vector& z()
+        {
+            return z_;
+        }
+
+        //- Update coefficients
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace incompressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletVelocity/ABLInletVelocityFvPatchVectorField.C

@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "ABLInletVelocityFvPatchVectorField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+ABLInletVelocityFvPatchVectorField::ABLInletVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchVectorField(p, iF),
+    Ustar_(0),
+    n_(pTraits<vector>::zero),
+    z_(pTraits<vector>::zero),
+    z0_(0),
+    kappa_(0.41)
+{}
+
+
+ABLInletVelocityFvPatchVectorField::ABLInletVelocityFvPatchVectorField
+(
+    const ABLInletVelocityFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchVectorField(ptf, p, iF, mapper),
+    Ustar_(ptf.Ustar_),
+    n_(ptf.n_),
+    z_(ptf.z_),
+    z0_(ptf.z0_),
+    kappa_(ptf.kappa_)
+{}
+
+
+ABLInletVelocityFvPatchVectorField::ABLInletVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchVectorField(p, iF),
+    Ustar_(readScalar(dict.lookup("Ustar"))),
+    n_(dict.lookup("n")),
+    z_(dict.lookup("z")),
+    z0_(readScalar(dict.lookup("z0"))),
+    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41))
+{
+    if (mag(n_) < SMALL || mag(z_) < SMALL)
+    {
+        FatalErrorIn("ABLInletVelocityFvPatchVectorField(dict)")
+            << "n or z given with zero size not correct"
+            << abort(FatalError);
+    }
+
+    n_ /= mag(n_);
+    z_ /= mag(z_);
+
+    evaluate();
+}
+
+
+ABLInletVelocityFvPatchVectorField::ABLInletVelocityFvPatchVectorField
+(
+    const ABLInletVelocityFvPatchVectorField& fcvpvf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchVectorField(fcvpvf, iF),
+    Ustar_(fcvpvf.Ustar_),
+    n_(fcvpvf.n_),
+    z_(fcvpvf.z_),
+    z0_(fcvpvf.z0_),
+    kappa_(fcvpvf.kappa_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void ABLInletVelocityFvPatchVectorField::updateCoeffs()
+{
+    const vectorField& c = patch().Cf();
+    scalarField coord = (c & z_);
+    vectorField::operator=(n_*(Ustar_/kappa_)*log((coord + z0_)/z0_));
+}
+
+
+// Write
+void ABLInletVelocityFvPatchVectorField::write(Ostream& os) const
+{
+    fvPatchVectorField::write(os);
+    os.writeKeyword("Ustar")
+        << Ustar_ << token::END_STATEMENT << nl;
+    os.writeKeyword("z0")
+        << z0_ << token::END_STATEMENT << nl;
+    os.writeKeyword("n")
+        << n_ << token::END_STATEMENT << nl;
+    os.writeKeyword("z")
+        << z_ << token::END_STATEMENT << nl;
+     os.writeKeyword("kappa")
+        << kappa_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField(fvPatchVectorField, ABLInletVelocityFvPatchVectorField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletVelocity/ABLInletVelocityFvPatchVectorField.H

@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    ABLInletVelocityFvPatchVectorField
+
+Description
+    Boundary condition specifies a atmospheric boundary layer (ABL)
+    velocity inlet profile given the friction velocity value,
+    flow direction n and direction of the parabolic coordinate z.
+
+    @verbatim
+        U = (Ustar/K) ln((z + z0)/z0)
+
+    where:
+
+        Ustar is the frictional velocity
+        K is karman's constant
+        z0 is the surface roughness lenght
+        z is the verical coordinate
+
+    and:
+
+        Ustar = K Uref/ln((Zref + z0)/z0)
+
+    where:
+
+        Uref is the reference velocity at Zref
+        Zref is the reference height.
+
+    @endverbatim
+
+    Reference:
+    D.M. Hargreaves and N.G. Wright
+    "On the use of the k-epsilon model in commercial CFD software to model the
+     neutral atmospheric boundary layer"
+    Journal of Wind Engineering and Industrial Aerodynamics 95(2007) 355-369.
+
+NOTE: D.M. Hargreaves and N.G. Wright recommend Gamma epsilon in the k-epsilon
+      model should be changed from 1.3 to 1.11 for consistency.
+      The roughness height (Er) is given by Er = 20 z0 following the same
+      reference
+
+SourceFiles
+    ABLInletVelocityFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ABLInletVelocityFvPatchVectorField_H
+#define ABLInletVelocityFvPatchVectorField_H
+
+#include "fvPatchFields.H"
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+
+/*---------------------------------------------------------------------------*\
+             Class ABLInletVelocityFvPatchVectorField Declaration
+\*---------------------------------------------------------------------------*/
+
+class ABLInletVelocityFvPatchVectorField
+:
+    public fixedValueFvPatchVectorField
+{
+    // Private data
+
+        //- Friction velocity
+        scalar Ustar_;
+
+        //- Flow direction
+        vector n_;
+
+        //- Direction of the z-coordinate
+        vector z_;
+
+        //- Surface roughness lenght
+        scalar z0_;
+
+        //- Von Karman constant
+        scalar kappa_;
+
+public:
+
+    //- Runtime type information
+    TypeName("ABLInletVelocity");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        ABLInletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        ABLInletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given ABLInletVelocityFvPatchVectorField
+        //  onto a new patch
+        ABLInletVelocityFvPatchVectorField
+        (
+            const ABLInletVelocityFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new ABLInletVelocityFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        ABLInletVelocityFvPatchVectorField
+        (
+            const ABLInletVelocityFvPatchVectorField&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchVectorField> clone
+        (
+            const DimensionedField<vector, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new ABLInletVelocityFvPatchVectorField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Return max value
+        scalar& Ustar()
+        {
+            return Ustar_;
+        }
+
+        //- Return flow direction
+        vector& n()
+        {
+            return n_;
+        }
+
+        //- Return y direction
+        vector& z()
+        {
+            return z_;
+        }
+
+        //- Update coefficients
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace incompressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -10,7 +10,7 @@ License
 
     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
+    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

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -10,7 +10,7 @@ License
 
     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
+    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

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.H

@@ -152,6 +152,21 @@ public:
 
     // Member functions
 
+        // Acces functions
+
+            // Return Ks
+             scalarField& Ks()
+            {
+                return Ks_;
+            }
+
+            // Return Cs
+            scalarField& Cs()
+            {
+                return Cs_;
+            }
+
+
         // Mapping functions
 
             //- Map (and resize as needed) from self given a mapping object