From aece1c1c4d65444ec9c17da8ac4983ccb6650917 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 13 Oct 2015 16:32:52 +0100
Subject: [PATCH] ENH: FaceCellWave, localPointRegion: support for local
 baffles

---
 .../algorithms/MeshWave/FaceCellWave.C        | 162 +++++++++++++++++-
 .../algorithms/MeshWave/FaceCellWave.H        |  27 ++-
 .../indexedOctree/treeDataPrimitivePatch.C    |  23 ++-
 src/meshTools/meshTools/meshTools.C           |  12 +-
 src/meshTools/meshTools/meshTools.H           |   6 +-
 src/meshTools/regionSplit/localPointRegion.C  |  46 ++++-
 src/meshTools/regionSplit/localPointRegion.H  |  12 +-
 7 files changed, 253 insertions(+), 35 deletions(-)

diff --git a/src/meshTools/algorithms/MeshWave/FaceCellWave.C b/src/meshTools/algorithms/MeshWave/FaceCellWave.C
index f254302a7f..08ff51e3ed 100644
--- a/src/meshTools/algorithms/MeshWave/FaceCellWave.C
+++ b/src/meshTools/algorithms/MeshWave/FaceCellWave.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -810,6 +810,79 @@ void Foam::FaceCellWave<Type, TrackingData>::handleAMICyclicPatches()
 }
 
 
+template<class Type, class TrackingData>
+void Foam::FaceCellWave<Type, TrackingData>::handleExplicitConnections()
+{
+    // Collect changed information
+
+    DynamicList<label> f0Baffle(explicitConnections_.size());
+    DynamicList<Type> f0Info(explicitConnections_.size());
+
+    DynamicList<label> f1Baffle(explicitConnections_.size());
+    DynamicList<Type> f1Info(explicitConnections_.size());
+
+    forAll(explicitConnections_, connI)
+    {
+        const labelPair& baffle = explicitConnections_[connI];
+
+        label f0 = baffle[0];
+        if (changedFace_[f0])
+        {
+            f0Baffle.append(connI);
+            f0Info.append(allFaceInfo_[f0]);
+        }
+
+        label f1 = baffle[1];
+        if (changedFace_[f1])
+        {
+            f1Baffle.append(connI);
+            f1Info.append(allFaceInfo_[f1]);
+        }
+    }
+
+
+    // Update other side with changed information
+
+    forAll(f1Info, index)
+    {
+        const labelPair& baffle = explicitConnections_[f1Baffle[index]];
+
+        label f0 = baffle[0];
+        Type& currentWallInfo = allFaceInfo_[f0];
+
+        if (!currentWallInfo.equal(f1Info[index], td_))
+        {
+            updateFace
+            (
+                f0,
+                f1Info[index],
+                propagationTol_,
+                currentWallInfo
+            );
+        }
+    }
+
+    forAll(f0Info, index)
+    {
+        const labelPair& baffle = explicitConnections_[f0Baffle[index]];
+
+        label f1 = baffle[1];
+        Type& currentWallInfo = allFaceInfo_[f1];
+
+        if (!currentWallInfo.equal(f0Info[index], td_))
+        {
+            updateFace
+            (
+                f1,
+                f0Info[index],
+                propagationTol_,
+                currentWallInfo
+            );
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Set up only. Use setFaceInfo and iterate() to do actual calculation.
@@ -823,6 +896,7 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
 )
 :
     mesh_(mesh),
+    explicitConnections_(0),
     allFaceInfo_(allFaceInfo),
     allCellInfo_(allCellInfo),
     td_(td),
@@ -878,6 +952,7 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
 )
 :
     mesh_(mesh),
+    explicitConnections_(0),
     allFaceInfo_(allFaceInfo),
     allCellInfo_(allCellInfo),
     td_(td),
@@ -939,6 +1014,87 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
 }
 
 
