ENH: overset: various improvements to trackingInverseDistance. See #736.

src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C

@@ -48,13 +48,13 @@ namespace cellCellStencils
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::cellCellStencils::trackingInverseDistance::markBoundaries
+bool Foam::cellCellStencils::trackingInverseDistance::markBoundaries
 (
     const fvMesh& mesh,
     const vector& smallVec,
 
     const boundBox& bb,
-    const labelVector& nDivs,
+    labelVector& nDivs,
     PackedList<2>& patchTypes,
 
     const labelList& cellMap,
@@ -63,6 +63,8 @@ void Foam::cellCellStencils::trackingInverseDistance::markBoundaries
 {
     // Mark all voxels that overlap the bounding box of any patch
 
+    static bool hasWarned = false;
+
     const fvBoundaryMesh& pbm = mesh.boundary();
 
     patchTypes = patchCellType::OTHER;
@@ -113,6 +115,7 @@ void Foam::cellCellStencils::trackingInverseDistance::markBoundaries
         {
             //Info<< "Marking cells on overset patch " << fvp.name() << endl;
             const polyPatch& pp = fvp.patch();
+            const vectorField::subField faceCentres(pp.faceCentres());
             forAll(pp, i)
             {
                 // Mark in overall patch types
@@ -122,8 +125,54 @@ void Foam::cellCellStencils::trackingInverseDistance::markBoundaries
                 boundBox faceBb(pp.points(), pp[i], false);
                 faceBb.min() -= smallVec;
                 faceBb.max() += smallVec;
-                if (bb.overlaps(faceBb))
+                if (!bb.contains(faceCentres[i]))
                 {
+                    if (!hasWarned)
+                    {
+                        hasWarned = true;
+                        WarningInFunction << "Patch " << pp.name()
+                            << " : face at " << faceCentres[i]
+                            << " is not inside search bounding box of"
+                            << " voxel mesh " << bb << endl
+                            << "    Is your 'searchBox' specification correct?"
+                            << endl;
+                    }
+                }
+                else
+                {
+                    // Test for having voxels already marked as patch
+                    // -> voxel mesh is too coarse
+                    if
+                    (
+                        voxelMeshSearch::overlaps
+                        (
+                            bb,
+                            nDivs,
+                            faceBb,
+                            patchTypes,
+                            static_cast<unsigned int>(patchCellType::PATCH)
+                        )
+                    )
+                    {
+                        // Determine number of voxels from number of cells
+                        // in mesh
+                        const labelVector& dim = mesh.geometricD();
+                        forAll(dim, i)
+                        {
+                            if (dim[i] != -1)
+                            {
+                                nDivs[i] *= 2;
+                            }
+                        }
+
+                        Pout<< "Voxel mesh too coarse. Bounding box "
+                            << faceBb
+                            << " contains both non-overset and overset patches"
+                            << ". Refining voxel mesh to " << nDivs << endl;
+
+                        return false;
+                    }
+
                     voxelMeshSearch::fill
                     (
                         patchTypes,
@@ -136,6 +185,8 @@ void Foam::cellCellStencils::trackingInverseDistance::markBoundaries
             }
         }
     }
+
+    return true;
 }
 
 
@@ -522,64 +573,15 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
 
     DebugInfo<< FUNCTION_NAME << " : Moved zone sub-meshes" << endl;
 
-
     // Calculate fast search structure for each zone
-    PtrList<voxelMeshSearch> meshSearches(nZones);
+    List<labelVector> searchBoxDivisions(dict_["searchBoxDivisions"]);
-
+    if (searchBoxDivisions.size() != nZones)
-    List<labelVector> searchBoxDivisions;
-    if (dict_.readIfPresent("searchBoxDivisions", searchBoxDivisions))
     {
-        forAll(meshParts_, zonei)
+        FatalIOErrorInFunction(dict_) << "Number of searchBoxDivisions "
-        {
+            << searchBoxDivisions.size()
-            meshSearches.set
+            << " should  equal the number of zones " << nZones
-            (
+            << exit(FatalIOError);
-                zonei,
-                new voxelMeshSearch
-                (
-                    meshParts_[zonei].subMesh(),
-                    searchBoxDivisions[zonei],
-                    true
-                )
-            );
     }
-    }
-    else
-    {
-        forAll(meshParts_, zonei)
-        {
-            meshSearches.set
-            (
-                zonei,
-                new voxelMeshSearch
-                (
-                    meshParts_[zonei].subMesh(),
-                    true
-                )
-            );
-        }
-    }
-    DebugInfo<< FUNCTION_NAME << " : Constructed cell voxel meshes" << endl;
-
-    PtrList<fvMesh> voxelMeshes;
-    if (debug)
-    {
-        voxelMeshes.setSize(meshSearches.size());
-        forAll(meshSearches, zonei)
-        {
-            IOobject meshIO
-            (
-                word("voxelMesh" + Foam::name(zonei)),
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ
-            );
-
-            Pout<< "Constructing voxel mesh " << meshIO.objectPath() << endl;
-            voxelMeshes.set(zonei, meshSearches[zonei].makeMesh(meshIO));
-        }
-    }
-
-
 
     const boundBox& allBb(mesh_.bounds());
 
@@ -652,21 +654,18 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
             patchBb = meshBb;
         }
     }
