From 7e735715a4189e02ea322cf507892f623474c018 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Tue, 3 Mar 2009 17:09:01 +0000
Subject: [PATCH] more flexible mrfzones

	modified:   finiteVolume/cfdTools/general/MRF/MRFZone.C
modified:   finiteVolume/cfdTools/general/MRF/MRFZone.H
	modified:   ../tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/constant/MRFZones
---
 .../cfdTools/general/MRF/MRFZone.C            | 441 +++++++++++++-----
 .../cfdTools/general/MRF/MRFZone.H            |  36 +-
 .../mixerVessel2D/constant/MRFZones           |   4 +-
 3 files changed, 348 insertions(+), 133 deletions(-)

diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
index 6250de7e24..8b689d0d66 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.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) 1991-2008 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,8 +29,264 @@ License
 #include "volFields.H"
 #include "surfaceFields.H"
 #include "fvMatrices.H"
-#include "PackedList.H"
 #include "syncTools.H"
+#include "faceSet.H"
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::MRFZone::setMRFFaces
+(
+    labelList& faceType,
+    const labelList& excludedPatchIDs
+)
+{
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    // Knock out coupled patches
+    forAll(patches, patchi)
+    {
+        if (patches[patchi].coupled())
+        {
+            const polyPatch& pp = patches[patchi];
+
+            forAll(pp, j)
+            {
+                label faceI = pp.start()+j;
+
+                if (faceType[faceI] == 1)
+                {
+                    faceType[faceI] = 2;
+                }
+            }
+        }
+    }
+
+    // All explicitly provided exclusions
+    forAll(excludedPatchIDs, i)
+    {
+        const polyPatch& pp = patches[excludedPatchIDs[i]];
+
+        forAll(pp, j)
+        {
+            label faceI = pp.start()+j;
+
+            if (faceType[faceI] == 1)
+            {
+                faceType[faceI] = 2;
+            }
+        }
+    }
+
+    // Collect into lists per patch.
+    internalFaces_.setSize(mesh_.nFaces());
+    label nInternal = 0;
+
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        if (faceType[faceI] == 1)
+        {
+            internalFaces_[nInternal++] = faceI;
+        }
+    }
+    internalFaces_.setSize(nInternal);
+
+    labelList nIncludedFaces(patches.size(), 0);
+    labelList nExcludedFaces(patches.size(), 0);
+
+    forAll(patches, patchi)
+    {
+        const polyPatch& pp = patches[patchi];
+
+        forAll(pp, patchFacei)
+        {
+            label faceI = pp.start()+patchFacei;
+
+            if (faceType[faceI] == 1)
+            {
+                nIncludedFaces[patchi]++;
+            }
+            else if (faceType[faceI] == 2)
+            {
+                nExcludedFaces[patchi]++;
+            }
+        }
+    }
+
+    includedFaces_.setSize(patches.size());
+    excludedFaces_.setSize(patches.size());
+    forAll(nIncludedFaces, patchi)
+    {
+        includedFaces_[patchi].setSize(nIncludedFaces[patchi]);
+        excludedFaces_[patchi].setSize(nExcludedFaces[patchi]);
+    }
+    nIncludedFaces = 0;
+    nExcludedFaces = 0;
+
+    forAll(patches, patchi)
+    {
+        const polyPatch& pp = patches[patchi];
+
+        forAll(pp, patchFacei)
+        {
+            label faceI = pp.start()+patchFacei;
+
+            if (faceType[faceI] == 1)
+            {
+                includedFaces_[patchi][nIncludedFaces[patchi]++] = patchFacei;
+            }
+            else if (faceType[faceI] == 2)
+            {
+                excludedFaces_[patchi][nExcludedFaces[patchi]++] = patchFacei;
+            }
+        }
+    }
+
+    //if (debug)
+    //{
+    //    faceSet internalFaces(mesh_, "internalFaces", internalFaces_);
+    //    Pout<< "Writing internal faces in MRF zone to faceSet "
+    //        << internalFaces.name() << endl;
+    //    internalFaces.write();
+    //}
+    //{
+    //    faceSet MRFFaces(mesh_, "includedFaces", 100);
+    //    forAll(includedFaces_, patchi)
+    //    {
+    //        forAll(includedFaces_[patchi], i)
+    //        {
+    //            label patchFacei = includedFaces_[patchi][i];
+    //            MRFFaces.insert(patches[patchi].start()+patchFacei);
+    //        }
+    //    }
+    //    Pout<< "Writing patch faces in MRF zone to faceSet "
+    //        << MRFFaces.name() << endl;
+    //    MRFFaces.write();
+    //}
+    //{
+    //    faceSet excludedFaces(mesh_, "excludedFaces", 100);
+    //    forAll(excludedFaces_, patchi)
+    //    {
+    //        forAll(excludedFaces_[patchi], i)
+    //        {
+    //            label patchFacei = excludedFaces_[patchi][i];
+    //            excludedFaces.insert(patches[patchi].start()+patchFacei);
+    //        }
+    //    }
+    //    Pout<< "Writing faces in MRF zone with special handling to faceSet "
+    //        << excludedFaces.name() << endl;
+    //    excludedFaces.write();
+    //}
+}
+
+
+void Foam::MRFZone::setMRFFaces()
+{
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    // Type per face:
+    //  0:not in zone
+    //  1:moving with frame
+    //  2:other
+    labelList faceType(mesh_.nFaces(), 0);
+
+    bool faceZoneFound = (faceZoneID_ != -1);
+    reduce(faceZoneFound, orOp<bool>());
+
+    if (faceZoneFound)
+    {
+        // Explicitly provided faces.
+        if (faceZoneID_ != -1)
+        {
+            const labelList& zoneFaces = mesh_.faceZones()[faceZoneID_];
+
+            forAll(zoneFaces, i)
+            {
+                faceType[zoneFaces[i]] = 1;
+            }
+
+            if (allPatchesMove_)
+            {
+                // Explicitly provided patches that do not move
+                setMRFFaces(faceType, patchLabels_);
+            }
+            else
+            {
+                setMRFFaces(faceType, labelList(0));
+            }
+        }
+    }
+    else
+    {
+        // Determine faces in cell zone
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // (without constructing cells)
+
+        const labelList& own = mesh_.faceOwner();
+        const labelList& nei = mesh_.faceNeighbour();
+
+        // Cells in zone
+        boolList zoneCell(mesh_.nCells(), false);
+
+        if (cellZoneID_ != -1)
+        {
+            const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
+            forAll(cellLabels, i)
+            {
+                zoneCell[cellLabels[i]] = true;
+            }
+        }
+
+
+        label nZoneFaces = 0;
+
+        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+        {
+            if (zoneCell[own[faceI]] || zoneCell[nei[faceI]])
+            {
+                faceType[faceI] = 1;
+                nZoneFaces++;
+            }
+        }
+        forAll(patches, patchI)
+        {
+            const polyPatch& pp = patches[patchI];
+
+            if (!isA<emptyPolyPatch>(pp))
+            {
+                forAll(pp, i)
+                {
+                    label faceI = pp.start()+i;
+
+                    if (zoneCell[own[faceI]])
+                    {
+                        faceType[faceI] = 1;
+                        nZoneFaces++;
+                    }
+                }
+            }
+        }
+        syncTools::syncFaceList(mesh_, faceType, maxEqOp<label>(), false);
+
+
+        Info<< nl
+            << "MRFZone " << name_ << " : did not find a faceZone; using "
+            << returnReduce(nZoneFaces, sumOp<label>())
+            << " faces from the cellZone instead." << endl;
+
+
+        if (allPatchesMove_)
+        {
+            // Explicitly provided excluded patches
+            setMRFFaces(faceType, patchLabels_);
+        }
+        else
+        {
+            setMRFFaces(faceType, labelList(0));
+        }
+    }    
+}
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -41,8 +297,13 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
     dict_(is),
     cellZoneID_(mesh_.cellZones().findZoneID(name_)),
     faceZoneID_(mesh_.faceZones().findZoneID(name_)),
