diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C index 57bcc5f761..7e857a1e27 100644 --- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C +++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C @@ -1807,23 +1807,25 @@ Foam::conformalVoronoiMesh::nearestFeatureEdgeLocations const Foam::point& pt ) const { + const auto& tree = edgeLocationTreePtr_(); + const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt); - labelList elems - = edgeLocationTreePtr_().findSphere(pt, exclusionRangeSqr); + labelList elems = tree.findSphere(pt, exclusionRangeSqr); - DynamicList dynPointHit; + DynamicList dynPointHit(elems.size()); - forAll(elems, elemI) + for (const label index : elems) { - label index = elems[elemI]; - - const Foam::point& pointi - = edgeLocationTreePtr_().shapes().shapePoints()[index]; - - pointIndexHit nearHit(true, pointi, index); - - dynPointHit.append(nearHit); + dynPointHit.append + ( + pointIndexHit + ( + true, + tree.shapes().centre(index), + index + ) + ); } return dynPointHit; diff --git a/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C b/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C index 955680be35..d5bc7d3e11 100644 --- a/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C +++ b/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C @@ -235,16 +235,14 @@ public: for (const label index : indices) { - const label edgeIndex = shape.edgeLabels()[index]; + const label edgeIndex = shape.objectIndex(index); if (shapeMask_.found(edgeIndex)) { continue; } - const edge& e = shape.edges()[edgeIndex]; - - pointHit nearHit = e.line(shape.points()).nearestDist(sample); + pointHit nearHit = shape.line(index).nearestDist(sample); // Only register hit if closest point is not an edge point if (nearHit.hit()) diff --git a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.C b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.C index a5914b35a6..ad48ca2c56 100644 --- a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.C +++ b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.C @@ -27,14 +27,13 @@ License \*---------------------------------------------------------------------------*/ #include "dynamicTreeDataPoint.H" -#include "treeBoundBox.H" #include "dynamicIndexedOctree.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { -defineTypeNameAndDebug(dynamicTreeDataPoint, 0); + defineTypeName(dynamicTreeDataPoint); } @@ -51,10 +50,10 @@ Foam::dynamicTreeDataPoint::dynamicTreeDataPoint // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -const Foam::DynamicList& -Foam::dynamicTreeDataPoint::shapePoints() const +Foam::treeBoundBox +Foam::dynamicTreeDataPoint::bounds(const labelUList& indices) const { - return points_; + return treeBoundBox(points_, indices); } @@ -71,10 +70,10 @@ Foam::volumeType Foam::dynamicTreeDataPoint::getVolumeType bool Foam::dynamicTreeDataPoint::overlaps ( const label index, - const treeBoundBox& cubeBb + const treeBoundBox& searchBox ) const { - return cubeBb.contains(points_[index]); + return searchBox.contains(centre(index)); } @@ -85,7 +84,7 @@ bool Foam::dynamicTreeDataPoint::overlaps const scalar radiusSqr ) const { - return (centre.distSqr(points_[index]) <= radiusSqr); + return (centre.distSqr(this->centre(index)) <= radiusSqr); } @@ -99,11 +98,9 @@ void Foam::dynamicTreeDataPoint::findNearest point& nearestPoint ) const { - forAll(indices, i) + for (const label index : indices) { - const label index = indices[i]; - - const point& pt = points_[index]; + const point& pt = centre(index); const scalar distSqr = sample.distSqr(pt); @@ -128,42 +125,30 @@ void Foam::dynamicTreeDataPoint::findNearest point& nearestPoint ) const { + const treeBoundBox lnBb(ln.box()); + // Best so far scalar nearestDistSqr = linePoint.distSqr(nearestPoint); - forAll(indices, i) + for (const label index : indices) { - const label index = indices[i]; + const point& pt = centre(index); - const point& shapePt = points_[index]; - - if (tightest.contains(shapePt)) + if (tightest.contains(pt)) { // Nearest point on line - pointHit pHit = ln.nearestDist(shapePt); - scalar distSqr = sqr(pHit.distance()); + pointHit pHit = ln.nearestDist(pt); + const scalar distSqr = sqr(pHit.distance()); if (distSqr < nearestDistSqr) { nearestDistSqr = distSqr; minIndex = index; linePoint = pHit.point(); - nearestPoint = shapePt; + nearestPoint = pt; - { - point& minPt = tightest.min(); - minPt = min(ln.start(), ln.end()); - minPt.x() -= pHit.distance(); - minPt.y() -= pHit.distance(); - minPt.z() -= pHit.distance(); - } - { - point& maxPt = tightest.max(); - maxPt = max(ln.start(), ln.end()); - maxPt.x() += pHit.distance(); - maxPt.y() += pHit.distance(); - maxPt.z() += pHit.distance(); - } + tightest = lnBb; + tightest.grow(pHit.distance()); } } } diff --git a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.H b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.H index 673f9e14e4..76ef8fca75 100644 --- a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.H +++ b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.H @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2012-2015 OpenFOAM Foundation + Copyright (C) 2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -61,100 +62,126 @@ template class dynamicIndexedOctree; class dynamicTreeDataPoint { - // Private data + // Private Data + //- The point field const DynamicList& points_; + public: - // Declare name of the class and its debug switch - ClassName("dynamicTreeDataPoint"); + // Declare name of the class + ClassNameNoDebug("dynamicTreeDataPoint"); - // Constructors + // Constructors (non-caching) //- Construct from List. Holds reference! - dynamicTreeDataPoint(const DynamicList& points); + explicit dynamicTreeDataPoint(const DynamicList& points); // Member Functions - // Access + //- Object dimension == 0 (point element) + int nDim() const noexcept { return 0; } - inline label size() const - { - return points_.size(); - } - - //- Get representative point cloud for all shapes inside - // (one point per shape) - const DynamicList& shapePoints() const; + //- Return bounding box for the specified point indices + treeBoundBox bounds(const labelUList& indices) const; - // Search + // Access - //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. - // Only makes sense for closed surfaces. - volumeType getVolumeType - ( - const dynamicIndexedOctree&, - const point& - ) const; + //- The original point field + const DynamicList& points() const noexcept { return points_; } - //- Does (bb of) shape at index overlap bb - bool overlaps - ( - const label index, - const treeBoundBox& sampleBb - ) const; + //TDB //- The subset of point ids to use + ///const labelList& pointLabels() const noexcept { labelList::null(); } - //- Check if any point on shape is inside sphere. - bool overlaps - ( - const label index, - const point& centre, - const scalar radiusSqr - ) const; + //- Use a subset of points + bool useSubset() const noexcept { return false; } - //- Calculates nearest (to sample) point in shape. - // Returns actual point and distance (squared) - void findNearest - ( - const labelUList& indices, - const point& sample, + //- Is the effective point field empty? + bool empty() const noexcept { return points_.empty(); } - scalar& nearestDistSqr, - label& nearestIndex, - point& nearestPoint - ) const; + //- The size of the effective point field + label size() const noexcept { return points_.size(); } - //- Calculates nearest (to line) point in shape. - // Returns point and distance (squared) - void findNearest - ( - const labelUList& indices, - const linePointRef& ln, + //- Map from shape index to original (non-subset) point label + label objectIndex(label index) const noexcept { return index; } - treeBoundBox& tightest, - label& minIndex, - point& linePoint, - point& nearestPoint - ) const; + //- Point at specified shape index + const point& operator[](label index) const { return points_[index]; } - //- Calculate intersection of shape with ray. Sets result - // accordingly - bool intersects - ( - const label index, - const point& start, - const point& end, - point& result - ) const - { - NotImplemented; - return false; - } + //- Point at specified shape index + const point& centre(label index) const { return points_[index]; } + //- Representative point cloud + const DynamicList& centres() const noexcept { return points_; } + + + // Search + + //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. + // Only makes sense for closed surfaces. + volumeType getVolumeType + ( + const dynamicIndexedOctree&, + const point& + ) const; + + //- Does (bb of) shape at index overlap bb + bool overlaps + ( + const label index, + const treeBoundBox& searchBox + ) const; + + //- Check if any point on shape is inside sphere. + bool overlaps + ( + const label index, + const point& centre, + const scalar radiusSqr + ) const; + + //- Calculates nearest (to sample) point in shape. + // Returns actual point and distance (squared) + void findNearest + ( + const labelUList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& nearestIndex, + point& nearestPoint + ) const; + + //- Calculates nearest (to line) point in shape. + // Returns point and distance (squared) + void findNearest + ( + const labelUList& indices, + const linePointRef& ln, + + treeBoundBox& tightest, + label& minIndex, + point& linePoint, + point& nearestPoint + ) const; + + //- Calculate intersection of shape with ray. + // Sets result accordingly + bool intersects + ( + const label index, + const point& start, + const point& end, + point& result + ) const + { + NotImplemented; + return false; + } }; diff --git a/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C index d2752db4a9..f405624f6a 100644 --- a/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C +++ b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C @@ -29,44 +29,119 @@ License #include "treeDataCell.H" #include "indexedOctree.H" #include "polyMesh.H" +#include // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { - defineTypeNameAndDebug(treeDataCell, 0); + defineTypeName(treeDataCell); +} + + +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Bound boxes corresponding to specified cells +template +static treeBoundBoxList boxesImpl +( + const primitiveMesh& mesh, + const ElementIds& elemIds +) +{ + treeBoundBoxList bbs(elemIds.size()); + + if (mesh.hasCellPoints()) + { + std::transform + ( + elemIds.cbegin(), + elemIds.cend(), + bbs.begin(), + [&](label celli) + { + return treeBoundBox(mesh.points(), mesh.cellPoints(celli)); + } + ); + } + else + { + std::transform + ( + elemIds.cbegin(), + elemIds.cend(), + bbs.begin(), + [&](label celli) + { + return treeBoundBox(mesh.cells()[celli].box(mesh)); + } + ); + } + + return bbs; +} + +} // End namespace Foam + + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +Foam::treeBoundBoxList +Foam::treeDataCell::boxes(const primitiveMesh& mesh) +{ + // All cells + return boxesImpl(mesh, labelRange(mesh.nCells())); +} + + +Foam::treeBoundBoxList +Foam::treeDataCell::boxes +( + const primitiveMesh& mesh, + const labelUList& cellIds +) +{ + return boxesImpl(mesh, cellIds); +} + + +Foam::treeBoundBox +Foam::treeDataCell::bounds +( + const primitiveMesh& mesh, + const labelUList& cellIds +) +{ + treeBoundBox bb; + + if (mesh.hasCellPoints()) + { + for (const label celli : cellIds) + { + bb.add(mesh.points(), mesh.cellPoints(celli)); + } + } + else + { + for (const label celli : cellIds) + { + bb.add(mesh.cells()[celli].box(mesh)); + } + } + + return bb; } // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -Foam::treeBoundBox Foam::treeDataCell::calcCellBb(const label celli) const +inline Foam::treeBoundBox +Foam::treeDataCell::getBounds(const label index) const { - const cellList& cells = mesh_.cells(); - const faceList& faces = mesh_.faces(); - const pointField& points = mesh_.points(); - - treeBoundBox cellBb - ( - vector(GREAT, GREAT, GREAT), - vector(-GREAT, -GREAT, -GREAT) - ); - - const cell& cFaces = cells[celli]; - - forAll(cFaces, cFacei) - { - const face& f = faces[cFaces[cFacei]]; - - forAll(f, fp) - { - const point& p = points[f[fp]]; - - cellBb.min() = min(cellBb.min(), p); - cellBb.max() = max(cellBb.max(), p); - } - } - return cellBb; + return treeBoundBox(mesh_.cellBb(objectIndex(index))); } @@ -74,11 +149,13 @@ void Foam::treeDataCell::update() { if (cacheBb_) { - bbs_.setSize(cellLabels_.size()); - - forAll(cellLabels_, i) + if (useSubset_) { - bbs_[i] = calcCellBb(cellLabels_[i]); + bbs_ = treeDataCell::boxes(mesh_, cellLabels_); + } + else + { + bbs_ = treeDataCell::boxes(mesh_); } } } @@ -86,6 +163,23 @@ void Foam::treeDataCell::update() // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // +Foam::treeDataCell::treeDataCell +( + const bool cacheBb, + const polyMesh& mesh, + const polyMesh::cellDecomposition decompMode +) +: + mesh_(mesh), + cellLabels_(), + useSubset_(false), + cacheBb_(cacheBb), + decompMode_(decompMode) +{ + update(); +} + + Foam::treeDataCell::treeDataCell ( const bool cacheBb, @@ -96,6 +190,7 @@ Foam::treeDataCell::treeDataCell : mesh_(mesh), cellLabels_(cellLabels), + useSubset_(true), cacheBb_(cacheBb), decompMode_(decompMode) { @@ -113,6 +208,7 @@ Foam::treeDataCell::treeDataCell : mesh_(mesh), cellLabels_(std::move(cellLabels)), + useSubset_(true), cacheBb_(cacheBb), decompMode_(decompMode) { @@ -120,22 +216,79 @@ Foam::treeDataCell::treeDataCell } -Foam::treeDataCell::treeDataCell +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::treeBoundBox +Foam::treeDataCell::bounds(const labelUList& indices) const +{ + if (useSubset_) + { + treeBoundBox bb; + + if (mesh_.hasCellPoints()) + { + for (const label index : indices) + { + const label celli = cellLabels_[index]; + + bb.add(mesh_.points(), mesh_.cellPoints(celli)); + } + } + else + { + for (const label index : indices) + { + const label celli = cellLabels_[index]; + + bb.add(mesh_.cells()[celli].box(mesh_)); + } + } + + return bb; + } + + return treeDataCell::bounds(mesh_, indices); +} + + +Foam::tmp Foam::treeDataCell::centres() const +{ + if (useSubset_) + { + return tmp::New(mesh_.cellCentres(), cellLabels_); + } + + return mesh_.cellCentres(); +} + + +bool Foam::treeDataCell::overlaps ( - const bool cacheBb, - const polyMesh& mesh, - const polyMesh::cellDecomposition decompMode -) -: - mesh_(mesh), - cellLabels_(identity(mesh_.nCells())), - cacheBb_(cacheBb), - decompMode_(decompMode) + const label index, + const treeBoundBox& searchBox +) const { - update(); + return + ( + cacheBb_ + ? searchBox.overlaps(bbs_[index]) + : searchBox.overlaps(getBounds(index)) + ); } +bool Foam::treeDataCell::contains +( + const label index, + const point& sample +) const +{ + return mesh_.pointInCell(sample, objectIndex(index), decompMode_); +} + + +// * * * * * * * * * * * * * * * * Searching * * * * * * * * * * * * * * * * // + Foam::treeDataCell::findNearestOp::findNearestOp ( const indexedOctree& tree @@ -154,43 +307,29 @@ Foam::treeDataCell::findIntersectOp::findIntersectOp {} -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -Foam::pointField Foam::treeDataCell::shapePoints() const -{ - pointField cc(cellLabels_.size()); - - forAll(cellLabels_, i) - { - cc[i] = mesh_.cellCentres()[cellLabels_[i]]; - } - - return cc; -} - - -bool Foam::treeDataCell::overlaps +void Foam::treeDataCell::findNearest ( - const label index, - const treeBoundBox& cubeBb + const labelUList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& minIndex, + point& nearestPoint ) const { - if (cacheBb_) + for (const label index : indices) { - return cubeBb.overlaps(bbs_[index]); + const point& pt = centre(index); + + const scalar distSqr = sample.distSqr(pt); + + if (distSqr < nearestDistSqr) + { + nearestDistSqr = distSqr; + minIndex = index; + nearestPoint = pt; + } } - - return cubeBb.overlaps(calcCellBb(cellLabels_[index])); -} - - -bool Foam::treeDataCell::contains -( - const label index, - const point& sample -) const -{ - return mesh_.pointInCell(sample, cellLabels_[index], decompMode_); } @@ -204,23 +343,14 @@ void Foam::treeDataCell::findNearestOp::operator() point& nearestPoint ) const { - const treeDataCell& shape = tree_.shapes(); - - forAll(indices, i) - { - label index = indices[i]; - label celli = shape.cellLabels()[index]; - const point& pt = shape.mesh().cellCentres()[celli]; - - scalar distSqr = sample.distSqr(pt); - - if (distSqr < nearestDistSqr) - { - nearestDistSqr = distSqr; - minIndex = index; - nearestPoint = pt; - } - } + tree_.shapes().findNearest + ( + indices, + sample, + nearestDistSqr, + minIndex, + nearestPoint + ); } @@ -262,7 +392,7 @@ bool Foam::treeDataCell::findIntersectOp::operator() } else { - const treeBoundBox cellBb = shape.calcCellBb(shape.cellLabels_[index]); + const treeBoundBox cellBb = shape.getBounds(index); if ((cellBb.posBits(start) & cellBb.posBits(end)) != 0) { @@ -276,23 +406,23 @@ bool Foam::treeDataCell::findIntersectOp::operator() // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // Disable picking up intersections behind us. - scalar oldTol = intersection::setPlanarTol(0.0); - - const cell& cFaces = shape.mesh_.cells()[shape.cellLabels_[index]]; + const scalar oldTol = intersection::setPlanarTol(0); const vector dir(end - start); scalar minDistSqr = magSqr(dir); bool hasMin = false; - forAll(cFaces, i) + const label celli = shape.objectIndex(index); + + for (const label facei : shape.mesh().cells()[celli]) { - const face& f = shape.mesh_.faces()[cFaces[i]]; + const face& f = shape.mesh().faces()[facei]; pointHit inter = f.ray ( start, dir, - shape.mesh_.points(), + shape.mesh().points(), intersection::HALF_RAY ); diff --git a/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H index a726408c60..e09bebdf62 100644 --- a/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H +++ b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation + Copyright (C) 2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -35,8 +36,8 @@ SourceFiles \*---------------------------------------------------------------------------*/ -#ifndef treeDataCell_H -#define treeDataCell_H +#ifndef Foam_treeDataCell_H +#define Foam_treeDataCell_H #include "polyMesh.H" #include "treeBoundBoxList.H" @@ -47,7 +48,7 @@ SourceFiles namespace Foam { -// Forward declaration of classes +// Forward Declarations template class indexedOctree; /*---------------------------------------------------------------------------*\ @@ -56,13 +57,17 @@ template class indexedOctree; class treeDataCell { - // Private data + // Private Data + //- Reference to the mesh const polyMesh& mesh_; //- Subset of cells to work on const labelList cellLabels_; + //- Use subset of cells (cellLabels) + const bool useSubset_; + //- Whether to precalculate and store cell bounding box const bool cacheBb_; @@ -75,8 +80,8 @@ class treeDataCell // Private Member Functions - //- Calculate cell bounding box - treeBoundBox calcCellBb(const label celli) const; + //- Get cell bounding box at specified index + inline treeBoundBox getBounds(const label index) const; //- Initialise all member data void update(); @@ -133,11 +138,19 @@ public: }; - // Declare name of the class and its debug switch - ClassName("treeDataCell"); + // Declare name of the class + ClassNameNoDebug("treeDataCell"); - // Constructors + // Constructors (cachable versions) + + //- Construct from mesh, using all cells in mesh. + treeDataCell + ( + const bool cacheBb, + const polyMesh& mesh, + polyMesh::cellDecomposition decompMode + ); //- Construct from mesh, copying subset of cells. treeDataCell @@ -145,7 +158,7 @@ public: const bool cacheBb, const polyMesh& mesh, const labelUList& cellLabels, - const polyMesh::cellDecomposition decompMode + polyMesh::cellDecomposition decompMode ); //- Construct from mesh, moving subset of cells @@ -154,74 +167,166 @@ public: const bool cacheBb, const polyMesh& mesh, labelList&& cellLabels, - const polyMesh::cellDecomposition decompMode + polyMesh::cellDecomposition decompMode ); - //- Construct from mesh. Uses all cells in mesh. + + // Constructors (non-caching) + + //- Construct from mesh, using all cells in mesh. treeDataCell ( - const bool cacheBb, const polyMesh& mesh, - const polyMesh::cellDecomposition decompMode + polyMesh::cellDecomposition decompMode + ) + : + treeDataCell(false, mesh, decompMode) + {} + + //- Construct from mesh, copying subset of cells. + treeDataCell + ( + const polyMesh& mesh, + const labelUList& cellLabels, + polyMesh::cellDecomposition decompMode + ) + : + treeDataCell(false, mesh, cellLabels, decompMode) + {} + + //- Construct from mesh, moving subset of cells + treeDataCell + ( + const polyMesh& mesh, + labelList&& cellLabels, + polyMesh::cellDecomposition decompMode + ) + : + treeDataCell(false, mesh, std::move(cellLabels), decompMode) + {} + + + // Static Functions + + //- Calculate and return bounding boxes for all mesh cells + static treeBoundBoxList boxes + ( + const primitiveMesh& mesh + ); + + //- Calculate and return bounding boxes for specified mesh cells + static treeBoundBoxList boxes + ( + const primitiveMesh& mesh, + const labelUList& cellIds + ); + + //- Return bounding box of specified mesh cells + static treeBoundBox bounds + ( + const primitiveMesh& mesh, + const labelUList& cellIds ); // Member Functions - // Access + //- Object dimension == 3 (volume element) + int nDim() const noexcept { return 3; } - inline const labelList& cellLabels() const - { - return cellLabels_; - } - - inline const polyMesh& mesh() const - { - return mesh_; - } - - inline polyMesh::cellDecomposition decompMode() const - { - return decompMode_; - } - - inline label size() const - { - return cellLabels_.size(); - } - - //- Get representative point cloud for all shapes inside - // (one point per shape) - pointField shapePoints() const; + //- Return bounding box for the specified cell indices + treeBoundBox bounds(const labelUList& indices) const; - // Search + // Access - //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. - // Only makes sense for closed surfaces. - volumeType getVolumeType - ( - const indexedOctree&, - const point& - ) const - { - NotImplemented; - return volumeType::UNKNOWN; - } + //- Reference to the supporting mesh + const polyMesh& mesh() const noexcept { return mesh_; } - //- Does (bb of) shape at index overlap bb - bool overlaps - ( - const label index, - const treeBoundBox& sampleBb - ) const; + //- The cell decomposition mode used + polyMesh::cellDecomposition decompMode() const noexcept + { + return decompMode_; + } - //- Does shape at index contain sample - bool contains - ( - const label index, - const point& sample - ) const; + //- The subset of cell ids to use + const labelList& cellLabels() const noexcept { return cellLabels_; } + + //- Use a subset of cells + bool useSubset() const noexcept { return useSubset_; } + + //- Is the effective cell selection empty? + bool empty() const noexcept + { + return useSubset_ ? cellLabels_.empty() : !mesh_.nCells(); + } + + //- The size of the cell selection + label size() const noexcept + { + return useSubset_ ? cellLabels_.size() : mesh_.nCells(); + } + + //- Map from shape index to original (non-subset) cell label + label objectIndex(const label index) const + { + return useSubset_ && index >= 0 ? cellLabels_[index] : index; + } + + //TBD //- Cell at specified shape index + ///const cell& operator[](const label index) const + ///{ return mesh_.cells()[objectIndex(index)]; } + + //- Representative point (cell centre) at shape index + const point& centre(const label index) const + { + return mesh_.cellCentres()[objectIndex(index)]; + } + + //- Representative point cloud for contained shapes. + //- One point per shape, corresponding to the cell centres. + tmp centres() const; + + + // Search + + //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. + // Only makes sense for closed surfaces. + volumeType getVolumeType + ( + const indexedOctree&, + const point& + ) const + { + NotImplemented; + return volumeType::UNKNOWN; + } + + //- Does (bb of) shape at index overlap searchBox + bool overlaps + ( + const label index, + const treeBoundBox& searchBox + ) const; + + //- Does shape at index contain sample + bool contains + ( + const label index, + const point& sample + ) const; + + //- Calculates nearest (to sample) point in shape. + // Returns actual point and distance (squared) + void findNearest + ( + const labelUList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& nearestIndex, + point& nearestPoint + ) const; }; diff --git a/src/OpenFOAM/algorithms/indexedOctree/treeDataEdge.C b/src/OpenFOAM/algorithms/indexedOctree/treeDataEdge.C index fb069c7d4b..64474327b3 100644 --- a/src/OpenFOAM/algorithms/indexedOctree/treeDataEdge.C +++ b/src/OpenFOAM/algorithms/indexedOctree/treeDataEdge.C @@ -28,36 +28,173 @@ License #include "treeDataEdge.H" #include "indexedOctree.H" +#include // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { - defineTypeNameAndDebug(treeDataEdge, 0); + defineTypeName(treeDataEdge); +} + + +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Bound boxes corresponding to specified edges +template +static treeBoundBoxList boxesImpl +( + const edgeList& edges, + const pointField& points, + const ElementIds& elemIds +) +{ + treeBoundBoxList bbs(elemIds.size()); + + std::transform + ( + elemIds.cbegin(), + elemIds.cend(), + bbs.begin(), + [&](label edgei) + { + return treeBoundBox(edges[edgei].box(points)); + } + ); + + return bbs; +} + + +// Overall bound box for specified edges +template +static treeBoundBox boundsImpl +( + const edgeList& edges, + const pointField& points, + const ElementIds& elemIds +) +{ + treeBoundBox bb; + + for (const label edgei : elemIds) + { + const edge& e = edges[edgei]; + + bb.add(points[e.first()], points[e.second()]); + } + + return bb; +} + +} // End namespace Foam + + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +Foam::treeBoundBoxList +Foam::treeDataEdge::boxes +( + const edgeList& edges, + const pointField& points +) +{ + treeBoundBoxList bbs(edges.size()); + + std::transform + ( + edges.cbegin(), + edges.cend(), + bbs.begin(), + [&](const edge& e) { return treeBoundBox(e.box(points)); } + ); + + return bbs; +} + + +Foam::treeBoundBoxList +Foam::treeDataEdge::boxes +( + const edgeList& edges, + const pointField& points, + const labelRange& range +) +{ + return boxesImpl(edges, points, range); +} + + +Foam::treeBoundBoxList +Foam::treeDataEdge::boxes +( + const edgeList& edges, + const pointField& points, + const labelUList& edgeIds +) +{ + return boxesImpl(edges, points, edgeIds); +} + + +Foam::treeBoundBox +Foam::treeDataEdge::bounds +( + const edgeList& edges, + const pointField& points +) +{ + treeBoundBox bb; + + for (const edge& e : edges) + { + bb.add(points[e.first()], points[e.second()]); + } + + return bb; +} + + +Foam::treeBoundBox +Foam::treeDataEdge::bounds +( + const edgeList& edges, + const pointField& points, + const labelRange& range +) +{ + return boundsImpl(edges, points, range); +} + + +Foam::treeBoundBox +Foam::treeDataEdge::bounds +( + const edgeList& edges, + const pointField& points, + const labelUList& edgeIds +) +{ + return boundsImpl(edges, points, edgeIds); } // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -Foam::treeBoundBox Foam::treeDataEdge::calcBb(const label edgeI) const -{ - const edge& e = edges_[edgeI]; - const point& p0 = points_[e[0]]; - const point& p1 = points_[e[1]]; - - return treeBoundBox(min(p0, p1), max(p0, p1)); -} - - void Foam::treeDataEdge::update() { if (cacheBb_) { - bbs_.setSize(edgeLabels_.size()); - - forAll(edgeLabels_, i) + if (useSubset_) { - bbs_[i] = calcBb(edgeLabels_[i]); + bbs_ = treeDataEdge::boxes(edges_, points_, edgeLabels_); + } + else + { + bbs_ = treeDataEdge::boxes(edges_, points_); } } } @@ -65,6 +202,41 @@ void Foam::treeDataEdge::update() // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // +Foam::treeDataEdge::treeDataEdge +( + const bool cacheBb, + const edgeList& edges, + const pointField& points +) +: + points_(points), + edges_(edges), + edgeLabels_(), + useSubset_(false), + cacheBb_(cacheBb) +{ + update(); +} + + +Foam::treeDataEdge::treeDataEdge +( + const bool cacheBb, + const edgeList& edges, + const pointField& points, + const labelRange& range +) +: + points_(points), + edges_(edges), + edgeLabels_(identity(range)), + useSubset_(true), + cacheBb_(cacheBb) +{ + update(); +} + + Foam::treeDataEdge::treeDataEdge ( const bool cacheBb, @@ -73,9 +245,10 @@ Foam::treeDataEdge::treeDataEdge const labelUList& edgeLabels ) : - edges_(edges), points_(points), + edges_(edges), edgeLabels_(edgeLabels), + useSubset_(true), cacheBb_(cacheBb) { update(); @@ -90,15 +263,107 @@ Foam::treeDataEdge::treeDataEdge labelList&& edgeLabels ) : - edges_(edges), points_(points), + edges_(edges), edgeLabels_(std::move(edgeLabels)), + useSubset_(true), cacheBb_(cacheBb) { update(); } +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::treeBoundBox Foam::treeDataEdge::bounds(const labelUList& indices) const +{ + if (useSubset_) + { + treeBoundBox bb; + + for (const label index : indices) + { + const edge& e = edges_[edgeLabels_[index]]; + + bb.add(points_[e.first()], points_[e.second()]); + } + + return bb; + } + + return treeDataEdge::bounds(edges_, points_, indices); +} + + +Foam::tmp Foam::treeDataEdge::centres() const +{ + tmp tpts; + + if (useSubset_) + { + tpts = tmp::New(edgeLabels_.size()); + + std::transform + ( + edgeLabels_.cbegin(), + edgeLabels_.cend(), + tpts.ref().begin(), + [&](label edgei) { return edges_[edgei].centre(points_); } + ); + } + else + { + tpts = tmp::New(edges_.size()); + + std::transform + ( + edges_.cbegin(), + edges_.cend(), + tpts.ref().begin(), + [&](const edge& e) { return e.centre(points_); } + ); + } + + return tpts; +} + + +Foam::volumeType Foam::treeDataEdge::getVolumeType +( + const indexedOctree& oc, + const point& sample +) const +{ + return volumeType::UNKNOWN; +} + + +bool Foam::treeDataEdge::overlaps +( + const label index, + const treeBoundBox& searchBox +) const +{ + point intersect; + return searchBox.intersects(this->line(index), intersect); +} + + +bool Foam::treeDataEdge::overlaps +( + const label index, + const point& centre, + const scalar radiusSqr +) const +{ + const pointHit nearHit = this->line(index).nearestDist(centre); + + return (sqr(nearHit.distance()) <= radiusSqr); +} + + +// * * * * * * * * * * * * * * * * Searching * * * * * * * * * * * * * * * * // + Foam::treeDataEdge::findNearestOp::findNearestOp ( const indexedOctree& tree @@ -115,67 +380,7 @@ Foam::treeDataEdge::findIntersectOp::findIntersectOp {} -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -Foam::pointField Foam::treeDataEdge::shapePoints() const -{ - pointField eMids(edgeLabels_.size()); - - forAll(edgeLabels_, i) - { - const edge& e = edges_[edgeLabels_[i]]; - - eMids[i] = e.centre(points_); - } - return eMids; -} - - -Foam::volumeType Foam::treeDataEdge::getVolumeType -( - const indexedOctree& oc, - const point& sample -) const -{ - return volumeType::UNKNOWN; -} - - -bool Foam::treeDataEdge::overlaps -( - const label index, - const treeBoundBox& cubeBb -) const -{ - const edge& e = edges_[edgeLabels_[index]]; - - const point& start = points_[e.start()]; - const point& end = points_[e.end()]; - - point intersect; - - return cubeBb.intersects(start, end, intersect); -} - - -bool Foam::treeDataEdge::overlaps -( - const label index, - const point& centre, - const scalar radiusSqr -) const -{ - const edge& e = edges_[edgeLabels_[index]]; - - const pointHit nearHit = e.line(points_).nearestDist(centre); - - const scalar distSqr = sqr(nearHit.distance()); - - return (distSqr <= radiusSqr); -} - - -void Foam::treeDataEdge::findNearestOp::operator() +void Foam::treeDataEdge::findNearest ( const labelUList& indices, const point& sample, @@ -185,13 +390,9 @@ void Foam::treeDataEdge::findNearestOp::operator() point& nearestPoint ) const { - const treeDataEdge& shape = tree_.shapes(); - for (const label index : indices) { - const edge& e = shape.edges()[shape.edgeLabels()[index]]; - - pointHit nearHit = e.line(shape.points()).nearestDist(sample); + pointHit nearHit = this->line(index).nearestDist(sample); const scalar distSqr = sqr(nearHit.distance()); @@ -205,6 +406,27 @@ void Foam::treeDataEdge::findNearestOp::operator() } +void Foam::treeDataEdge::findNearestOp::operator() +( + const labelUList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& minIndex, + point& nearestPoint +) const +{ + tree_.shapes().findNearest + ( + indices, + sample, + nearestDistSqr, + minIndex, + nearestPoint + ); +} + + void Foam::treeDataEdge::findNearestOp::operator() ( const labelUList& indices, @@ -218,19 +440,20 @@ void Foam::treeDataEdge::findNearestOp::operator() { const treeDataEdge& shape = tree_.shapes(); + const treeBoundBox lnBb(ln.box()); + // Best so far scalar nearestDistSqr = linePoint.distSqr(nearestPoint); for (const label index : indices) { - const edge& e = shape.edges()[shape.edgeLabels()[index]]; - // Note: could do bb test ? Worthwhile? // Nearest point on line point ePoint, lnPt; - scalar dist = e.line(shape.points()).nearestDist(ln, ePoint, lnPt); - scalar distSqr = sqr(dist); + const scalar dist = shape.line(index).nearestDist(ln, ePoint, lnPt); + + const scalar distSqr = sqr(dist); if (distSqr < nearestDistSqr) { @@ -239,20 +462,8 @@ void Foam::treeDataEdge::findNearestOp::operator() linePoint = lnPt; nearestPoint = ePoint; - { - point& minPt = tightest.min(); - minPt = min(ln.start(), ln.end()); - minPt.x() -= dist; - minPt.y() -= dist; - minPt.z() -= dist; - } - { - point& maxPt = tightest.max(); - maxPt = max(ln.start(), ln.end()); - maxPt.x() += dist; - maxPt.y() += dist; - maxPt.z() += dist; - } + tightest = lnBb; + tightest.grow(dist); } } } diff --git a/src/OpenFOAM/algorithms/indexedOctree/treeDataEdge.H b/src/OpenFOAM/algorithms/indexedOctree/treeDataEdge.H index 04cbe30612..1f490df6d8 100644 --- a/src/OpenFOAM/algorithms/indexedOctree/treeDataEdge.H +++ b/src/OpenFOAM/algorithms/indexedOctree/treeDataEdge.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2015 OpenFOAM Foundation - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -39,6 +39,7 @@ SourceFiles #define Foam_treeDataEdge_H #include "treeBoundBoxList.H" +#include "labelRange.H" #include "line.H" #include "volumeType.H" @@ -56,41 +57,37 @@ template class indexedOctree; class treeDataEdge { - // Static Data - - //- Tolerance on linear dimensions - static scalar tol; - - // Private Data + //- Reference to the underlying support points + const pointField& points_; + //- Reference to edgeList const edgeList& edges_; - //- Reference to points - const pointField& points_; - - //- Labels of edges + //- Subset of edges to work on const labelList edgeLabels_; - //- Whether to precalculate and store face bounding box + //- Use subset of edges (edgeLabels) + const bool useSubset_; + + //- Whether to precalculate and store edge bounding box const bool cacheBb_; - //- Bbs for all above edges (valid only if cacheBb_) + //- Edge bounding boxes (valid only if cacheBb_) treeBoundBoxList bbs_; // Private Member Functions - //- Calculate edge bounding box - treeBoundBox calcBb(const label edgeI) const; - //- Initialise all member data void update(); + public: + //- Forward to treeDataEdge findNearest operations class findNearestOp { const indexedOctree& tree_; @@ -122,6 +119,7 @@ public: }; + //- Forward to treeDataEdge findIntersect operations class findIntersectOp { public: @@ -141,10 +139,29 @@ public: // Declare name of the class and its debug switch - ClassName("treeDataEdge"); + ClassNameNoDebug("treeDataEdge"); - // Constructors + // Constructors (cachable versions) + + //- Construct from all edges. + // \note Holds references to edges and points! + treeDataEdge + ( + const bool cacheBb, + const edgeList& edges, + const pointField& points + ); + + //- Construct from range of edges. + // \note Holds references to edges and points! + treeDataEdge + ( + const bool cacheBb, + const edgeList& edges, + const pointField& points, + const labelRange& range + ); //- Construct from selected edges. // \note Holds references to edges and points! @@ -167,59 +184,201 @@ public: ); + // Constructors (non-caching) + + //- Construct from all edges. + // \note Holds references to edges and points! + treeDataEdge(const edgeList& edges, const pointField& points) + : + treeDataEdge(false, edges, points) + {} + + //- Construct from range of edges. + // \note Holds references to edges and points! + treeDataEdge + ( + const edgeList& edges, + const pointField& points, + const labelRange& range + ) + : + treeDataEdge(false, edges, points, range) + {} + + //- Construct from selected edges. + // \note Holds references to edges and points! + treeDataEdge + ( + const edgeList& edges, + const pointField& points, + const labelUList& edgeLabels + ) + : + treeDataEdge(false, edges, points, edgeLabels) + {} + + //- Construct from selected edges, transferring contents. + // \note Holds references to edges and points! + treeDataEdge + ( + const edgeList& edges, + const pointField& points, + labelList&& edgeLabels + ) + : + treeDataEdge(false, edges, points, std::move(edgeLabels)) + {} + + + // Static Functions + + //- Calculate and return bounding boxes for all edges + static treeBoundBoxList boxes + ( + const edgeList& edges, + const pointField& points + ); + + //- Calculate and return bounding boxes for specified range of edges + static treeBoundBoxList boxes + ( + const edgeList& edges, + const pointField& points, + const labelRange& range + ); + + //- Calculate and return bounding boxes for specified edges + static treeBoundBoxList boxes + ( + const edgeList& edges, + const pointField& points, + const labelUList& edgeIds + ); + + //- Return bounding box of all edges + static treeBoundBox bounds + ( + const edgeList& edges, + const pointField& points + ); + + //- Return bounding box of specified range of edges + static treeBoundBox bounds + ( + const edgeList& edges, + const pointField& points, + const labelRange& range + ); + + //- Return bounding box of specified edges + static treeBoundBox bounds + ( + const edgeList& edges, + const pointField& points, + const labelUList& edgeIds + ); + + // Member Functions - // Access + //- Object dimension == 1 (line element) + int nDim() const noexcept { return 1; } - const edgeList& edges() const - { - return edges_; - } - - const pointField& points() const - { - return points_; - } - - const labelList& edgeLabels() const - { - return edgeLabels_; - } - - label size() const - { - return edgeLabels_.size(); - } - - //- Representative point cloud for all shapes inside - //- (one point per shape) - pointField shapePoints() const; + //- Return bounding box for the specified edge indices + treeBoundBox bounds(const labelUList& indices) const; - // Search + // Access - //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. - // Only makes sense for closed surfaces. - volumeType getVolumeType - ( - const indexedOctree&, - const point& - ) const; + //- The reference point field + const pointField& points() const noexcept { return points_; } - //- Does (bb of) shape at index overlap bb - bool overlaps - ( - const label index, - const treeBoundBox& sampleBb - ) const; + //- The original list of edges + const edgeList& edges() const noexcept { return edges_; } - //- Does (bb of) shape at index overlap bb - bool overlaps - ( - const label index, - const point& centre, - const scalar radiusSqr - ) const; + //- The subset of edge ids to use + const labelList& edgeLabels() const noexcept { return edgeLabels_; } + + //- Use a subset of edges + bool useSubset() const noexcept { return useSubset_; } + + //- Is the effective edge selection empty? + bool empty() const noexcept + { + return useSubset_ ? edgeLabels_.empty() : edges_.empty(); + } + + //- The size of edge selection + label size() const noexcept + { + return useSubset_ ? edgeLabels_.size() : edges_.size(); + } + + //- Map from shape index to original (non-subset) edge label + label objectIndex(const label index) const + { + return useSubset_ && index >= 0 ? edgeLabels_[index] : index; + } + + //- Edge at specified shape index + const edge& operator[](const label index) const + { + return edges_[objectIndex(index)]; + } + + //- Geometric line for edge at specified shape index. Frequently used + const linePointRef line(const label index) const + { + return edges_[objectIndex(index)].line(points_); + } + + //- Representative point (edge centre) at shape index + point centre(const label index) const + { + return edges_[objectIndex(index)].centre(points_); + } + + //- Representative point cloud for contained shapes. + //- One point per shape, corresponding to the edge centres. + tmp centres() const; + + + // Search + + //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. + // Only makes sense for closed surfaces. + volumeType getVolumeType + ( + const indexedOctree&, + const point& + ) const; + + //- Does (bb of) shape at index searchBox + bool overlaps + ( + const label index, + const treeBoundBox& searchBox + ) const; + + //- Does (bb of) shape at index overlap bb + bool overlaps + ( + const label index, + const point& centre, + const scalar radiusSqr + ) const; + + //- Calculates nearest (to sample) point in shape. + // Returns actual point and distance (squared) + void findNearest + ( + const labelUList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& nearestIndex, + point& nearestPoint + ) const; }; diff --git a/src/OpenFOAM/algorithms/indexedOctree/treeDataPoint.C b/src/OpenFOAM/algorithms/indexedOctree/treeDataPoint.C index d2329f6058..991c560f23 100644 --- a/src/OpenFOAM/algorithms/indexedOctree/treeDataPoint.C +++ b/src/OpenFOAM/algorithms/indexedOctree/treeDataPoint.C @@ -33,7 +33,7 @@ License namespace Foam { - defineTypeNameAndDebug(treeDataPoint, 0); + defineTypeName(treeDataPoint); } @@ -72,29 +72,29 @@ Foam::treeDataPoint::treeDataPoint {} -Foam::treeDataPoint::findNearestOp::findNearestOp -( - const indexedOctree& tree -) -: - tree_(tree) -{} - - -Foam::treeDataPoint::findIntersectOp::findIntersectOp -( - const indexedOctree& tree -) -{} - - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::pointField Foam::treeDataPoint::shapePoints() const +Foam::treeBoundBox Foam::treeDataPoint::bounds(const labelUList& indices) const { if (useSubset_) { - return pointField(points_, pointLabels_); + treeBoundBox bb; + for (const label index : indices) + { + bb.add(points_[pointLabels_[index]]); + } + return bb; + } + + return treeBoundBox(points_, indices); +} + + +Foam::tmp Foam::treeDataPoint::centres() const +{ + if (useSubset_) + { + return tmp::New(points_, pointLabels_); } return points_; @@ -114,10 +114,10 @@ Foam::volumeType Foam::treeDataPoint::getVolumeType bool Foam::treeDataPoint::overlaps ( const label index, - const treeBoundBox& cubeBb + const treeBoundBox& searchBox ) const { - return cubeBb.contains(shapePoint(index)); + return searchBox.contains(centre(index)); } @@ -128,11 +128,29 @@ bool Foam::treeDataPoint::overlaps const scalar radiusSqr ) const { - return (centre.distSqr(shapePoint(index)) <= radiusSqr); + return (centre.distSqr(this->centre(index)) <= radiusSqr); } -void Foam::treeDataPoint::findNearestOp::operator() +// * * * * * * * * * * * * * * * * Searching * * * * * * * * * * * * * * * * // + +Foam::treeDataPoint::findNearestOp::findNearestOp +( + const indexedOctree& tree +) +: + tree_(tree) +{} + + +Foam::treeDataPoint::findIntersectOp::findIntersectOp +( + const indexedOctree& tree +) +{} + + +void Foam::treeDataPoint::findNearest ( const labelUList& indices, const point& sample, @@ -142,11 +160,9 @@ void Foam::treeDataPoint::findNearestOp::operator() point& nearestPoint ) const { - const treeDataPoint& shape = tree_.shapes(); - for (const label index : indices) { - const point& pt = shape.shapePoint(index); + const point& pt = centre(index); const scalar distSqr = sample.distSqr(pt); @@ -160,6 +176,27 @@ void Foam::treeDataPoint::findNearestOp::operator() } +void Foam::treeDataPoint::findNearestOp::operator() +( + const labelUList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& minIndex, + point& nearestPoint +) const +{ + tree_.shapes().findNearest + ( + indices, + sample, + nearestDistSqr, + minIndex, + nearestPoint + ); +} + + void Foam::treeDataPoint::findNearestOp::operator() ( const labelUList& indices, @@ -173,6 +210,8 @@ void Foam::treeDataPoint::findNearestOp::operator() { const treeDataPoint& shape = tree_.shapes(); + const treeBoundBox lnBb(ln.box()); + // Best so far scalar nearestDistSqr = GREAT; if (minIndex >= 0) @@ -182,12 +221,12 @@ void Foam::treeDataPoint::findNearestOp::operator() for (const label index : indices) { - const point& shapePt = shape.shapePoint(index); + const point& pt = shape.centre(index); - if (tightest.contains(shapePt)) + if (tightest.contains(pt)) { // Nearest point on line - pointHit pHit = ln.nearestDist(shapePt); + pointHit pHit = ln.nearestDist(pt); const scalar distSqr = sqr(pHit.distance()); if (distSqr < nearestDistSqr) @@ -195,22 +234,10 @@ void Foam::treeDataPoint::findNearestOp::operator() nearestDistSqr = distSqr; minIndex = index; linePoint = pHit.point(); - nearestPoint = shapePt; + nearestPoint = pt; - { - point& minPt = tightest.min(); - minPt = min(ln.start(), ln.end()); - minPt.x() -= pHit.distance(); - minPt.y() -= pHit.distance(); - minPt.z() -= pHit.distance(); - } - { - point& maxPt = tightest.max(); - maxPt = max(ln.start(), ln.end()); - maxPt.x() += pHit.distance(); - maxPt.y() += pHit.distance(); - maxPt.z() += pHit.distance(); - } + tightest = lnBb; + tightest.grow(pHit.distance()); } } } diff --git a/src/OpenFOAM/algorithms/indexedOctree/treeDataPoint.H b/src/OpenFOAM/algorithms/indexedOctree/treeDataPoint.H index e368cc84fa..391e770239 100644 --- a/src/OpenFOAM/algorithms/indexedOctree/treeDataPoint.H +++ b/src/OpenFOAM/algorithms/indexedOctree/treeDataPoint.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2013 OpenFOAM Foundation - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -64,16 +64,19 @@ class treeDataPoint { // Private Data + //- Reference to the underlying point field const pointField& points_; //- Subset of points to work on (or empty) const labelList pointLabels_; + //- Use subset of points (pointLabels) const bool useSubset_; + public: - + //- Forward to treeDataPoint findNearest operations class findNearestOp { const indexedOctree& tree_; @@ -105,6 +108,7 @@ public: }; + //- Forward to treeDataPoint findIntersect operations (not possible) class findIntersectOp { public: @@ -123,11 +127,11 @@ public: }; - // Declare name of the class and its debug switch - ClassName("treeDataPoint"); + // Declare name of the class + ClassNameNoDebug("treeDataPoint"); - // Constructors + // Constructors (non-caching) //- Construct from pointField // \note Holds reference to the points! @@ -154,108 +158,101 @@ public: // Member Functions - // Access + //- Object dimension == 0 (point element) + int nDim() const noexcept { return 0; } - //- An empty effective point field? - inline bool empty() const - { - return - ( - useSubset_ - ? pointLabels_.empty() - : points_.empty() - ); - } - - //- The effective point field size - inline label size() const - { - return - ( - useSubset_ - ? pointLabels_.size() - : points_.size() - ); - } - - //- The original point field - inline const pointField& points() const - { - return points_; - } - - //- The original point ids - inline const labelList& pointLabels() const - { - return pointLabels_; - } - - //- Use a subset of points - inline bool useSubset() const - { - return useSubset_; - } - - //- The original (non-subset) point label - inline label pointLabel(const label index) const - { - return - ( - useSubset_ && index >= 0 - ? pointLabels_[index] - : index - ); - } - - //- Point at specified index - inline const point& shapePoint(const label index) const - { - return - ( - useSubset_ - ? points_[pointLabels_[index]] - : points_[index] - ); - } - - //- Representative point cloud for all shapes inside - //- (one point per shape) - pointField shapePoints() const; + //- Return bounding box for the specified point indices + treeBoundBox bounds(const labelUList& indices) const; - // Search + // Access - //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. - // Only makes sense for closed surfaces. - volumeType getVolumeType - ( - const indexedOctree& os, - const point& sample - ) const; + //- The original point field + const pointField& points() const noexcept { return points_; } - //- Does (bb of) shape at index overlap bb - bool overlaps - ( - const label index, - const treeBoundBox& sampleBb - ) const; + //- The subset of point ids to use + const labelList& pointLabels() const noexcept { return pointLabels_; } - //- Does shape at index overlap the sphere - bool overlaps - ( - const label index, - const point& centre, - const scalar radiusSqr - ) const; + //- Use a subset of points + bool useSubset() const noexcept { return useSubset_; } - - // Member Operators - - //- The point at the specified index - inline const point& operator[](const label index) const + //- Is the effective point field empty? + bool empty() const noexcept { - return shapePoint(index); + return useSubset_ ? pointLabels_.empty() : points_.empty(); } + + //- The size of the effective point field + label size() const noexcept + { + return useSubset_ ? pointLabels_.size() : points_.size(); + } + + //- Map to the original (non-subset) point label + label objectIndex(const label index) const + { + return useSubset_ && index >= 0 ? pointLabels_[index] : index; + } + + //- Point at specified shape index + const point& operator[](const label index) const + { + return points_[objectIndex(index)]; + } + + //- Point at specified shape index + const point& centre(const label index) const + { + return points_[objectIndex(index)]; + } + + //- Point cloud + tmp centres() const; + + + // Search + + //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. + // Only makes sense for closed surfaces. + volumeType getVolumeType + ( + const indexedOctree& os, + const point& sample + ) const; + + //- Does (bb of) shape at index searchBox + bool overlaps + ( + const label index, + const treeBoundBox& searchBox + ) const; + + //- Does shape at index overlap the sphere + bool overlaps + ( + const label index, + const point& centre, + const scalar radiusSqr + ) const; + + //- Calculates nearest (to sample) point in shape. + // Returns actual point and distance (squared) + void findNearest + ( + const labelUList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& nearestIndex, + point& nearestPoint + ) const; + + + // Housekeeping + + //- Map to the original (non-subset) point label + ///FOAM_DEPRECATED_FOR(2022-10, "objectIndex()") + label pointLabel(label index) const { return objectIndex(index); } }; diff --git a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C index 0e96c24080..477a492518 100644 --- a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C +++ b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C @@ -1039,7 +1039,7 @@ void Foam::faceCoupleInfo::findSlavesCoveringMaster if (nearInfo.hit()) { - label mesh0Facei = tree.shapes().faceLabels()[nearInfo.index()]; + label mesh0Facei = tree.shapes().objectIndex(nearInfo.index()); // Check if points of f1 actually lie on top of mesh0 face // This is the bit that might fail since if f0 is severely warped diff --git a/src/functionObjects/forces/propellerInfo/propellerInfo.C b/src/functionObjects/forces/propellerInfo/propellerInfo.C index 8ad3146a3b..7f233b3f25 100644 --- a/src/functionObjects/forces/propellerInfo/propellerInfo.C +++ b/src/functionObjects/forces/propellerInfo/propellerInfo.C @@ -435,7 +435,7 @@ void Foam::functionObjects::propellerInfo::updateSampleDiskCells() // Kick the tet base points calculation to avoid parallel comms later (void)mesh_.tetBasePtIs(); - const auto& cellLabels = tree.shapes().cellLabels(); + const auto& treeData = tree.shapes(); forAll(points_, pointi) { @@ -453,9 +453,11 @@ void Foam::functionObjects::propellerInfo::updateSampleDiskCells() if (proci >= 0) { cellIds_[pointi] = + ( proci == Pstream::myProcNo() - ? cellLabels[treeCelli] - : -1; + ? treeData.objectIndex(treeCelli) + : -1 + ); } else { diff --git a/src/lagrangian/basic/InteractionLists/InteractionLists.C b/src/lagrangian/basic/InteractionLists/InteractionLists.C index 352f51224c..5084b4956e 100644 --- a/src/lagrangian/basic/InteractionLists/InteractionLists.C +++ b/src/lagrangian/basic/InteractionLists/InteractionLists.C @@ -67,19 +67,7 @@ void Foam::InteractionLists::buildInteractionLists() extendedProcBbsOrigProc ); - treeBoundBoxList cellBbs(mesh_.nCells()); - - forAll(cellBbs, celli) - { - cellBbs[celli] = treeBoundBox - ( - mesh_.cells()[celli].points - ( - mesh_.faces(), - mesh_.points() - ) - ); - } + treeBoundBoxList cellBbs(treeDataCell::boxes(mesh_)); const globalIndexAndTransform& globalTransforms = mesh_.globalData().globalTransforms(); @@ -214,7 +202,7 @@ void Foam::InteractionLists::buildInteractionLists() // i.e. a more accurate bounding volume like a OBB or // convex hull or an exact geometrical test. - label c = coupledPatchRangeTree.shapes().cellLabels()[elemI]; + label c = coupledPatchRangeTree.shapes().objectIndex(elemI); ril_[bbI][i] = c; } @@ -237,7 +225,7 @@ void Foam::InteractionLists::buildInteractionLists() // At this point, cellBbsToExchange does not need to be maintained // or distributed as it is not longer needed. - cellBbsToExchange.setSize(0); + cellBbsToExchange.clearStorage(); cellMap().reverseDistribute ( @@ -434,7 +422,7 @@ void Foam::InteractionLists::buildInteractionLists() // i.e. a more accurate bounding volume like a OBB or // convex hull or an exact geometrical test. - label c = coupledPatchRangeTree.shapes().cellLabels()[elemI]; + label c = coupledPatchRangeTree.shapes().objectIndex(elemI); rwfil_[bbI][i] = c; } @@ -586,7 +574,7 @@ void Foam::InteractionLists::buildInteractionLists() for (const label elemi : interactingElems) { - const label c = allCellsTree.shapes().cellLabels()[elemi]; + const label c = allCellsTree.shapes().objectIndex(elemi); // Here, a more detailed geometric test could be applied, // i.e. a more accurate bounding volume like a OBB or @@ -611,7 +599,7 @@ void Foam::InteractionLists::buildInteractionLists() { const label elemi = interactingElems[i]; - const label f = wallFacesTree.shapes().faceLabels()[elemi]; + const label f = wallFacesTree.shapes().objectIndex(elemi); dwfil_[celli][i] = f; } diff --git a/src/lumpedPointMotion/movement/lumpedPointMovement.C b/src/lumpedPointMotion/movement/lumpedPointMovement.C index 5c23a55c3a..75679203d2 100644 --- a/src/lumpedPointMotion/movement/lumpedPointMovement.C +++ b/src/lumpedPointMotion/movement/lumpedPointMovement.C @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2016-2021 OpenCFD Ltd. + Copyright (C) 2016-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -435,25 +435,23 @@ void Foam::lumpedPointMovement::setPatchControl } - treeDataPoint treePoints - ( - lumpedCentres0, - subsetPointIds.sortedToc(), - subsetPointIds.size() - ); - - treeBoundBox bb(lumpedCentres0); bb.inflate(0.01); indexedOctree ppTree ( - treePoints, + treeDataPoint + ( + lumpedCentres0, + subsetPointIds.sortedToc(), + subsetPointIds.size() // subset is optional/required + ), bb, // overall search domain 8, // maxLevel 10, // leafsize 3.0 // duplicity ); + const auto& treePoints = ppTree.shapes(); const scalar searchDistSqr(sqr(GREAT)); @@ -464,7 +462,7 @@ void Foam::lumpedPointMovement::setPatchControl // Store the original pointId, not subset id faceToPoint[patchFacei] = - treePoints.pointLabel + treePoints.objectIndex ( ppTree.findNearest(fc, searchDistSqr).index() ); @@ -655,24 +653,23 @@ void Foam::lumpedPointMovement::setInterpolator } } - treeDataPoint treePoints - ( - lumpedCentres0, - subsetPointIds.sortedToc(), - subsetPointIds.size() - ); - treeBoundBox bb(lumpedCentres0); bb.inflate(0.01); indexedOctree ppTree ( - treePoints, + treeDataPoint + ( + lumpedCentres0, + subsetPointIds.sortedToc(), + subsetPointIds.size() // subset is optional/required + ), bb, // overall search domain 8, // maxLevel 10, // leafsize 3.0 // duplicity ); + const auto& treePoints = ppTree.shapes(); // Searching @@ -692,7 +689,7 @@ void Foam::lumpedPointMovement::setInterpolator // Nearest (original) point id const label nearest = - treePoints.pointLabel + treePoints.objectIndex ( ppTree.findNearest(ptOnMesh, searchDistSqr).index() ); diff --git a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C index 763d488eb9..02ace61b73 100644 --- a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C +++ b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C @@ -651,8 +651,9 @@ void Foam::refinementFeatures::findNearestEdge forAll(edgeTrees_, featI) { const indexedOctree& tree = edgeTrees_[featI]; + const treeDataEdge& treeData = tree.shapes(); - if (tree.shapes().size() > 0) + if (!treeData.empty()) { forAll(samples, sampleI) { @@ -677,13 +678,9 @@ void Foam::refinementFeatures::findNearestEdge ( info.hit(), info.point(), - tree.shapes().edgeLabels()[info.index()] + treeData.objectIndex(info.index()) ); - - const treeDataEdge& td = tree.shapes(); - const edge& e = td.edges()[nearInfo[sampleI].index()]; - - nearNormal[sampleI] = e.unitVec(td.points()); + nearNormal[sampleI] = treeData.line(info.index()).unitVec(); } } } @@ -714,6 +711,7 @@ void Foam::refinementFeatures::findNearestRegionEdge forAll(regionTrees, featI) { const indexedOctree& regionTree = regionTrees[featI]; + const treeDataEdge& treeData = regionTree.shapes(); forAll(samples, sampleI) { @@ -734,19 +732,14 @@ void Foam::refinementFeatures::findNearestRegionEdge if (info.hit()) { - const treeDataEdge& td = regionTree.shapes(); - nearFeature[sampleI] = featI; nearInfo[sampleI] = pointIndexHit ( info.hit(), info.point(), - regionTree.shapes().edgeLabels()[info.index()] + treeData.objectIndex(info.index()) ); - - const edge& e = td.edges()[nearInfo[sampleI].index()]; - - nearNormal[sampleI] = e.unitVec(td.points()); + nearNormal[sampleI] = treeData.line(info.index()).unitVec(); } } } @@ -770,7 +763,7 @@ void Foam::refinementFeatures::findNearestRegionEdge // { // const indexedOctree& tree = pointTrees_[featI]; // -// if (tree.shapes().pointLabels().size() > 0) +// if (!tree.shapes().empty()) // { // forAll(samples, sampleI) // { @@ -779,14 +772,11 @@ void Foam::refinementFeatures::findNearestRegionEdge // scalar distSqr; // if (nearFeature[sampleI] != -1) // { -// label nearFeatI = nearFeature[sampleI]; -// const indexedOctree& nearTree = -// pointTrees_[nearFeatI]; -// label featPointi = -// nearTree.shapes().pointLabels()[nearIndex[sampleI]]; -// const point& featPt = -// operator[](nearFeatI).points()[featPointi]; -// distSqr = magSqr(featPt-sample); +// const nearTree = pointTrees_[nearFeature[sampleI]]; +// distSqr = sample.distSqr +// ( +// nearTree.shapes()[nearIndex[sampleI]] +// ); // } // else // { @@ -822,8 +812,9 @@ void Foam::refinementFeatures::findNearestPoint forAll(pointTrees_, featI) { const indexedOctree& tree = pointTrees_[featI]; + const auto& treeData = tree.shapes(); - if (tree.shapes().pointLabels().size() > 0) + if (!treeData.empty()) { forAll(samples, sampleI) { @@ -848,7 +839,7 @@ void Foam::refinementFeatures::findNearestPoint ( info.hit(), info.point(), - tree.shapes().pointLabels()[info.index()] + treeData.objectIndex(info.index()) ); } } diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C index 1e8d1f287c..2e3b373f4f 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C @@ -3296,8 +3296,13 @@ void Foam::snappySnapDriver::reverseAttractMeshPoints if (nearInfo.hit()) { label pointi = - ppTree.shapes().pointLabels()[nearInfo.index()]; - const point attraction = featPt-pp.localPoints()[pointi]; + ppTree.shapes().objectIndex(nearInfo.index()); + + const point attraction + ( + featPt + - ppTree.shapes()[nearInfo.index()] + ); // Check if this point is already being attracted. If so // override it only if nearer. @@ -3368,10 +3373,13 @@ void Foam::snappySnapDriver::reverseAttractMeshPoints if (nearInfo.hit()) { label pointi = - ppTree.shapes().pointLabels()[nearInfo.index()]; + ppTree.shapes().objectIndex(nearInfo.index()); - const point& pt = pp.localPoints()[pointi]; - const point attraction = featPt-pt; + const point attraction + ( + featPt + - ppTree.shapes()[nearInfo.index()] + ); // - already attracted to feature edge : point always wins // - already attracted to feature point: nearest wins diff --git a/src/meshTools/cellClassification/cellClassification.C b/src/meshTools/cellClassification/cellClassification.C index b8de84255c..af53b27d9a 100644 --- a/src/meshTools/cellClassification/cellClassification.C +++ b/src/meshTools/cellClassification/cellClassification.C @@ -211,7 +211,7 @@ Foam::boolList Foam::cellClassification::markFaces } else { - label facei = faceTree.shapes().faceLabels()[pHit.index()]; + label facei = faceTree.shapes().objectIndex(pHit.index()); if (!cutFace[facei]) { diff --git a/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C b/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C index ecf0abbf89..ec49e31603 100644 --- a/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C +++ b/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C @@ -709,29 +709,20 @@ void Foam::extendedEdgeMesh::nearestFeatureEdgeByType { const PtrList>& edgeTrees = edgeTreesByType(); - info.setSize(edgeTrees.size()); - - labelList sliceStarts(edgeTrees.size()); - - sliceStarts[0] = externalStart_; - sliceStarts[1] = internalStart_; - sliceStarts[2] = flatStart_; - sliceStarts[3] = openStart_; - sliceStarts[4] = multipleStart_; + info.resize(edgeTrees.size()); forAll(edgeTrees, i) { - info[i] = edgeTrees[i].findNearest - ( - sample, - searchDistSqr[i] - ); + const auto& tree = edgeTrees[i]; + const auto& treeData = edgeTrees[i].shapes(); + + info[i] = tree.findNearest(sample, searchDistSqr[i]); // The index returned by the indexedOctree is local to the slice of // edges it was supplied with, return the index to the value in the // complete edge list - info[i].setIndex(info[i].index() + sliceStarts[i]); + info[i].setIndex(treeData.objectIndex(info[i].index())); } } @@ -752,11 +743,11 @@ void Foam::extendedEdgeMesh::allNearestFeaturePoints DynamicList dynPointHit(elems.size()); - forAll(elems, elemI) + const auto& treeData = pointTree().shapes(); + + for (const label index : elems) { - label index = elems[elemI]; - label ptI = pointTree().shapes().pointLabels()[index]; - const point& pt = points()[ptI]; + const point& pt = treeData.centre(index); pointIndexHit nearHit(true, pt, index); @@ -776,37 +767,25 @@ void Foam::extendedEdgeMesh::allNearestFeatureEdges { const PtrList>& edgeTrees = edgeTreesByType(); - info.setSize(edgeTrees.size()); - - labelList sliceStarts(edgeTrees.size()); - - sliceStarts[0] = externalStart_; - sliceStarts[1] = internalStart_; - sliceStarts[2] = flatStart_; - sliceStarts[3] = openStart_; - sliceStarts[4] = multipleStart_; - DynamicList dynEdgeHit(edgeTrees.size()*3); // Loop over all the feature edge types forAll(edgeTrees, i) { + const auto& tree = edgeTrees[i]; + const auto& treeData = tree.shapes(); + // Pick up all the edges that intersect the search sphere - labelList elems = edgeTrees[i].findSphere - ( - sample, - searchRadiusSqr - ); + labelList elems = tree.findSphere(sample, searchRadiusSqr); - forAll(elems, elemI) + for (const label index : elems) { - label index = elems[elemI]; - label edgeI = edgeTrees[i].shapes().edgeLabels()[index]; - const edge& e = edges()[edgeI]; + pointHit hitPoint = treeData.line(index).nearestDist(sample); - pointHit hitPoint = e.line(points()).nearestDist(sample); + // The index returned by indexedOctree is local its slice of + // edges. Return the index into the complete edge list - label hitIndex = index + sliceStarts[i]; + const label hitIndex = treeData.objectIndex(index); dynEdgeHit.append(pointIndexHit(hitPoint, hitIndex)); } diff --git a/src/meshTools/indexedOctree/treeDataFace.C b/src/meshTools/indexedOctree/treeDataFace.C index a5bc221160..f6cec8259c 100644 --- a/src/meshTools/indexedOctree/treeDataFace.C +++ b/src/meshTools/indexedOctree/treeDataFace.C @@ -29,6 +29,7 @@ License #include "treeDataFace.H" #include "polyMesh.H" #include "triangleFuncs.H" +#include // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -40,25 +41,139 @@ namespace Foam } +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Bound boxes corresponding to specified faces +template +static treeBoundBoxList boxesImpl +( + const primitiveMesh& mesh, + const ElementIds& elemIds +) +{ + treeBoundBoxList bbs(elemIds.size()); + + std::transform + ( + elemIds.cbegin(), + elemIds.cend(), + bbs.begin(), + [&](label facei) + { + return treeBoundBox(mesh.points(), mesh.faces()[facei]); + } + ); + + return bbs; +} + + +// Overall bound box for specified faces +template +static treeBoundBox boundsImpl +( + const primitiveMesh& mesh, + const ElementIds& elemIds +) +{ + treeBoundBox bb; + + for (const label facei : elemIds) + { + bb.add(mesh.points(), mesh.faces()[facei]); + } + + return bb; +} + +} // End namespace Foam + + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +Foam::treeBoundBoxList +Foam::treeDataFace::boxes(const primitiveMesh& mesh) +{ + // All faces + return boxesImpl(mesh, labelRange(mesh.nFaces())); +} + + +Foam::treeBoundBoxList +Foam::treeDataFace::boxes +( + const primitiveMesh& mesh, + const labelRange& range +) +{ + return boxesImpl(mesh, range); +} + + +Foam::treeBoundBoxList +Foam::treeDataFace::boxes +( + const primitiveMesh& mesh, + const labelUList& faceIds +) +{ + return boxesImpl(mesh, faceIds); +} + + +Foam::treeBoundBox +Foam::treeDataFace::bounds +( + const primitiveMesh& mesh, + const labelRange& range +) +{ + return boundsImpl(mesh, range); +} + + +Foam::treeBoundBox +Foam::treeDataFace::bounds +( + const primitiveMesh& mesh, + const labelUList& faceIds +) +{ + return boundsImpl(mesh, faceIds); +} + + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -Foam::treeBoundBox Foam::treeDataFace::calcBb(const label facei) const +inline bool +Foam::treeDataFace::usesFace(const label facei) const { + return (!useSubset_ || isTreeFace_.test(facei)); +} + + +inline Foam::treeBoundBox +Foam::treeDataFace::getBounds(const label index) const +{ + const label facei = objectIndex(index); return treeBoundBox(mesh_.points(), mesh_.faces()[facei]); } void Foam::treeDataFace::update() { - isTreeFace_.set(faceLabels_); - if (cacheBb_) { - bbs_.setSize(faceLabels_.size()); - - forAll(faceLabels_, i) + if (useSubset_) { - bbs_[i] = calcBb(faceLabels_[i]); + bbs_ = treeDataFace::boxes(mesh_, faceLabels_); + } + else + { + bbs_ = treeDataFace::boxes(mesh_); } } } @@ -66,6 +181,49 @@ void Foam::treeDataFace::update() // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // +Foam::treeDataFace::treeDataFace +( + const bool cacheBb, + const primitiveMesh& mesh +) +: + mesh_(mesh), + faceLabels_(), + isTreeFace_(), + useSubset_(false), + cacheBb_(cacheBb) +{ + update(); +} + + +Foam::treeDataFace::treeDataFace +( + const bool cacheBb, + const primitiveMesh& mesh, + const labelRange& range +) +: + mesh_(mesh), + faceLabels_(identity(range)), + isTreeFace_(range), + useSubset_(true), + cacheBb_(cacheBb) +{ + update(); +} + + +Foam::treeDataFace::treeDataFace +( + const bool cacheBb, + const polyPatch& patch +) +: + treeDataFace(cacheBb, patch.boundaryMesh().mesh(), patch.range()) +{} + + Foam::treeDataFace::treeDataFace ( const bool cacheBb, @@ -75,7 +233,8 @@ Foam::treeDataFace::treeDataFace : mesh_(mesh), faceLabels_(faceLabels), - isTreeFace_(mesh.nFaces(), false), + isTreeFace_(mesh_.nFaces(), faceLabels_), + useSubset_(true), cacheBb_(cacheBb) { update(); @@ -91,73 +250,45 @@ Foam::treeDataFace::treeDataFace : mesh_(mesh), faceLabels_(std::move(faceLabels)), - isTreeFace_(mesh.nFaces(), false), + isTreeFace_(mesh_.nFaces(), faceLabels_), + useSubset_(true), cacheBb_(cacheBb) { update(); } -Foam::treeDataFace::treeDataFace -( - const bool cacheBb, - const primitiveMesh& mesh -) -: - mesh_(mesh), - faceLabels_(identity(mesh_.nFaces())), - isTreeFace_(mesh.nFaces(), false), - cacheBb_(cacheBb) -{ - update(); -} - - -Foam::treeDataFace::treeDataFace -( - const bool cacheBb, - const polyPatch& patch -) -: - mesh_(patch.boundaryMesh().mesh()), - faceLabels_(identity(patch.range())), - isTreeFace_(mesh_.nFaces(), false), - cacheBb_(cacheBb) -{ - update(); -} - - -Foam::treeDataFace::findNearestOp::findNearestOp -( - const indexedOctree& tree -) -: - tree_(tree) -{} - - -Foam::treeDataFace::findIntersectOp::findIntersectOp -( - const indexedOctree& tree -) -: - tree_(tree) -{} - - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::pointField Foam::treeDataFace::shapePoints() const +Foam::treeBoundBox +Foam::treeDataFace::bounds(const labelUList& indices) const { - pointField cc(faceLabels_.size()); - - forAll(faceLabels_, i) + if (useSubset_) { - cc[i] = mesh_.faceCentres()[faceLabels_[i]]; + treeBoundBox bb; + + for (const label index : indices) + { + const label facei = faceLabels_[index]; + + bb.add(mesh_.points(), mesh_.faces()[facei]); + } + + return bb; } - return cc; + return treeDataFace::bounds(mesh_, indices); +} + + +Foam::tmp Foam::treeDataFace::centres() const +{ + if (useSubset_) + { + return tmp::New(mesh_.faceCentres(), faceLabels_); + } + + return mesh_.faceCentres(); } @@ -181,7 +312,9 @@ Foam::volumeType Foam::treeDataFace::getVolumeType // Find nearest face to sample pointIndexHit info = oc.findNearest(sample, sqr(GREAT)); - if (info.index() == -1) + const label index = info.index(); + + if (index == -1) { FatalErrorInFunction << "Could not find " << sample << " in octree." @@ -190,7 +323,7 @@ Foam::volumeType Foam::treeDataFace::getVolumeType // Get actual intersection point on face - label facei = faceLabels_[info.index()]; + const label facei = objectIndex(index); if (debug & 2) { @@ -238,32 +371,32 @@ Foam::volumeType Foam::treeDataFace::getVolumeType const scalar typDimSqr = mag(area) + VSMALL; - forAll(f, fp) + for (const label fp : f) { - if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < tolSqr) + const scalar relDistSqr = (magSqr(points[fp] - curPt)/typDimSqr); + + if (relDistSqr < tolSqr) { - // Face intersection point equals face vertex fp + // Face intersection point equals face vertex // Calculate point normal (wrong: uses face normals instead of // triangle normals) - const labelList& pFaces = mesh_.pointFaces()[f[fp]]; vector pointNormal(Zero); - for (const label facei : pFaces) + for (const label ptFacei : mesh_.pointFaces()[fp]) { - if (isTreeFace_.test(facei)) + if (usesFace(ptFacei)) { - pointNormal += normalised(mesh_.faceAreas()[facei]); + pointNormal += normalised(mesh_.faceAreas()[ptFacei]); } } if (debug & 2) { - Pout<< " -> face point hit :" << points[f[fp]] + Pout<< " -> face point hit :" << points[fp] << " point normal:" << pointNormal - << " distance:" - << magSqr(points[f[fp]] - curPt)/typDimSqr << endl; + << " distance:" << relDistSqr << endl; } return indexedOctree::getSide ( @@ -272,7 +405,10 @@ Foam::volumeType Foam::treeDataFace::getVolumeType ); } } - if ((magSqr(fc - curPt)/typDimSqr) < tolSqr) + + const scalar relDistSqr = (magSqr(fc - curPt)/typDimSqr); + + if (relDistSqr < tolSqr) { // Face intersection point equals face centre. Normal at face centre // is already average of face normals @@ -280,47 +416,38 @@ Foam::volumeType Foam::treeDataFace::getVolumeType if (debug & 2) { Pout<< " -> centre hit:" << fc - << " distance:" << magSqr(fc - curPt)/typDimSqr << endl; + << " distance:" << relDistSqr << endl; } return indexedOctree::getSide(area, sample - curPt); } - // // 3] Get the 'real' edge the face intersection is on // - const labelList& myEdges = mesh_.faceEdges()[facei]; - - forAll(myEdges, myEdgeI) + for (const label edgei : mesh_.faceEdges()[facei]) { - const edge& e = mesh_.edges()[myEdges[myEdgeI]]; - pointHit edgeHit = - line - ( - points[e.start()], - points[e.end()] - ).nearestDist(sample); + mesh_.edges()[edgei].line(points).nearestDist(sample); + const scalar relDistSqr = edgeHit.point().distSqr(curPt)/typDimSqr; - if ((edgeHit.point().distSqr(curPt)/typDimSqr) < tolSqr) + if (relDistSqr < tolSqr) { // Face intersection point lies on edge e // Calculate edge normal (wrong: uses face normals instead of // triangle normals) - const labelList& eFaces = mesh_.edgeFaces()[myEdges[myEdgeI]]; vector edgeNormal(Zero); - for (const label facei : eFaces) + for (const label eFacei : mesh_.edgeFaces()[edgei]) { - if (isTreeFace_.test(facei)) + if (usesFace(eFacei)) { - edgeNormal += normalised(mesh_.faceAreas()[facei]); + edgeNormal += normalised(mesh_.faceAreas()[eFacei]); } } @@ -348,13 +475,11 @@ Foam::volumeType Foam::treeDataFace::getVolumeType forAll(f, fp) { - pointHit edgeHit = line - ( - points[f[fp]], - fc - ).nearestDist(sample); + pointHit edgeHit = linePointRef(points[f[fp]], fc).nearestDist(sample); - if ((edgeHit.point().distSqr(curPt)/typDimSqr) < tolSqr) + const scalar relDistSqr = edgeHit.point().distSqr(curPt)/typDimSqr; + + if (relDistSqr < tolSqr) { // Face intersection point lies on edge between two face triangles @@ -406,37 +531,32 @@ Foam::volumeType Foam::treeDataFace::getVolumeType } -// Check if any point on shape is inside cubeBb. +// Check if any point on shape is inside searchBox. bool Foam::treeDataFace::overlaps ( const label index, - const treeBoundBox& cubeBb + const treeBoundBox& searchBox ) const { // 1. Quick rejection: bb does not intersect face bb at all - if (cacheBb_) + if + ( + cacheBb_ + ? !searchBox.overlaps(bbs_[index]) + : !searchBox.overlaps(getBounds(index)) + ) { - if (!cubeBb.overlaps(bbs_[index])) - { - return false; - } - } - else - { - if (!cubeBb.overlaps(calcBb(faceLabels_[index]))) - { - return false; - } + return false; } const pointField& points = mesh_.points(); // 2. Check if one or more face points inside - label facei = faceLabels_[index]; + const label facei = objectIndex(index); const face& f = mesh_.faces()[facei]; - if (cubeBb.containsAny(points, f)) + if (searchBox.containsAny(points, f)) { return true; } @@ -452,7 +572,7 @@ bool Foam::treeDataFace::overlaps points[f[fp]], points[f[f.fcIndex(fp)]], fc, - cubeBb + searchBox ); if (triIntersects) @@ -465,6 +585,89 @@ bool Foam::treeDataFace::overlaps } +// Check if any point on shape is inside sphere. +bool Foam::treeDataFace::overlaps +( + const label index, + const point& centre, + const scalar radiusSqr +) const +{ + // 1. Quick rejection: sphere does not intersect face bb at all + if + ( + cacheBb_ + ? !bbs_[index].overlaps(centre, radiusSqr) + : !getBounds(index).overlaps(centre, radiusSqr) + ) + { + return false; + } + + const label facei = objectIndex(index); + + const face& f = mesh().faces()[facei]; + + pointHit nearHit = f.nearestPoint(centre, mesh().points()); + + if (sqr(nearHit.distance()) < radiusSqr) + { + return true; + } + + return false; +} + + +// * * * * * * * * * * * * * * * * Searching * * * * * * * * * * * * * * * * // + +Foam::treeDataFace::findNearestOp::findNearestOp +( + const indexedOctree& tree +) +: + tree_(tree) +{} + + +Foam::treeDataFace::findIntersectOp::findIntersectOp +( + const indexedOctree& tree +) +: + tree_(tree) +{} + + +void Foam::treeDataFace::findNearest +( + const labelUList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& minIndex, + point& nearestPoint +) const +{ + for (const label index : indices) + { + const label facei = objectIndex(index); + + const face& f = mesh().faces()[facei]; + + pointHit nearHit = f.nearestPoint(sample, mesh().points()); + scalar distSqr = sqr(nearHit.distance()); + + if (distSqr < nearestDistSqr) + { + nearestDistSqr = distSqr; + minIndex = index; + nearestPoint = nearHit.point(); + } + } +} + + void Foam::treeDataFace::findNearestOp::operator() ( const labelUList& indices, @@ -475,22 +678,14 @@ void Foam::treeDataFace::findNearestOp::operator() point& nearestPoint ) const { - const treeDataFace& shape = tree_.shapes(); - - for (const label index : indices) - { - const face& f = shape.mesh().faces()[shape.faceLabels()[index]]; - - pointHit nearHit = f.nearestPoint(sample, shape.mesh().points()); - scalar distSqr = sqr(nearHit.distance()); - - if (distSqr < nearestDistSqr) - { - nearestDistSqr = distSqr; - minIndex = index; - nearestPoint = nearHit.point(); - } - } + tree_.shapes().findNearest + ( + indices, + sample, + nearestDistSqr, + minIndex, + nearestPoint + ); } @@ -531,7 +726,7 @@ bool Foam::treeDataFace::findIntersectOp::operator() } } - const label facei = shape.faceLabels_[index]; + const label facei = shape.objectIndex(index); const vector dir(end - start); diff --git a/src/meshTools/indexedOctree/treeDataFace.H b/src/meshTools/indexedOctree/treeDataFace.H index 495143dd6b..fa223fc859 100644 --- a/src/meshTools/indexedOctree/treeDataFace.H +++ b/src/meshTools/indexedOctree/treeDataFace.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -28,15 +28,15 @@ Class Foam::treeDataFace Description - Encapsulation of data needed to search for faces. + Encapsulation of data for searching on faces. SourceFiles treeDataFace.C \*---------------------------------------------------------------------------*/ -#ifndef treeDataFace_H -#define treeDataFace_H +#ifndef Foam_treeDataFace_H +#define Foam_treeDataFace_H #include "face.H" #include "indexedOctree.H" @@ -50,7 +50,7 @@ SourceFiles namespace Foam { -// Forward declaration of classes +// Forward Declarations class polyPatch; /*---------------------------------------------------------------------------*\ @@ -67,6 +67,7 @@ class treeDataFace // Private Data + //- Reference to the mesh const primitiveMesh& mesh_; //- Subset of faces to work on @@ -75,6 +76,9 @@ class treeDataFace //- Inverse of faceLabels. For every mesh whether face is in faceLabels. bitSet isTreeFace_; + //- Use subset of faces (faceLabels) + const bool useSubset_; + //- Whether to precalculate and store face bounding box const bool cacheBb_; @@ -84,8 +88,11 @@ class treeDataFace // Private Member Functions - //- Calculate face bounding box - treeBoundBox calcBb(const label celli) const; + //- True if specified mesh facei is being used + inline bool usesFace(const label facei) const; + + //- Get face bounding box at specified index + inline treeBoundBox getBounds(const label index) const; //- Initialise all member data void update(); @@ -148,7 +155,25 @@ public: ClassName("treeDataFace"); - // Constructors + // Constructors (cachable versions) + + //- Construct from mesh, using all faces + treeDataFace + ( + const bool cacheBb, + const primitiveMesh& mesh + ); + + //- Construct from mesh, using polyPatch faces + treeDataFace(const bool cacheBb, const polyPatch& patch); + + //- Construct from mesh, using specified range of faces + treeDataFace + ( + const bool cacheBb, + const primitiveMesh& mesh, + const labelRange& range + ); //- Construct from mesh, copying subset of faces treeDataFace @@ -166,53 +191,168 @@ public: labelList&& faceLabels ); - //- Construct from mesh. Uses all faces in mesh. - treeDataFace(const bool cacheBb, const primitiveMesh& mesh); - //- Construct from mesh. Uses all faces in patch. - treeDataFace(const bool cacheBb, const polyPatch& patch); + // Constructors (non-caching) + + //- Construct from mesh, using all faces + explicit treeDataFace(const primitiveMesh& mesh) + : + treeDataFace(false, mesh) + {} + + //- Construct from mesh, using polyPatch faces + explicit treeDataFace(const polyPatch& patch) + : + treeDataFace(false, patch) + {} + + //- Construct from mesh, using specified range of faces + treeDataFace(const primitiveMesh& mesh, const labelRange& range) + : + treeDataFace(false, mesh, range) + {} + + //- Construct from mesh, copying subset of faces + treeDataFace(const primitiveMesh& mesh, const labelUList& faceLabels) + : + treeDataFace(false, mesh, faceLabels) + {} + + //- Construct from mesh, moving subset of faces + treeDataFace(const primitiveMesh& mesh, labelList&& faceLabels) + : + treeDataFace(false, mesh, std::move(faceLabels)) + {} + + + // Static Functions + + //- Calculate and return bounding boxes for all mesh faces + static treeBoundBoxList boxes + ( + const primitiveMesh& mesh + ); + + //- Calculate and return bounding boxes for specified range of faces + static treeBoundBoxList boxes + ( + const primitiveMesh& mesh, + const labelRange& range + ); + + //- Calculate and return bounding boxes for specified mesh faces + static treeBoundBoxList boxes + ( + const primitiveMesh& mesh, + const labelUList& faceIds + ); + + //- Return bounding box of specified range of faces + static treeBoundBox bounds + ( + const primitiveMesh& mesh, + const labelRange& range + ); + + //- Return bounding box of specified mesh faces + static treeBoundBox bounds + ( + const primitiveMesh& mesh, + const labelUList& faceIds + ); // Member Functions - // Access + //- Object dimension == 2 (face element) + int nDim() const noexcept { return 2; } - inline const labelList& faceLabels() const - { - return faceLabels_; - } - - inline const primitiveMesh& mesh() const - { - return mesh_; - } - - inline label size() const - { - return faceLabels_.size(); - } - - //- Representative point cloud for all shapes inside - //- (one point per shape) - pointField shapePoints() const; + //- Return bounding box for the specified face indices + treeBoundBox bounds(const labelUList& indices) const; - // Search + // Access - //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. - // Only makes sense for closed surfaces. - volumeType getVolumeType - ( - const indexedOctree&, - const point& - ) const; + //- Reference to the supporting mesh + const primitiveMesh& mesh() const noexcept { return mesh_; } - //- Does (bb of) shape at index overlap bb - bool overlaps - ( - const label index, - const treeBoundBox& sampleBb - ) const; + //- The subset of face ids to use + const labelList& faceLabels() const noexcept { return faceLabels_; } + + //- Use a subset of faces + bool useSubset() const noexcept { return useSubset_; } + + //- Is the effective face selection empty? + bool empty() const noexcept + { + return useSubset_ ? faceLabels_.empty() : !mesh_.nFaces(); + } + + //- The size of the face selection + label size() const noexcept + { + return useSubset_ ? faceLabels_.size() : mesh_.nFaces(); + } + + //- Map from shape index to original (non-subset) face label + label objectIndex(const label index) const + { + return useSubset_ && index >= 0 ? faceLabels_[index] : index; + } + + //- Face at specified shape index + const face& operator[](const label index) const + { + return mesh_.faces()[objectIndex(index)]; + } + + //- Representative point (face centre) at shape index + const point& centre(const label index) const + { + return mesh_.faceCentres()[objectIndex(index)]; + } + + //- Representative point cloud for contained shapes. + //- One point per shape, corresponding to the face centres. + tmp centres() const; + + + // Search + + //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. + // Only makes sense for closed surfaces. + volumeType getVolumeType + ( + const indexedOctree&, + const point& + ) const; + + //- Does (bb of) shape at index overlap searchBox + bool overlaps + ( + const label index, + const treeBoundBox& searchBox + ) const; + + //- Does shape at index overlap sphere + bool overlaps + ( + const label index, + const point& centre, + const scalar radiusSqr + ) const; + + //- Calculates nearest (to sample) point in shape. + // Returns actual point and distance (squared) + void findNearest + ( + const labelUList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& nearestIndex, + point& nearestPoint + ) const; }; diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C index ca60b6525b..51b39a820b 100644 --- a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C +++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C @@ -31,20 +31,51 @@ License #include "triangleFuncs.H" #include "triSurfaceTools.H" #include "triFace.H" +#include + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +template +Foam::treeBoundBoxList +Foam::treeDataPrimitivePatch::boxes(const PatchType& pp) +{ + const auto& points = pp.points(); + + treeBoundBoxList bbs(pp.size()); + + // Like std::transform with [&](const auto& f) + // which may not work with older C++ versions + + { + auto iter = bbs.begin(); + + for (const auto& f : pp) + { + *iter = treeBoundBox(points, f); + ++iter; + } + } + + return bbs; +} + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +template +inline Foam::treeBoundBox +Foam::treeDataPrimitivePatch::getBounds(const label patchFacei) const +{ + return treeBoundBox(patch_.points(), patch_[patchFacei]); +} + + template void Foam::treeDataPrimitivePatch::update() { if (cacheBb_) { - bbs_.setSize(patch_.size()); - - forAll(patch_, i) - { - bbs_[i] = treeBoundBox(patch_.points(), patch_[i]); - } + bbs_ = treeDataPrimitivePatch::boxes(patch_); } } @@ -115,16 +146,20 @@ findSelfIntersectOp::findSelfIntersectOp // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -Foam::pointField Foam::treeDataPrimitivePatch::shapePoints() const +Foam::treeBoundBox +Foam::treeDataPrimitivePatch::bounds +( + const labelUList& indices +) const { - pointField cc(patch_.size()); + treeBoundBox bb; - forAll(patch_, i) + for (const label patchFacei : indices) { - cc[i] = patch_[i].centre(patch_.points()); + bb.add(patch_.points(), patch_[patchFacei]); } - return cc; + return bb; } @@ -168,7 +203,6 @@ Foam::volumeType Foam::treeDataPrimitivePatch::getVolumeType const typename PatchType::face_type& localF = patch_.localFaces()[facei]; const typename PatchType::face_type& 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. @@ -212,7 +246,9 @@ Foam::volumeType Foam::treeDataPrimitivePatch::getVolumeType forAll(f, fp) { - if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < planarTol_) + const scalar relDistSqr = (magSqr(points[f[fp]] - curPt)/typDimSqr); + + if (relDistSqr < planarTol_) { // Face intersection point equals face vertex fp @@ -229,7 +265,9 @@ Foam::volumeType Foam::treeDataPrimitivePatch::getVolumeType const point fc(f.centre(points)); - if ((magSqr(fc - curPt)/typDimSqr) < planarTol_) + const scalar relDistSqr = (magSqr(fc - curPt)/typDimSqr); + + if (relDistSqr < planarTol_) { // Face intersection point equals face centre. Normal at face centre // is already average of face normals @@ -237,7 +275,7 @@ Foam::volumeType Foam::treeDataPrimitivePatch::getVolumeType if (debug & 2) { Pout<< " -> centre hit:" << fc - << " distance:" << magSqr(fc - curPt)/typDimSqr << endl; + << " distance:" << relDistSqr << endl; } return indexedOctree::getSide @@ -248,33 +286,28 @@ Foam::volumeType Foam::treeDataPrimitivePatch::getVolumeType } - // // 3] Get the 'real' edge the face intersection is on // - const labelList& fEdges = patch_.faceEdges()[facei]; - - forAll(fEdges, fEdgeI) + for (const label edgei : patch_.faceEdges()[facei]) { - label edgeI = fEdges[fEdgeI]; - const edge& e = patch_.edges()[edgeI]; - const linePointRef ln(points[mp[e.start()]], points[mp[e.end()]]); - pointHit edgeHit = ln.nearestDist(sample); + pointHit edgeHit = + patch_.meshEdge(edgei).line(points).nearestDist(sample); - if ((edgeHit.point().distSqr(curPt)/typDimSqr) < planarTol_) + const scalar relDistSqr = (edgeHit.point().distSqr(curPt)/typDimSqr); + + if (relDistSqr < planarTol_) { // Face intersection point lies on edge e // Calculate edge normal (wrong: uses face normals instead of // triangle normals) - const labelList& eFaces = patch_.edgeFaces()[edgeI]; - vector edgeNormal(Zero); - forAll(eFaces, i) + for (const label eFacei : patch_.edgeFaces()[edgei]) { - edgeNormal += patch_.faceNormals()[eFaces[i]]; + edgeNormal += patch_.faceNormals()[eFacei]; } if (debug & 2) @@ -303,7 +336,9 @@ Foam::volumeType Foam::treeDataPrimitivePatch::getVolumeType { pointHit edgeHit = linePointRef(points[f[fp]], fc).nearestDist(sample); - if ((edgeHit.point().distSqr(curPt)/typDimSqr) < planarTol_) + const scalar relDistSqr = (edgeHit.point().distSqr(curPt)/typDimSqr); + + if (relDistSqr < planarTol_) { // Face intersection point lies on edge between two face triangles @@ -356,20 +391,20 @@ Foam::volumeType Foam::treeDataPrimitivePatch::getVolumeType } -// Check if any point on shape is inside cubeBb. +// Check if any point on shape is inside searchBox. template bool Foam::treeDataPrimitivePatch::overlaps ( const label index, - const treeBoundBox& cubeBb + const treeBoundBox& searchBox ) const { // 1. Quick rejection: bb does not intersect face bb at all if ( cacheBb_ - ? !cubeBb.overlaps(bbs_[index]) - : !cubeBb.overlaps(treeBoundBox(patch_.points(), patch_[index])) + ? !searchBox.overlaps(bbs_[index]) + : !searchBox.overlaps(getBounds(index)) ) { return false; @@ -381,7 +416,7 @@ bool Foam::treeDataPrimitivePatch::overlaps const pointField& points = patch_.points(); const typename PatchType::face_type& f = patch_[index]; - if (cubeBb.containsAny(points, f)) + if (searchBox.containsAny(points, f)) { return true; } @@ -397,7 +432,7 @@ bool Foam::treeDataPrimitivePatch::overlaps points[f[0]], points[f[1]], points[f[2]], - cubeBb + searchBox ); } else @@ -409,7 +444,7 @@ bool Foam::treeDataPrimitivePatch::overlaps points[f[fp]], points[f[f.fcIndex(fp)]], fc, - cubeBb + searchBox ); if (triIntersects) @@ -437,7 +472,7 @@ bool Foam::treeDataPrimitivePatch::overlaps ( cacheBb_ ? !bbs_[index].overlaps(centre, radiusSqr) - : !treeBoundBox(patch_.points(),patch_[index]).overlaps(centre, radiusSqr) + : !getBounds(index).overlaps(centre, radiusSqr) ) { return false; @@ -459,8 +494,10 @@ bool Foam::treeDataPrimitivePatch::overlaps } +// * * * * * * * * * * * * * * * * Searching * * * * * * * * * * * * * * * * // + template -void Foam::treeDataPrimitivePatch::findNearestOp::operator() +void Foam::treeDataPrimitivePatch::findNearest ( const labelUList& indices, const point& sample, @@ -470,14 +507,11 @@ void Foam::treeDataPrimitivePatch::findNearestOp::operator() point& nearestPoint ) const { - const treeDataPrimitivePatch& shape = tree_.shapes(); - const PatchType& patch = shape.patch(); - - const pointField& points = patch.points(); + const pointField& points = patch_.points(); for (const label index : indices) { - const typename PatchType::face_type& f = patch[index]; + const typename PatchType::face_type& f = patch_[index]; const pointHit nearHit = f.nearestPoint(sample, points); const scalar distSqr = sqr(nearHit.distance()); @@ -492,6 +526,28 @@ void Foam::treeDataPrimitivePatch::findNearestOp::operator() } +template +void Foam::treeDataPrimitivePatch::findNearestOp::operator() +( + const labelUList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& minIndex, + point& nearestPoint +) const +{ + tree_.shapes().findNearest + ( + indices, + sample, + nearestDistSqr, + minIndex, + nearestPoint + ); +} + + template void Foam::treeDataPrimitivePatch::findNearestOp::operator() ( diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.H b/src/meshTools/indexedOctree/treeDataPrimitivePatch.H index 87f0a3dd41..d6a6fda20d 100644 --- a/src/meshTools/indexedOctree/treeDataPrimitivePatch.H +++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.H @@ -82,6 +82,9 @@ class treeDataPrimitivePatch // Private Member Functions + //- Get face bounding box + inline treeBoundBox getBounds(const label patchFacei) const; + //- Initialise all member data void update(); @@ -194,71 +197,126 @@ public: }; - // Constructors + + // Constructors (cachable versions) //- Construct from patch. treeDataPrimitivePatch ( const bool cacheBb, - const PatchType&, + const PatchType& patch, const scalar planarTol ); + // Constructors (non-caching) + + //- Construct from patch. + treeDataPrimitivePatch(const PatchType& patch, const scalar planarTol) + : + treeDataPrimitivePatch(false, patch, planarTol) + {} + + + // Static Functions + + //- Calculate and return bounding boxes for each patch face + static treeBoundBoxList boxes(const PatchType& patch); + + // Member Functions - // Access + //- Object dimension == 2 (face element) + int nDim() const noexcept { return 2; } - label size() const - { - return patch_.size(); - } - - //- Representative point cloud for all shapes inside - //- (one point per shape) - pointField shapePoints() const; - - //- Return access to the underlying patch - const PatchType& patch() const - { - return patch_; - } + //- Return bounding box for the specified face indices + treeBoundBox bounds(const labelUList& indices) const; - // Search + // Access - //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. - // Only makes sense for closed surfaces. - volumeType getVolumeType - ( - const indexedOctree>&, - const point& - ) const; + //- The underlying patch + const PatchType& patch() const noexcept { return patch_; } - //- Does shape at index overlap bb - bool overlaps - ( - const label index, - const treeBoundBox& sampleBb - ) const; + //TDB //- The subset of face ids to use + ///const labelList& faceLabels() const noexcept { labelList::null(); } - //- Does shape at index overlap sphere - bool overlaps - ( - const label index, - const point& centre, - const scalar radiusSqr - ) const; + //- Use a subset of the patch + bool useSubset() const noexcept { return false; } - //- Helper: find intersection of line with shapes - static bool findIntersection - ( - const indexedOctree>& tree, - const label index, - const point& start, - const point& end, - point& intersectionPoint - ); + //- Is the patch empty (no faces)? + bool empty() const { return !patch_.size(); } + + //- The patch size + label size() const { return patch_.size(); } + + //- Map from shape index to original (non-subset) face label + label objectIndex(const label index) const noexcept { return index; } + + //TBD //- Face at specified shape index + ///const typename PatchType::face_type& + ///operator[](label index) const { return patch_[index]; } + + //- Representative point (face centre) at shape index + const point& centre(const label index) const + { + return patch_.faceCentres()[index]; + } + + //- Representative point cloud for contained shapes. + //- One point per shape, corresponding to the face centres. + tmp centres() const + { + return patch_.faceCentres(); + } + + + // Search + + //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. + // Only makes sense for closed surfaces. + volumeType getVolumeType + ( + const indexedOctree>&, + const point& + ) const; + + //- Does shape at index overlap searchBox + bool overlaps + ( + const label index, + const treeBoundBox& searchBox + ) const; + + //- Does shape at index overlap sphere + bool overlaps + ( + const label index, + const point& centre, + const scalar radiusSqr + ) const; + + //- Helper: find intersection of line with shapes + static bool findIntersection + ( + const indexedOctree>& tree, + const label index, + const point& start, + const point& end, + point& intersectionPoint + ); + + //- Calculates nearest (to sample) point in shape. + // Returns actual point and distance (squared) + void findNearest + ( + const labelUList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& nearestIndex, + point& nearestPoint + ) const; }; diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C index 383f6f7932..1db1558ce6 100644 --- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C +++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C @@ -414,6 +414,7 @@ void Foam::mappedPatchBase::findLocalSamples 10, // leafsize 3.0 // duplicity ); + const auto& treeData = boundaryTree.shapes(); forAll(samples, sampleI) { @@ -435,7 +436,8 @@ void Foam::mappedPatchBase::findLocalSamples } else { - point fc(pp[nearInfo.index()].centre(pp.points())); + const point& fc = treeData.centre(nearInfo.index()); + nearInfo.setPoint(fc); near.first().second().first() = sample.distSqr(fc); near.first().second().second() = myRank; diff --git a/src/meshTools/meshSearch/meshSearch.C b/src/meshTools/meshSearch/meshSearch.C index 8a2240d06e..347532c5dd 100644 --- a/src/meshTools/meshSearch/meshSearch.C +++ b/src/meshTools/meshSearch/meshSearch.C @@ -84,7 +84,7 @@ namespace Foam // the new face shares a point with an existing hit face and the // line passes through both faces in the same direction, then this // is also assumed to be a duplicate hit. - const label newFacei = tree_.shapes().faceLabels()[index]; + const label newFacei = tree_.shapes().objectIndex(index); const face& newFace = mesh.faces()[newFacei]; const scalar newDot = mesh.faceAreas()[newFacei] & (end - start); forAll(hits_, hiti) @@ -161,10 +161,8 @@ bool Foam::meshSearch::findNearer { bool nearer = false; - forAll(indices, i) + for (const label pointi : indices) { - label pointi = indices[i]; - scalar distSqr = sample.distSqr(points[pointi]); if (distSqr < nearestDistSqr) @@ -620,9 +618,10 @@ Foam::meshSearch::boundaryTree() const if (!boundaryTreePtr_) { // All boundary faces (not just walls) - labelList bndFaces + const labelRange bndFaces ( - identity(mesh_.nBoundaryFaces(), mesh_.nInternalFaces()) + mesh_.nInternalFaces(), + mesh_.nBoundaryFaces() ); boundaryTreePtr_.reset @@ -818,7 +817,7 @@ Foam::label Foam::meshSearch::findNearestBoundaryFace ); } - return tree.shapes().faceLabels()[info.index()]; + return tree.shapes().objectIndex(info.index()); } else { @@ -867,7 +866,7 @@ Foam::pointIndexHit Foam::meshSearch::intersection if (curHit.hit()) { // Change index into octreeData into face label - curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]); + curHit.setIndex(boundaryTree().shapes().objectIndex(curHit.index())); } return curHit; } @@ -889,7 +888,7 @@ Foam::List Foam::meshSearch::intersections if (!curHit.hit()) break; // Change index into octreeData into face label - curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]); + curHit.setIndex(boundaryTree().shapes().objectIndex(curHit.index())); hits.append(curHit); } diff --git a/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C b/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C index c58a76fdf0..052d1c2838 100644 --- a/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C +++ b/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C @@ -223,6 +223,7 @@ void Foam::searchableExtrudedCircle::findParametricNearest { const edgeMesh& mesh = eMeshPtr_(); const indexedOctree& tree = edgeTree_(); + const auto& treeData = tree.shapes(); const edgeList& edges = mesh.edges(); const pointField& points = mesh.points(); const labelListList& pointEdges = mesh.pointEdges(); @@ -243,11 +244,11 @@ void Foam::searchableExtrudedCircle::findParametricNearest const pointIndexHit startInfo = tree.findNearest(start, maxDistSqr); curvePoints[0] = startInfo.hitPoint(); - axialVecs[0] = edges[startInfo.index()].vec(points); + axialVecs[0] = treeData.line(startInfo.index()).vec(); const pointIndexHit endInfo = tree.findNearest(end, maxDistSqr); curvePoints.last() = endInfo.hitPoint(); - axialVecs.last() = edges[endInfo.index()].vec(points); + axialVecs.last() = treeData.line(endInfo.index()).vec(); diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C index 5e13516813..c140fb39b6 100644 --- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C +++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C @@ -73,12 +73,11 @@ static bool onLine(const Foam::point& p, const linePointRef& line) Foam::pointIndexHit Foam::surfaceFeatures::edgeNearest ( - const point& start, - const point& end, + const linePointRef& line, const point& sample ) { - pointHit eHit = linePointRef(start, end).nearestDist(sample); + pointHit eHit = line.nearestDist(sample); // Classification of position on edge. label endPoint; @@ -96,8 +95,8 @@ Foam::pointIndexHit Foam::surfaceFeatures::edgeNearest // which one. if ( - eHit.point().distSqr(start) - < eHit.point().distSqr(end) + eHit.point().distSqr(line.start()) + < eHit.point().distSqr(line.end()) ) { endPoint = 0; @@ -1255,6 +1254,7 @@ Foam::Map Foam::surfaceFeatures::nearestSamples 10, // leafsize 3.0 // duplicity ); + const auto& treeData = ppTree.shapes(); // From patch point to surface point Map