+// Iterate, propagating changedFacesInfo across mesh, until no change (or
+// maxIter reached). Initial cell values specified.
+template<class Type, class TrackingData>
+Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
+(
+    const polyMesh& mesh,
+    const List<labelPair>& explicitConnections,
+    const bool handleCyclicAMI,
+    const labelList& changedFaces,
+    const List<Type>& changedFacesInfo,
+    UList<Type>& allFaceInfo,
+    UList<Type>& allCellInfo,
+    const label maxIter,
+    TrackingData& td
+)
+:
+    mesh_(mesh),
+    explicitConnections_(explicitConnections),
+    allFaceInfo_(allFaceInfo),
+    allCellInfo_(allCellInfo),
+    td_(td),
+    changedFace_(mesh_.nFaces(), false),
+    changedFaces_(mesh_.nFaces()),
+    nChangedFaces_(0),
+    changedCell_(mesh_.nCells(), false),
+    changedCells_(mesh_.nCells()),
+    nChangedCells_(0),
+    hasCyclicPatches_(hasPatch<cyclicPolyPatch>()),
+    hasCyclicAMIPatches_
+    (
+        handleCyclicAMI
+     && returnReduce(hasPatch<cyclicAMIPolyPatch>(), orOp<bool>())
+    ),
+    nEvals_(0),
+    nUnvisitedCells_(mesh_.nCells()),
+    nUnvisitedFaces_(mesh_.nFaces())
+{
+    if
+    (
+        allFaceInfo.size() != mesh_.nFaces()
+     || allCellInfo.size() != mesh_.nCells()
+    )
+    {
+        FatalErrorIn
+        (
+            "FaceCellWave<Type, TrackingData>::FaceCellWave"
+            "(const polyMesh&, const List<labelPair>&, const labelList&"
+            ", const List<Type>,"
+            " UList<Type>&, UList<Type>&, const label maxIter)"
+        )   << "face and cell storage not the size of mesh faces, cells:"
+            << endl
+            << "    allFaceInfo   :" << allFaceInfo.size() << endl
+            << "    mesh_.nFaces():" << mesh_.nFaces() << endl
+            << "    allCellInfo   :" << allCellInfo.size() << endl
+            << "    mesh_.nCells():" << mesh_.nCells()
+            << exit(FatalError);
+    }
+
+    // Copy initial changed faces data
+    setFaceInfo(changedFaces, changedFacesInfo);
+
+    // Iterate until nothing changes
+    label iter = iterate(maxIter);
+
+    if ((maxIter > 0) && (iter >= maxIter))
+    {
+        FatalErrorIn
+        (
+            "FaceCellWave<Type, TrackingData>::FaceCellWave"
+            "(const polyMesh&, const List<labelPair>&, const labelList&"
+            ", const List<Type>, UList<Type>&, UList<Type>&"
+            ", const label maxIter)"
+        )   << "Maximum number of iterations reached. Increase maxIter." << endl
+            << "    maxIter:" << maxIter << endl
+            << "    nChangedCells:" << nChangedCells_ << endl
+            << "    nChangedFaces:" << nChangedFaces_ << endl
+            << exit(FatalError);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 
@@ -1093,6 +1249,10 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::cellToFace()
     // Handled all changed cells by now
     nChangedCells_ = 0;
 
+
+    // Transfer across any explicitly provided internal connections
+    handleExplicitConnections();
+
     if (hasCyclicPatches_)
     {
         // Transfer changed faces across cyclic halves
diff --git a/src/meshTools/algorithms/MeshWave/FaceCellWave.H b/src/meshTools/algorithms/MeshWave/FaceCellWave.H
index 93e2277daa..f8d43e6865 100644
--- a/src/meshTools/algorithms/MeshWave/FaceCellWave.H
+++ b/src/meshTools/algorithms/MeshWave/FaceCellWave.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,6 +50,7 @@ SourceFiles
 #include "boolList.H"
 #include "labelList.H"
 #include "primitiveFieldsFwd.H"
+#include "labelPair.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -81,6 +82,9 @@ class FaceCellWave
         //- Reference to mesh
         const polyMesh& mesh_;
 
+        //- Optional boundary faces that information should travel through
+        const List<labelPair> explicitConnections_;
+
         //- Information for all faces
         UList<Type>& allFaceInfo_;
 
@@ -232,6 +236,10 @@ class FaceCellWave
             //- Merge data from across AMI cyclics
             void handleAMICyclicPatches();
 
+            //- Merge data across explicitly provided local connections (usually
+            //  baffles)
+            void handleExplicitConnections();
+
 
       // Private static data
 
@@ -286,6 +294,23 @@ public:
             TrackingData& td = dummyTrackData_
         );
 
+        //- Construct from mesh and explicitly connected boundary faces
+        //  and list of changed faces with the Type
+        //  for these faces. Iterates until nothing changes or maxIter reached.
+        //  (maxIter can be 0)
+        FaceCellWave
+        (
+            const polyMesh&,
+            const List<labelPair>& explicitConnections,
+            const bool handleCyclicAMI,
+            const labelList& initialChangedFaces,
+            const List<Type>& changedFacesInfo,
+            UList<Type>& allFaceInfo,
+            UList<Type>& allCellInfo,
+            const label maxIter,
+            TrackingData& td = dummyTrackData_
+        );
+
 
     // Member Functions
 
diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
index 51e5252121..c5771463c5 100644
--- a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
+++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -188,8 +188,10 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
             << " nearest face:" << faceI;
     }
 
-    const pointField& points = patch_.localPoints();
-    const typename PatchType::FaceType& f = patch_.localFaces()[faceI];
+    const typename PatchType::FaceType& localF = patch_.localFaces()[faceI];
+    const typename PatchType::FaceType& f = patch_[faceI];
+    const pointField& points = patch_.points();
+    const labelList& mp = patch_.meshPoints();
 
     // Retest to classify where on face info is. Note: could be improved. We
     // already have point.
@@ -242,7 +244,7 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
 
             return indexedOctree<treeDataPrimitivePatch>::getSide
             (
-                patch_.pointNormals()[f[fp]],
+                patch_.pointNormals()[localF[fp]],
                 sample - curPt
             );
         }