-    outsideFaces_(0),
-    patchNames_(dict_.lookup("patches")),
+    allPatchesMove_(dict_.found("nonRotatingPatches")),
+    patchNames_
+    (
+        allPatchesMove_
+      ? dict_.lookup("nonRotatingPatches")
+      : dict_.lookup("patches")
+    ),
     origin_(dict_.lookup("origin")),
     axis_(dict_.lookup("axis")),
     omega_(dict_.lookup("omega")),
@@ -53,100 +314,6 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
     axis_ = axis_/mag(axis_);
     Omega_ = omega_*axis_;
 
-    bool cellZoneFound = (cellZoneID_ != -1);
-    reduce(cellZoneFound, orOp<bool>());
-
-    if (!cellZoneFound)
-    {
-        FatalErrorIn
-        (
-            "Foam::MRFZone::MRFZone(const fvMesh& , const dictionary&)"
-        )   << "cannot find MRF cellZone " << name_
-            << exit(FatalError);
-    }
-
-    bool faceZoneFound = (faceZoneID_ != -1);
-    reduce(faceZoneFound, orOp<bool>());
-
-    if (!faceZoneFound)
-    {
-        // Determine faces in cell zone
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        // (does not construct cells)
-
-        const labelList& own = mesh_.faceOwner();
-        const labelList& nei = mesh_.faceNeighbour();
-
-        // Cells in zone
-        PackedBoolList zoneCell(mesh_.nCells());
-
-        if (cellZoneID_ != -1)
-        {
-            const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
-            forAll(cellLabels, i)
-            {
-                zoneCell[cellLabels[i]] = 1u;
-            }
-        }
-
-
-        // Faces in zone
-        PackedBoolList zoneFacesSet(mesh_.nFaces());
-
-        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-        {
-            if
-            (
-                zoneCell.get(own[faceI]) == 1u
-             || zoneCell.get(nei[faceI]) == 1u
-            )
-            {
-                zoneFacesSet[faceI] = 1u;
-            }
-        }
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-
-            if (pp.coupled())
-            {
-                forAll(pp, i)
-                {
-                    label faceI = pp.start()+i;
-
-                    if (zoneCell.get(own[faceI]) == 1u)
-                    {
-                        zoneFacesSet[faceI] = 1u;
-                    }
-                }
-            }
-        }
-        syncTools::syncFaceList(mesh_, zoneFacesSet, orEqOp<unsigned int>());
-
-
-        // Transfer to labelList
-        label n = zoneFacesSet.count();
-        outsideFaces_.setSize(n);
-        n = 0;
-        forAll(zoneFacesSet, faceI)
-        {
-            if (zoneFacesSet.get(faceI) == 1u)
-            {
-                outsideFaces_[n++] = faceI;
-            }
-        }
-
-        Info<< nl
-            << "MRFZone " << name_ << " : did not find a faceZone; using "
-            << returnReduce(outsideFaces_.size(), sumOp<label>())
-            << " faces internal to cellZone instead." << endl;
-
-
-        // Flag use of outsideFaces
-        faceZoneID_ = -2;
-    }
-
-
     patchLabels_.setSize(patchNames_.size());
 
     forAll(patchNames_, i)