-    if (patchDivisions.empty())
-    {
-        patchDivisions.setSize(nZones);
-        forAll(patchDivisions, zonei)
-        {
-            patchDivisions[zonei] = meshSearches[zonei].nDivs();
-        }
-    }
 
     forAll(patchParts, zonei)
     {
-        const labelVector& g = patchDivisions[zonei];
+        while (true)
-        patchParts.set(zonei, new PackedList<2>(cmptProduct(g)));
+        {
+            patchParts.set
+            (
+                zonei,
+                new PackedList<2>(cmptProduct(patchDivisions[zonei]))
+            );
 
-        markBoundaries
+            bool ok = markBoundaries
             (
                 meshParts_[zonei].subMesh(),
                 smallVec_,
@@ -678,9 +677,55 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
                 meshParts_[zonei].cellMap(),
                 allPatchTypes
             );
+
+            //if (returnReduce(ok, andOp<bool>()))
+            if (ok)
+            {
+                break;
+            }
+        }
     }
     DebugInfo<< FUNCTION_NAME << " : Calculated boundary voxel meshes" << endl;
 
+
+    PtrList<voxelMeshSearch> meshSearches(meshParts_.size());
+    forAll(meshParts_, zonei)
+    {
+        meshSearches.set
+        (
+            zonei,
+            new voxelMeshSearch
+            (
+                meshParts_[zonei].subMesh(),
+                patchBb[zonei][Pstream::myProcNo()],
+                patchDivisions[zonei],
+                true
+            )
+        );
+    }
+
+    DebugInfo<< FUNCTION_NAME << " : Constructed cell voxel meshes" << endl;
+
+    PtrList<fvMesh> voxelMeshes;
+    if (debug)
+    {
+        voxelMeshes.setSize(meshSearches.size());
+        forAll(meshSearches, zonei)
+        {
+            IOobject meshIO
+            (
+                word("voxelMesh" + Foam::name(zonei)),
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ
+            );
+
+            Pout<< "Constructing voxel mesh " << meshIO.objectPath() << endl;
+            voxelMeshes.set(zonei, meshSearches[zonei].makeMesh(meshIO));
+        }
+    }
+
+
     if (debug)
     {
         forAll(patchParts, zonei)

src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2017-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,13 +74,13 @@ protected:
     // Protected Member Functions
 
         //- Mark voxels of patchTypes with type of patch face
-        static void markBoundaries
+        static bool markBoundaries
         (
             const fvMesh& mesh,
             const vector& smallVec,
 
             const boundBox& bb,
-            const labelVector& nDivs,
+            labelVector& nDivs,
             PackedList<2>& patchTypes,
 
             const labelList& cellMap,

src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.C

@@ -281,6 +281,14 @@ Foam::voxelMeshSearch::voxelMeshSearch
         }
     }
 