@@ -280,8 +282,8 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
     {
         label edgeI = fEdges[fEdgeI];
         const edge& e = patch_.edges()[edgeI];
-
-        pointHit edgeHit = e.line(points).nearestDist(sample);
+        const linePointRef ln(points[mp[e.start()]], points[mp[e.end()]]);
+        pointHit edgeHit = ln.nearestDist(sample);
 
         if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < planarTol_)
         {
@@ -322,11 +324,7 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
 
     forAll(f, fp)
     {
-        pointHit edgeHit = linePointRef
-        (
-            points[f[fp]],
-            fc
-        ).nearestDist(sample);
+        pointHit edgeHit = linePointRef(points[f[fp]], fc).nearestDist(sample);
 
         if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < planarTol_)
         {
@@ -369,7 +367,8 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
 
         forAll(f, fp)
         {
-            Pout<< "    vertex:" << f[fp] << "  coord:" << points[f[fp]]
+            Pout<< "    vertex:" << f[fp]
+                << "  coord:" << points[f[fp]]
                 << endl;
         }
     }
diff --git a/src/meshTools/meshTools/meshTools.C b/src/meshTools/meshTools/meshTools.C
index 5c98fe087b..cf208822a4 100644
--- a/src/meshTools/meshTools/meshTools.C
+++ b/src/meshTools/meshTools/meshTools.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,7 +29,6 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-// Check if n is in same direction as normals of all faceLabels
 bool Foam::meshTools::visNormal
 (
     const vector& n,
@@ -309,7 +308,6 @@ bool Foam::meshTools::edgeOnFace
 }
 
 
-// Return true if faceI part of cellI
 bool Foam::meshTools::faceOnCell
 (
     const primitiveMesh& mesh,
@@ -465,7 +463,6 @@ Foam::label Foam::meshTools::getSharedFace
 }
 
 
-// Get the two faces on cellI using edgeI.
 void Foam::meshTools::getEdgeFaces
 (
     const primitiveMesh& mesh,
@@ -511,7 +508,6 @@ void Foam::meshTools::getEdgeFaces
 }
 
 
-// Return label of other edge connected to vertex
 Foam::label Foam::meshTools::otherEdge
 (
     const primitiveMesh& mesh,
@@ -549,7 +545,6 @@ Foam::label Foam::meshTools::otherEdge
 }
 
 
-// Return face on other side of edgeI
 Foam::label Foam::meshTools::otherFace
 (
     const primitiveMesh& mesh,
@@ -574,7 +569,6 @@ Foam::label Foam::meshTools::otherFace
 }
 
 
-// Return face on other side of edgeI
 Foam::label Foam::meshTools::otherCell
 (
     const primitiveMesh& mesh,
@@ -602,8 +596,6 @@ Foam::label Foam::meshTools::otherCell
 }
 
 
-// Returns label of edge nEdges away from startEdge (in the direction of
-// startVertI)
 Foam::label Foam::meshTools::walkFace
 (
     const primitiveMesh& mesh,
@@ -688,7 +680,6 @@ void Foam::meshTools::constrainToMeshCentre
 }
 
 
-//- Set the constrained components of directions/velocity to zero
 void Foam::meshTools::constrainDirection
 (
     const polyMesh& mesh,
@@ -814,7 +805,6 @@ Foam::vector Foam::meshTools::edgeToCutDir
 }
 
 
-// Find edges most aligned with cutDir
 Foam::label Foam::meshTools::cutDirToEdge
 (
     const primitiveMesh& mesh,
diff --git a/src/meshTools/meshTools/meshTools.H b/src/meshTools/meshTools/meshTools.H
index 08d36f3208..74bf2fbfc3 100644
--- a/src/meshTools/meshTools/meshTools.H
+++ b/src/meshTools/meshTools/meshTools.H
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -200,7 +200,7 @@ namespace meshTools
             const label v1
         );
 
-        //- Return edge between two vertices. Returns -1 if no edge.
+        //- Return edge between two mesh vertices. Returns -1 if no edge.
         label findEdge
         (
             const primitiveMesh&,
diff --git a/src/meshTools/regionSplit/localPointRegion.C b/src/meshTools/regionSplit/localPointRegion.C
index bad9ca8d57..5eb8d10c52 100644
--- a/src/meshTools/regionSplit/localPointRegion.C
+++ b/src/meshTools/regionSplit/localPointRegion.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -242,6 +242,7 @@ void Foam::localPointRegion::countPointRegions
 void Foam::localPointRegion::calcPointRegions
 (
     const polyMesh& mesh,
+    const labelPairList& baffles,
     boolList& candidatePoint
 )
 {
@@ -423,6 +424,13 @@ void Foam::localPointRegion::calcPointRegions
             minEqOpFace(),
             Foam::dummyTransform()  // dummy transformation
         );
+        forAll(baffles, i)
+        {
+            label f0 = baffles[i].first();
+            label f1 = baffles[i].second();
+            minEqOpFace()(minRegion[f0], minRegion[f1]);
+            minRegion[f1] = minRegion[f0];
+        }
     }
 
 
@@ -469,7 +477,7 @@ Foam::localPointRegion::localPointRegion(const polyMesh& mesh)
         }
     }
 
