Merge remote-tracking branch 'origin/master' into feature/procAgglom

applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    uncoupledKinematicParcelFoam
+    icoUncoupledKinematicParcelFoam
 
 Description
     Transient solver for the passive transport of a single kinematic

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,7 @@ Application
 
 Description
     Transient solver for the passive transport of a single kinematic
-    particle could.
+    particle cloud.
 
     Uses a pre- calculated velocity field to evolve the cloud.
 

applications/test/dictionary/testDict

@@ -47,6 +47,10 @@ x 5;
 varName x;
 
 
+//Indirection for keys
+key inlet_9;
+
+
 boundaryField
 {
     Default_Boundary_Region
@@ -63,6 +67,7 @@ boundaryField
     inlet_6  { $.inactive }     // Test scoping
     inlet_7  { ${inactive}}     // Test variable expansion
     inlet_8  { $inactive }
+    ${key}   { $inactive }
 
     #include "testDictInc"
 

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -402,9 +402,10 @@ addLayersControls
         nLayerIter 50;
 
         // Max number of iterations after which relaxed meshQuality controls
-        // get used. Up to nRelaxIter it uses the settings in
+        // get used. Up to nRelaxedIter it uses the settings in
         // meshQualityControls,
-        // after nRelaxIter it uses the values in meshQualityControls::relaxed.
+        // after nRelaxedIter it uses the values in
+        // meshQualityControls::relaxed.
         nRelaxedIter 20;
 
         // Additional reporting: if there are just a few faces where there

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C

@@ -52,8 +52,8 @@ Foam::porosityModels::DarcyForchheimer::DarcyForchheimer
 )
 :
     porosityModel(name, modelType, mesh, dict, cellZoneName),
-    D_(cellZoneIds_.size()),
+    D_(cellZoneIDs_.size()),
-    F_(cellZoneIds_.size()),
+    F_(cellZoneIDs_.size()),
     rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho")),
     muName_(coeffs_.lookupOrDefault<word>("mu", "thermo:mu")),
     nuName_(coeffs_.lookupOrDefault<word>("nu", "nu"))