@@ -157,11 +324,26 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
         {
             FatalErrorIn
             (
-                "Foam::MRFZone::MRFZone(const fvMesh& , const dictionary&)"
+                "Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
             )   << "cannot find MRF patch " << patchNames_[i]
                 << exit(FatalError);
         }
     }
+
+
+    bool cellZoneFound = (cellZoneID_ != -1);
+    reduce(cellZoneFound, orOp<bool>());
+
+    if (!cellZoneFound)
+    {
+        FatalErrorIn
+        (
+            "Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
+        )   << "cannot find MRF cellZone " << name_
+            << exit(FatalError);
+    }
+
+    setMRFFaces();
 }
 
 
@@ -189,48 +371,42 @@ void Foam::MRFZone::addCoriolis(fvVectorMatrix& UEqn) const
 
 void Foam::MRFZone::relativeFlux(surfaceScalarField& phi) const
 {
-    if (faceZoneID_ == -1)
-    {
-        return;
-    }
-
     const surfaceVectorField& Cf = mesh_.Cf();
     const surfaceVectorField& Sf = mesh_.Sf();
 
     const vector& origin = origin_.value();
     const vector& Omega = Omega_.value();
 
-    // Use either zone faces or outsideFaces_
-    const labelList& faces =
-    (
-        faceZoneID_ == -2
-      ? outsideFaces_
-      : mesh_.faceZones()[faceZoneID_]
-    );
-
-    forAll(faces, i)
+    // Internal faces
+    forAll(internalFaces_, i)
     {
-        label facei = faces[i];
+        label facei = internalFaces_[i];
+        phi[facei] -= (Omega ^ (Cf[facei] - origin)) & Sf[facei];
+    }
 
-        if (facei < mesh_.nInternalFaces())
+    // Included patches
+    forAll(includedFaces_, patchi)
+    {
+        forAll(includedFaces_[patchi], i)
         {
-            phi[facei] -= (Omega ^ (Cf[facei] - origin)) & Sf[facei];
+            label patchFacei = includedFaces_[patchi][i];
+
+            phi.boundaryField()[patchi][patchFacei] = 0.0;
         }
-        else
+    }
+
+    // Excluded patches
+    forAll(excludedFaces_, patchi)
+    {
+        forAll(excludedFaces_[patchi], i)
         {
-            label patchi = mesh_.boundaryMesh().whichPatch(facei);
-            label patchFacei = mesh_.boundaryMesh()[patchi].whichFace(facei);
+            label patchFacei = excludedFaces_[patchi][i];
 
             phi.boundaryField()[patchi][patchFacei] -=
                 (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
               & Sf.boundaryField()[patchi][patchFacei];
         }
     }
-
-    forAll(patchLabels_, i)
-    {
-        phi.boundaryField()[patchLabels_[i]] = 0.0;
-    }
 }
 
 
@@ -239,10 +415,21 @@ void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
     const vector& origin = origin_.value();
     const vector& Omega = Omega_.value();
 
-    forAll(patchLabels_, i)
+    // Included patches
+    forAll(includedFaces_, patchi)
     {
-        U.boundaryField()[patchLabels_[i]] ==
-            (Omega ^ (mesh_.Cf().boundaryField()[patchLabels_[i]] - origin));
+        const vectorField& patchC = mesh_.Cf().boundaryField()[patchi];
+
+        vectorField pfld(U.boundaryField()[patchi]);
+
+        forAll(includedFaces_[patchi], i)
+        {
+            label patchFacei = includedFaces_[patchi][i];
+
+            pfld[patchFacei] = (Omega ^ (patchC[patchFacei] - origin));
+        }
+
+        U.boundaryField()[patchi] == pfld;
     }
 }
 
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
index 7c489ed25b..5926e87f10 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.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) 1991-2008 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -75,16 +75,24 @@ class MRFZone
         label cellZoneID_;
 
         //- label of face zone with faces on outside of cell zone.
