ENH: add primitiveMesh cellBb()

- the boundBox for a given cell, using the cheapest calculation:

  - cellPoints if already available, since this will involve the
    fewest number of min/max comparisions.

  - otherwise walk the cell faces: via the cell box() method
    to avoid creating demand-driven cellPoints etc.

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C

@@ -167,7 +167,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
             {
                 if (volumeStatus[celli] == volumeType::UNKNOWN)
                 {
-                    treeBoundBox cellBb(mesh_.cells()[celli].box(mesh_));
+                    treeBoundBox cellBb(mesh_.cellBb(celli));
 
                     if (geometry.overlaps(cellBb))
                     {
@@ -279,7 +279,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
             {
                 if (volumeStatus[celli] == volumeType::UNKNOWN)
                 {
-                    treeBoundBox cellBb(mesh_.cells()[celli].box(mesh_));
+                    treeBoundBox cellBb(mesh_.cellBb(celli));
 
                     if (geometry.overlaps(cellBb))
                     {
@@ -498,7 +498,7 @@ bool Foam::backgroundMeshDecomposition::refineCell
 
 //    const conformationSurfaces& geometry = geometryToConformTo_;
 
-    treeBoundBox cellBb(mesh_.cells()[celli].box(mesh_));
+    treeBoundBox cellBb(mesh_.cellBb(celli));
 
     weightEstimate = 1.0;
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -71,6 +71,7 @@ namespace Foam
 
 // Forward Declarations
 class bitSet;
+class boundBox;
 
 /*---------------------------------------------------------------------------*\
                       Class primitiveMesh Declaration
@@ -769,6 +770,9 @@ public:
 
         // Useful derived info
 
+            //- The bounding box for given cell index
+            boundBox cellBb(const label celli) const;
+
             //- Return true if the point in the cell bounding box.
             //  The bounding box may be isotropically inflated by the fraction
             //  inflationFraction

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,6 +32,18 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::boundBox Foam::primitiveMesh::cellBb(const label celli) const
+{
+    if (hasCellPoints())
+    {
+        // No reduction!
+        return boundBox(points(), cellPoints()[celli], false);
+    }
+
+    return boundBox(cells()[celli].box(points(), faces()));
+}
+
+
 bool Foam::primitiveMesh::pointInCellBB
 (
     const point& p,
@@ -39,15 +51,7 @@ bool Foam::primitiveMesh::pointInCellBB
     scalar inflationFraction
 ) const
 {
-    boundBox bb
+    boundBox bb(cellBb(celli));
-    (
-        cells()[celli].points
-        (
-            faces(),
-            points()
-        ),
-        false
-    );
 
     if (inflationFraction > SMALL)
     {

src/functionObjects/randomProcesses/energySpectrum/energySpectrum.C

@@ -150,13 +150,7 @@ bool Foam::functionObjects::energySpectrum::read(const dictionary& dict)
     const boundBox meshBb(mesh_.bounds());
 
     // Assume all cells are the same size...
-    const cell& c = mesh_.cells()[0];
+    boundBox cellBb(mesh_.cellBb(0));
-    boundBox cellBb(boundBox::invertedBox);
-    forAll(c, facei)
-    {
-        const face& f = mesh_.faces()[c[facei]];
-        cellBb.add(mesh_.points(), f);
-    }
 
     const vector L(meshBb.max() - meshBb.min());
     const vector nCellXYZ(cmptDivide(L, cellBb.max() - cellBb.min()));

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.C

@@ -133,9 +133,9 @@ void Foam::snappyVoxelMeshDriver::isInside
     }
     else
     {
-        for (const cell& c : mesh.cells())
+        for (label celli = 0; celli < mesh.nCells(); ++celli)
         {
-            const boundBox cellBb(c.box(mesh));
+            const boundBox cellBb(mesh.cellBb(celli));
 
             voxelMeshSearch::fill
             (

src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C

@@ -236,20 +236,6 @@ void Foam::cellCellStencils::inverseDistance::markBoundaries
 }
 
 
-Foam::treeBoundBox Foam::cellCellStencils::inverseDistance::cellBb
-(
-    const primitiveMesh& mesh,
-    const label celli
-)
-{
-    if (mesh.hasCellPoints())
-    {
-        return treeBoundBox(mesh.points(), mesh.cellPoints(celli));
-    }
-    return treeBoundBox(mesh.cells()[celli].box(mesh));
-}
-
-
 bool Foam::cellCellStencils::inverseDistance::overlaps
 (
     const boundBox& bb,
@@ -326,7 +312,7 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
             forAll(tgtCellMap, tgtCelli)
             {
                 label celli = tgtCellMap[tgtCelli];
-                treeBoundBox cBb(cellBb(mesh_, celli));
+                treeBoundBox cBb(mesh_.cellBb(celli));
                 cBb.min() -= smallVec_;
                 cBb.max() += smallVec_;
 
@@ -396,9 +382,10 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
                 forAll(tgtCellMap, tgtCelli)
                 {
                     label celli = tgtCellMap[tgtCelli];
-                    treeBoundBox cBb(cellBb(mesh_, celli));
+                    treeBoundBox cBb(mesh_.cellBb(celli));
                     cBb.min() -= smallVec_;
                     cBb.max() += smallVec_;
+
                     if
                     (
                         overlaps
@@ -554,7 +541,7 @@ void Foam::cellCellStencils::inverseDistance::markDonors
         label celli = tgtCellMap[tgtCelli];
         if (srcOverlapProcs.size())
         {
-            treeBoundBox subBb(cellBb(mesh_, celli));
+            treeBoundBox subBb(mesh_.cellBb(celli));
             subBb.min() -= smallVec_;
             subBb.max() += smallVec_;
 

src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H

@@ -179,13 +179,6 @@ protected:
             labelList& patchCellTypes
         );
 
-        //- Calculate bounding box of cell
-        static treeBoundBox cellBb
-        (
-            const primitiveMesh& mesh,
-            const label celli
-        );
-
         //- Mark all cells overlapping (a voxel covered by) a src patch
         //  with type HOLE
         void markPatchesAsHoles

src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C

@@ -206,9 +206,6 @@ void Foam::cellCellStencils::trackingInverseDistance::markPatchesAsHoles
     labelList& allCellTypes
 ) const
 {
-    const pointField& allPoints = mesh_.points();
-    const labelListList& allCellPoints = mesh_.cellPoints();
-
     const treeBoundBoxList& srcPatchBbs = patchBb[srcI];
     const treeBoundBoxList& tgtPatchBbs = patchBb[tgtI];
     const labelList& tgtCellMap = meshParts_[tgtI].cellMap();
@@ -226,7 +223,7 @@ void Foam::cellCellStencils::trackingInverseDistance::markPatchesAsHoles
             forAll(tgtCellMap, tgtCelli)
             {
                 label celli = tgtCellMap[tgtCelli];
-                boundBox cBb(allPoints, allCellPoints[celli], false);
+                boundBox cBb(mesh_.cellBb(celli));
                 cBb.min() -= smallVec_;
                 cBb.max() += smallVec_;
 
@@ -295,7 +292,7 @@ void Foam::cellCellStencils::trackingInverseDistance::markPatchesAsHoles
                 forAll(tgtCellMap, tgtCelli)
                 {
                     label celli = tgtCellMap[tgtCelli];
-                    boundBox cBb(allPoints, allCellPoints[celli], false);
+                    boundBox cBb(mesh_.cellBb(celli));
                     cBb.min() -= smallVec_;
                     cBb.max() += smallVec_;
 
@@ -412,7 +409,7 @@ void Foam::cellCellStencils::trackingInverseDistance::markDonors
         label celli = tgtCellMap[tgtCelli];
         if (srcOverlapProcs.size())
         {
-            treeBoundBox subBb(cellBb(mesh_, celli));
+            treeBoundBox subBb(mesh_.cellBb(celli));
             subBb.min() -= smallVec_;
             subBb.max() += smallVec_;
 

src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C

@@ -54,23 +54,20 @@ Foam::labelList Foam::meshToMeshMethod::maskCells() const
 
     intersectBb.inflate(0.01);
 
-    const cellList& srcCells = src_.cells();
-    const faceList& srcFaces = src_.faces();
-    const pointField& srcPts = src_.points();
-
     DynamicList<label> cells(src_.nCells());
-    forAll(srcCells, srcI)
+
+    for (label srcCelli = 0; srcCelli < src_.nCells(); ++srcCelli)
     {
-        boundBox cellBb(srcCells[srcI].points(srcFaces, srcPts), false);
+        boundBox cellBb(src_.cellBb(srcCelli));
         if (intersectBb.overlaps(cellBb))
         {
-            cells.append(srcI);
+            cells.append(srcCelli);
         }
     }
 
     if (debug)
     {
-        Pout<< "participating source mesh cells: " << cells.size() << endl;
+        Pout<< "participating source mesh cells: " << src_.nCells() << endl;
     }
 
     return cells;
@@ -87,14 +84,7 @@ bool Foam::meshToMeshMethod::intersect
 
     tetOverlapVolume overlapEngine;
 
-    // Note: avoid demand-driven construction of cellPoints
+    treeBoundBox bbTgtCell(tgt_.cellBb(tgtCelli));
-    // treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
-    const UList<label>& cellFaces = tgt_.cells()[tgtCelli];
-    treeBoundBox bbTgtCell(tgt_.points(), tgt_.faces()[cellFaces[0]]);
-    for (label i = 1; i < cellFaces.size(); ++i)
-    {
-        bbTgtCell.add(tgt_.points(), tgt_.faces()[cellFaces[i]]);
-    }
 
     return overlapEngine.cellCellOverlapMinDecomp
     (
@@ -116,14 +106,7 @@ Foam::scalar Foam::meshToMeshMethod::interVol
 {
     tetOverlapVolume overlapEngine;
 
-    // Note: avoid demand-driven construction of cellPoints
+    treeBoundBox bbTgtCell(tgt_.cellBb(tgtCelli));
-    // treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
-    const UList<label>& cellFaces = tgt_.cells()[tgtCelli];
-    treeBoundBox bbTgtCell(tgt_.points(), tgt_.faces()[cellFaces[0]]);
-    for (label i = 1; i < cellFaces.size(); ++i)
-    {
-        bbTgtCell.add(tgt_.points(), tgt_.faces()[cellFaces[i]]);
-    }
 
     scalar vol = overlapEngine.cellCellOverlapVolumeMinDecomp
     (
@@ -147,14 +130,7 @@ Foam::meshToMeshMethod::interVolAndCentroid
 {
     tetOverlapVolume overlapEngine;
 
-    // Note: avoid demand-driven construction of cellPoints
+    treeBoundBox bbTgtCell(tgt_.cellBb(tgtCelli));
-    // treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
-    const UList<label>& cellFaces = tgt_.cells()[tgtCelli];
-    treeBoundBox bbTgtCell(tgt_.points(), tgt_.faces()[cellFaces[0]]);
-    for (label i = 1; i < cellFaces.size(); ++i)
-    {
-        bbTgtCell.add(tgt_.points(), tgt_.faces()[cellFaces[i]]);
-    }
 
     Tuple2<scalar, point> volAndInertia =
     overlapEngine.cellCellOverlapMomentMinDecomp

src/sampling/meshToMesh/meshToMeshParallelOps.C

@@ -178,12 +178,8 @@ Foam::autoPtr<Foam::mapDistribute> Foam::meshToMesh::calcProcMap
             }
 
 
-            // determine which cells of tgt mesh overlaps src mesh per proc
+            // Determine which cells of tgt mesh overlaps src mesh per proc
-            const cellList& cells = tgt.cells();
+            labelListList sendMap(Pstream::nProcs());
-            const faceList& faces = tgt.faces();
-            const pointField& points = tgt.points();
-
-            labelListList sendMap;
 
             {
                 // per processor indices into all segments to send
@@ -198,20 +194,11 @@ Foam::autoPtr<Foam::mapDistribute> Foam::meshToMesh::calcProcMap
                 // work array - whether src processor bb overlaps the tgt cell
                 // bounds
                 boolList procBbOverlaps(Pstream::nProcs());
-                forAll(cells, celli)
-                {
-                    const cell& c = cells[celli];
 
-                    // determine bounding box of tgt cell
+                for (label celli = 0; celli < tgt.nCells(); ++celli)
-                    boundBox cellBb(boundBox::invertedBox);
-                    forAll(c, facei)
                 {
-                        const face& f = faces[c[facei]];
+                    // Bounding box of tgt cell
-                        forAll(f, fp)
+                    boundBox cellBb(tgt.cellBb(celli));
-                        {
-                            cellBb.add(points, f);
-                        }
-                    }
 
                     // find the overlapping tgt cells on each src processor
                     (void)calcOverlappingProcs(procBb, cellBb, procBbOverlaps);
@@ -226,23 +213,22 @@ Foam::autoPtr<Foam::mapDistribute> Foam::meshToMesh::calcProcMap
                 }
 
                 // convert dynamicList to labelList
-                sendMap.setSize(Pstream::nProcs());
                 forAll(sendMap, proci)
                 {
                     sendMap[proci].transfer(dynSendMap[proci]);
                 }
-            }
 
                 // debug printing
                 if (debug)
                 {
-                Pout<< "Of my " << cells.size()
+                    Pout<< "Of my " << tgt.nCells()
                         << " target cells I need to send to:" << nl
                         << "\tproc\tcells" << endl;
                     forAll(sendMap, proci)
                     {
-                    Pout<< '\t' << proci << '\t' << sendMap[proci].size()
+                        Pout<< '\t' << proci << '\t'
-                        << endl;
+                            << sendMap[proci].size() << endl;
+                    }
                 }
             }
 
@@ -250,7 +236,7 @@ Foam::autoPtr<Foam::mapDistribute> Foam::meshToMesh::calcProcMap
             // send over how many tgt cells I need to receive from each
             // processor
             labelListList sendSizes(Pstream::nProcs());
-            sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
+            sendSizes[Pstream::myProcNo()].resize(Pstream::nProcs());
             forAll(sendMap, proci)
             {
                 sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size();

src/sampling/surface/cutting/cuttingSurfaceCuts.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -89,20 +89,14 @@ void Foam::cuttingSurface::calcCellCuts
     );
 
 
-    boundBox cellBb;
-
     for (const label celli : cellCuts)
     {
-        cellBb.clear();
+        if (!fvm.cellBb(celli).contains(nearest[celli].point()))
-        cellBb.add(pts, fvm.cellPoints(celli));
-
-        if (!cellBb.contains(nearest[celli].hitPoint()))
         {
             cellCuts.unset(celli);
         }
     }
 
-
     if (debug)
     {
         volScalarField cCuts

src/sampling/surface/distanceSurface/distanceSurface.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -204,23 +204,18 @@ static bitSet simpleGeometricFilter
     // A deny filter. Initially false (accept everything)
     ignoreCells.resize(mesh.nCells());
 
-    bitSet pointFilter;
+    bitSet ignorePoints;
     if (WantPointFilter)
     {
-        // Create as accept filter. Initially false (deny everything)
+        // Create deny filter
-        pointFilter.resize(mesh.nPoints());
+        ignorePoints.resize(mesh.nPoints(), true);
     }
 
-    boundBox cellBb;
-
     forAll(nearest, celli)
     {
         const point& pt = nearest[celli].point();
 
-        const labelList& cPoints = mesh.cellPoints(celli);
+        boundBox cellBb(mesh.cellBb(celli));
-
-        cellBb.clear();
-        cellBb.add(mesh.points(), cPoints);
         cellBb.inflate(boundBoxInflate);
 
         if (!cellBb.contains(pt))
@@ -229,15 +224,12 @@ static bitSet simpleGeometricFilter
         }
         else if (WantPointFilter)
         {
-            // Good cell candidate, accept its points
+            // Good cell candidate, do not ignore its points
-            pointFilter.set(cPoints);
+            ignorePoints.unset(mesh.cellPoints(celli));
         }
     }
 
-    // Flip from accept to deny filter
+    return ignorePoints;
-    pointFilter.flip();
-
-    return pointFilter;
 }