@@ -67,7 +67,7 @@ Foam::porosityModels::DarcyForchheimer::DarcyForchheimer
 
     if (coordSys_.R().uniform())
     {
-        forAll (cellZoneIds_, zoneI)
+        forAll (cellZoneIDs_, zoneI)
         {
             D_[zoneI].setSize(1, tensor::zero);
             F_[zoneI].setSize(1, tensor::zero);
@@ -89,9 +89,9 @@ Foam::porosityModels::DarcyForchheimer::DarcyForchheimer
     }
     else
     {
-        forAll(cellZoneIds_, zoneI)
+        forAll(cellZoneIDs_, zoneI)
         {
-            const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+            const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
 
             D_[zoneI].setSize(cells.size(), tensor::zero);
             F_[zoneI].setSize(cells.size(), tensor::zero);
@@ -122,6 +122,24 @@ Foam::porosityModels::DarcyForchheimer::~DarcyForchheimer()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::porosityModels::DarcyForchheimer::calcForce
+(
+    const volVectorField& U,
+    const volScalarField& rho,
+    const volScalarField& mu,
+    vectorField& force
+) const
+{
+    scalarField Udiag(U.size());
+    vectorField Usource(U.size());
+    const scalarField& V = mesh_.V();
+
+    apply(Udiag, Usource, V, rho, mu, U);
+
+    force = Udiag*U - Usource;
+}
+
+
 void Foam::porosityModels::DarcyForchheimer::correct
 (
     fvVectorMatrix& UEqn
@@ -194,10 +212,12 @@ void Foam::porosityModels::DarcyForchheimer::correct
 }
 
 
-void Foam::porosityModels::DarcyForchheimer::writeData(Ostream& os) const
+bool Foam::porosityModels::DarcyForchheimer::writeData(Ostream& os) const
 {
     os  << indent << name_ << endl;
     dict_.write(os);
+
+    return true;
 }
 
 

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.H

@@ -143,6 +143,15 @@ public:
 
     // Member Functions
 
+        //- Calculate the porosity force
+        virtual void calcForce
+        (
+            const volVectorField& U,
+            const volScalarField& rho,
+            const volScalarField& mu,
+            vectorField& force
+        ) const;
+
         //- Add resistance
         virtual void correct(fvVectorMatrix& UEqn) const;
 
@@ -165,7 +174,7 @@ public:
     // I-O
 
         //- Write
-        void writeData(Ostream& os) const;
+        bool writeData(Ostream& os) const;
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimerTemplates.C

@@ -36,12 +36,12 @@ void Foam::porosityModels::DarcyForchheimer::apply
     const vectorField& U
 ) const
 {
-    forAll(cellZoneIds_, zoneI)
+    forAll(cellZoneIDs_, zoneI)
     {
         const tensorField& dZones = D_[zoneI];
         const tensorField& fZones = F_[zoneI];
 
-        const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+        const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
 
         forAll(cells, i)
         {
@@ -68,12 +68,12 @@ void Foam::porosityModels::DarcyForchheimer::apply
     const vectorField& U
 ) const
 {
-    forAll(cellZoneIds_, zoneI)
+    forAll(cellZoneIDs_, zoneI)
     {
         const tensorField& dZones = D_[zoneI];
         const tensorField& fZones = F_[zoneI];
 
-        const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+        const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
 
         forAll(cells, i)
         {

src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.C

@@ -50,12 +50,12 @@ void Foam::porosityModels::fixedCoeff::apply
     const scalar rho
 ) const
 {
-    forAll(cellZoneIds_, zoneI)
+    forAll(cellZoneIDs_, zoneI)
     {
         const tensorField& alphaZones = alpha_[zoneI];
         const tensorField& betaZones = beta_[zoneI];
 
-        const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+        const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
 
         forAll(cells, i)
         {
@@ -79,12 +79,12 @@ void Foam::porosityModels::fixedCoeff::apply
 ) const
 {
 
-    forAll(cellZoneIds_, zoneI)
+    forAll(cellZoneIDs_, zoneI)
     {
         const tensorField& alphaZones = alpha_[zoneI];
         const tensorField& betaZones = beta_[zoneI];
 
-        const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+        const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
 
         forAll(cells, i)
         {
@@ -111,8 +111,8 @@ Foam::porosityModels::fixedCoeff::fixedCoeff
 )
 :
     porosityModel(name, modelType, mesh, dict, cellZoneName),
-    alpha_(cellZoneIds_.size()),
+    alpha_(cellZoneIDs_.size()),
-    beta_(cellZoneIds_.size())
+    beta_(cellZoneIDs_.size())
 {
     dimensionedVector alpha(coeffs_.lookup("alpha"));
     dimensionedVector beta(coeffs_.lookup("beta"));
@@ -122,7 +122,7 @@ Foam::porosityModels::fixedCoeff::fixedCoeff
 
     if (coordSys_.R().uniform())
     {
-        forAll (cellZoneIds_, zoneI)
+        forAll (cellZoneIDs_, zoneI)
         {
             alpha_[zoneI].setSize(1, tensor::zero);
             beta_[zoneI].setSize(1, tensor::zero);
@@ -140,9 +140,9 @@ Foam::porosityModels::fixedCoeff::fixedCoeff
     }
     else
     {
-        forAll(cellZoneIds_, zoneI)
+        forAll(cellZoneIDs_, zoneI)
         {
-            const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+            const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
 
             alpha_[zoneI].setSize(cells.size(), tensor::zero);
             beta_[zoneI].setSize(cells.size(), tensor::zero);
@@ -175,6 +175,25 @@ Foam::porosityModels::fixedCoeff::~fixedCoeff()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::porosityModels::fixedCoeff::calcForce
+(
+    const volVectorField& U,
+    const volScalarField& rho,
+    const volScalarField& mu,
+    vectorField& force
+) const
+{
+    scalarField Udiag(U.size());
+    vectorField Usource(U.size());
+    const scalarField& V = mesh_.V();
+    scalar rhoRef = readScalar(coeffs_.lookup("rhoRef"));
+
+    apply(Udiag, Usource, V, U, rhoRef);
+
+    force = Udiag*U - Usource;
+}
+
+
 void Foam::porosityModels::fixedCoeff::correct
 (
     fvVectorMatrix& UEqn
@@ -235,10 +254,12 @@ void Foam::porosityModels::fixedCoeff::correct
 }
 
 
-void Foam::porosityModels::fixedCoeff::writeData(Ostream& os) const
+bool Foam::porosityModels::fixedCoeff::writeData(Ostream& os) const
 {
     os  << indent << name_ << endl;
     dict_.write(os);
+
+    return true;
 }
 
 

src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.H

@@ -119,6 +119,15 @@ public:
 
     // Member Functions
 
+        //- Calculate the porosity force
+        virtual void calcForce
+        (
+            const volVectorField& U,
+            const volScalarField& rho,
+            const volScalarField& mu,
+            vectorField& force
+        ) const;
+
         //- Add resistance
         virtual void correct(fvVectorMatrix& UEqn) const;
 
@@ -141,7 +150,7 @@ public:
     // I-O
 
         //- Write
-        void writeData(Ostream& os) const;
+        bool writeData(Ostream& os) const;
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C

@@ -88,13 +88,14 @@ Foam::porosityModel::porosityModel
     const word& cellZoneName
 )
 :
+    MeshObject<fvMesh, Foam::UpdateableMeshObject, porosityModel>(mesh),
     name_(name),
     mesh_(mesh),
     dict_(dict),
     coeffs_(dict.subDict(modelType + "Coeffs")),
     active_(true),
     zoneName_(cellZoneName),
-    cellZoneIds_(),
+    cellZoneIDs_(),
     coordSys_(coordinateSystem::New(mesh, coeffs_))
 {
     if (zoneName_ == word::null)
@@ -103,11 +104,11 @@ Foam::porosityModel::porosityModel
         dict_.lookup("cellZone") >> zoneName_;
     }
 
-    cellZoneIds_ = mesh_.cellZones().findIndices(zoneName_);
+    cellZoneIDs_ = mesh_.cellZones().findIndices(zoneName_);
 
     Info<< "    creating porous zone: " << zoneName_ << endl;
 
-    bool foundZone = !cellZoneIds_.empty();
+    bool foundZone = !cellZoneIDs_.empty();
     reduce(foundZone, orOp<bool>());
 
     if (!foundZone && Pstream::master())
@@ -136,12 +137,30 @@ Foam::porosityModel::~porosityModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::tmp<Foam::vectorField> Foam::porosityModel::porosityModel::force
+(
+    const volVectorField& U,
+    const volScalarField& rho,
+    const volScalarField& mu
+) const
+{
+    tmp<vectorField> tforce(new vectorField(U.size(), vector::zero));
+
+    if (!cellZoneIDs_.empty())
+    {
+        this->calcForce(U, rho, mu, tforce());
+    }
+
+    return tforce;
+}
+
+
 void Foam::porosityModel::porosityModel::addResistance
 (
     fvVectorMatrix& UEqn
 ) const
 {
-    if (cellZoneIds_.empty())
+    if (cellZoneIDs_.empty())
     {
         return;
     }
@@ -157,7 +176,7 @@ void Foam::porosityModel::porosityModel::addResistance
     const volScalarField& mu
 ) const
 {
-    if (cellZoneIds_.empty())
+    if (cellZoneIDs_.empty())
     {
         return;
     }
@@ -173,7 +192,7 @@ void Foam::porosityModel::porosityModel::addResistance
     bool correctAUprocBC
 ) const
 {
-    if (cellZoneIds_.empty())
+    if (cellZoneIDs_.empty())
     {
         return;
     }
@@ -190,12 +209,30 @@ void Foam::porosityModel::porosityModel::addResistance
 }
 
 
+bool Foam::porosityModel::porosityModel::movePoints()
+{
+    // no updates necessary; all member data independent of mesh
+    return true;
+}
+
+
+void Foam::porosityModel::porosityModel::updateMesh(const mapPolyMesh& mpm)
+{
+    // no updates necessary; all member data independent of mesh
+}
+
+
+bool Foam::porosityModel::writeData(Ostream& os) const
+{
+    return true;
+}
+
 bool Foam::porosityModel::read(const dictionary& dict)
 {
     active_ = readBool(dict.lookup("active"));
     coeffs_ = dict.subDict(type() + "Coeffs");
     dict.lookup("cellZone") >> zoneName_;
-    cellZoneIds_ = mesh_.cellZones().findIndices(zoneName_);
+    cellZoneIDs_ = mesh_.cellZones().findIndices(zoneName_);
 
     return true;
 }

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H

@@ -36,6 +36,7 @@ SourceFiles
 #ifndef porosityModel_H
 #define porosityModel_H
 
+#include "MeshObject.H"
 #include "fvMesh.H"
 #include "dictionary.H"
 #include "fvMatricesFwd.H"
@@ -54,6 +55,8 @@ namespace Foam
 \*---------------------------------------------------------------------------*/
 
 class porosityModel
+:
+    public MeshObject<fvMesh, UpdateableMeshObject, porosityModel>
 {
 private:
 
@@ -88,8 +91,8 @@ protected:
         //- Name(s) of cell-zone
         keyType zoneName_;
 
-        //- Cell zone Ids
+        //- Cell zone IDs
-        labelList cellZoneIds_;
+        labelList cellZoneIDs_;
 
         //- Local co-ordinate system
         coordinateSystem coordSys_;
@@ -100,6 +103,15 @@ protected:
         //- Adjust negative resistance values to be multiplier of max value
         void adjustNegativeResistance(dimensionedVector& resist);
 
+        //- Calculate the porosity force
+        virtual void calcForce
+        (
+            const volVectorField& U,
+            const volScalarField& rho,
+            const volScalarField& mu,
+            vectorField& force
+        ) const = 0;
+
         virtual void correct(fvVectorMatrix& UEqn) const = 0;
 
         virtual void correct
@@ -207,6 +219,17 @@ public:
         //- Return const access to the porosity active flag
         inline bool active() const;
 
+        //- Return const access to the cell zone IDs
+        inline const labelList& cellZoneIDs() const;
+
+        //- Return the force over the cell zone(s)
+        virtual tmp<vectorField> force
+        (
+            const volVectorField& U,
+            const volScalarField& rho,
+            const volScalarField& mu
+        ) const;
+
         //- Add resistance
         virtual void addResistance(fvVectorMatrix& UEqn) const;
 
@@ -227,10 +250,19 @@ public:
         ) const;
 
 
+    // Topology change
+
+        //- Move points
+        virtual bool movePoints();
+
+        //- Update on meshUpdate
+        virtual void updateMesh(const mapPolyMesh& mpm);
+
+
     // I-O
 
         //- Write
-        virtual void writeData(Ostream& os) const = 0;
+        virtual bool writeData(Ostream& os) const;
 
         //- Read porosity dictionary
         virtual bool read(const dictionary& dict);

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,4 +35,10 @@ inline bool Foam::porosityModel::active() const
 }
 
 
+inline const Foam::labelList& Foam::porosityModel::cellZoneIDs() const
+{
+    return cellZoneIDs_;
+}
+
+
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,6 +66,23 @@ Foam::porosityModels::powerLaw::~powerLaw()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::porosityModels::powerLaw::calcForce
+(
+    const volVectorField& U,
+    const volScalarField& rho,
+    const volScalarField& mu,
+    vectorField& force
+) const
+{
+    scalarField Udiag(U.size());
+    const scalarField& V = mesh_.V();
+
+    apply(Udiag, V, rho, U);
+
+    force = Udiag*U;
+}
+
+
 void Foam::porosityModels::powerLaw::correct
 (
     fvVectorMatrix& UEqn
@@ -126,10 +143,12 @@ void Foam::porosityModels::powerLaw::correct
 }
 
 
-void Foam::porosityModels::powerLaw::writeData(Ostream& os) const
+bool Foam::porosityModels::powerLaw::writeData(Ostream& os) const
 {
     os  << indent << name_ << endl;
     dict_.write(os);
+
+    return true;
 }
 
 

src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -127,6 +127,15 @@ public:
 
     // Member Functions
 
+        //- Calculate the porosity force
+        virtual void calcForce
+        (
+            const volVectorField& U,
+            const volScalarField& rho,
+            const volScalarField& mu,
+            vectorField& force
+        ) const;
+
         //- Add resistance
         virtual void correct(fvVectorMatrix& UEqn) const;
 
@@ -149,7 +158,7 @@ public:
     // I-O
 
         //- Write
-        void writeData(Ostream& os) const;
+        bool writeData(Ostream& os) const;
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLawTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,9 +37,9 @@ void Foam::porosityModels::powerLaw::apply
     const scalar C0 = C0_;
     const scalar C1m1b2 = (C1_ - 1.0)/2.0;
 
-    forAll(cellZoneIds_, zoneI)
+    forAll(cellZoneIDs_, zoneI)
     {
-        const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+        const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
 
         forAll(cells, i)
         {
@@ -63,9 +63,9 @@ void Foam::porosityModels::powerLaw::apply
     const scalar C0 = C0_;
     const scalar C1m1b2 = (C1_ - 1.0)/2.0;
 
-    forAll(cellZoneIds_, zoneI)
+    forAll(cellZoneIDs_, zoneI)
     {
-        const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+        const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
 
         forAll(cells, i)
         {

src/postProcessing/functionObjects/forces/forceCoeffs/forceCoeffs.C

@@ -120,8 +120,8 @@ void Foam::forceCoeffs::write()
 
             scalar pDyn = 0.5*rhoRef_*magUInf_*magUInf_;
 
-            Field<vector> totForce(force_[0] + force_[1]);
+            Field<vector> totForce(force_[0] + force_[1] + force_[2]);
-            Field<vector> totMoment(moment_[0] + moment_[1]);
+            Field<vector> totMoment(moment_[0] + moment_[1] + moment_[2]);
 
             List<Field<scalar> > coeffs(3);
             coeffs[0].setSize(nBin_);

src/postProcessing/functionObjects/forces/forces/forces.C

@@ -37,6 +37,7 @@ License
 #include "compressible/RAS/RASModel/RASModel.H"
 #include "compressible/LES/LESModel/LESModel.H"
 
+#include "porosityModel.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -52,14 +53,15 @@ void Foam::forces::writeFileHeader(const label i)
 {
     file()
         << "# Time" << tab
-        << "forces(pressure, viscous) moment(pressure, viscous)";
+        << "forces(pressure, viscous, porous) "
+        << "moment(pressure, viscous, porous)";
 
     if (localSystem_)
     {
         file()
             << tab
-            << "local forces(pressure, viscous) "
+            << "local forces(pressure, viscous, porous) "
-            << "local moment(pressure, viscous)";
+            << "local moment(pressure, viscous, porous)";
     }
 
     file()<< endl;
@@ -132,7 +134,7 @@ Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
     else
     {
         FatalErrorIn("forces::devRhoReff()")
-            << "No valid model for viscous stress calculation."
+            << "No valid model for viscous stress calculation"
             << exit(FatalError);
 
         return volSymmTensorField::null();
@@ -140,6 +142,46 @@ Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
 }
 
 
+Foam::tmp<Foam::volScalarField> Foam::forces::mu() const
+{
+    if (obr_.foundObject<fluidThermo>("thermophysicalProperties"))
+    {
+        const fluidThermo& thermo =
+             obr_.lookupObject<fluidThermo>("thermophysicalProperties");
+
+        return thermo.mu();
+    }
+    else if
+    (
+        obr_.foundObject<singlePhaseTransportModel>("transportProperties")
+    )
+    {
+        const singlePhaseTransportModel& laminarT =
+            obr_.lookupObject<singlePhaseTransportModel>
+            ("transportProperties");
+
+        return rho()*laminarT.nu();
+    }
+    else if (obr_.foundObject<dictionary>("transportProperties"))
+    {
+        const dictionary& transportProperties =
+             obr_.lookupObject<dictionary>("transportProperties");
+
+        dimensionedScalar nu(transportProperties.lookup("nu"));
+
+        return rho()*nu;
+    }
+    else
+    {
+        FatalErrorIn("forces::mu()")
+            << "No valid model for dynamic viscosity calculation"
+            << exit(FatalError);
+
+        return volScalarField::null();
+    }
+}
+
+
 Foam::tmp<Foam::volScalarField> Foam::forces::rho() const
 {
     if (rhoName_ == "rhoInf")
@@ -190,32 +232,35 @@ Foam::scalar Foam::forces::rho(const volScalarField& p) const
 
 void Foam::forces::applyBins
 (
-    const label patchI,
+    const vectorField& Md,
-    const vectorField fN,
+    const vectorField& fN,
-    const vectorField Md,
+    const vectorField& fT,
-    const vectorField fT
+    const vectorField& fP,
+    const vectorField& d
 )
 {
     if (nBin_ == 1)
     {
         force_[0][0] += sum(fN);
         force_[1][0] += sum(fT);
-        moment_[0][0] += sum(Md ^ fN);
+        force_[2][0] += sum(fP);
-        moment_[1][0] += sum(Md ^ fT);
+        moment_[0][0] += sum(Md^fN);
+        moment_[1][0] += sum(Md^fT);
+        moment_[2][0] += sum(Md^fP);
     }
     else
     {
-        const fvMesh& mesh = refCast<const fvMesh>(obr_);
+        scalarField dd((d & binDir_) - binMin_);
-        const vectorField& Cf = mesh.C().boundaryField()[patchI];
-        scalarField d((Cf & binDir_) - binMin_);
 
-        forAll(d, i)
+        forAll(dd, i)
         {
-            label binI = floor(d[i]/binDx_);
+            label binI = floor(dd[i]/binDx_);
             force_[0][binI] += fN[i];
             force_[1][binI] += fT[i];
+            force_[2][binI] += fP[i];
             moment_[0][binI] += Md[i]^fN[i];
             moment_[1][binI] += Md[i]^fT[i];
+            moment_[2][binI] += Md[i]^fP[i];
         }
     }
 }
@@ -239,7 +284,7 @@ void Foam::forces::writeBins() const
         Info<< "    Writing bins to " << forcesDir << endl;
     }
 
-    wordList fieldNames(IStringStream("(pressure viscous)")());
+    wordList fieldNames(IStringStream("(pressure viscous porous)")());
 
     List<Field<vector> > f(force_);
     List<Field<vector> > m(moment_);
@@ -250,8 +295,11 @@ void Foam::forces::writeBins() const
         {
             f[0][i] += f[0][i-1];
             f[1][i] += f[1][i-1];
+            f[2][i] += f[2][i-1];
+
             m[0][i] += m[0][i-1];
             m[1][i] += m[1][i-1];
+            m[2][i] += m[2][i-1];
         }
     }
 
@@ -264,12 +312,14 @@ void Foam::forces::writeBins() const
 
     if (localSystem_)
     {
-        List<Field<vector> > localForce(2);
+        List<Field<vector> > localForce(3);
-        List<Field<vector> > localMoment(2);
+        List<Field<vector> > localMoment(3);
         localForce[0] = coordSys_.localVector(force_[0]);
         localForce[1] = coordSys_.localVector(force_[1]);
+        localForce[2] = coordSys_.localVector(force_[2]);
         localMoment[0] = coordSys_.localVector(moment_[0]);
         localMoment[1] = coordSys_.localVector(moment_[1]);
+        localMoment[2] = coordSys_.localVector(moment_[2]);
 
         if (binCumulative_)
         {
@@ -277,8 +327,10 @@ void Foam::forces::writeBins() const
             {
                 localForce[0][i] += localForce[0][i-1];
                 localForce[1][i] += localForce[1][i-1];
+                localForce[2][i] += localForce[2][i-1];
                 localMoment[0][i] += localMoment[0][i-1];
                 localMoment[1][i] += localMoment[1][i-1];
+                localMoment[2][i] += localMoment[2][i-1];
             }
         }
 
@@ -306,8 +358,8 @@ Foam::forces::forces
     obr_(obr),
     active_(true),
     log_(false),
-    force_(2),
+    force_(3),
-    moment_(2),
+    moment_(3),
     patchSet_(),
     pName_(word::null),
     UName_(word::null),
@@ -318,6 +370,7 @@ Foam::forces::forces
     pRef_(0),
     coordSys_(),
     localSystem_(false),
+    porosity_(false),
     nBin_(1),
     binDir_(vector::zero),
     binDx_(0.0),
@@ -365,8 +418,8 @@ Foam::forces::forces
     obr_(obr),
     active_(true),
     log_(false),
-    force_(2),
+    force_(3),
-    moment_(2),
+    moment_(3),
     patchSet_(patchSet),
     pName_(pName),
     UName_(UName),
@@ -377,6 +430,7 @@ Foam::forces::forces
     pRef_(pRef),
     coordSys_(coordSys),
     localSystem_(false),
+    porosity_(false),
     nBin_(1),
     binDir_(vector::zero),
     binDx_(0.0),
@@ -475,6 +529,18 @@ void Foam::forces::read(const dictionary& dict)
             localSystem_ = true;
         }
 
+        if (dict.readIfPresent("porosity", porosity_))
+        {
+            if (porosity_)
+            {
+                Info<< "Including porosity effects" << endl;
+            }
+            else
+            {
+                Info<< "Not including porosity effects" << endl;
+            }
+        }
+
         if (dict.found("binData"))
         {
             const dictionary& binDict(dict.subDict("binData"));
@@ -491,10 +557,11 @@ void Foam::forces::read(const dictionary& dict)
             else if ((nBin_ == 0) || (nBin_ == 1))
             {
                 nBin_ = 1;
-                force_[0].setSize(1);
+                forAll(force_, i)
-                force_[1].setSize(1);
+                {
-                moment_[0].setSize(1);
+                    force_[i].setSize(1);
-                moment_[1].setSize(1);
+                    moment_[i].setSize(1);
+                }
             }
 
             if (nBin_ > 1)
@@ -546,8 +613,10 @@ void Foam::forces::read(const dictionary& dict)
             // allocate storage for forces and moments
             force_[0].setSize(1);
             force_[1].setSize(1);
+            force_[2].setSize(1);
             moment_[0].setSize(1);
             moment_[1].setSize(1);
+            moment_[2].setSize(1);
         }
     }
 }
@@ -581,28 +650,50 @@ void Foam::forces::write()
         if (log_)
         {
             Info<< type() << " output:" << nl
-                << "    forces(pressure,viscous)"
+                << "    forces(pressure,viscous,porous) = ("
-                << "(" << sum(force_[0]) << "," << sum(force_[1]) << ")" << nl
+                << sum(force_[0]) << ","
-                << "    moment(pressure,viscous)"
+                << sum(force_[1]) << ","
-                << "(" << sum(moment_[0]) << "," << sum(moment_[1]) << ")"
+                << sum(force_[2]) << ")" << nl
+                << "    moment(pressure,viscous,porous) = ("
+                << sum(moment_[0]) << ","
+                << sum(moment_[1]) << ","
+                << sum(moment_[2]) << ")"
                 << nl;
         }
 
         file() << obr_.time().value() << tab
-            << "(" << sum(force_[0]) << "," << sum(force_[1]) << ") "
+            << "("
-            << "(" << sum(moment_[0]) << "," << sum(moment_[1]) << ")"
+            << sum(force_[0]) << ","
+            << sum(force_[1]) << ","
+            << sum(force_[2])
+            << ") "
+            << "("
+            << sum(moment_[0]) << ","
+            << sum(moment_[1]) << ","
+            << sum(moment_[2])
+            << ")"
             << endl;
 
         if (localSystem_)
         {
-            vectorField localForceP(coordSys_.localVector(force_[0]));
+            vectorField localForceN(coordSys_.localVector(force_[0]));
-            vectorField localForceV(coordSys_.localVector(force_[1]));
+            vectorField localForceT(coordSys_.localVector(force_[1]));
-            vectorField localMomentP(coordSys_.localVector(moment_[0]));
+            vectorField localForceP(coordSys_.localVector(force_[2]));
-            vectorField localMomentV(coordSys_.localVector(moment_[1]));
+            vectorField localMomentN(coordSys_.localVector(moment_[0]));
+            vectorField localMomentT(coordSys_.localVector(moment_[1]));
+            vectorField localMomentP(coordSys_.localVector(moment_[2]));
 
             file() << obr_.time().value() << tab
-                << "(" << sum(localForceP) << "," << sum(localForceV) << ") "
+                << "("
-                << "(" << sum(localMomentP) << "," << sum(localMomentV) << ")"
+                << sum(localForceN) << ","
+                << sum(localForceT) << ","
+                << sum(localForceP)
+                << ") "
+                << "("
+                << sum(localMomentN) << ","
+                << sum(localMomentT) << ","
+                << sum(localMomentP)
+                << ")"
                 << endl;
         }
 
@@ -620,9 +711,11 @@ void Foam::forces::calcForcesMoment()
 {
     force_[0] = vector::zero;
     force_[1] = vector::zero;
+    force_[2] = vector::zero;
 
     moment_[0] = vector::zero;
     moment_[1] = vector::zero;
+    moment_[2] = vector::zero;
 
     if (directForceDensity_)
     {
@@ -656,7 +749,10 @@ void Foam::forces::calcForcesMoment()
             // Tangential force (total force minus normal fN)
             vectorField fT(sA*fD.boundaryField()[patchI] - fN);
 
-            applyBins(patchI, fN, Md, fT);
+            //- Porous force
+            vectorField fP(Md.size(), vector::zero);
+
+            applyBins(Md, fN, fT, fP, mesh.C().boundaryField()[patchI]);
         }
     }
     else
@@ -685,14 +781,51 @@ void Foam::forces::calcForcesMoment()
                 mesh.C().boundaryField()[patchI] - coordSys_.origin()
             );
 
-            vectorField pf
+            vectorField fN
             (
                 rho(p)*Sfb[patchI]*(p.boundaryField()[patchI] - pRef)
             );
 
-            vectorField vf(Sfb[patchI] & devRhoReffb[patchI]);
+            vectorField fT(Sfb[patchI] & devRhoReffb[patchI]);
 
-            applyBins(patchI, pf, Md, vf);
+            vectorField fP(Md.size(), vector::zero);
+
+            applyBins(Md, fN, fT, fP, mesh.C().boundaryField()[patchI]);
+        }
+    }
+
+    if (porosity_)
+    {
+        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
+        const volScalarField rho(this->rho());
+        const volScalarField mu(this->mu());
+
+        const fvMesh& mesh = U.mesh();
+
+        const HashTable<const porosityModel*> models =
+            obr_.lookupClass<porosityModel>();
+
+        forAllConstIter(HashTable<const porosityModel*>, models, iter)
+        {
+            const porosityModel& pm = *iter();
+
+            vectorField fPTot(pm.force(U, rho, mu));
+
+            const labelList& cellZoneIDs = pm.cellZoneIDs();
+
+            forAll(cellZoneIDs, i)
+            {
+                label zoneI = cellZoneIDs[i];
+                const cellZone& cZone = mesh.cellZones()[zoneI];
+
+                const vectorField d(mesh.C(), cZone);
+                const vectorField fP(fPTot, cZone);
+                const vectorField Md(d - coordSys_.origin());
+
+                const vectorField fDummy(Md.size(), vector::zero);
+
+                applyBins(Md, fDummy, fDummy, fP, d);
+            }
         }
     }
 
@@ -703,13 +836,13 @@ void Foam::forces::calcForcesMoment()
 
 Foam::vector Foam::forces::forceEff() const
 {
-    return sum(force_[0]) + sum(force_[1]);
+    return sum(force_[0]) + sum(force_[1]) + sum(force_[2]);
 }
 
 
 Foam::vector Foam::forces::momentEff() const
 {
-    return sum(moment_[0]) + sum(moment_[1]);
+    return sum(moment_[0]) + sum(moment_[1]) + sum(moment_[2]);
 }
 
 

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

@@ -149,10 +149,10 @@ protected:
         //- Switch to send output to Info as well as to file
         Switch log_;
 
-        //- Pressure and viscous force per bin
+        //- Pressure, viscous and porous force per bin
         List<Field<vector> > force_;
 
-        //- Pressure and viscous pressure per bin
+        //- Pressure, viscous and porous moment per bin
         List<Field<vector> > moment_;
 
 
@@ -188,6 +188,9 @@ protected:
             //- Flag to indicate whether we are using a local co-ordinate sys
             bool localSystem_;
 
+            //- Flag to include porosity effects
+            bool porosity_;
+
 
             // Bin information
 
@@ -221,6 +224,9 @@ protected:
         //- Return the effective viscous stress (laminar + turbulent).
         tmp<volSymmTensorField> devRhoReff() const;
 
+        //- Dynamic viscosity field
+        tmp<volScalarField> mu() const;
+
         //- Return rho if rhoName is specified otherwise rhoRef
         tmp<volScalarField> rho() const;
 
@@ -231,10 +237,11 @@ protected:
         //- Accumulate bin data
         void applyBins
         (
-            const label patchI,
+            const vectorField& Md,
-            const vectorField fN,
+            const vectorField& fN,
-            const vectorField Md,
+            const vectorField& fT,
-            const vectorField fT
+            const vectorField& fP,
+            const vectorField& d
         );
 
         //- Helper function to write bin data