+    // Redo the local bounding box
+    localBb_ = boundBox(mesh_.points(), false);
+
+    const point eps(1e-10, 1e-10, 1e-10);
+
+    localBb_.min() = localBb_.min()-eps;
+    localBb_.max() = localBb_.max()+eps;
+
     if (debug)
     {
         Pout<< "voxelMeshSearch : mesh:" << mesh_.name()
@@ -297,11 +305,13 @@ Foam::voxelMeshSearch::voxelMeshSearch
 Foam::voxelMeshSearch::voxelMeshSearch
 (
     const polyMesh& mesh,
+    const boundBox& localBb,
     const labelVector& nDivs,
     const bool doUpdate
 )
 :
     mesh_(mesh),
+    localBb_(localBb),
     nDivs_(nDivs)
 {
     if (doUpdate)
@@ -315,14 +325,6 @@ Foam::voxelMeshSearch::voxelMeshSearch
 
 bool Foam::voxelMeshSearch::update()
 {
-    // Redo the local bounding box
-    localBb_ = boundBox(mesh_.points(), false);
-
-    const point eps(1e-10, 1e-10, 1e-10);
-
-    localBb_.min() = localBb_.min()-eps;
-    localBb_.max() = localBb_.max()+eps;
-
     // Initialise seed cell array
 
     seedCell_.setSize(cmptProduct(nDivs_));
@@ -373,6 +375,7 @@ Foam::label Foam::voxelMeshSearch::findCell(const point& p) const
         return -1;
     }
 
+
     // Locate the voxel index for this point. Do not clip.
     label voxeli = index(localBb_, nDivs_, p, false);
 
@@ -397,10 +400,14 @@ Foam::label Foam::voxelMeshSearch::findCell(const point& p) const
             // found does not have to be the absolute 'correct' one as
             // long as at least one of the processors finds a cell.
 
-            label nextCellOld = -1;
+            track_.clear();
-
             while (true)
             {
+                if (track_.size() < 5)
+                {
+                    track_.append(celli);
+                }
+
                 // I am in celli now. How many faces do I have ?
                 label facei = findIntersectedFace(celli, p);
 
@@ -409,6 +416,7 @@ Foam::label Foam::voxelMeshSearch::findCell(const point& p) const
                     return celli;
                 }
 
+                const label startOfTrack(max(0, track_.size()-5));
 
                 label nextCell;
                 if (mesh_.isInternalFace(facei))
@@ -416,21 +424,24 @@ Foam::label Foam::voxelMeshSearch::findCell(const point& p) const
                     label own = mesh_.faceOwner()[facei];
                     label nei = mesh_.faceNeighbour()[facei];
                     nextCell = (own == celli ? nei : own);
+
+                    if (findIndex(track_, nextCell, startOfTrack) != -1)
+                    {
+                        return celli;
+                    }
                 }
                 else
                 {
                     nextCell = searchProcPatch(facei, p);
 
-                    if (nextCell == nextCellOld)
-                    {
-                        return -1; // point is really out
-                    }
-
                     if (nextCell == -1 || nextCell == celli)
                     {
                         return nextCell;
                     }
-                    nextCellOld = nextCell;
+                    else if (findIndex(track_, nextCell, startOfTrack) != -1)
+                    {
+                        return -1;  // point is really out
+                    }
                 }
 
                 celli = nextCell;

src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.H

@@ -68,6 +68,9 @@ class voxelMeshSearch
         //- Voxel to seed cell
         labelList seedCell_;
 
+        //- Cells visited
+        mutable DynamicList<label> track_;
+
 
    // Private Member Functions
 
@@ -93,7 +96,8 @@ public:
         voxelMeshSearch
         (
             const polyMesh&,
-            const labelVector&,
+            const boundBox& bb,
+            const labelVector& nDivs,
             const bool doUpdate = true
         );
 

tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/fvSchemes

@@ -54,6 +54,11 @@ snGradSchemes
 oversetInterpolation
 {
     method          cellVolumeWeight;
+
+    // Faster but less accurate
+    //method              trackingInverseDistance;
+    //searchBox           (0 0 0)(0.02 0.01 0.01);
+    //searchBoxDivisions  3{(64 64 1)};
 }
 
 fluxRequired