-    calcPointRegions(mesh, candidatePoint);
+    calcPointRegions(mesh, labelPairList(0), candidatePoint);
 }
 
 
@@ -492,7 +500,31 @@ Foam::localPointRegion::localPointRegion
         candidatePoint[candidatePoints[i]] = true;
     }
 
-    calcPointRegions(mesh, candidatePoint);
+    calcPointRegions(mesh, labelPairList(0), candidatePoint);
+}
+
+
+Foam::localPointRegion::localPointRegion
+(
+    const polyMesh& mesh,
+    const labelPairList& baffles,
+    const labelList& candidatePoints
+)
+:
+    //nRegions_(0),
+    meshPointMap_(0),
+    pointRegions_(0),
+    meshFaceMap_(0),
+    faceRegions_(0)
+{
+    boolList candidatePoint(mesh.nPoints(), false);
+
+    forAll(candidatePoints, i)
+    {
+        candidatePoint[candidatePoints[i]] = true;
+    }
+
+    calcPointRegions(mesh, baffles, candidatePoint);
 }
 
 
@@ -630,14 +662,18 @@ Foam::List<Foam::labelPair> Foam::localPointRegion::findDuplicateFacePairs
                     << " processorPolyPatch."
                     << "This is not allowed." << nl
                     << "Face:" << meshFace0
+                    << " fc:" << mesh.faceCentres()[meshFace0]
                     << " is on patch:" << patches[patch0].name()
                     << nl
                     << "Face:" << meshFace1
+                    << " fc:" << mesh.faceCentres()[meshFace1]
                     << " is on patch:" << patches[patch1].name()
                     << abort(FatalError);
             }
-
-            baffles.append(labelPair(meshFace0, meshFace1));
+            else
+            {
+                baffles.append(labelPair(meshFace0, meshFace1));
+            }
         }
     }
     return baffles.shrink();
diff --git a/src/meshTools/regionSplit/localPointRegion.H b/src/meshTools/regionSplit/localPointRegion.H
index 83bbce1e51..aadb28d904 100644
--- a/src/meshTools/regionSplit/localPointRegion.H
+++ b/src/meshTools/regionSplit/localPointRegion.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -101,10 +101,10 @@ class localPointRegion
         void calcPointRegions
         (
             const polyMesh& mesh,
+            const labelPairList& baffles,
             boolList& candidatePoint
         );
 
-
         //- Check if two faces are equal. If forward = false checks f1 in
         //  reverse order.
         static bool isDuplicate
@@ -133,6 +133,14 @@ public:
             const labelList& candidatePoints
         );
 
+        //- Construct from mesh and candidate points for duplication
+        localPointRegion
+        (
+            const polyMesh& mesh,
+            const labelPairList& baffles,
+            const labelList& candidatePoints
+        );
+
 
     // Member Functions