ENH: Adding tensor fields option for porosity models DarcyForchheimer

and fixedCoeff

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

@@ -52,66 +52,65 @@ Foam::porosityModels::DarcyForchheimer::DarcyForchheimer
 )
 :
     porosityModel(name, modelType, mesh, dict, cellZoneName),
-    coordSys_(coeffs_, mesh),
-    D_("D", dimless/sqr(dimLength), tensor::zero),
-    F_("F", dimless/dimLength, tensor::zero),
+    D_(cellZoneIds_.size()),
+    F_(cellZoneIds_.size()),
     rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho")),
     muName_(coeffs_.lookupOrDefault<word>("mu", "thermo:mu")),
     nuName_(coeffs_.lookupOrDefault<word>("nu", "nu"))
 {
-    // local-to-global transformation tensor
-    const tensor& E = coordSys_.R();
 
     dimensionedVector d(coeffs_.lookup("d"));
-    if (D_.dimensions() != d.dimensions())
-    {
-        FatalIOErrorIn
-        (
-            "Foam::porosityModels::DarcyForchheimer::DarcyForchheimer"
-            "("
-                "const word&, "
-                "const word&, "
-                "const fvMesh&, "
-                "const dictionary&"
-            ")",
-            coeffs_
-        )   << "incorrect dimensions for d: " << d.dimensions()
-            << " should be " << D_.dimensions()
-            << exit(FatalIOError);
-    }
+    dimensionedVector f(coeffs_.lookup("f"));
 
     adjustNegativeResistance(d);
-
-    D_.value().xx() = d.value().x();
-    D_.value().yy() = d.value().y();
-    D_.value().zz() = d.value().z();
-    D_.value() = (E & D_ & E.T()).value();
-
-    dimensionedVector f(coeffs_.lookup("f"));
-    if (F_.dimensions() != f.dimensions())
-    {
-        FatalIOErrorIn
-        (
-            "Foam::porosityModels::DarcyForchheimer::DarcyForchheimer"
-            "("
-                "const word&, "
-                "const word&, "
-                "const fvMesh&, "
-                "const dictionary&"
-            ")",
-            coeffs_
-        )   << "incorrect dimensions for f: " << f.dimensions()
-            << " should be " << F_.dimensions()
-            << exit(FatalIOError);
-    }
-
     adjustNegativeResistance(f);
 
-    // leading 0.5 is from 1/2*rho
-    F_.value().xx() = 0.5*f.value().x();
-    F_.value().yy() = 0.5*f.value().y();
-    F_.value().zz() = 0.5*f.value().z();
-    F_.value() = (E & F_ & E.T()).value();
+    if (coordSys_.R().uniform())
+    {
+        forAll (cellZoneIds_, zoneI)
+        {
+            D_[zoneI].setSize(1, tensor::zero);
+            F_[zoneI].setSize(1, tensor::zero);
+
+            D_[zoneI][0].xx() = d.value().x();
+            D_[zoneI][0].yy() = d.value().y();
+            D_[zoneI][0].zz() = d.value().z();
+
+            D_[zoneI][0] = coordSys_.R().transformTensor(D_[zoneI][0]);
+
+            // leading 0.5 is from 1/2*rho
+            F_[zoneI][0].xx() = 0.5*f.value().x();
+            F_[zoneI][0].yy() = 0.5*f.value().y();
+            F_[zoneI][0].zz() = 0.5*f.value().z();
+
+            F_[zoneI][0] = coordSys_.R().transformTensor(F_[zoneI][0]);
+        }
+
+    }
+    else
+    {
+        forAll(cellZoneIds_, zoneI)
+        {
+            const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+
+            D_[zoneI].setSize(cells.size(), tensor::zero);
+            F_[zoneI].setSize(cells.size(), tensor::zero);
+
+            forAll(cells, i)
+            {
+                D_[zoneI][i].xx() = d.value().x();
+                D_[zoneI][i].yy() = d.value().y();
+                D_[zoneI][i].zz() = d.value().z();
+
+                F_[zoneI][i].xx() = f.value().x();
+                F_[zoneI][i].yy() = f.value().y();
+                F_[zoneI][i].zz() = f.value().z();
+            }
+
+            D_[zoneI] = coordSys_.R().transformTensor(D_[zoneI], cells);
+            F_[zoneI] = coordSys_.R().transformTensor(F_[zoneI], cells);
+        }
+    }
 }
 
 
