ENH: Updated porosity models

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)
         {