-        //  If -2 determines outside faces itself
         label faceZoneID_;
 
-        //- outside faces (only if faceZoneID = -2)
-        labelList outsideFaces_;
-
+        //- Do patches move with frame (true) or are explicitly provided (false,
+        //  old behaviour)
+        Switch allPatchesMove_;
 
         const wordList patchNames_;
         labelList patchLabels_;
 
+        //- Internal faces that are part of MRF
+        labelList internalFaces_;
+
+        //- Outside faces (per patch) that move with the MRF
+        labelListList includedFaces_;
+
+        //- Excluded faces (per patch) that do not move with the MRF
+        labelListList excludedFaces_;
+
         const dimensionedVector origin_;
         dimensionedVector axis_;
         const dimensionedScalar omega_;
@@ -93,6 +101,17 @@ class MRFZone
 
     // Private Member Functions
 
+        //- Divide faces in frame according to patch
+        void setMRFFaces
+        (
+            labelList& faceType,
+            const labelList& excludedPatchIDs
+        );
+
+        //- Divide faces in frame according to patch
+        void setMRFFaces();
+
+
         //- Disallow default bitwise copy construct
         MRFZone(const MRFZone&);
 
@@ -136,6 +155,13 @@ public:
 
     // Member Functions
 
+        //- Update the mesh corresponding to given map
+        void updateMesh(const mapPolyMesh& mpm)
+        {
+            // Only updates face addressing
+            setMRFFaces();
+        }
+
         //- Add the Coriolis force contribution to the momentum equation
         void addCoriolis(fvVectorMatrix& UEqn) const;
 
diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/constant/MRFZones b/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/constant/MRFZones
index 7c8989f567..e7b7f0aa63 100644
--- a/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/constant/MRFZones
+++ b/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/constant/MRFZones
@@ -19,7 +19,9 @@ FoamFile
 (
     rotor
     {
-        patches   (rotor);
+        // Fixed patches (by default they 'move' with the MRF zone)
+        nonRotatingPatches ();
+
         origin    origin [0 1 0 0 0 0 0]  (0 0 0);
         axis      axis   [0 0 0 0 0 0 0]  (0 0 1);
         omega     omega  [0 0 -1 0 0 0 0] 104.72;