@@ -132,7 +131,7 @@ void Foam::porosityModels::DarcyForchheimer::correct
     const scalarField& V = mesh_.V();
     scalarField& Udiag = UEqn.diag();
     vectorField& Usource = UEqn.source();
- 
+
     if (UEqn.dimensions() == dimForce)
     {
         const volScalarField& rho =
@@ -163,7 +162,7 @@ void Foam::porosityModels::DarcyForchheimer::correct
     const scalarField& V = mesh_.V();
     scalarField& Udiag = UEqn.diag();
     vectorField& Usource = UEqn.source();
- 
+
     apply(Udiag, Usource, V, rho, mu, U);
 }
 

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

@@ -53,7 +53,6 @@ SourceFiles
 #define DarcyForchheimer_H
 
 #include "porosityModel.H"
-#include "coordinateSystem.H"
 #include "dimensionedTensor.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -75,14 +74,12 @@ private:
 
     // Private data
 
-        //- Local co-ordinate system
-        coordinateSystem coordSys_;
 
         //- Darcy coefficient [1/m2]
-        dimensionedTensor D_;
+        List<tensorField> D_;
 
         //- Forchheimer coefficient [1/m]
-        dimensionedTensor F_;
+        List<tensorField> F_;
 
         //- Name of density field
         word rhoName_;
@@ -161,7 +158,7 @@ public:
         virtual void correct
         (
             const fvVectorMatrix& UEqn,
-            volTensorField& AU            
+            volTensorField& AU
         ) const;
 
 

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

@@ -36,18 +36,19 @@ void Foam::porosityModels::DarcyForchheimer::apply
     const vectorField& U
 ) const
 {
-    const tensor& D = D_.value();
-    const tensor& F = F_.value();
-
     forAll(cellZoneIds_, zoneI)
     {
+        const tensorField& dZones = D_[zoneI];
+        const tensorField& fZones = F_[zoneI];
+
         const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
 
         forAll(cells, i)
         {
             const label cellI = cells[i];
-
-            const tensor Cd = mu[cellI]*D + (rho[cellI]*mag(U[cellI]))*F;
+            const label j = this->fieldIndex(i);
+            const tensor Cd =
+                mu[cellI]*dZones[j] + (rho[cellI]*mag(U[cellI]))*fZones[j];
 
             const scalar isoCd = tr(Cd);
 
@@ -67,16 +68,20 @@ void Foam::porosityModels::DarcyForchheimer::apply
     const vectorField& U
 ) const
 {
-    const tensor& D = D_.value();
-    const tensor& F = F_.value();
-
     forAll(cellZoneIds_, zoneI)
     {
+        const tensorField& dZones = D_[zoneI];
+        const tensorField& fZones = F_[zoneI];
+
         const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
 
         forAll(cells, i)
         {
             const label cellI = cells[i];
+            const label j = this->fieldIndex(i);
+            const tensor D = dZones[j];
+            const tensor F = fZones[j];
+
             AU[cellI] += mu[cellI]*D + (rho[cellI]*mag(U[cellI]))*F;
         }
     }

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

@@ -50,19 +50,18 @@ void Foam::porosityModels::fixedCoeff::apply
     const scalar rho
 ) const
 {
-    const tensor& alpha = alpha_.value();
-    const tensor& beta = beta_.value();
-
     forAll(cellZoneIds_, zoneI)
     {
+        const tensorField& alphaZones = alpha_[zoneI];
+        const tensorField& betaZones = beta_[zoneI];
+
         const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
 
         forAll(cells, i)
         {
             const label cellI = cells[i];
-
-            const tensor Cd = rho*(alpha + beta*mag(U[cellI]));
-
+            const label j = fieldIndex(i);
+            const tensor Cd = rho*(alphaZones[j] + betaZones[j]*mag(U[cellI]));
             const scalar isoCd = tr(Cd);
 
             Udiag[cellI] += V[cellI]*isoCd;
@@ -79,16 +78,21 @@ void Foam::porosityModels::fixedCoeff::apply
     const scalar rho
 ) const
 {
-    const tensor& alpha = alpha_.value();
-    const tensor& beta = beta_.value();
 
     forAll(cellZoneIds_, zoneI)
     {
+        const tensorField& alphaZones = alpha_[zoneI];
+        const tensorField& betaZones = beta_[zoneI];
+
         const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
 
         forAll(cells, i)
         {
             const label cellI = cells[i];
+            const label j = fieldIndex(i);
+            const tensor alpha = alphaZones[j];
+            const tensor beta = betaZones[j];
+
             AU[cellI] += rho*(alpha + beta*mag(U[cellI]));
         }
     }
@@ -107,62 +111,59 @@ Foam::porosityModels::fixedCoeff::fixedCoeff
 )
 :
     porosityModel(name, modelType, mesh, dict, cellZoneName),
-    coordSys_(coeffs_, mesh),
-    alpha_("alpha", dimless/dimTime, tensor::zero),
-    beta_("beta", dimless/dimLength, tensor::zero)
+    alpha_(cellZoneIds_.size()),
+    beta_(cellZoneIds_.size())
 {
-    // local-to-global transformation tensor
-    const tensor& E = coordSys_.R();
-
     dimensionedVector alpha(coeffs_.lookup("alpha"));
-    if (alpha_.dimensions() != alpha.dimensions())
-    {
-        FatalIOErrorIn
-        (
-            "Foam::porosityModels::fixedCoeff::fixedCoeff"
-            "("
-                "const word&, "
-                "const word&, "
-                "const fvMesh&, "
-                "const dictionary&"
-            ")",
-            coeffs_
-        )   << "incorrect dimensions for alpha: " << alpha.dimensions()
-            << " should be " << alpha_.dimensions()
-            << exit(FatalIOError);
-    }
+    dimensionedVector beta(coeffs_.lookup("beta"));
 
     adjustNegativeResistance(alpha);
-
-    alpha_.value().xx() = alpha.value().x();
-    alpha_.value().yy() = alpha.value().y();
-    alpha_.value().zz() = alpha.value().z();
-    alpha_.value() = (E & alpha_ & E.T()).value();
-
-    dimensionedVector beta(coeffs_.lookup("beta"));
-    if (beta_.dimensions() != beta.dimensions())
-    {
-        FatalIOErrorIn
-        (
-            "Foam::porosityModels::fixedCoeff::fixedCoeff"
-            "("
-                "const word&, "
-                "const word&, "
-                "const fvMesh&, "
-                "const dictionary&"
-            ")",
-            coeffs_
-        )   << "incorrect dimensions for beta: " << beta.dimensions()
-            << " should be " << beta_.dimensions()
-            << exit(FatalIOError);
-    }
-
     adjustNegativeResistance(beta);
 
-    beta_.value().xx() = beta.value().x();
-    beta_.value().yy() = beta.value().y();
-    beta_.value().zz() = beta.value().z();
-    beta_.value() = (E & beta_ & E.T()).value();
+    if (coordSys_.R().uniform())
+    {
+        forAll (cellZoneIds_, zoneI)
+        {
+            alpha_[zoneI].setSize(1, tensor::zero);
+            beta_[zoneI].setSize(1, tensor::zero);
+
+            alpha_[zoneI][0].xx() = alpha.value().x();
+            alpha_[zoneI][0].yy() = alpha.value().y();
+            alpha_[zoneI][0].zz() = alpha.value().z();
+            alpha_[zoneI][0] = coordSys_.R().transformTensor(alpha_[zoneI][0]);
+
+            beta_[zoneI][0].xx() = beta.value().x();
+            beta_[zoneI][0].yy() = beta.value().y();
+            beta_[zoneI][0].zz() = beta.value().z();
+            beta_[zoneI][0] = coordSys_.R().transformTensor(beta_[zoneI][0]);
+        }
+    }
+    else
+    {
+        forAll(cellZoneIds_, zoneI)
+        {
+            const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+
+            alpha_[zoneI].setSize(cells.size(), tensor::zero);
+            beta_[zoneI].setSize(cells.size(), tensor::zero);
+
+            forAll(cells, i)
+            {
+                alpha_[zoneI][i].xx() = alpha.value().x();
+                alpha_[zoneI][i].yy() = alpha.value().y();
+                alpha_[zoneI][i].zz() = alpha.value().z();
+
+                beta_[zoneI][i].xx() = beta.value().x();
+                beta_[zoneI][i].yy() = beta.value().y();
+                beta_[zoneI][i].zz() = beta.value().z();
+            }
+
+            alpha_[zoneI] =
+                coordSys_.R().transformTensor(alpha_[zoneI], cells);
+
+            beta_[zoneI] = coordSys_.R().transformTensor(beta_[zoneI], cells);
+        }
+    }
 }
 
 
@@ -183,7 +184,7 @@ void Foam::porosityModels::fixedCoeff::correct
     const scalarField& V = mesh_.V();
     scalarField& Udiag = UEqn.diag();
     vectorField& Usource = UEqn.source();
- 
+
     scalar rho = 1.0;
     if (UEqn.dimensions() == dimForce)
     {
@@ -211,7 +212,7 @@ void Foam::porosityModels::fixedCoeff::correct
     {
         coeffs_.lookup("rhoRef") >> rho;
     }
- 
+
     apply(Udiag, Usource, V, U, rho);
 }
 

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

@@ -43,7 +43,6 @@ SourceFiles
 #define fixedCoeff_H
 
 #include "porosityModel.H"
-#include "coordinateSystem.H"
 #include "dimensionedTensor.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -65,14 +64,11 @@ private:
 
     // Private data
 
-        //- Local co-ordinate system
-        coordinateSystem coordSys_;
-
         //- Model alpha coefficient [1/s]
-        dimensionedTensor alpha_;
+        List<tensorField> alpha_;
 
         //- Model beta coefficient [1/m]
-        dimensionedTensor beta_;
+        List<tensorField> beta_;
 
 
     // Private Member Functions
@@ -138,7 +134,7 @@ public:
         virtual void correct
         (
             const fvVectorMatrix& UEqn,
-            volTensorField& AU            
+            volTensorField& AU
         ) const;
 
 

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

@@ -66,6 +66,17 @@ void Foam::porosityModel::adjustNegativeResistance(dimensionedVector& resist)
 }
 
 
+Foam::label Foam::porosityModel::fieldIndex(const label i) const
+{
+    label index = 0;
+    if (!coordSys_.R().uniform())
+    {
+        index = i;
+    }
+    return index;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::porosityModel::porosityModel
@@ -83,7 +94,8 @@ Foam::porosityModel::porosityModel
     coeffs_(dict.subDict(modelType + "Coeffs")),
     active_(true),
     zoneName_(cellZoneName),
-    cellZoneIds_()
+    cellZoneIds_(),
+    coordSys_(coordinateSystem::New(mesh, coeffs_))
 {
     if (zoneName_ == word::null)
     {
@@ -158,7 +170,7 @@ void Foam::porosityModel::porosityModel::addResistance
 (
     const fvVectorMatrix& UEqn,
     volTensorField& AU,
-    bool correctAUprocBC         
+    bool correctAUprocBC
 ) const
 {
     if (cellZoneIds_.empty())

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

@@ -40,6 +40,7 @@ SourceFiles
 #include "dictionary.H"
 #include "fvMatricesFwd.H"
 #include "runTimeSelectionTables.H"
+#include "coordinateSystem.H"
 #include "dimensionedVector.H"
 #include "keyType.H"
 
@@ -90,6 +91,9 @@ protected:
         //- Cell zone Ids
         labelList cellZoneIds_;
 
+        //- Local co-ordinate system
+        coordinateSystem coordSys_;
+
 
     // Protected Member Functions
 
@@ -111,6 +115,9 @@ protected:
             volTensorField& AU
         ) const = 0;
 
+        //- Return label index
+        label fieldIndex(const label index) const;
+
 
 public: